C * * * * * * * * * * * * * * * * * * * * * * * * * * *
C --- DRIVER FOR DMV10 FOR THE MOTION OF A RIGID BODY
C * * * * * * * * * * * * * * * * * * * * * * * * * * *
        include 'dmv10.f'
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        DIMENSION Q(4),AM(3),RPAR(20),IPAR(20)
        REAL TIME0,TIME1
        EXTERNAL POTENP
C --- INITIAL VALUES
C ANGULAR MOMENTUM
        AM(1)=1.8D0
        AM(2)=0.4D0
        AM(3)=-0.9D0
C QUATERNION
        Q(1)=1.0D0
        Q(2)=0.0D0
        Q(3)=0.0D0
        Q(4)=0.0D0
C --- MOMENTS OF INERTIA
        AI1=0.6D0
        AI2=0.8D0
        AI3=1.0D0
        RPAR(11)=AI1
        RPAR(12)=AI2
        RPAR(13)=AI3
C ---
        H=0.01
        XEND=10.0D0
        NSTEP=XEND/H
        H=XEND/NSTEP
        WRITE (6,*) 'XEND=',XEND,' H=',H,' NSTEP=',NSTEP
C ---
        DO I=1,10
          RPAR(I)=0.0D0
          IPAR(I)=0
        END DO
C
        WRITE (6,*) '--- INITIAL CONDITION'
        WRITE (6,*) '  AM    ',AM(1),AM(2),AM(3)
        WRITE (6,*) '  QQ    ',Q(1),Q(2),Q(3),Q(4)
        CALL HAMIL(Q,AM,HAM0,RPAR,IPAR)
        CALL CPU_TIME(TIME0)
C ---
        CALL DMV10 (AM,Q,POTENP,H,NSTEP,RPAR,IPAR)
C ---
        CALL CPU_TIME(TIME1)
        CALL HAMIL(Q,AM,HAM1,RPAR,IPAR)
        WRITE (6,*) '--- SOLUTION AT ENDPOINT'
        WRITE (6,*) '  AM    ',AM(1),AM(2),AM(3)
        WRITE (6,*) '  QQ    ',Q(1),Q(2),Q(3),Q(4)
        WRITE (6,*) 'ERR HAM=',(HAM1-HAM0)/HAM0,' TIME=',TIME1-TIME0
        STOP
        END
C
        SUBROUTINE HAMIL (Q,AM,HAM,RPAR,IPAR)
        IMPLICIT REAL*8 (A-H,O-Z) 
        DIMENSION Q(4),AM(3)
        DIMENSION IPAR(*),RPAR(*)
        CALL POTEN(Q,POT,RPAR,IPAR)
        HAM=AM(1)**2/RPAR(11)+AM(2)**2/RPAR(12)+AM(3)**2/RPAR(13)
        HAM=HAM/2.0D0+POT
        RETURN
        END 
c
        SUBROUTINE POTEN(Q,POT,RPAR,IPAR)
        IMPLICIT REAL*8 (A-H,O-Z) 
        DIMENSION Q(4)
        DIMENSION IPAR(*),RPAR(*)
        POT=Q(1)**2-Q(2)**2-Q(3)**2+Q(4)**2
        POT=0.0D0
        RETURN
        END 
c
        SUBROUTINE POTENP(Q,POTP,RPAR,IPAR)
        IMPLICIT REAL*8 (A-H,O-Z) 
        DIMENSION Q(4),POTP(3)
        DIMENSION IPAR(*),RPAR(*)
        POTP(1)=-2*(Q(1)*Q(2)+Q(3)*Q(4))
        POTP(2)=-2*(Q(1)*Q(3)-Q(2)*Q(4))
        POTP(1)=0.0d0
        POTP(2)=0.0d0
        POTP(3)=0.0d0
        RETURN
        END