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

' Initialize Window
  nomainwin
  button#1, "-2", [m2], UL, 26, 10, 20, 20
  button#1, "-1", [m1], UL, 56, 10, 20, 20
  button#1, "+1", [p1], UL, 86, 10, 20, 20
  button#1, "Orbit", [p2], UL, 140, 10, 50, 20
  WindowWidth=442     ' pixel scale 0-432
  WindowHeight=389    ' pixel scale 0-360
  UpperLeftX=100: UpperLeftY=100
  open "Apsidal-Angle Integral" for graphics_nsb as #1   ' no slide bars
  #1 "trapclose [quit] ; down"
  #1 "place 10 45": #1 "\Force Power Law"
  #1 "line 216 180 432 180 ; line 216 0 216 360"
  dim dtheta(1001)

[waitHere]
  wait

[m2]
  n=-2: norbit=1: goto [apint]

[m1]
  n=-1: norbit=16: goto [apint]

[p1]
  n=1: norbit=2: goto [apint]

[p2]
  if orbit=1 then
    n=2: goto [apint]
  else
    goto [waitHere]
  end if

[apint]
  #1 "cls": #1 "color black"
  np=1000: cf=5: pi=4*atn(1)
  if n=2 then
    #1 "backcolor green"
    #1 "place 216 180": #1 "circlefilled 10"
    #1 "place 216 180": #1 "circle 10"
    #1 "backcolor white"
    #1 "color red"
    #1 "place 216 180": #1 "circle 180"
    #1 "place 216 180": #1 "circle "; int(180*c+.5)
    #1 "color blue"
    r=1: theta=0: #1 "place 396 180"
    for j=1 to norbit
      for i=1 to np
        r=r-dr: theta=theta+dtheta(i)
        #1 "goto "; int(180*(1.2+r*cos(theta))); " "; int(180*(1-r*sin(theta)))
      next i
      for i=1 to np
        r=r+dr: theta=theta+dtheta(np+1-i)
        #1 "goto "; int(180*(1.2+r*cos(theta))); " "; int(180*(1-r*sin(theta)))
      next i
    next j
    orbit=0: goto [waitHere]
  end if
  #1 "place 10 45": #1 "\Force Power Law"
  #1 "line 216 180 432 180"
  #1 "line 216 0 216 360"
  #1 "color green"
  L=.488: L2=L*L/2: x=1.2: y=L2/x/x
  #1 "place 432 "; int(180*(1-y/cf))
  while y<cf
    x=x-.002: y=L2/x/x
    #1 "goto "; int(180*(1.2+x)); " "; int(180*(1-y/cf))
  wend
  #1 "color darkpink"
  x=1.2: gosub [pot]
  #1 "place 432 "; int(180*(1-y/cf))
  while y>0-cf and x>0
    x=x-.002: gosub [pot]
    #1 "goto "; int(180*(1.2+x)); " "; int(180*(1-y/cf))
  wend
  #1 "color blue"
  x=1.2: gosub [pot]: y=y+L2/x/x
  #1 "place 432 "; int(180*(1-y/cf))
  while y<cf
    x=x-.002: gosub [pot]: y=y+L2/x/x
    #1 "goto "; int(180*(1.2+x)); " "; int(180*(1-y/cf))
  wend
  #1 "color red"
  x=1: gosub [pot]: en=L2+y: y1=int(180*(1-en/cf))
  #1 "line 216 "; y1; " 432 "; y1
  e=.00000001: a=.1: b=.5
  x=a: gosub [pot]: ya=L2/a/a+y: x=b: gosub [pot]: yb=L2/b/b+y
  while (b-a)>e
    x=(a+b)/2: gosub [pot]: yx=L2/x/x+y
    if yx>en then a=x else b=x
  wend
  x1=int(180*(1.2+x)): #1 "line "; x1; " 180 "; x1; " "; y1
  #1 "line 396 180 396 "; y1
  c=x: theta=0: z1=0: dr=(1-c)/np
  for i=1 to np-1
    x=1-i*dr: gosub [pot]: z2=x*x/L*(2*(en-L2/x/x-y))^.5
    dtheta(i)=2*dr/(z1+z2): theta=theta+dtheta(i): z1=z2
  next i
  dtheta(np)=2*dr/z1: theta=theta+dtheta(np)
  #1 "place 0 60": #1 "\Apsidal Angle"
  #1 "\"; left$(str$(theta*180/pi), 8); " deg"
  orbit=1
goto [waitHere]

[pot]
  if n=0-2 then y=0-1/x
  if n=0-1 then y=log(x)
  if n=1 then y=x*x/2
return

[quit]
  close #1
end