REM ftjnbbc.bas - jackord@kw.igs.net - 26 Feb 10 - BBC Basic 5.50b
      REM Restricted 3-body problem: the Trojan asteroid group
      REM   (1) Locating equilibrium sites in the rotating frame using a colour
      REM       plot of the effective potential
      REM   (2) Testing the stability of orbits at the Trojan and saddlepoint sites
      
      *FLOAT 64
      
      REM Install libraries
      INSTALL @lib$+"WINLIB5"
      
      REM Initialize window
      WindowWidth=400                       : REM Pixel Scale 0-400
      WindowHeight=300                      : REM Pixel Scale 0-300
      VDU 23,22,WindowWidth;WindowHeight;8,15,16,128
      *FONT Arial,10
      b1% = FN_button("Tro", 5, 275, 40, 20, FN_setproc(PROCq1), 0)
      b2% = FN_button("Sad", 50, 275, 40, 20, FN_setproc(PROCq2), 0)
      c1% = FN_combobox("", 5, 5, 80, 20, 0, 3)
      SYS "SetWindowText", @hwnd%, "Trojan Asteroid Group"
      SYS "SendMessage", c1%, &143, 0, "M/m=9"
      SYS "SendMessage", c1%, &143, 1, "M/m=19"
      SYS "SendMessage", c1%, &143, 1, "M/m=29"
      SYS "SendMessage", c1%, &143, 1, "M/m=99"
      SYS "SendMessage", c1%, &143, 1, "Jupiter"
      SYS "SendMessage", c1%, &143, 1, "Earth"
      SYS "SendMessage", c1%, &14E, 0, 0
      OFF: ON CLOSE QUIT
      DIM m(6): DIM cf(6): DIM aa(3): DIM bb(3): DIM c%(6)
      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
      c%(0)=12: c%(1)=14: c%(2)=9
      c%(3)=11: c%(4)=10: c%(5)=13
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      
      DEF PROCq1
      opt=1: PROCcplot
      ENDPROC
      
      DEF PROCq2
      opt=2: PROCcplot
      ENDPROC
      
      DEF PROCcplot
      SYS "SendMessage", c1%, &147, 0, 0 TO ii
      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
      CLS : VDU 5: GCOL 0
      CIRCLE FILL 2*INT(x1+.5), 2*y0, 28
      CIRCLE FILL 2*INT(x2+.5), 2*y0, 20
      FOR i=1 TO 399
        FOR j=0 TO 149
          rr=(i-x0)*(i-x0)+j*j
          r1=SQR((i-x1)*(i-x1)+j*j)
          r2=SQR((i-x2)*(i-x2)+j*j)
          IF r1<=13 OR r2<=9 THEN
            GCOL 0
          ELSE
            u=w2*rr/2+m1/r1+m2/r2
            plt=INT(cf(ii)/u): plt=plt-6*INT(plt/6)
            GCOL c%(plt)
          ENDIF
          PLOT 2*i, 2*(y0+j): PLOT 2*i, 2*(y0-j)
        ENDIF
      NEXT j
      NEXT i
      IF opt=1 THEN
      dt=5: x=x1+d/2: y=y0+d*SQR(3)/2
      imax=INT(32*PI/w/dt): IF ii=5 THEN imax=6*imax
      GCOL 15: MOVE 2*INT(x1), 2*y0: DRAW 2*INT(x2), 2*y0
      DRAW 2*INT(x), 2*INT(y): DRAW 2*INT(x1), 2*y0
      y=y+d/400: PROCorbit
      ELSE
      GCOL 15: LINE 0, 300, 800, 300
      dt=1: imax=INT(16*PI/w/dt): xm=400: jm=2: ya=0
      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 jm=0: imax=2*imax
      IF ii>3 THEN
        xm=180: jm=0: dt=5: imax=5*imax
        IF ii=5 THEN dt=50
      ENDIF
      FOR i=0 TO xm
        IF ABS(i-x1)<4 THEN i=i+8: ya=0
        IF ABS(i-x2)<2 THEN i=i+4: ya=0
        r13=ABS((i-x1)^3): r23=ABS((i-x2)^3)
        yb=INT(y0-200000*(w2*(x0-i)+m1*(i-x1)/r13+m2*(i-x2)/r23))
        IF i>0 AND ya>0 AND yb<300 THEN
          LINE 2*(i-1), 2*ya, 2*i, 2*yb
        ENDIF
        ya=yb
      NEXT i
      FOR j=0 TO jm
        e=aa(j): f=bb(j)
        WHILE (f-e)>.001
          h=(e+f)/2: PROCforce
          IF ff<0 THEN e=h ELSE f=h
        ENDWHILE
        x=(e+f)/2: y=y0+d/400
        PROCorbit
      NEXT j
      ENDIF
      ENDPROC
      
      DEF PROCforce
      r13=ABS(h-x1)^3: r23=ABS(h-x2)^3
      ff=w2*(h-x0)+m1*(x1-h)/r13+m2*(x2-h)/r23
      ENDPROC
      
      DEF PROCorbit
      vx=0: vy=0: dvx=0: dvy=0: PROCnewton
      GCOL 0: MOVE 2*INT(x), 2*INT(y)
      vx=vx+dvx/2: vy=vy+dvy/2
      FOR i=1 TO imax
      x=x+vx*dt: y=y+vy*dt
      DRAW 2*INT(x), 2*INT(y)
      PROCnewton
      vx=vx+dvx: vy=vy+dvy
      NEXT i
      ENDPROC
      
      DEF PROCnewton
      r13=((x-x1)*(x-x1)+(y-y0)*(y-y0))^1.5
      r23=((x-x2)*(x-x2)+(y-y0)*(y-y0))^1.5: dvt=dvx
      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+dvt/2)-m1*(y-y0)/r13-m2*(y-y0)/r23)*dt
      IF r13<2197 OR r23<729 THEN i=imax
      ENDPROC