CLS
REM This program does solution of a set of four equations
pi = 3.141593
alpha = 30 * (pi / 180)   ' half angle of bow
INPUT "Half angle of bow (degrees) <30>"; x$
IF x$ <> "" THEN alpha = VAL(x$) * (pi / 180)
beta = 20 * (pi / 180)  'half angle of stern
INPUT "Half angle of stern (degrees) <20>"; x$
IF x$ <> "" THEN beta = VAL(x$) * (pi / 180)
b = 1.625 'b is in inches and is the maximum half width of the boat hull
PRINT "Note: the beam of the boat is defined as the maximum width of the boat"
INPUT "Beam half width in inches <1-5/8>"; x$
IF x$ <> "" THEN b = VAL(x$)
l1 = 7  'length where hull is at maximum width
INPUT "Distance in inches from bow where beam lies <7>"; x$
IF x$ <> "" THEN l1 = VAL(x$)
s = .5
INPUT "Input stern half width in inches <.5> (any value ok - including zero)"; x$
IF x$ <> "" THEN s = VAL(x$)' s is the stern half width in inches
l = 14  'total length of boat          
INPUT "Total length of boat in inches <14>"; x$
IF x$ <> "" THEN l = VAL(x$)
a1 = TAN(alpha)
c11 = l1 ^ 2: c12 = l1 ^ 3: c13 = l1 ^ 4: c14 = l1 ^ 5
c21 = 2 * l1: c22 = 3 * (l1 ^ 2): c23 = 4 * (l1 ^ 3): c24 = 5 * (l1 ^ 4)
c31 = 2 * l: c32 = 3 * (l ^ 2): c33 = 4 * (l ^ 3): c34 = 5 * (l ^ 4)
c41 = l ^ 2: c42 = l ^ 3: c43 = l ^ 4: c44 = l ^ 5
GOSUB mydet
olddet = det
'PRINT det
'INPUT a$

'now reset the first column
c11 = b - l1 * TAN(alpha): c21 = -TAN(alpha): c31 = -TAN(beta) - TAN(alpha)
c41 = s - l * TAN(alpha)
GOSUB mydet

a2 = det / olddet

c11 = l1 ^ 2: c12 = l1 ^ 3: c13 = l1 ^ 4: c14 = l1 ^ 5
c21 = 2 * l1: c22 = 3 * l1 ^ 2: c23 = 4 * l1 ^ 3: c24 = 5 * l1 ^ 4
c31 = 2 * l: c32 = 3 * l ^ 2: c33 = 4 * l ^ 3: c34 = 5 * l ^ 4
c41 = l ^ 2: c42 = l ^ 3: c43 = l ^ 4: c44 = l ^ 5
'now reset the second column
c12 = b - l1 * TAN(alpha): c22 = -TAN(alpha): c32 = -TAN(beta) - TAN(alpha)
c42 = s - l * TAN(alpha)
GOSUB mydet
a3 = det / olddet

c11 = l1 ^ 2: c12 = l1 ^ 3: c13 = l1 ^ 4: c14 = l1 ^ 5
c21 = 2 * l1: c22 = 3 * l1 ^ 2: c23 = 4 * l1 ^ 3: c24 = 5 * l1 ^ 4
c31 = 2 * l: c32 = 3 * l ^ 2: c33 = 4 * l ^ 3: c34 = 5 * l ^ 4
c41 = l ^ 2: c42 = l ^ 3: c43 = l ^ 4: c44 = l ^ 5
'now reset the third column
c13 = b - l1 * TAN(alpha): c23 = -TAN(alpha): c33 = -TAN(beta) - TAN(alpha)
c43 = s - l * TAN(alpha)
GOSUB mydet
a4 = det / olddet

c11 = l1 ^ 2: c12 = l1 ^ 3: c13 = l1 ^ 4: c14 = l1 ^ 5
c21 = 2 * l1: c22 = 3 * l1 ^ 2: c23 = 4 * l1 ^ 3: c24 = 5 * l1 ^ 4
c31 = 2 * l: c32 = 3 * l ^ 2: c33 = 4 * l ^ 3: c34 = 5 * l ^ 4
c41 = l ^ 2: c42 = l ^ 3: c43 = l ^ 4: c44 = l ^ 5
'now reset the fourth column
c14 = b - l1 * TAN(alpha): c24 = -TAN(alpha): c34 = -TAN(beta) - TAN(alpha)
c44 = s - l * TAN(alpha)
GOSUB mydet
a5 = det / olddet


'now draw the result
sw = 600: sh = 480
SCREEN 12
oldx% = 0: oldy% = sh / 2  'draw the top line
FOR ii = 0 TO 100
x = l * ii / 100
wid = a1 * x + a2 * x * x + a3 * x ^ 3 + a4 * x ^ 4 + a5 * x ^ 5
xout% = 2.7 * ii
yout% = sh / 2 - wid * 20
'PSET (xout%, yout%)
LINE (oldx%, oldy%)-(xout%, yout%)
oldx% = xout%: oldy% = yout%
NEXT ii
xout% = 2.7 * 100: yout1 = sh / 2 - s * 20: yout2 = sh / 2 + s * 20
LINE (xout%, yout1)-(xout%, yout2)  'draw the stern
LINE (0, sh / 2)-(270, sh / 2)      'draw a center line

oldx% = 0: oldy% = sh / 2  'draw the bottom line
FOR ii = 0 TO 100
x = l * ii / 100
wid = a1 * x + a2 * x * x + a3 * x ^ 3 + a4 * x ^ 4 + a5 * x ^ 5
xout% = 2.7 * ii
yout% = sh / 2 + wid * 20
'PSET (xout%, yout%)
LINE (oldx%, oldy%)-(xout%, yout%)
oldx% = xout%: oldy% = yout%
NEXT ii

INPUT a$
CLS
PRINT "X pos.(in.) from bow", "half beam(in.)"
FOR ii = 0 TO l
x = ii: wid = a1 * x + a2 * x * x + a3 * x ^ 3 + a4 * x ^ 4 + a5 * x ^ 5
'convert the width to inches and 32nds

my2 = 0: my4 = 0: my8 = 0: my16 = 0
myin = INT(wid + 1 / 64): my32 = INT((wid - myin) * 32 + .5): my16 = 0
IF my32 MOD 2 = 0 THEN my16 = INT(my32 / 2)
IF my16 MOD 2 = 0 THEN my8 = INT(my16 / 2)
IF my8 MOD 2 = 0 THEN my4 = INT(my8 / 2)
IF my4 MOD 2 = 0 THEN my2 = INT(my4 / 2)
PRINT x, , myin; "-";
IF my2 <> 0 THEN PRINT my2; "/2 in.": GOTO skipit
IF my4 <> 0 THEN PRINT my4; "/4 in.": GOTO skipit
IF my8 <> 0 THEN PRINT my8; "/8 in.": GOTO skipit
IF my16 <> 0 THEN PRINT my16; "/16 in.": GOTO skipit
IF my32 = 0 THEN PRINT " in.": GOTO skipit
PRINT my32; "/32 in."
skipit:
NEXT ii
INPUT a$
CLS

INPUT "Print values out to a text file called boat.txt (y/n)"; a$
IF a$ = "y" OR a$ = "Y" THEN
OPEN "boat.txt" FOR OUTPUT AS #1
PRINT #1, "X pos.(in.) from bow", "half beam(in.)"
FOR ii = 0 TO l
x = ii: wid = a1 * x + a2 * x * x + a3 * x ^ 3 + a4 * x ^ 4 + a5 * x ^ 5
'convert the width to inches and 32nds

my2 = 0: my4 = 0: my8 = 0: my16 = 0
myin = INT(wid + 1 / 64): my32 = INT((wid - myin) * 32 + .5): my16 = 0
IF my32 MOD 2 = 0 THEN my16 = INT(my32 / 2)
IF my16 MOD 2 = 0 THEN my8 = INT(my16 / 2)
IF my8 MOD 2 = 0 THEN my4 = INT(my8 / 2)
IF my4 MOD 2 = 0 THEN my2 = INT(my4 / 2)
PRINT #1, x, , myin; "-";
IF my2 <> 0 THEN PRINT #1, my2; "/2 in.": GOTO skipit1
IF my4 <> 0 THEN PRINT #1, my4; "/4 in.": GOTO skipit1
IF my8 <> 0 THEN PRINT #1, my8; "/8 in.": GOTO skipit1
IF my16 <> 0 THEN PRINT #1, my16; "/16 in.": GOTO skipit1
IF my32 = 0 THEN PRINT #1, " in.": GOTO skipit1
PRINT #1, my32; "/32 in."
skipit1:
NEXT ii
CLOSE #1
END IF

INPUT "Create a PostScript file (y/n)"; a$
IF a$ = "y" OR a$ = "Y" THEN
OPEN "boat.eps" FOR OUTPUT AS #1
PRINT #1, "%!PS-Adobe-2.0 EPSF-1.2"
PRINT #1, "%%Title: Boat Design"
PRINT #1, "%%Creator: God"
PRINT #1, "%%BoundingBox: 0 0 540 720"
PRINT #1, "%%EndComments"
PRINT #1, "270 36 translate"
PRINT #1, "0.0 setgray"
PRINT #1, "/Helvetica findfont 11 scalefont setfont"
PRINT #1, "72 36 moveto (Measurements are the half-beam widths) show"
PRINT #1, "/l "; l; " def"
PRINT #1, "/nn l cvi def %To be used for looping control"
PRINT #1, "/mm l def nn 9 gt { /mm 9 def } if"
PRINT #1, "/xold 0 def   /yold 0 def"
PRINT #1, "/a1 "; a1; " def "; "/a2 "; a2; " def "; "/a3 "; a3; " def "
PRINT #1, "/a4 "; a4; " def "; "/a5 "; a5; " def"
PRINT #1, "1 1 mm { /x exch def"
PRINT #1, "/wid a1 x mul a2 x x mul mul add a3 x x x mul mul mul add"
PRINT #1, "a4 x x x x mul mul mul mul add a5 x x x x x mul mul mul"
PRINT #1, "mul mul add def %calculate the half beam using polynomial"
PRINT #1, "/yp x 72 mul def   /xp wid 72 mul def  %convert coordinates to points"
PRINT #1, "xold yold moveto xp yp lineto stroke"
PRINT #1, "xold neg yold moveto xp neg yp lineto stroke"
PRINT #1, "%Also draw a section line"
PRINT #1, "xp yp moveto xp 2 mul neg 0 rlineto stroke"
PRINT #1, "/xold xp def   /yold yp def"
PRINT #1, "%Now label each section w. half-beam measurement"
PRINT #1, "/mystr (0-00/16 in.) def"
PRINT #1, "/myin wid cvi 48 add def mystr 0 myin put"
PRINT #1, "/my16 wid wid cvi cvr sub 16 mul .5 add cvi def"
PRINT #1, "/my16ms my16 10 idiv 48 add def mystr 2 my16ms put"
PRINT #1, "/my16ls my16 10 mod 48 add def mystr 3 my16ls put"
PRINT #1, "xp 24 add yp moveto mystr show"
PRINT #1, "} for"
PRINT #1, "0 0 moveto 0 yp lineto stroke %draw centerline"
'If L>9 inches, we must finish the design on a second piece of paper
IF l <= 9 THEN
PRINT #1, "%If L is noninteger, we must finish the job"
PRINT #1, "/x l def"
PRINT #1, "/wid a1 x mul a2 x x mul mul add a3 x x x mul mul mul add"
PRINT #1, "a4 x x x x mul mul mul mul add a5 x x x x x mul mul mul"
PRINT #1, "mul mul add def %calculate the half beam using polynomial"
PRINT #1, "/yp x 72 mul def   /xp wid 72 mul def  %convert coordinates to points"
PRINT #1, "xold yold moveto xp yp lineto stroke"
PRINT #1, "xold neg yold moveto xp neg yp lineto stroke"
PRINT #1, "%Now draw the stern"
PRINT #1, "xp yp moveto xp 2 mul neg yp lineto stroke"
PRINT #1, "0 yold moveto 0 yp lineto stroke %draw centerline"
PRINT #1, "%Now label stern w. half-beam measurement"
PRINT #1, "/mystr (0-00/16 in.) def"
PRINT #1, "/myin wid cvi 48 add def mystr 0 myin put"
PRINT #1, "/my16 wid wid cvi cvr sub 16 mul .5 add cvi def"
PRINT #1, "/my16ms my16 10 idiv 48 add def mystr 2 my16ms put"
PRINT #1, "/my16ls my16 10 mod 48 add def mystr 3 my16ls put"
PRINT #1, "xp 24 add yp moveto mystr show"

END IF
CLOSE #1

'if L>9 then open a second EPS file and finish the job
IF l > 9 THEN

OPEN "boat2.eps" FOR OUTPUT AS #1
PRINT #1, "%!PS-Adobe-2.0 EPSF-1.2"
PRINT #1, "%%Title: Boat Design - 2nd Page"
PRINT #1, "%%Creator: God"
PRINT #1, "%%BoundingBox: 0 0 540 720"
PRINT #1, "%%EndComments"
PRINT #1, "270 36 translate"
PRINT #1, "0.0 setgray"
PRINT #1, "/Helvetica findfont 11 scalefont setfont"
PRINT #1, "/l "; l; " def"
PRINT #1, "/nn l cvi def %To be used for looping control"
PRINT #1, "/xold 0 def   /yold 0 def"
PRINT #1, "/a1 "; a1; " def "; "/a2 "; a2; " def "; "/a3 "; a3; " def "
PRINT #1, "/a4 "; a4; " def "; "/a5 "; a5; " def"
PRINT #1, "9 1 nn { /x exch def"
PRINT #1, "/wid a1 x mul a2 x x mul mul add a3 x x x mul mul mul add"
PRINT #1, "a4 x x x x mul mul mul mul add a5 x x x x x mul mul mul"
PRINT #1, "mul mul add def %calculate the half beam using polynomial"
PRINT #1, "/yp x 9 sub 72 mul def   /xp wid 72 mul def  %convert coordinates to points"
PRINT #1, "xold yold moveto xp yp lineto stroke"
PRINT #1, "xold neg yold moveto xp neg yp lineto stroke"
PRINT #1, "xp yp moveto xp 2 mul neg 0 rlineto stroke %draw section line"
PRINT #1, "/xold xp def   /yold yp def"
PRINT #1, "%Now label each section w. half-beam measurement"
PRINT #1, "/mystr (0-00/16 in.) def"
PRINT #1, "/myin wid cvi 48 add def mystr 0 myin put"
PRINT #1, "/my16 wid wid cvi cvr sub 16 mul .5 add cvi def"
PRINT #1, "/my16ms my16 10 idiv 48 add def mystr 2 my16ms put"
PRINT #1, "/my16ls my16 10 mod 48 add def mystr 3 my16ls put"
PRINT #1, "xp 24 add yp moveto mystr show"
PRINT #1, "} for"
PRINT #1, "%If L is noninteger, we must finish the job"
PRINT #1, "/x l def"
PRINT #1, "/wid a1 x mul a2 x x mul mul add a3 x x x mul mul mul add"
PRINT #1, "a4 x x x x mul mul mul mul add a5 x x x x x mul mul mul"
PRINT #1, "mul mul add def %calculate the half beam using polynomial"
PRINT #1, "/yp x 9 sub 72 mul def   /xp wid 72 mul def  %convert coordinates to points"
PRINT #1, "xold yold moveto xp yp lineto stroke"
PRINT #1, "xold neg yold moveto xp neg yp lineto stroke"
PRINT #1, "%Now draw the stern"
PRINT #1, "xp yp moveto xp neg yp lineto stroke"
PRINT #1, "0 0 moveto 0 yp rlineto stroke %draw centerline"
PRINT #1, "%Now label stern w. half-beam measurement"
PRINT #1, "/mystr (0-00/16 in.) def"
PRINT #1, "/myin wid cvi 48 add def mystr 0 myin put"
PRINT #1, "/my16 wid wid cvi cvr sub 16 mul .5 add cvi def"
PRINT #1, "/my16ms my16 10 idiv 48 add def mystr 2 my16ms put"
PRINT #1, "/my16ls my16 10 mod 48 add def mystr 3 my16ls put"
PRINT #1, "xp 24 add yp moveto mystr show"

CLOSE #1
END IF


END IF

END



mydet:  'subroutine to evaluate a 4x4 determinant
det = c11 * (c22 * c33 * c44 + c23 * c34 * c42 + c24 * c32 * c43)
det = det - c11 * (c42 * c33 * c24 + c43 * c34 * c22 + c44 * c32 * c23)
det = det - c21 * (c12 * c33 * c44 + c13 * c34 * c42 + c14 * c32 * c43)
det = det + c21 * (c42 * c33 * c14 + c43 * c34 * c12 + c44 * c32 * c13)
det = det + c31 * (c12 * c23 * c44 + c13 * c24 * c42 + c14 * c22 * c43)
det = det - c31 * (c42 * c23 * c14 + c43 * c24 * c12 + c44 * c22 * c13)
det = det - c41 * (c12 * c23 * c34 + c13 * c24 * c32 + c14 * c22 * c33)
det = det + c41 * (c32 * c23 * c14 + c33 * c24 * c12 + c34 * c22 * c13)
RETURN