      PROGRAM planet2
*     大小2つの惑星のある2次元太陽系
      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                      ! 重力定数
      dt = 0.001D0                          ! 時間幅 (年)
      nshow = 20                            ! 各出力間の計算回数
*     惑星1の初期条件
      x(1) = 2.52D0
      y(1) = 0
      vx(1) = 0
      vy(1) = sqrt(GM/x(1))
*     惑星2の質量は太陽の 0.04 倍
      ratio(1) = 0.04D0*GM
*     惑星2の初期条件
      x(2) = 5.24D0
      y(2) = 0
      vx(2) = 0
      vy(2) = sqrt(GM/x(2))
*     惑星1の質量は太陽の 0.001 倍
      ratio(2) = -0.001D0*GM
*     ウィンドウの設定
      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)
*     惑星1と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)  ! 太陽を原点とする
         CALL GWputmrk(IR_, 0.0, 0.0)
      END IF
      CALL GWsetpxl(IR_, REAL(x(1)), REAL(y(1)), 16) ! 惑星1
      CALL GWsetpxl(IR_, REAL(x(2)), REAL(y(2)), 14) ! 惑星2
      END
