PROGRAM planet                    ! planetary motion
LIBRARY "csgraphics"
DIM pos(2),vel(2)                 ! define dimension of arrays
CALL initial(pos(),vel(),t,GM,dt,nshow)
CALL output(pos(),t)
LET counter = 0
DO
   CALL Euler(pos(),vel(),t,GM,dt)
   LET counter = counter + 1
   IF mod(counter,nshow) = 0 then
      CALL output(pos(),t)        ! draw position of "earth"
   END IF
LOOP until key input
END

SUB initial(pos(),vel(),t,GM,dt,nshow)
    LET GM = 4.0*pi^2             ! GM constant in astronomical units
    LET t = 0
    INPUT prompt "time step = ": dt
    INPUT prompt "number of time steps between plots = ": nshow
    INPUT prompt "initial x position = ": pos(1)
    LET pos(2) = 0                ! initial y-position
    LET vel(1) = 0                ! initial x-velocity
    INPUT prompt "initial y-velocity = ": vel(2)
    ! assumed maximum value of semi-major axis
    LET r = 2*pos(1)
    CALL compute_aspect_ratio(r,xwin,ywin)  ! in library
    SET WINDOW -xwin,xwin,-ywin,ywin
END SUB

SUB Euler(pos(),vel(),t,GM,dt)    ! Euler-Cromer algorithm
    DIM accel(2)
    LET r2 = pos(1)*pos(1) + pos(2)*pos(2)
    LET r3 = r2*sqr(r2)
    FOR i = 1 to 2
        LET accel(i) = -GM*pos(i)/r3
        LET vel(i) = vel(i) + accel(i)*dt
        LET pos(i) = pos(i) + vel(i)*dt
    NEXT i
    LET t = t + dt
END SUB

SUB output(pos(),t)
    IF t = 0 then
       CLEAR
       LET radius = 0.05          ! radius of "sun" (not to scale)
       ASK MAX COLOR mc
       IF mc > 2 then SET COLOR "red"
       BOX CIRCLE -radius,radius,-radius,radius  ! "sun" at origin
       FLOOD 0,0                  ! paint sun in foreground color
       IF mc > 2 then SET COLOR "green"
    END IF
    PLOT pos(1),pos(2)
END SUB