' ftjnlib.bas - jackord@kw.igs.net - revised 7 Dec 04 - Liberty Basic v4.0
' restricted 3-body problem:
'   (1) locating equilibrium sites in the rotating frame using a color
'       plot of effective potential
'   (2) testing the stability of orbits at the Trojan and saddle point sites

  nomainwin
  WindowWidth=410     ' pixel scale 0-400
  WindowHeight=329    ' pixel scale 0-300
  UpperLeftX=100: UpperLeftY=50
  dim a$(6): a$(0)="M/m=9": a$(1)="M/m=19": a$(2)="M/m=29": a$(3)="M/m=99"
  a$(4)="Jupiter": a$(5)="Earth"
  dim c$(6): c$(0)="blue": c$(1)="cyan": c$(2)="red"      ' Colors
  c$(3)="yellow": c$(4)="green": c$(5)="pink"
  dim m(6), cf(6), aa(3), bb(3)
  m(0)=.1: m(1)=.05: m(2)=1/30: m(3)=.01: m(4)=.001: m(5)=.000003
  cf(0)=.209: cf(1)=.2116: cf(2)=.2123: cf(3)=.2137: cf(4)=.21427: cf(5)=.21429
  button#1, "Plot", [plot], UL, 3, 277, 40, 20
  combobox#1.c1, a$(, [waitHere], 3, 3, 80, 150
  radiobutton#1.rb1, "Trojan", [trojan], [waitHere], 50, 277, 68, 20
  radiobutton#1.rb2, "Saddle", [saddle], [waitHere], 120, 277, 68, 20
  open "Trojan Asteroid Group" for graphics_nsb as #1
  #1.c1 "selectindex 1"
  #1.rb1 "Set": opt=1
  #1 "trapclose [quit]": #1 "down"

[waitHere]
  wait

[trojan]
  #1.rb1 "Set": opt=1: goto [waitHere]

[saddle]
  #1.rb2 "Set": opt=2: goto [waitHere]

[plot]
  #1.c1 "selectionindex?": input#1.c1, ii: ii=ii-1        ' Set M/m ratio
  m2=m(ii): m1=1-m2
  x0=200: y0=150: d=140: w2=1/(d*d*d): w=sqr(w2)
  x1=x0-m2*d: x2=x0+m1*d                                  ' Site M and m
  #1 "cls ; color darkgray ; backcolor darkgray"          ' Initial screen
  #1 "place "; int(x1+.5); " "; y0: #1 "circlefilled 13"
  #1 "place "; int(x2+.5); " "; y0: #1 "circlefilled 9 ; backcolor white"
  for i=1 to 399                                          ' Scan the screen
    #1 "up"
    for j=1 to 299
      rr=(i-x0)*(i-x0)+(j-y0)*(j-y0)
      r1=sqr((i-x1)*(i-x1)+(j-y0)*(j-y0))
      r2=sqr((i-x2)*(i-x2)+(j-y0)*(j-y0))
      if r1>13 and r2>9 then
        u=w2*rr/2+m1/r1+m2/r2                             ' Effective potential
        plt=int(cf(ii)/u): plt=plt-6*int(plt/6)           ' Color pixel
        #1 "color "; c$(plt); " ; goto "; i; " "; j; " ; down"
      else
        #1 "up"
      end if
    next j
    scan
  next i
  if opt=1 then                                           ' Trojan site
    dt=5: x=x1+d/2: y=y0-d*sqr(3)/2
    imax=int(32*4*atn(1)/w/dt): if ii=5 then imax=10*imax
    #1 "color white ; place "; int(x1); " "; y0           ' Equilateral triangle
    #1 "goto "; int(x2); " "; y0
    #1 "goto "; int(x); " "; int(y)
    #1 "goto "; int(x1); " "; y0
    y=y-d/400: gosub [orbit]                              ' Plot orbit
  else                                                    ' Saddle points
    #1 "color white ; line 0 150 400 150"
    dt=1: imax=int(8*4*atn(1)/w/dt): xm=400: jm=2
    aa(0)=0: bb(0)=x1-13: aa(1)=x0: bb(1)=x2-9: aa(2)=x2+9: bb(2)=400
    if ii>3 then
      xm=180: jm=0: dt=5: imax=10*imax
      if ii=5 then dt=50
    end if
    for i=0 to xm                                         ' Plot the force
      if abs(i-x1)<13 then i=i+26
      if abs(i-x2)<9 then i=i+18
      r13=abs((i-x1)^3): r23=abs((i-x2)^3)
      yb=int(150.5+200000*(w2*(x0-i)+m1*(i-x1)/r13+m2*(i-x2)/r23))
      if i>0 and ya>0 and yb<300 then
        #1 "line "; i-1; " "; ya; " "; i; " "; yb
      end if
      ya=yb
    next i
    for j=0 to jm                                         ' Three saddle points
      e=aa(j): f=bb(j)                                    '   if ii<4
      while (f-e)>.001                                    ' Bisection algorithm
        h=(e+f)/2: gosub [force]
        if ff<0 then e=h else f=h
      wend
      x=(e+f)/2: y=y0-d/400
      gosub [orbit]                                       ' Plot the orbits
    next j
  end if
goto [waitHere]

[force]                                                   ' Routine used in force plot
  r13=abs(h-x1)^3: r23=abs(h-x2)^3
  ff=w2*(h-x0)+m1*(x1-h)/r13+m2*(x2-h)/r23
return

[orbit]                                                   ' Routine that plots orbits
  vx=0: vy=0: gosub [newton]
  #1 "color black ; place "; int(x); " "; int(y)
  vx=vx+dvx/2: vy=vy+dvy/2
  for i=1 to imax
    x=x+vx*dt: y=y+vy*dt
    #1 "goto "; int(x); " "; int(y)
    gosub [newton]
    vx=vx+dvx: vy=vy+dvy
  next i
return

[newton]                                                  ' Newton's law
  r13=((x-x1)*(x-x1)+(y-y0)*(y-y0))^1.5
  r23=((x-x2)*(x-x2)+(y-y0)*(y-y0))^1.5
  dvx=(w2*(x-x0)-2*w*(vy+dvy/2)-m1*(x-x1)/r13-m2*(x-x2)/r23)*dt
  dvy=(w2*(y-y0)+2*w*(vx+dvx/2)-m1*(y-y0)/r13-m2*(y-y0)/r23)*dt
  if r13<2197 or r23<729 then i=imax
return

[quit]
  close#1
end