' fgunlib.bas - jackord@kw.igs.net - revised 6 Apr 05 - Liberty Basic v4.01
' a version of fhomer.bas modified to deal with the variation of
'   density with elevation in the atmosphere
' Range finds the initial velocity required to give 40 cm diameter iron
'   sphere a horizontal range of 40 km if there is
'      (a) no atmosphere
'      (b) a constant density atmosphere
'      (c) a variable density atmosphere
' Motion shows the motion for cases (a) and (c)

' 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
  open "Naval Warfare" for graphics_nsb as #1
  #1 "trapclose [quit]"
  kk=0: rr=0: gosub [grid]

[waitHere]
  wait

[rng]
  kk=1: goto [plt]

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

[plt]
  gosub [grid]
  r=.2: pi=4*atn(1): m=4*pi/3*7870*r^3
  gr=9.8: k=.87*r*r/m: range=40000
  if kk=1 then
    dens=0: vo=(gr*range)^.5: phi=45                    ' No Air Resistance
    #1 "color red"
    show=1: vi=vo: gosub [hit]
    #1 "color darkpink"
    for dens=1 to 2
      ntry=0: show=0: xo=0: phi=45
      dvi=vo: if dens=1 then dvi=12*dvi: phi=20
      vi=0: dphi=2
      while xo<range                                    ' Bracket the range
        vi=vi+dvi: gosub [ang]
      wend
      va=vi-dvi: vb=vi
      while vb/va>1.001                                 ' Bisection Algorithm
        vi=(va+vb)/2: gosub [ang]
        if xo>range then vb=vi else va=vi
      wend
      vi=(va+vb)/2: show=1: gosub [hit]
      #1 "color blue"
    next dens
    rr=1
  else
    gosub [hit]
    #1 "color blue ; backcolor white ; place 50 266"
    #1 "\T = "; int(ta); " sec without air      T = "; int(t); " sec with air"
  end if
goto [waitHere]

[grid]
  #1 "cls ; down ; color black"
  for i=1 to 8
    #1 "line "; 50*i; " 0 "; 50*i; " 250"
  next i
  for i=1 to 5
    #1 "line 0 "; 50*i; " 450 "; 50*i
  next i
return

[ang]                                               ' Find optimum angle
  dph=dphi: gosub [hit]
  while abs(dph)>.1
    phi=phi+dph: xo=x: gosub [hit]
    if x<xo then dph=0-dph/2
  wend
  phi=phi+2*dph: ntry=ntry+1
  if ntry<5 then show=1: gosub [hit]: show=0
return

[hit]                                               ' Feynman Algorithm
  dt=.1: t=0: x=0: y=0: q=phi*pi/180: h=100000/9.8/1.21
  #1 "place 0 250": x1=0: y1=250: xx1=x1: yy1=y1
  vx=vi*cos(q): vy=vi*sin(q)
  if kk=1 then
    if dens=0 then                                    ' No Air Resistance
      vy=vy-gr*dt/2
      while y>=0
        x=x+vx*dt: y=y+vy*dt: vy=vy-gr*dt
        #1 "goto "; int(x/100); " "; int(250-y/100)
      wend
    end if
    if dens>0 then                                    ' Sea-level or variable
      dvx=0-k*vx*vi*dt: vx=vx+dvx/2                   '   density of air...
      dvy=0-(gr+k*vy*vi)*dt: vy=vy+dvy/2
      while y>=0
        x=x+vx*dt: y=y+vy*dt
        if show=1 then print#1, "goto "; int(x/100); " "; int(250-y/100)
        ux=vx+dvx/2: uy=vy+dvy/2: u=(ux*ux+uy*uy)^.5
        if dens=1 then cc=1 else cc=exp(0-y/h)        ' ...choose which
        dvx=0-cc*k*ux*u*dt: vx=vx+dvx
        dvy=0-(gr+cc*k*uy*u)*dt: vy=vy+dvy
      wend
    end if
    if show=1 then
      #1 "place "; 152*dens; " 266"
      #1 "\"; int(vi); " m/s at "; int(10*phi)/10; " deg  "
    end if
  else
    #1 "rule xor ; color red ; backcolor red"
    xx=0: yy=0: vxx=vo*cos(pi/4): vyy=vxx: ta=2*vyy/gr
    dvx=0-k*vx*vi*dt: vx=vx+dvx/2
    dvy=0-(gr+k*vy*vi)*dt: vy=vy+dvy/2
    tt=time$("ms"): del=4
    while y>=0
      t=t+dt: x=x+vx*dt: y=y+vy*dt
      xx=vxx*t: yy=vyy*t-gr*t*t/2
      x2=int(x/100): y2=int(250-y/100): xx2=int(xx/100): yy2=int(250-yy/100)
      if y2=250 then y2=250: x2=400: #1 "rule over"
      if yy2>250 then yy2=250: xx2=400
      #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
      xx1=xx2: yy1=yy2
      #1 "rule xor ; color red ; discard"
      x1=x2: y1=y2
      ux=vx+dvx/2: uy=vy+dvy/2: u=(ux*ux+uy*uy)^.5
      cc=exp(0-y/h)
      dvx=0-cc*k*ux*u*dt: vx=vx+dvx
      dvy=0-(gr+cc*k*uy*u)*dt: vy=vy+dvy
      while time$("ms")<tt: wend: tt=tt+del
    wend
  end if
  x=x-vx*y/vy
  scan
return

[quit]
  close #1
end