' fhmrlib.bas - jackord@kw.igs.net - revised 6 Apr 05 - Liberty Basic v4.01
' find Vinit then show motion with and without air resistance

' Initialize Window
  nomainwin
  button#1, "Range", [rng], UL, 5, 5, 50, 20
  button#1, "Motion", [mot], UL, 60, 5, 50, 20
  WindowWidth=460     ' pixel scale 0-450
  WindowHeight=299    ' pixel scale 0-270
  UpperLeftX=100: UpperLeftY=100: kk=0: rr=0
  open "Hit a Homer" for graphics_nsb as #1
  #1 "trapclose [quit]"
  gosub [grid]

[waitHere]
  wait

[rng]
  kk=1: goto [plt]

[mot]
  kk=2: if rr=1 then goto [plt] else goto [waitHere]

[plt]
  gosub [grid]
  if kk=1 then
    m=5.125/16/2.2: pi=4*atn(1): r=9.125*.0254/2/pi
    gr=9.8: k=.87*r*r/m: cmph=3600/5280/.0254/12
    range=450*12*.0254: cf=100/2.54/12: dt=.02: xo=0
    dvi=(gr*range)^.5: vi=0: phi=30: dphi=2.56
    #1 "color blue"
    while xo<range                                    ' Bracket the range
      vi=vi+dvi: gosub [ang]
    wend
    va=vi-dvi: vb=vi
    while vb-va>.02                                   ' Bisection Algorithm
      vi=(va+vb)/2: gosub [ang]
      if xo>range then vb=vi else va=vi
    wend
    #1 "color red": vi=(va+vb)/2: gosub [hit]: rr=1
  else
    dt=.005: #1 "color blue": gosub [hit]
  end if
  range2=int(vi*vi*sin(2*phi*pi/180)/gr/.0254/12): vim=vi*cmph
  #1 "backcolor white"
  if kk=2 then
    #1 "place 120 35"
    #1 "\Flight time "; int(100*t+.5)/100; " sec     Range "; range2; " ft without air"
  end if
  #1 "place 120 20"
  #1 "\Vinit "; int(10*vim+.5)/10; " mph at "; phi; " deg     Range "; int(10*cf*x+.5)/10; " ft    "
goto [waitHere]

[grid]
  #1 "backcolor white ; cls ; down ; color black"
  for i=1 to 9
    #1 "line "; 45*i; " 45 "; 45*i; " 270"
  next i
  for i=1 to 5
    #1 "line 0 "; 45*i; " 450 "; 45*i
  next i
return

[ang]                                               ' Find optimum angle
  dph=dphi: show=0: gosub [hit]
  while abs(dph)>.02
    phi=phi+dph: xo=x: gosub [hit]
    if x<xo then dph=0-dph/2
  wend
  phi=phi+2*dph: show=1
  gosub [hit]: vim=vi*cmph
  #1 "place 120 20"
  #1 "\Vinit "; int(10*vim+.5)/10; " mph at "; phi; " deg     Range "; int(10*cf*x+.5)/10; " ft     "
return

[hit]                                               ' Feynman Algorithm
  t=0: x=0: y=0: q=phi*pi/180
  #1 "place 0 270": x1=0: y1=270: xx1=x1: yy1=y1
  vx=vi*cos(q): vy=vi*sin(q): dvx=0: dvy=0: t=0: vxx=vx: vyy=vy
  gosub [acc]: dvx=ax*dt: dvy=ay*dt: vx=vx+dvx/2: vy=vy+dvy/2
  tt=time$("ms"): del=10
  while y>=0
    t=t+dt: x=x+vx*dt: y=y+vy*dt: x2=int(cf*x): y2=int(270-cf*y)
    if kk=1 then
      if show=1 then #1 "goto "; x2; " "; y2
    else
      #1 "rule xor ; color red ; backcolor red"
      xx=vxx*t: yy=vyy*t-gr*t*t/2: xx2=int(cf*xx): yy2=int(270-cf*yy)
      if y2>270 then y2=270: x2=450
      #1 "place "; x1; " "; y1: #1 "circlefilled 6 ; circle 6"
      #1 "place "; x2; " "; y2: #1 "circlefilled 6 ; circle 6"
      #1 "place "; xx1; " "; yy1: #1 "circlefilled 6 ; circle 6"
      #1 "place "; xx2; " "; yy2: #1 "circlefilled 6 ; circle 6"
      #1 "rule over ; color green ; line "; x1; " "; y1; " "; x2; " "; y2
      #1 "line "; xx1; " "; yy1; " "; xx2; " "; yy2
      #1 "rule xor ; color red"
      x1=x2: y1=y2: xx1=xx2: yy1=yy2
      while time$("ms")<tt: wend: tt=tt+del
    end if
    gosub [acc]: dvx=ax*dt: dvy=ay*dt: vx=vx+dvx: vy=vy+dvy
  wend
  #1 "rule over"
  x=x-vx*y/vy
  scan
return

[acc]                                               ' Newton's Law
  ux=vx+dvx/2: uy=vy+dvy/2: u=(ux*ux+uy*uy)^.5
  ax=0-k*ux*u: ay=0-k*uy*u-gr
return

[quit]
  close #1
end