' forblib.bas - jackord@kw.igs.net - revised 8 Sep 02 - Liberty Basic v3.01

' Initialize Window
  nomainwin
  button#1, "Plot", [plt], UL, 5, 5, 40, 20
  WindowWidth=410     ' pixel scale 0-400
  WindowHeight=329    ' pixel scale 0-300
  UpperLeftX=100: UpperLeftY=100
  open "Projectile launched from north pole" for graphics_nsb as #1
  #1 "trapclose [quit] ; backcolor green"

[waitHere]
  wait

[plt]
  xmin=0-2*10^7: ymax=15*10^6: sc=10^(0-5): pi=4*atn(1)
  R=6.37*10^6: M=5.98*10^24: G=6.672*10^(0-11): GM=G*M: d=R: dmax=R
  #1 "cls ; color black ; down"
  #1 "place 200 150 ; circlefilled "; int(sc*R)
  t=0: x=0: y=R: dt=.1: v=9000                         ' Feynman algorithm
  #1 "color blue ; place "; int(sc*(x-xmin)); " "; int(sc*(ymax-y))
  vx=v/2^.5: vy=vx: vy=vy-GM/(R*R)*dt/2
  while d>=R
    t=t+dt: x=x+vx*dt: y=y+vy*dt
    #1 "goto "; int(sc*(x-xmin)); " "; int(sc*(ymax-y))
    d=(x*x+y*y)^.5: dcube=d*d*d
    if d>dmax then dmax=d
    vx=vx-GM*x/dcube*dt: vy=vy-GM*y/dcube*dt
  wend
  v=(vx*vx+vy*vy)^.5
  phi=atn(y/x)*180/pi: phiv=atn(vy/vx)*180/pi-phi
  if phi<0 then ns$="south" else ns$="north"
  #1 "color black ; backcolor white ; place 10 45"
  #1 "\Lands at "; ns$; " latitude "; left$(str$(abs(phi)), 6); " deg"
  #1 "\with velocity "; left$(str$(v), 6); " m/s"
  #1 "\at "; left$(str$(phiv), 6); " deg to the vertical"
  #1 "\Flight time "; left$(str$(t/60), 6); " min"
  #1 "\dmax = "; left$(str$(dmax/R), 6); " R"
goto [waitHere]

[quit]
  close #1
end