REM fplwbbc.bas - jackord@kw.igs.net - revised 26 Apr 07 - BBC Basic 5.50b
      REM central-force motion in a power-law force field
      
      *FLOAT 64
      INSTALL @lib$+ "WINLIB5"
      
      ba% = FN_button("-2", 6, 4, 26, 18, FN_setproc(PROCGravity), 0)
      bb% = FN_button("-1", 42, 4, 26, 18, FN_setproc(PROCInverseR), 0)
      bc% = FN_button("+1", 78, 4, 26, 18, FN_setproc(PROCSpring), 0)
      WindowWidth=400: WindowHeight=400
      VDU 23, 22, WindowWidth; WindowHeight; 8, 16, 16, 128
      SYS "SetWindowText", @hwnd%, "Central-Force Motion"
      VDU 5: OFF
      PROCinit
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      QUIT
      
      DEF PROCGravity
      dt=.001: test=1: lp=1
      CLS: PROCinit: MOVE -400, -214: PRINT "Gravity: n=-2"
      t=0: x=1: y=0: vx=0: vy=.488
      GCOL 12: MOVE 360, 0: tt=TIME
      rm=(x*x+y*y)^1.5: vx=vx-x/rm*dt/2: vy=vy-y/rm*dt/2
      WHILE lp=1
        xo=x: yo=y: t=t+dt: x=x+vx*dt: y=y+vy*dt
        IF test=1 AND y<=0 AND yo>0 THEN rmin=-x: test=0
        DRAW INT(360*x), INT(360*y)
        rm=(x*x+y*y)^1.5: vx=vx-x/rm*dt: vy=vy-y/rm*dt
        IF y>=0 AND yo<0 THEN lp=0
        WHILE TIME-tt<300*t: WAIT 1: ENDWHILE
      ENDWHILE
      t=t-y/vy
      GCOL 0: MOVE -400, -244: PRINT "Closure Test: Rmax="; x
      MOVE -400, -274: PRINT "Rmin="; rmin
      MOVE -400, -304: PRINT "Apsidal Angle=180 degrees"
      MOVE -400, -334: PRINT "Torbit/Tcircle= "; t/2/PI
      MOVE -400, -364: PRINT "By Kepler's Law "; ((x+rmin)/2)^1.5;
      ENDPROC
      
      DEF PROCSpring
      dt=.001: test=1: lp=1
      CLS: PROCinit: MOVE -400, -244: PRINT "Spring: n=1"
      t=0: x=1: y=0: vx=0: vy=.488
      GCOL 12: MOVE 360, 0: tt=TIME
      vx=vx-x*dt/2: vy=vy-y*dt/2
      WHILE lp=1
        xo=x: yo=y: t=t+dt: x=x+vx*dt: y=y+vy*dt
        IF test=1 AND x<=0 AND xo>0 THEN rmin=y: test=0
        DRAW INT(360*x), INT(360*y)
        vx=vx-x*dt: vy=vy-y*dt
        IF y>=0 AND yo<0 THEN lp=0
        WHILE TIME-tt<130*t: WAIT 1: ENDWHILE
      ENDWHILE
      t=t-y/vy
      GCOL 0: MOVE -400, -274: PRINT "Closure Test: Rmax="; x
      MOVE -400, -304: PRINT "Rmin="; rmin
      MOVE -400, -334: PRINT "Apsidal Angle=90 degrees"
      MOVE -400, -364: PRINT "Torbit/Tcircle="; t/2/PI;
      ENDPROC
      
      DEF PROCInverseR
      dt=.001: test=1
      CLS: PROCinit: MOVE -400, -334: PRINT "InverseR: n=-1": TIME=0
      t=0: x=1: y=0: vx=0: vy=.488
      GCOL 12: MOVE 360, 0: tt=TIME
      rm=x*x+y*y: vx=vx-x/rm*dt/2: vy=vy-y/rm*dt/2
      WHILE t<47.5
        IF test=1 THEN
          xo=x: yo=y: rmin=rm
        ENDIF
        t=t+dt: x=x+vx*dt: y=y+vy*dt
        DRAW INT(360*x), INT(360*y)
        rm=x*x+y*y: vx=vx-x/rm*dt: vy=vy-y/rm*dt
        IF rm>rmin THEN test=0
        WHILE TIME-tt<50*t: WAIT 1: ENDWHILE
      ENDWHILE
      rmin=INT(10000*SQR(rmin)+.5)/10000
      ang=INT(18000*(1+ATN(yo/xo)/PI))/100
      GCOL 0: MOVE -400, -364
      PRINT "Rmin="; rmin; "  Apsidal Angle="; ang; "=360*11/32 degrees";
      ENDPROC
      
      DEF PROCinit
      ORIGIN 400, 400
      GCOL 0: MOVE -386, 350: PRINT "Select Power";
      GCOL 2: CIRCLE FILL 2, -2, 38
      GCOL 1: CIRCLE FILL 362, -2, 14
      GCOL 0
      LINE 354, 48, 360, 60: LINE 366, 48, 360, 60: DRAW 360, 0
      ENDPROC