C
      IMPLICIT NONE
      CHARACTER buff*60
      REAL th(2),om(2),t,dt,dt2,energy
      REAL a,b,e0,m1,m2,r1,r2,s1,s2
      REAL x1,y1,x2,y2
      INTEGER IR,n,nl,MS,IT,IDT,IVK,IS           ! modified
      COMMON /PRMS/a,b,e0,m1,m2,r1,r2,s1,s2,MS
      DATA IT,IVK,IS/3*0/                        ! modified
C
      CALL initial(th,om,t,dt,dt2,nl)
      IDT = dt*nl * 1000
C
      CALL GWSETMSG(IR,'click to quit')
      CALL GWEVENT(IR, x1,x2)
      DO WHILE(IVK .EQ. 0)
        CALL GWSETOGN(IR, 0, -2)
        x1 =  r1*sin(th(1))
        y1 = -r1*cos(th(1))
        CALL GWLINE(IR,0.0,0.0,x1,y1)
        x2 = x1 + r2*sin(th(2))
        y2 = y1 - r2*cos(th(2))
        CALL GWLINE2(IR,x2,y2)
        CALL GWPUTCMB(IR,MS,x1,y1,s1,s1,0)
        CALL GWPUTCMB(IR,MS,x2,y2,s2,s2,0)
        CALL GWFLUSH(IR,-2)
        DO n = 1,nl
          CALL UPDATE(th,om,t,dt,dt2)
        ENDDO
        WRITE(buff,'(3(A,F7.2),A,1P,E8.1)') 
     +     't =',t,',  x1 =',th(1),',  x2 =',th(2),
     +     ',  dE =',(energy(t,th,om)-e0)/abs(e0)
        CALL GWSETOGN(IR, 0, -3)
        CALL GWPUTTXT(IR, -1.0, -1.07, buff)
        IT = IT + IDT
        CALL GWSLEEP2(IR,IT)
        CALL GWFLUSH(IR,-3)
        CALL GWEVENT(IVK, x1,x2)
        IS = IS + 1                              ! added
        WRITE(buff,'(A,I6,A)') 'dp',IS,'.emf'    ! added
        n = 3                                    ! added
        DO WHILE(buff(n:n) .EQ. ' ')             ! added
          buff(n:n) = '0'                        ! added
          n = n + 1                              ! added
        ENDDO                                    ! added
        CALL GWSAVEAS(IR,2,2,8000,8000,buff)     ! added
      ENDDO
      CALL GWQUIT(IR)
      END
C
      SUBROUTINE initial(th,om,t,dt,dt2,nl)
      IMPLICIT NONE
      CHARACTER buff*40
      REAL th(2),om(2),t,dt,dt2,energy
      REAL a,b,e0,m1,m2,r1,r2,s1,s2
      REAL x1,y1,x2,y2
      INTEGER IR,nl,lb,MS
      COMMON /PRMS/a,b,e0,m1,m2,r1,r2,s1,s2,MS
C
      CALL GWINIT(IR)
      CALL GWOPEN(IR,0)
C     CALL GWMODE(IR, 2, 1)
      CALL GWVPORT(IR,0.0,0.0,1.0,1.0)
      CALL GWINDOW(IR,-1.1,-1.1,1.1,1.1)
      CALL GWRECT(IR,-1.0,-1.0,1.0,1.0)
      CALL GWSETTXT(IR, 0.06, 0.0, -1, -1, -1, '')
C
      t = 0
      dt = 0.001
      dt2 = dt/2
      nl = 50
C
      m1 = 0
      m2 = 0
      DO WHILE(m1.LT.0.01d0 .OR. m2.LE.0.01d0)
        CALL GWINPUT(lb, 'm1, m2:', buff)
        IF(lb .GT. 2) THEN
          READ(buff,*) m1,m2
        ELSE
          m1 = 1
          m2 = 2
        ENDIF
      ENDDO
      WRITE(buff,'(2(A,F4.2))') 'm1 = ', m1, ',  m2 = ', m2
      CALL GWPUTTXT(IR, -1.0, 1.01, buff)
      CALL GWANCHOR(IR,1)
      CALL GWSETOGN(IR,0,1)
C
      s1 = 0.05*m1**(1.0/3.0)
      s2 = 0.05*m2**(1.0/3.0)
C
      CALL GWCMBMRK(MS,0,6,0.1,13,-1,4)
C
      CALL GWCAPVEC(IR,0.0,0.0,x1,y1,'')
      CALL GWLINE(IR,0.0,0.0,x1,y1)
      CALL GWPUTCMB(IR,MS,x1,y1,s1,s1,0)
      CALL GWCAPVEC(IR,x1,y1,x2,y2,'')
      CALL GWLINE2(IR,x2,y2)
      CALL GWPUTCMB(IR,MS,x2,y2,s2,s2,0)
      r1 = SQRT(x1**2+y1**2)
      r2 = SQRT((x2-x1)**2+(y2-y1)**2)
      th(1) = atan2(x1,-y1)
      th(2) = atan2(x2-x1,y1-y2)
      om(1) = 0
      om(2) = 0
c
c     g == r1
      a = (m1+m2)*r1/m2/r2
      b = r2/r1
      e0 = energy(t,th,om)
      WRITE(*,'(4(A,F6.2))') 
     +   't =',t,', x1 =',th(1),', x2 =',th(2),', E =',e0
c
      CALL GWERASE(IR,1,-1)
      END
C
      SUBROUTINE UPDATE(th,om,t,dt,dt2)
      IMPLICIT NONE
      REAL th(2),om(2),t,dt,dt2
      REAL x(2),v(2),acc(2)
      REAL k1v(2),k1x(2),k2x(2),k2v(2),k3x(2),k3v(2),k4x(2),k4v(2)
      INTEGER I
      CALL accel(t,th,om,acc)
      DO I = 1,2
        k1v(I) = acc(I)*dt
        k1x(I) = om(I)*dt
        x(I) = th(I)+k1x(I)/2
        v(I) = om(I)+k1v(I)/2
      ENDDO
      CALL accel(t+dt2,x,v,acc)
      DO I = 1,2
        k2v(I) = acc(I)*dt
        k2x(I) = (om(I)+k1v(I)/2)*dt
        x(I) = th(I)+k2x(I)/2
        v(I) = om(I)+k2v(I)/2
      ENDDO
      CALL accel(t+dt2,x,v,acc)
      DO I = 1,2
        k3v(I) = acc(I)*dt
        k3x(I) = (om(I)+k2v(I)/2)*dt
        x(I) = th(I)+k3x(I)/2
        v(I) = om(I)+k3v(I)/2
      ENDDO
      CALL accel(t+dt,x,v,acc)
      DO I = 1,2
        k4v(I) = acc(I)*dt
        k4x(I) = (om(I)+k3v(I))*dt
        th(I) = th(I)+(k1x(I)+2*k2x(I)+2*k3x(I)+k4x(I))/6
        om(I) = om(I)+(k1v(I)+2*k2v(I)+2*k3v(I)+k4v(I))/6
      ENDDO
      t = t + dt
      END
C
      SUBROUTINE accel(t,x,v,acc)
      IMPLICIT NONE
      REAL t,x(2),v(2),acc(2),den
      REAL c12,s12,f1,f2
      REAL a,b,e0,m1,m2,r1,r2,s1,s2
      INTEGER MS
      COMMON /PRMS/a,b,e0,m1,m2,r1,r2,s1,s2,MS
      c12 = cos(x(1)-x(2))
      s12 = sin(x(1)-x(2))
      f1 = -s12*v(2)**2-a*sin(x(1))
      f2 =  s12*v(1)**2-  sin(x(2))
      den = a*b - c12*c12
      acc(1) = (b*f1-c12*f2)/den
      acc(2) = (a*f2-c12*f1)/den
      END
C
      REAL FUNCTION energy(t,x,v)
      IMPLICIT NONE
      REAL t,x(2),v(2)
      REAL a,b,e0,m1,m2,r1,r2,s1,s2
      INTEGER MS
      COMMON /PRMS/a,b,e0,m1,m2,r1,r2,s1,s2,MS
      energy = (a*v(1)**2+b*v(2)**2)/2
     +    +cos(x(1)-x(2))*v(1)*v(2)
     +    -a*cos(x(1))-cos(x(2))
      END