      PROGRAM planet2
*     d = 2 solar system with major and minor planet
      IMPLICIT NONE
      REAL*8 GM, dt, ratio(2), t, vx(2), vy(2), x(2), y(2)
      INTEGER IR_, counter, nshow
      LOGICAL Ltmp1_, TBkeyinput
      CALL GWopen(IR_, 0)
      CALL initial(x,y,vx,vy,t,GM,ratio,dt,nshow)
      CALL output(x,y,t)
      counter = 0
      Ltmp1_ = .TRUE.
      DO WHILE(Ltmp1_)
         CALL Euler(x,y,vx,vy,t,GM,ratio,dt)
         counter = counter + 1
         IF(mod(counter,nshow) .EQ. 0) THEN
            CALL output(x,y,t)
         END IF
         Ltmp1_ = .NOT.TBkeyinput()
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(x,y,vx,vy,t,GM,ratio,dt,nshow)
      IMPLICIT NONE
      REAL*8 pi, x(2), y(2), vx(2), vy(2), t, GM, ratio(2), dt
      REAL xwin, ywin
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER IR_, nshow
      t = 0
      GM = 4.0D0*pi*pi                      ! gravitational constant
      dt = 0.001D0                          ! time step (yrs)
      nshow = 20                            ! # of iterations between output
*     planet one initial conditions
      x(1) = 2.52D0
      y(1) = 0
      vx(1) = 0
      vy(1) = sqrt(GM/x(1))
*     mass of planet 2 is 0.04 mass of sun 
      ratio(1) = 0.04D0*GM
*     planet two initial conditions
      x(2) = 5.24D0
      y(2) = 0
      vx(2) = 0
      vy(2) = sqrt(GM/x(2))
*     mass of planet 1 is 0.001 mass of sun 
      ratio(2) = -0.001D0*GM
*     set up window   
      CALL compute_aspect_ratio(REAL(2*x(2)),xwin,ywin)
      CALL GWindow(IR_, -xwin, -ywin, xwin, ywin)
      END
C
      SUBROUTINE Euler(x,y,vx,vy,t,GM,ratio,dt)
*     compute distance between planets 1 and 2
      IMPLICIT NONE
      REAL*8 ax, ay, dr, dr2, dr3, dx, dy, r, r2, r3, x(2), y(2),
     +       vx(2), vy(2), t, GM, ratio(2), dt
      INTEGER i
      dx = x(2) - x(1)
      dy = y(2) - y(1)
      dr2 = dx*dx + dy*dy
      dr = sqrt(dr2)
      dr3 = dr2*dr
      DO i = 1, 2
         r2 = x(i)*x(i) + y(i)*y(i)
         r = sqrt(r2)
         r3 = r2*r
         ax = -GM*x(i)/r3 + ratio(i)*dx/dr3
         ay = -GM*y(i)/r3 + ratio(i)*dy/dr3
         vx(i) = vx(i) + ax*dt
         vy(i) = vy(i) + ay*dt
         x(i) = x(i) + vx(i)*dt
         y(i) = y(i) + vy(i)*dt
      END DO
      t = t + dt
      END
C
      SUBROUTINE output(x,y,t)
      IMPLICIT NONE
      REAL*8 x(2), y(2), t
      INTEGER IR_
      IF(t .EQ. 0) THEN
         CALL GWsetpen(IR_, 13, -1, -1, -1)
         CALL GWsetmrk(IR_, 6, 0.2, 13, -1, 4)  ! sun at origin
         CALL GWputmrk(IR_, 0.0, 0.0)
      END IF
      CALL GWsetpxl(IR_, REAL(x(1)), REAL(y(1)), 16) ! planet one
      CALL GWsetpxl(IR_, REAL(x(2)), REAL(y(2)), 14) ! planet two
      END
