SUBROUTINE DMV10 (AM,Q,POTENP,H,NSTEP,RPAR,IPAR) C----------------------------------------------------------------------- C PREPROCESSED DISCRETE MOSER-VESELOV ALGORITHM C----------------------------------------------------------------------- C PREPROCESSED DISCRETE MOSER-VESELOV ALGORITHM OF ORDER 10 FOR THE C NUMERICAL SOLUTION OF THE EQUATIONS OF MOTION OF THE FREE RIGID BODY. C ORTHOGONAL MATRICES ARE REPRESENTED BY QUATERNIONS. C THE CODE IS READY TO INCLUDE AN EXERNAL POTENTIAL (SYMMETRIC STRANG C SPLITTING OF ORDER 2). C----------------------------------------------------------------------- C AUTHORS: ERNST HAIRER (1) AND GILLES VILMART (1)(2) C (1) UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES C 2-4 RUE DU LIEVRE, CASE POSTALE 64 C CH-1211 GENEVE 4, SWITZERLAND C (2) IRISA/INRIA RENNES, PROJET IPSO C CAMPUS DE BEAULIEU, F-35042 RENNES CEDEX, FRANCE C E-MAILS: Ernst.Hairer@math.unige.ch C Gilles.Vilmart@math.unige.ch C C THIS CODE IS DESCRIBED IN: C E. HAIRER AND G. VILMART, PREPROCESSED DISCRETE MOSER-VESELOV C ALGORITHM FOR THE FULL DYNAMICS OF A RIGID BODY C J. Phys. A: Math. Gen. 39 (2006) 13225-13235. C http://stacks.iop.org/0305-4470/39/13225 C C VERSION: AUGUST 29, 2006 C (latest correction of a small bug: December 28, 2007) C----------------------------------------------------------------------- C INPUT PARAMETERS C ---------------- C AM(I) INITIAL ANGULAR MOMENTUM (I=1,2,3) C Q(I) INITIAL QUATERNION FOR ORTHOGONAL MATRIX (I=1,2,3,4) C H STEP SIZE C NSTEP NUMBER OF STEPS C POTENP NAME (EXTERNAL) OF SUBROUTINE FOR AN EXTERNAL POTENTIAL C SUBROUTINE POTENP(Q,POTP,RPAR,IPAR) C DIMENSION Q(4),POTP(3) C POTP(1)=... ETC. C RPAR,IPAR REAL AND INTEGER PARAMETERS (OR PARAMETER ARRAYS) WHICH C CAN BE USED FOR COMMUNICATION BETWEEN SUBROUTINES C RPAR(11), RPAR(12), RPAR(13) ARE THE THREE MOMENTS OF INERTIA OF C THE RIGID BODY C C OUTPUT PARAMETERS C ----------------- C AM(I) SOLUTION (ANGULAR MOMENTUM) AT ENDPOINT C Q(I) SOLUTION (QUATERNION) AT ENDPOINT C----------------------------------------------------------------------- PARAMETER (MSPLIT=1) C----------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) DOUBLE PRECISION Q(4),AM(3),POTP(3) DIMENSION IPAR(20),RPAR(20) EPS=ABS(1.D-15*H) HA=H/2.0D0 HB=H/MSPLIT C HD=HB/2.0D0 HC=4.0D0/HB AI1=RPAR(11) AI2=RPAR(12) AI3=RPAR(13) C--- CONSTANT COEFFICIENTS FOR THE MODIFIED MOMENTS OF INERTIA XI1=1.0D0/AI1+1.0D0/AI2+1.0D0/AI3 XI2=1.0D0/AI1**2+1.0D0/AI2**2+1.0D0/AI3**2 XI3=1.0D0/AI1**3+1.0D0/AI2**3+1.0D0/AI3**3 XI4=1.0D0/AI1**4+1.0D0/AI2**4+1.0D0/AI3**4 XDET=AI1*AI2*AI3 XDET2=XDET**2 XDET3=XDET**3 XDET4=XDET2**2 XS1=AI1+AI2+AI3 XS2=AI1**2+AI2**2+AI3**2 XS3=AI1**3+AI2**3+AI3**3 XS4=AI1**4+AI2**4+AI3**4 XSI=(AI1+AI2)/AI3+(AI2+AI3)/AI1+(AI3+AI1)/AI2 XSI12=(AI1+AI2)/AI3**2+(AI2+AI3)/AI1**2+(AI3+AI1)/AI2**2 XSI21=(AI1**2+AI2**2)/AI3+(AI2**2+AI3**2)/AI1+(AI3**2+AI1**2)/AI2 XSI13=(AI1+AI2)/AI3**3+(AI2+AI3)/AI1**3+(AI3+AI1)/AI2**3 XSI31=(AI1**3+AI2**3)/AI3+(AI2**3+AI3**3)/AI1+(AI3**3+AI1**3)/AI2 XSI22=(AI1**2+AI2**2)/AI3**2+(AI2**2+AI3**2)/AI1**2 & +(AI3**2+AI1**2)/AI2**2 C S3C1=XS1/(6.0D0*XDET) S3C2=-XI1/3.0D0 D3C1=-1.0D0/(3.0D0*XDET) D3C2=S3C1 C S5C1=(XS2/XDET-XI1)/(30.0D0*XDET) S5C2=(1.0D0-XSI)/(30.0D0*XDET) S5C3=(3.0D0*XS1/XDET+2.0D0*XI2)/60.0D0 D5C1=-XS1/XDET2/60.0D0 D5C2=XI1/XDET/10.0D0-XS2/XDET2/60.0D0 D5C3=-(9.0D0+XSI)/XDET/60.0D0 C S7C1=(4.0D0*XDET+17.0D0*XS3-15.0D0*XDET*XSI)/(2520.0D0*XDET3) S7C2=(9.0D0*XS1+10.0D0*XDET*XI2-6.0D0*XSI21)/(420.0D0*XDET2) S7C3=((6.0D0*XSI12-1.0D2*XI1)*XDET+53.0D0*XS2)/(2520.0D0*XDET2) S7C4=(15.0D0-XDET*XI3-2.0D0*XSI)/(630.0D0*XDET) D7C1=(34.0D0*XDET*XI1-19.0D0*XS2)/(2520.0D0*XDET3) D7C2=(XS3+2.0D0*XDET*XSI-85.0D0*XDET)/(1260.0D0*XDET3) D7C3=(47.0D0*XS1+13.0D0*XSI21-38.0D0*XDET*XI2)/(2520.0D0*XDET2) D7C4=(9.0D0*XDET*XI1+XDET*XSI12-11.0D0*XS2)/(1260.0D0*XDET2) C S9C1=(62.0D0*XS4-94.0D0*XDET*XSI21 & +66.0D0*XDET2*XI2+81.0D0*XDET*XS1)/(45360.0D0*XDET4) S9C2=(-77.0D0*XSI31+75.0D0*XDET*XSI12 & +214.0D0*XS2-240.0D0*XDET*XI1)/(22680.0D0*XDET3) S9C3=(26.0D0*XDET*XSI22+55.0D0*XS3+204.0D0*XDET & -50.0D0*XDET2*XI3-59.0D0*XDET*XSI)/(7560.0D0*XDET3) S9C4=(137.0D0*XDET*XI2-XDET*XSI13 & +3.0D0*XS1-69.0D0*XSI21)/(11340.0D0*XDET2) S9C5=(2.0D0*XDET2*XI4+5.0D0*XDET*XSI12 & -171.0D0*XDET*XI1+159.0D0*XS2)/(45360.0D0*XDET2) D9C1=(60.0D0*XSI*XDET-61.0D0*XS3-247.0D0*XDET)/(45360.0D0*XDET4) D9C2=(54.0D0*XDET*XS1-XS4+218.0D0*XDET*XSI21 & -426.0D0*XDET2*XI2)/(45360.0D0*XDET4) D9C3=(125.0D0*XDET*XI1-5.0D0*XSI31-130.0D0*XS2 & +4.0D0*XDET*XSI12)/(7560.0D0*XDET3) D9C4=(67.0D0*XS3-735.0D0*XDET-15.0D0*XDET*XSI22 & +87.0D0*XDET*XSI+34.0D0*XDET2*XI3)/(22680.0D0*XDET3) D9C5=(165.0D0*XS1-XDET*XSI13-9.0D0*XSI21 & -145.0D0*XDET*XI2)/(45360.0D0*XDET2) c--- CALL POTENP(Q,POTP,RPAR,IPAR) AM1=AM(1)-HA*POTP(1) AM2=AM(2)-HA*POTP(2) AM3=AM(3)-HA*POTP(3) c DO ISTEP=1,NSTEP C --- COMPUTATION OF THE MODIFIED MOMENTS OF INERTIA HAM0=0.5D0*HB**2*(AM1**2/AI1+AM2**2/AI2+AM3**2/AI3) ANOR0=0.5D0*HB**2*(AM1**2+AM2**2+AM3**2) HAM2=HAM0**2 ANOR2=ANOR0**2 HAM3=HAM0**3 ANOR3=ANOR0**3 HAM4=HAM2**2 ANOR4=ANOR2**2 ANORHAM=HAM0*ANOR0 ANOR2HAM=HAM0*ANOR2 ANORHAM2=HAM2*ANOR0 ANOR3HAM=HAM0*ANOR3 ANOR2HAM2=HAM2*ANOR2 ANORHAM3=HAM3*ANOR0 CSS=1.0D0+(S3C1*ANOR0+S3C2*HAM0) & +(S5C1*ANOR2+S5C2*ANORHAM+S5C3*HAM2) & +(S7C1*ANOR3+S7C2*ANOR2HAM+S7C3*ANORHAM2+S7C4*HAM3) & +(S9C1*ANOR4+S9C2*ANOR3HAM+S9C3*ANOR2HAM2 & +S9C4*ANORHAM3+S9C5*HAM4) CDD=(D3C1*ANOR0+D3C2*HAM0) & +(D5C1*ANOR2+D5C2*ANORHAM+D5C3*HAM2) & +(D7C1*ANOR3+D7C2*ANOR2HAM+D7C3*ANORHAM2+D7C4*HAM3) & +(D9C1*ANOR4+D9C2*ANOR3HAM+D9C3*ANOR2HAM2 & +D9C4*ANORHAM3+D9C5*HAM4) AI1MODI=CSS/AI1+CDD AI2MODI=CSS/AI2+CDD AI3MODI=CSS/AI3+CDD AI1MOD=1.0D0/AI1MODI AI2MOD=1.0D0/AI2MODI AI3MOD=1.0D0/AI3MODI FAD1=AI2MOD-AI3MOD FAD2=AI3MOD-AI1MOD FAD3=AI1MOD-AI2MOD FAC1=FAD1*AI1MODI FAC2=FAD2*AI2MODI FAC3=FAD3*AI3MODI DO ISPLIT=1,MSPLIT C --- SOLVE FOR INTERNAL STAGE AM1I=AM1*HD*AI1MODI AM2I=AM2*HD*AI2MODI AM3I=AM3*HD*AI3MODI CM1=AM1I+FAC1*AM2I*AM3I CM2=AM2I+FAC2*CM1*AM3I CM3=AM3I+FAC3*CM1*CM2 DO I=1,50 CM1B=CM1 CM2B=CM2 CM3B=CM3 CALPHA=1+CM1**2+CM2**2+CM3**2 CM1=CALPHA*AM1I+FAC1*CM2*CM3 CM2=CALPHA*AM2I+FAC2*CM1*CM3 CM3=CALPHA*AM3I+FAC3*CM1*CM2 ERR=ABS(CM1B-CM1)+ABS(CM2B-CM2)+ABS(CM3B-CM3) IF (ERR.LT.EPS) GOTO 22 END DO 22 CONTINUE C --- UPDATE Q Q0=Q(1) Q1=Q(2) Q2=Q(3) Q3=Q(4) Q(1)=Q0-CM1*Q1-CM2*Q2-CM3*Q3 Q(2)=Q1+CM1*Q0+CM3*Q2-CM2*Q3 Q(3)=Q2+CM2*Q0+CM1*Q3-CM3*Q1 Q(4)=Q3+CM3*Q0+CM2*Q1-CM1*Q2 C --- UPDATE M CALPHA=HC/CALPHA AM1=AM1+FAD1*CM2*CM3*CALPHA AM2=AM2+FAD2*CM1*CM3*CALPHA AM3=AM3+FAD3*CM1*CM2*CALPHA END DO C --- PROJECTION QUAT=1.0D0/SQRT(Q(1)**2+Q(2)**2+Q(3)**2+Q(4)**2) Q(1)=Q(1)*QUAT Q(2)=Q(2)*QUAT Q(3)=Q(3)*QUAT Q(4)=Q(4)*QUAT C CALL POTENP(Q,POTP,RPAR,IPAR) AM1=AM1-H*POTP(1) AM2=AM2-H*POTP(2) AM3=AM3-H*POTP(3) END DO AM(1)=AM1+HA*POTP(1) AM(2)=AM2+HA*POTP(2) AM(3)=AM3+HA*POTP(3) RETURN END C C----------------------------------------------------------------------- C SUBROUTINE DMV2 (AM,Q,POTENP,H,NSTEP,RPAR,IPAR) C--- DISCRETE MOSER-VESELOV ALGORITHM OF ORDER 2 C--- ORTHOGONAL MATRICES ARE REPRESENTED BY QUATERNIONS IMPLICIT REAL*8 (A-H,O-Z) DOUBLE PRECISION Q(4),AM(3),POTP(3) DIMENSION IPAR(20),RPAR(20) EPS=ABS(1.D-15*H) HA=H/2.0D0 AI1=RPAR(11) AI2=RPAR(12) AI3=RPAR(13) HC=4.0D0/H FAD1=HC*(AI2-AI3) FAD2=HC*(AI3-AI1) FAD3=HC*(AI1-AI2) FAC1=(AI2-AI3)/AI1 FAC2=(AI3-AI1)/AI2 FAC3=(AI1-AI2)/AI3 c--- CALL POTENP(Q,POTP,RPAR,IPAR) AM1=AM(1)-HA*POTP(1) AM2=AM(2)-HA*POTP(2) AM3=AM(3)-HA*POTP(3) c DO ISTEP=1,NSTEP C --- SOLVE FOR INTERNAL STAGE AM1I=AM1*HA/AI1 AM2I=AM2*HA/AI2 AM3I=AM3*HA/AI3 CM1=AM1I+FAC1*AM2I*AM3I CM2=AM2I+FAC2*CM1*AM3I CM3=AM3I+FAC3*CM1*CM2 DO I=1,50 CM1B=CM1 CM2B=CM2 CM3B=CM3 CALPHA=1+CM1**2+CM2**2+CM3**2 CM1=CALPHA*AM1I+FAC1*CM2*CM3 CM2=CALPHA*AM2I+FAC2*CM1*CM3 CM3=CALPHA*AM3I+FAC3*CM1*CM2 ERR=ABS(CM1B-CM1)+ABS(CM2B-CM2)+ABS(CM3B-CM3) IF (ERR.LT.EPS) GOTO 22 END DO 22 CONTINUE C --- UPDATE Q Q0=Q(1) Q1=Q(2) Q2=Q(3) Q3=Q(4) Q(1)=Q0-CM1*Q1-CM2*Q2-CM3*Q3 Q(2)=Q1+CM1*Q0+CM3*Q2-CM2*Q3 Q(3)=Q2+CM2*Q0+CM1*Q3-CM3*Q1 Q(4)=Q3+CM3*Q0+CM2*Q1-CM1*Q2 C --- PROJECTION QUAT=1.0D0/SQRT(Q(1)**2+Q(2)**2+Q(3)**2+Q(4)**2) Q(1)=Q(1)*QUAT Q(2)=Q(2)*QUAT Q(3)=Q(3)*QUAT Q(4)=Q(4)*QUAT C --- UPDATE M CALL POTENP(Q,POTP,RPAR,IPAR) AM1=AM1-H*POTP(1)+FAD1*CM2*CM3/CALPHA AM2=AM2-H*POTP(2)+FAD2*CM1*CM3/CALPHA AM3=AM3-H*POTP(3)+FAD3*CM1*CM2/CALPHA END DO AM(1)=AM1+HA*POTP(1) AM(2)=AM2+HA*POTP(2) AM(3)=AM3+HA*POTP(3) RETURN END C SUBROUTINE DMV4 (AM,Q,POTENP,H,NSTEP,RPAR,IPAR) C--- PREPROCESSED DISCRETE MOSER-VESELOV ALGORITHM OF ORDER 4 C--- ORTHOGONAL MATRICES ARE REPRESENTED BY QUATERNIONS IMPLICIT REAL*8 (A-H,O-Z) DOUBLE PRECISION Q(4),AM(3),POTP(3) DIMENSION IPAR(20),RPAR(20) EPS=ABS(1.D-15*H) HA=H/2.0D0 AI1=RPAR(11) AI2=RPAR(12) AI3=RPAR(13) HC=4.0D0/H C XI1=1.0D0/AI1+1.0D0/AI2+1.0D0/AI3 XDET=AI1*AI2*AI3 XS1=AI1+AI2+AI3 S3C1=XS1/(6.0D0*XDET) S3C2=-XI1/3.0D0 D3C1=-1.0D0/(3.0D0*XDET) D3C2=S3C1 c--- CALL POTENP(Q,POTP,RPAR,IPAR) AM1=AM(1)-HA*POTP(1) AM2=AM(2)-HA*POTP(2) AM3=AM(3)-HA*POTP(3) c DO ISTEP=1,NSTEP HAM0=0.5D0*(AM1**2/AI1+AM2**2/AI2+AM3**2/AI3) ANOR0=0.5D0*(AM1**2+AM2**2+AM3**2) CSS=1.0D0+H**2*(S3C1*ANOR0+S3C2*HAM0) CDD=H**2*(D3C1*ANOR0+D3C2*HAM0) AI1MODI=CSS/AI1+CDD AI2MODI=CSS/AI2+CDD AI3MODI=CSS/AI3+CDD AI1MOD=1.0D0/AI1MODI AI2MOD=1.0D0/AI2MODI AI3MOD=1.0D0/AI3MODI FAD1=AI2MOD-AI3MOD FAD2=AI3MOD-AI1MOD FAD3=AI1MOD-AI2MOD FAC1=FAD1*AI1MODI FAC2=FAD2*AI2MODI FAC3=FAD3*AI3MODI C --- SOLVE FOR INTERNAL STAGE AM1I=AM1*HA*AI1MODI AM2I=AM2*HA*AI2MODI AM3I=AM3*HA*AI3MODI CM1=AM1I+FAC1*AM2I*AM3I CM2=AM2I+FAC2*CM1*AM3I CM3=AM3I+FAC3*CM1*CM2 DO I=1,50 CM1B=CM1 CM2B=CM2 CM3B=CM3 CALPHA=1+CM1**2+CM2**2+CM3**2 CM1=CALPHA*AM1I+FAC1*CM2*CM3 CM2=CALPHA*AM2I+FAC2*CM1*CM3 CM3=CALPHA*AM3I+FAC3*CM1*CM2 ERR=ABS(CM1B-CM1)+ABS(CM2B-CM2)+ABS(CM3B-CM3) IF (ERR.LT.EPS) GOTO 22 END DO 22 CONTINUE C --- UPDATE Q Q0=Q(1) Q1=Q(2) Q2=Q(3) Q3=Q(4) Q(1)=Q0-CM1*Q1-CM2*Q2-CM3*Q3 Q(2)=Q1+CM1*Q0+CM3*Q2-CM2*Q3 Q(3)=Q2+CM2*Q0+CM1*Q3-CM3*Q1 Q(4)=Q3+CM3*Q0+CM2*Q1-CM1*Q2 C --- PROJECTION QUAT=1.0D0/SQRT(Q(1)**2+Q(2)**2+Q(3)**2+Q(4)**2) Q(1)=Q(1)*QUAT Q(2)=Q(2)*QUAT Q(3)=Q(3)*QUAT Q(4)=Q(4)*QUAT C --- UPDATE M CALL POTENP(Q,POTP,RPAR,IPAR) CALPHA=HC/CALPHA AM1=AM1-H*POTP(1)+FAD1*CM2*CM3*CALPHA AM2=AM2-H*POTP(2)+FAD2*CM1*CM3*CALPHA AM3=AM3-H*POTP(3)+FAD3*CM1*CM2*CALPHA END DO AM(1)=AM1+HA*POTP(1) AM(2)=AM2+HA*POTP(2) AM(3)=AM3+HA*POTP(3) RETURN END C SUBROUTINE DMV6 (AM,Q,POTENP,H,NSTEP,RPAR,IPAR) C--- PREPROCESSED DISCRETE MOSER-VESELOV ALGORITHM OF ORDER 6 C--- ORTHOGONAL MATRICES ARE REPRESENTED BY QUATERNIONS IMPLICIT REAL*8 (A-H,O-Z) DOUBLE PRECISION Q(4),AM(3),POTP(3) DIMENSION IPAR(20),RPAR(20) EPS=ABS(1.D-15*H) HA=H/2.0D0 AI1=RPAR(11) AI2=RPAR(12) AI3=RPAR(13) HC=4.0D0/H C XI1=1.0D0/AI1+1.0D0/AI2+1.0D0/AI3 XI2=1.0D0/AI1**2+1.0D0/AI2**2+1.0D0/AI3**2 XDET=AI1*AI2*AI3 XS1=AI1+AI2+AI3 XS2=AI1**2+AI2**2+AI3**2 XSI=(AI1+AI2)/AI3+(AI2+AI3)/AI1+(AI3+AI1)/AI2 XSJ=AI1/AI2/AI3+AI2/AI3/AI1+AI3/AI1/AI2 C S3C1=XS1/(6.0D0*XDET) S3C2=-XI1/3.0D0 D3C1=-1.0D0/(3.0D0*XDET) D3C2=S3C1 C S5C1=(XS2/XDET-XI1)/(30.0D0*XDET) S5C2=(1.0D0-XSI)/(30.0D0*XDET) S5C3=(3.0D0*XS1/XDET+2.0D0*XI2)/60.0D0 D5C1=-XS1/XDET**2/60.0D0 D5C2=XI1/XDET/10.0D0-XS2/XDET**2/60.0D0 D5C3=-(9.0D0+XSI)/XDET/60.0D0 c--- CALL POTENP(Q,POTP,RPAR,IPAR) AM1=AM(1)-HA*POTP(1) AM2=AM(2)-HA*POTP(2) AM3=AM(3)-HA*POTP(3) c DO ISTEP=1,NSTEP HAM0=0.5D0*H**2*(AM1**2/AI1+AM2**2/AI2+AM3**2/AI3) ANOR0=0.5D0*H**2*(AM1**2+AM2**2+AM3**2) HAM2=HAM0**2 ANOR2=ANOR0**2 ANORHAM=HAM0*ANOR0 CSS=1.0D0+(S3C1*ANOR0+S3C2*HAM0) & +(S5C1*ANOR2+S5C2*ANORHAM+S5C3*HAM2) CDD=(D3C1*ANOR0+D3C2*HAM0) & +(D5C1*ANOR2+D5C2*ANORHAM+D5C3*HAM2) AI1MODI=CSS/AI1+CDD AI2MODI=CSS/AI2+CDD AI3MODI=CSS/AI3+CDD AI1MOD=1.0D0/AI1MODI AI2MOD=1.0D0/AI2MODI AI3MOD=1.0D0/AI3MODI FAD1=AI2MOD-AI3MOD FAD2=AI3MOD-AI1MOD FAD3=AI1MOD-AI2MOD FAC1=FAD1*AI1MODI FAC2=FAD2*AI2MODI FAC3=FAD3*AI3MODI C --- SOLVE FOR INTERNAL STAGE AM1I=AM1*HA*AI1MODI AM2I=AM2*HA*AI2MODI AM3I=AM3*HA*AI3MODI CM1=AM1I+FAC1*AM2I*AM3I CM2=AM2I+FAC2*CM1*AM3I CM3=AM3I+FAC3*CM1*CM2 DO I=1,50 CM1B=CM1 CM2B=CM2 CM3B=CM3 CALPHA=1+CM1**2+CM2**2+CM3**2 CM1=CALPHA*AM1I+FAC1*CM2*CM3 CM2=CALPHA*AM2I+FAC2*CM1*CM3 CM3=CALPHA*AM3I+FAC3*CM1*CM2 ERR=ABS(CM1B-CM1)+ABS(CM2B-CM2)+ABS(CM3B-CM3) IF (ERR.LT.EPS) GOTO 22 END DO 22 CONTINUE C --- UPDATE Q Q0=Q(1) Q1=Q(2) Q2=Q(3) Q3=Q(4) Q(1)=Q0-CM1*Q1-CM2*Q2-CM3*Q3 Q(2)=Q1+CM1*Q0+CM3*Q2-CM2*Q3 Q(3)=Q2+CM2*Q0+CM1*Q3-CM3*Q1 Q(4)=Q3+CM3*Q0+CM2*Q1-CM1*Q2 C --- PROJECTION QUAT=1.0D0/SQRT(Q(1)**2+Q(2)**2+Q(3)**2+Q(4)**2) Q(1)=Q(1)*QUAT Q(2)=Q(2)*QUAT Q(3)=Q(3)*QUAT Q(4)=Q(4)*QUAT C --- UPDATE M CALL POTENP(Q,POTP,RPAR,IPAR) CALPHA=HC/CALPHA AM1=AM1-H*POTP(1)+FAD1*CM2*CM3*CALPHA AM2=AM2-H*POTP(2)+FAD2*CM1*CM3*CALPHA AM3=AM3-H*POTP(3)+FAD3*CM1*CM2*CALPHA END DO AM(1)=AM1+HA*POTP(1) AM(2)=AM2+HA*POTP(2) AM(3)=AM3+HA*POTP(3) RETURN END C SUBROUTINE DMV8 (AM,Q,POTENP,H,NSTEP,RPAR,IPAR) C--- PREPROCESSED DISCRETE MOSER-VESELOV ALGORITHM OF ORDER 8 C--- ORTHOGONAL MATRICES ARE REPRESENTED BY QUATERNIONS IMPLICIT REAL*8 (A-H,O-Z) DOUBLE PRECISION Q(4),AM(3),POTP(3) DIMENSION IPAR(20),RPAR(20) EPS=ABS(1.D-15*H) HA=H/2.0D0 AI1=RPAR(11) AI2=RPAR(12) AI3=RPAR(13) HC=4.0D0/H C--- CONSTANT COEFFICIENTS FOR THE MODIFIED MOMENTS OF INERTIA XI1=1.0D0/AI1+1.0D0/AI2+1.0D0/AI3 XI2=1.0D0/AI1**2+1.0D0/AI2**2+1.0D0/AI3**2 XI3=1.0D0/AI1**3+1.0D0/AI2**3+1.0D0/AI3**3 XDET=AI1*AI2*AI3 XS1=AI1+AI2+AI3 XS2=AI1**2+AI2**2+AI3**2 XS3=AI1**3+AI2**3+AI3**3 XSI=(AI1+AI2)/AI3+(AI2+AI3)/AI1+(AI3+AI1)/AI2 XSI12=(AI1+AI2)/AI3**2+(AI2+AI3)/AI1**2+(AI3+AI1)/AI2**2 XSI21=(AI1**2+AI2**2)/AI3+(AI2**2+AI3**2)/AI1+(AI3**2+AI1**2)/AI2 C S3C1=XS1/(6.0D0*XDET) S3C2=-XI1/3.0D0 D3C1=-1.0D0/(3.0D0*XDET) D3C2=S3C1 C S5C1=(XS2/XDET-XI1)/(30.0D0*XDET) S5C2=(1.0D0-XSI)/(30.0D0*XDET) S5C3=(3.0D0*XS1/XDET+2.0D0*XI2)/60.0D0 D5C1=-XS1/XDET**2/60.0D0 D5C2=XI1/XDET/10.0D0-XS2/XDET**2/60.0D0 D5C3=-(9.0D0+XSI)/XDET/60.0D0 C S7C1=(4.0D0*XDET+17.0D0*XS3-15.0D0*XDET*XSI)/(2520.0D0*XDET**3) S7C2=(9.0D0*XS1+10.0D0*XDET*XI2-6.0D0*XSI21)/(420.0D0*XDET**2) S7C3=((6.0D0*XSI12-1.0D2*XI1)*XDET+53.0D0*XS2)/(2520.0D0*XDET**2) S7C4=(15.0D0-XDET*XI3-2.0D0*XSI)/(630.0D0*XDET) D7C1=(34.0D0*XDET*XI1-19.0D0*XS2)/(2520.0D0*XDET**3) D7C2=(XS3+2.0D0*XDET*XSI-85.0D0*XDET)/(1260.0D0*XDET**3) D7C3=(47.0D0*XS1+13.0D0*XSI21-38.0D0*XDET*XI2)/(2520.0D0*XDET**2) D7C4=(9.0D0*XDET*XI1+XDET*XSI12-11.0D0*XS2)/(1260.0D0*XDET**2) c--- CALL POTENP(Q,POTP,RPAR,IPAR) AM1=AM(1)-HA*POTP(1) AM2=AM(2)-HA*POTP(2) AM3=AM(3)-HA*POTP(3) c DO ISTEP=1,NSTEP C --- COMPUTATION OF THE MODIFIED MOMENTS OF INERTIA HAM0=0.5D0*H**2*(AM1**2/AI1+AM2**2/AI2+AM3**2/AI3) ANOR0=0.5D0*H**2*(AM1**2+AM2**2+AM3**2) HAM2=HAM0**2 ANOR2=ANOR0**2 HAM3=HAM0**3 ANOR3=ANOR0**3 ANORHAM=HAM0*ANOR0 ANOR2HAM=HAM0*ANOR2 ANORHAM2=HAM2*ANOR0 CSS=1.0D0+(S3C1*ANOR0+S3C2*HAM0) & +(S5C1*ANOR2+S5C2*ANORHAM+S5C3*HAM2) & +(S7C1*ANOR3+S7C2*ANOR2HAM+S7C3*ANORHAM2+S7C4*HAM3) CDD=(D3C1*ANOR0+D3C2*HAM0) & +(D5C1*ANOR2+D5C2*ANORHAM+D5C3*HAM2) & +(D7C1*ANOR3+D7C2*ANOR2HAM+D7C3*ANORHAM2+D7C4*HAM3) AI1MODI=CSS/AI1+CDD AI2MODI=CSS/AI2+CDD AI3MODI=CSS/AI3+CDD AI1MOD=1.0D0/AI1MODI AI2MOD=1.0D0/AI2MODI AI3MOD=1.0D0/AI3MODI FAD1=AI2MOD-AI3MOD FAD2=AI3MOD-AI1MOD FAD3=AI1MOD-AI2MOD FAC1=FAD1*AI1MODI FAC2=FAD2*AI2MODI FAC3=FAD3*AI3MODI C --- SOLVE FOR INTERNAL STAGE AM1I=AM1*HA*AI1MODI AM2I=AM2*HA*AI2MODI AM3I=AM3*HA*AI3MODI CM1=AM1I+FAC1*AM2I*AM3I CM2=AM2I+FAC2*CM1*AM3I CM3=AM3I+FAC3*CM1*CM2 DO I=1,50 CM1B=CM1 CM2B=CM2 CM3B=CM3 CALPHA=1+CM1**2+CM2**2+CM3**2 CM1=CALPHA*AM1I+FAC1*CM2*CM3 CM2=CALPHA*AM2I+FAC2*CM1*CM3 CM3=CALPHA*AM3I+FAC3*CM1*CM2 ERR=ABS(CM1B-CM1)+ABS(CM2B-CM2)+ABS(CM3B-CM3) IF (ERR.LT.EPS) GOTO 22 END DO 22 CONTINUE C --- UPDATE Q Q0=Q(1) Q1=Q(2) Q2=Q(3) Q3=Q(4) Q(1)=Q0-CM1*Q1-CM2*Q2-CM3*Q3 Q(2)=Q1+CM1*Q0+CM3*Q2-CM2*Q3 Q(3)=Q2+CM2*Q0+CM1*Q3-CM3*Q1 Q(4)=Q3+CM3*Q0+CM2*Q1-CM1*Q2 C --- PROJECTION QUAT=1.0D0/SQRT(Q(1)**2+Q(2)**2+Q(3)**2+Q(4)**2) Q(1)=Q(1)*QUAT Q(2)=Q(2)*QUAT Q(3)=Q(3)*QUAT Q(4)=Q(4)*QUAT C --- UPDATE M CALL POTENP(Q,POTP,RPAR,IPAR) CALPHA=HC/CALPHA AM1=AM1-H*POTP(1)+FAD1*CM2*CM3*CALPHA AM2=AM2-H*POTP(2)+FAD2*CM1*CM3*CALPHA AM3=AM3-H*POTP(3)+FAD3*CM1*CM2*CALPHA END DO AM(1)=AM1+HA*POTP(1) AM(2)=AM2+HA*POTP(2) AM(3)=AM3+HA*POTP(3) RETURN END