REM fplwbbc.bas - jackord@kw.igs.net - revised 12 Nov 2010 - BBC Basic 5.50b
      REM central-force motion in a power-law force field for n=1, 0, -1,-2, -3, and -5
      
      *FLOAT 64
      INSTALL @lib$+ "WINLIB5"
      
      b1% = FN_button("Plot", 70, 20, 50, 25, FN_setproc(PROCq1), 0)
      c1% = FN_combobox("", 5, 20, 60, 0, 0, 3)
      WindowWidth=400: WindowHeight=400
      VDU 23, 22, WindowWidth; WindowHeight; 8, 16, 16, 128
      SYS "SetWindowText", @hwnd%, "Central-Force Motion"
      SYS "SendMessage", c1%, &143, 0, "n=1"
      SYS "SendMessage", c1%, &143, 1, "n=0"
      SYS "SendMessage", c1%, &143, 1, "n=-1"
      SYS "SendMessage", c1%, &143, 1, "n=-2"
      SYS "SendMessage", c1%, &143, 1, "n=-3"
      SYS "SendMessage", c1%, &143, 1, "n=-5"
      SYS "SendMessage", c1%, &14E, 0, 0
      VDU 5: OFF
      PROCinit
      
      REPEAT
        WAIT 1
      UNTIL FALSE
      QUIT
      
      DEF PROCq1
      SYS "SendMessage", c1%, &147, 0, 0 TO n: n=1-n
      IF n=-4 THEN n=-5
      dt=.0004: test=1: rmin=0: xo=0: yo=0
      CLS: PROCinit
      t=0: x=1: y=0: vx=0: vy=.488
      IF n=-3 THEN vy=.999
      IF n=-5 THEN vy=1/SQR(2)
      GCOL 12: MOVE 360, 0
      r=SQR(x*x+y*y): vx=vx-dt/2
      WHILE t<47.5 AND 360*r>34
        IF test=1 THEN xo=x: yo=y: rmin=r
        t=t+dt: x=x+vx*dt: y=y+vy*dt
        DRAW INT(360*x), INT(360*y)
        r=SQR(x*x+y*y): rm=r^(n-1)
        vx=vx-x*rm*dt: vy=vy-y*rm*dt
        IF r>rmin THEN test=0
      ENDWHILE
      IF n>-2.5 THEN
        GCOL13: CIRCLE FILL INT(360*x)+2, INT(360*y)-2, 14
        GCOL 0
        MOVE -400, -364: PRINT " Apsidal Angle="; 180*(1+ATN(yo/xo)/PI)
        MOVE 100, -364: PRINT " Rmin="; rmin
      ENDIF
      IF n=-3 THEN
        GCOL 0: MOVE -400, -364: PRINT " Vinit=0.999"
      ENDIF
      IF n=-5 THEN
        MOVE 360,0
        x=1: y=0: vx=-dt/2: vy=-1/SQR(2): r=1
        WHILE 360*r>34
          x=x+vx*dt: y=y+vy*dt
          DRAW INT(360*x), INT(360*y)
          r=SQR(x*x+y*y): rm=r^(n-1)
          vx=vx-x*rm*dt: vy=vy-y*rm*dt
        ENDWHILE
        GCOL 0: MOVE -400, -364: PRINT " Vinit=1/sqrt(2)"
      ENDIF
      ENDPROC
      
      DEF PROCinit
      ORIGIN 400, 400
      GCOL 0: MOVE -390, 390: PRINT "Central Force F=-k*r^n"
      GCOL 10: CIRCLE FILL 2, -2, 38
      GCOL 9: CIRCLE FILL 362, -2, 14
      GCOL 0
      LINE 354, 48, 360, 60: LINE 366, 48, 360, 60: DRAW 360, 0
      ENDPROC