PROGRAM planet2
! d = 2 solar system with major and minor planet
DIM x(2),y(2),vx(2),vy(2),ratio(2)
LIBRARY "csgraphics"
CALL initial(x(),y(),vx(),vy(),t,GM,ratio(),dt,nshow)
CALL output(x(),y(),t)
LET counter = 0
DO
   CALL Euler(x(),y(),vx(),vy(),t,GM,ratio(),dt)
   LET counter = counter + 1
   IF mod(counter,nshow) = 0 then
      CALL output(x(),y(),t)
   END IF
LOOP until key input
END

SUB initial(x(),y(),vx(),vy(),t,GM,ratio(),dt,nshow)
    LET t = 0
    LET GM = 4.0*pi*pi            ! gravitational constant
    LET dt = 0.001                ! time step (yrs)
    LET nshow = 20                ! # of iterations between output
    ! planet one initial conditions
    LET x(1) = 2.52
    LET y(1) = 0
    LET vx(1) = 0
    LET vy(1) = sqr(GM/x(1))
    ! mass of planet 2 is 0.04 mass of sun 
    LET ratio(1) = 0.04*GM
    ! planet two initial conditions
    LET x(2) = 5.24
    LET y(2) = 0
    LET vx(2) = 0
    LET vy(2) = sqr(GM/x(2))
    ! mass of planet 1 is 0.001 mass of sun 
    LET ratio(2) = -0.001*GM
    ! set up window   
    CALL compute_aspect_ratio(2*x(2),xwin,ywin)
    SET WINDOW -xwin,xwin,-ywin,ywin
END SUB

SUB Euler(x(),y(),vx(),vy(),t,GM,ratio(),dt)
    ! compute distance between planets 1 and 2
    LET dx = x(2) - x(1)
    LET dy = y(2) - y(1)
    LET dr2 = dx*dx + dy*dy
    LET dr = sqr(dr2)
    LET dr3 = dr2*dr
    FOR i = 1 to 2                ! sum over planets
        LET r2 = x(i)*x(i) + y(i)*y(i)
        LET r = sqr(r2)
        LET r3 = r2*r
        LET ax = -GM*x(i)/r3 + ratio(i)*dx/dr3
        LET ay = -GM*y(i)/r3 + ratio(i)*dy/dr3
        LET vx(i) = vx(i) + ax*dt
        LET vy(i) = vy(i) + ay*dt
        LET x(i) = x(i) + vx(i)*dt
        LET y(i) = y(i) + vy(i)*dt
    NEXT i
    LET t = t + dt
END SUB

SUB output(x(),y(),t)             ! plot orbits
    IF t = 0 then
       SET COLOR "red"
       BOX CIRCLE -0.1,0.1,-0.1,0.1    ! sun at origin
       FLOOD 0,0
    END IF
    SET COLOR "blue"
    PLOT POINTS: x(1),y(1)        ! planet one
    SET COLOR "green"
    PLOT POINTS: x(2),y(2)        ! planet two
END SUB