' h6alib.bas - jackord@kw.igs.net - revised 24 Nov 2010 - Liberty Basic v4.03
' Central Force Motion: apsidal angle integral

' Initialize Window
  nomainwin
  button#1, "+1", [b1], UL, 4, 4, 20, 20
  button#1, "0", [b2], UL, 28, 4, 20, 20
  button#1, "-1", [b3], UL, 52, 4, 20, 20
  button#1, "-2", [b4], UL, 76, 4, 20, 20
  button#1, "-3", [b5], UL, 100, 4, 20, 20
  button#1, "-5", [b6], UL, 124, 4, 20, 20
  button#1, "Orbit", [b7], UL, 148, 4, 40, 20
  WindowWidth=450     ' pixel scale 0-432
  WindowHeight=438    ' pixel scale 0-400
  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"
  xo=216: yo=200: ro=180
  np=10000: cf=5: pi=4*atn(1)
  #1 "line "; xo; " "; yo; " "; 2*xo; " "; yo: #1 "line "; xo; " 0 "; xo; " "; 2*yo
  dim dtheta(10001)

[waitHere]
  wait

[b1]
  kk=1: n=1: norbit=2: L=.488: goto [apint]

[b2]
  kk=1: n=0: norbit=16: L=.488: goto [apint]

[b3]
  kk=1: n=-1: norbit=16: L=.488: goto [apint]

[b4]
  kk=1: n=-2: norbit=1: L=.488:goto [apint]

[b5]
  kk=1: n=-3: norbit=1: L=.999: goto [apint]

[b6]
  kk=1: n=-5: norbit=1: L=1/sqr(2): goto [apint]

[b7]
  if orbit=1 then kk=2
  goto [apint]

[apint]
  #1 "cls": #1 "color black"
  if kk=1 then
    #1 "line "; xo; " "; yo; " "; 2*xo; " "; yo
    #1 "line "; xo; " 0 "; xo; " "; 2*yo
    #1 "color green"
    L2=L*L/2: x=1.2: y=L2/x/x: x1=2*xo: y1=yo-ro*y/cf
    do
      x=x-.002: y=L2/x/x: x2=xo+ro*x: y2=yo-ro*y/cf
      #1 "line "; x1; " "; y1; " "; x2; " "; y2: x1=x2: y1=y2
    loop while y2>0
    #1 "color darkpink"
    x=1.2: gosub [pot]: x1=2*xo: y1=yo-ro*y/cf
    do
      x=x-.002: gosub [pot]: x2=xo+ro*x: y2=yo-ro*y/cf
      #1 "line "; x1; " "; y1; " "; x2; " "; y2: x1=x2: y1=y2
    loop while y2<2*yo and x2>xo
    #1 "color blue"
    x=1.2: gosub [pot]: y=y+L2/x/x: x1=2*xo: y1=yo-ro*y/cf
    do
      x=x-.002: gosub [pot]: y=y+L2/x/x: x2=xo+ro*x: y2=yo-ro*y/cf
      #1 "line "; x1; " "; y1; " "; x2; " "; y2: x1=x2: y1=y2
    loop while y2>0 and y2<2*yo
    #1 "color red"
    x=1: gosub [pot]: en=L2+y: y1=yo-ro*en/cf
    #1 "line "; xo; " "; y1; " "; 2*xo; " "; y1
    if n>-2.5 then
      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
      c=x
    else
      c=.1
    end if
    x1=xo+ro*c: #1 "line "; x1; " "; yo; " "; x1; " "; y1
    #1 "line "; xo+ro; " "; yo; " "; xo+ro; " "; y1
    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 10 50": #1 "\Force Law Power = "; n: #1 "place 10 70"
    if n>-2.5 then #1 "\Apsidal Angle = "; left$(str$(theta*180/pi), 8); " deg"
    if n=-3 then #1 "\Spirals Toward Force Centre"
    if n=-5 then #1 "\Circle Through Force Centre"
    orbit=1
  end if
  if kk=2 then
    #1 "backcolor green"
    #1 "place "; xo; " "; yo: #1 "circlefilled 18"
    #1 "place "; xo; " "; yo: #1 "circle 18"
    #1 "backcolor white"
    #1 "color red"
    #1 "place "; xo; " "; yo: #1 "circle "; ro
    if n>-2.5 then #1 "place "; xo; " "; yo: #1 "circle "; ro*c
    #1 "color blue"
    r=1: theta=0: #1 "place "; xo+ro; " "; yo
    for j=1 to norbit
      for i=1 to np
        r=r-dr: theta=theta+dtheta(i)
        #1 "goto "; xo+ro*r*cos(theta); " "; yo-ro*r*sin(theta)
      next i
      if n>-2.5 then
        for i=1 to np
          r=r+dr: theta=theta+dtheta(np+1-i)
          #1 "goto "; xo+ro*r*cos(theta); " "; yo-ro*r*sin(theta)
        next i
      end if
      if n=-5 then
        r=1: theta=0: #1 "place "; xo+ro; " "; yo
        for i=1 to np
          r=r-dr: theta=theta-dtheta(i)
          #1 "goto "; xo+ro*r*cos(theta); " "; yo-ro*r*sin(theta)
        next i
      end if
    next j
    orbit=0
  end if
  kk=0
goto [waitHere]

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

[quit]
  close #1
end