' fudplib.bas - jackord@kw.igs.net - revised 6 Sep 02 - Liberty Basic v3.01
'   ball thrown upward at its terminal velocity (30 m/s)
'   eps file option for output on a PostScript printer

' Initialize Window
  nomainwin
  button#1, "Plot", [plt], UL, 5, 5, 40, 20
  button#1, "File", [eps], UL, 50, 5, 40, 20
  WindowWidth=430     ' pixel scale 0-420
  WindowHeight=329    ' pixel scale 0-300
  UpperLeftX=100: UpperLeftY=100
  open "Air Resistance" for graphics_nsb as #1
  #1 "trapclose [quit]"
  kk=0

[waitHere]
  wait

[plt]
  #1 "cls ; down ; color red ; place 0 300"
  gr=9.8: vt=30: dt=.01: h=0
  t=0: y=0: v=vt                                       ' Feynman algorithm
  v=v-gr*dt/2
  while y>=0                                           ' test for return
    t=t+dt: y=y+v*dt: v=v-gr*dt
    #1 "goto "; int(60*t); " "; int(6*(50-y))
  wend
  t=0: y=0: v=vt: dv=0                                 ' Feynman algorithm
  #1 "color blue ; place 0 300"
  gosub [acc]
  dv=a*dt: v=v+dv/2
  while y>=0                                           ' test for return
    t=t+dt: y=y+v*dt: gosub [acc]: dv=a*dt: v=v+dv
    #1 "goto "; int(60*t); " "; int(6*(50-y))
    if v-dv>0 and v<=0 then h=y
  wend
  v=v-(.5+y/v/dt)*dv                                   ' correct v
  vret=vt/2^.5                                         ' analytic v return
  ym=vt*vt*log(2)/2/gr                                 ' analytic y max
  #1 "place 60 200"
  #1 "\Initial speed = "; str$(vt)
  #1 "\Return speed = "; str$(0-v)
  #1 "\Analytic ratio = "; str$(0-v/vret)
  #1 "\ "
  #1 "\Maximum height = "; str$(h)
  #1 "\Analytic ratio = "; str$(h/ym)
  kk=1
goto [waitHere]

[acc]                                                  ' Newton's Law
  u=v+dv/2: a=0-(1+u*abs(u)/(vt*vt))*gr
return

[quit]
  close #1
end

[eps] ' Write PostScript .eps File ------------------------------------------
  if kk=0 then goto [waitHere]
'Set user scales, axis labels, title
  xmax=7: nx=7: ymax=50: ny=5: wide=420: high=300
  xt$="TIME  seconds": yt$="HEIGHT  meters"
'Set origin of hardcopy graph in points (size in points = size in pixels)
  xorg=120: yorg=360
'Set axis tic option (0 for grid, 1 for tics)
  tic=1

  open "plot.eps" for output as #2
  #2 "%!PS-Adobe-3.0 EPSF-2.0"
  #2 "%%Creator: Jack Ord"
  #2 "%%Title: Air Resistance"
  #2 "%%CreationDate: 6 Sep 02 99"
  #2 "%%Pages: 0"
  #2 "%%BoundingBox: 60 200 580 720"
  #2 "/scaleHelv {/Helvetica findfont exch scalefont setfont} def"
  #2 "/centerx"
  #2 "{dup stringwidth pop 2 div neg 0 rmoveto show} def"
  #2 "/justify"
  #2 "{dup stringwidth pop 6 add neg 0 rmoveto show} def"
  #2 "/centery"
  #2 "{dup stringwidth pop 2 div neg 0 exch rmoveto"
  #2 "90 rotate show -90 rotate} def"
  #2 "14 scaleHelv .5 setlinewidth"
  if tic=0 then zt=high else zt=9
  for i=0 to nx: x=wide*i/nx+xorg                      ' x axis
    if i>0 and i<nx then
      #2 "newpath "; x; " "; yorg; " moveto 0 "; zt; " rlineto stroke"
    end if
    #2 x; " "; yorg; " 16 sub moveto (";
    #2 xmax*i/nx; : print#2, ") centerx"
  next i
  if tic=1 then                                        ' tic option
    for i=1 to nx-1: x=wide*i/nx+xorg
      #2 "newpath "; x; " "; yorg + high; " moveto 0 "; 0-zt; " rlineto stroke"
    next i
    for i=1 to ny-1: y=high*i/ny+yorg
      #2 "newpath "; xorg + wide; " "; y; " moveto "; 0-zt; " 0 rlineto stroke"
    next i
  end if
  for i=0 to ny: y=high*i/ny+yorg                      ' y axis
    if tic=0 then zt=wide
    if i>0 and i<ny then
      #2 "newpath "; xorg; " "; y; " moveto "; zt; " 0 rlineto stroke"
    end if
    #2 xorg; " "; y-5; " moveto (";
    #2 ymax*i/ny; : print#2, ") justify"
  next i
  x=xorg+wide/2: y=yorg-36                             ' x axis caption
  #2 x; " "; y; " moveto ("; xt$; ") centerx"
  x=xorg: y=yorg+high/2                                ' y axis caption
  #2 x; " "; y; " moveto ("; str$(ymax); ") stringwidth pop neg 14 sub 0 rmoveto (";
  #2 yt$; ") centery"
  #2 "newpath "; xorg; " "; yorg; " moveto "; wide; " 0 rlineto 0 "
  #2 high; " rlineto "; 0-wide; " 0 rlineto closepath 1.2 setlinewidth"
  #2 "stroke 12 scaleHelv .9 setlinewidth"

  h=0: t=0: y=0: v=vt                                  ' Feynman algorithm
  gosub [lstrt]
  v=v-gr*dt/2
  while y>=0                                           ' test for return
    t=t+dt: y=y+v*dt: v=v-gr*dt
    gosub [lcont]
  wend
  t=0: y=0: v=vt: dv=0                                 ' Feynman algorithm
  gosub [lstrt]
  gosub [acc]
  dv=a*dt: v=v+dv/2
  while y>=0                                           ' test for return
    t=t+dt: y=y+v*dt: gosub [acc]: dv=a*dt: v=v+dv
    gosub [lcont]
    if v-dv>0 and v<=0 then h=y
  wend
  v=v-(.5+y/v/dt)*dv                                   ' correct v
  z$="190 470": lbl$="Initial speed = "+str$(vt): gosub [label]
  z$="190 455": lbl$="Return speed = "+str$(0-v): gosub [label]
  z$="190 440": lbl$="Analytic ratio = "+str$(0-v/vret): gosub [label]
  z$="190 410": lbl$="Maximum height = "+str$(h): gosub [label]
  z$="190 395": lbl$="Analytic ratio = "+str$(h/ym): gosub [label]
  z$="120 300": lbl$="Resisted Motion in One Dimension": gosub [label]
  #2 "stroke showpage": #2 "%%EOF": CLOSE #2
  kk=0
goto [waitHere]

[label]
  #2 "stroke newpath "; z$; " moveto ("; lbl$; ") show"
return

[lstrt]
  gosub [sc]
  #2 "stroke newpath "; z$; " moveto"
return

[lcont]
  gosub [sc]
  #2 z$; " lineto"
return

[sc]
  zx=xorg+int(10*wide*t/xmax)/10
  zy=yorg+int(10*high*y/ymax)/10
  z$=str$(zx)+" "+str$(zy)
return