      SUBROUTINE DMV10CONJ (AM,Q,AMVAR,QVAR,H,NSTEP,X,ISCONJ,RPAR,IPAR)
C-----------------------------------------------------------------------
C            PREPROCESSED DISCRETE MOSER-VESELOV ALGORITHM OF ORDER 10
C            FOR THE COMPUTATION OF THE MOTION OF A FREE RIGID BODY
C            WITH SIMULTANEOUS COMPUTATION OF THE TANGENT MAP 
C            FOR THE COMPUTATION OF FIRST CONJUGATE POINTS
C-----------------------------------------------------------------------
C  ORTHOGONAL MATRICES ARE REPRESENTED BY QUATERNIONS. 
C-----------------------------------------------------------------------
C  THIS CODE IS DESCRIBED IN SECTION 3.6 OF:
C    G. VILMART, ÉTUDE D'INTÉGRATEURS GÉOMÉTRIQUES POUR DES ÉQUATIONS 
C    DIFFÉRENTIELLES
C    Ph.D. Thesis, Univ. Rennes 1 (2008) Sc. 3758, Univ. Genève (2008) Sc. 4038.
C    http://archive-ouverte.unige.ch/vital/access/manager/Repository/unige:1337
C    ftp://ftp.irisa.fr/techreports/theses/2008/vilmart.pdf
C   
C  THIS CODE IS A MODIFICATION OF THE CODE 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     AUTHORS: ERNST HAIRER AND GILLES VILMART
C                    UNIVERSITE DE GENEVE, DEPT. DE MATHEMATIQUES
C                    2-4 RUE DU LIEVRE, CASE POSTALE 64
C                    CH-1211 GENEVE 4, SWITZERLAND 
C     E-MAILS: Ernst.Hairer@unige.ch
C              Gilles.Vilmart@unige.ch
C
C     VERSION: MAY 17, 2007
C     LAST UPDATE: JULY 30, 2009 
C                  (added compensated summations for reducing roundoff errors)
C-----------------------------------------------------------------------
C  For given moments of inertia and angular momentum, this code
C  computes the the derivative AMVAR of the angular momentum with respect 
C  to the initial angular momentum, and the derivative QVAR of the attitude 
C  matrix with respect to the initial angular momentum (see details below). 
C  The first conjugate point is detected when matrix QVAR becomes singular 
C  using a bisection method.
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      MAXIMUM NUMBER OF STEPS
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     OUTPUT AT FIRST CONJUGATE POINT:
C     ISCONJ      BOOLEAN VARIABLE, TRUE IF A CONJUGATE POINT IS DETECTED
C     X           TIME OF FIRST CONJUGATE POINT
C     AM(I)       ANGULAR MOMENTUM
C     Q(I)        QUATERNION (ORTHOGONAL MATRIX)
C     AMVAR(I,I)  3x3 MATRIX OF DERIVATIVES OF ANGULAR MOMENTUM 
C                 WITH RESPECT TO THE INITIAL ANGULAR MOMENTUM
C     QVAR(I,I)   3x3 MATRIX OF DERIVATIVES FOR ATTITUDE MATRIX
C                 WITH RESPECT TO THE INITIAL ANGULAR MOMENTUM
C     QVAR is a 3x3 matrix that contains in the j-th column the vector a_j 
C     corresponding to a skew-symmetric matrix A_j = \widehat a_j, such that 
C     the derivative of the orthogonal matrix R (represented by the quaternion Q) 
C     with respect to the initial angular momentum equals the product  R A_j.
C----------------------------------------------------------------------- 
      IMPLICIT REAL*8 (A-H,O-Z) 
      DOUBLE PRECISION Q(4),AM(3),QVAR(3,3),AMVAR(3,3)
      DOUBLE PRECISION QCS(4),AMCS(3),QVARCS(3,3),AMVARCS(3,3),
     &    QCSOLD(4), AMCSOLD(3),QVARCSOLD(3,3),AMVARCSOLD(3,3)
      DOUBLE PRECISION AMAT(3,3),BMAT(3,3),CMAT(3,3)
      DOUBLE PRECISION QOLD(4),AMOLD(3),QVAROLD(3,3),AMVAROLD(3,3)
      DIMENSION IPAR(20),RPAR(20)
      LOGICAL ISCONJ
      ISCONJ=.FALSE.
      X=0.0D0
      EPS=1.D-16
      HB=H
C
      AI1=RPAR(11)
      AI2=RPAR(12)
      AI3=RPAR(13) 
c constants for variational eq.
      AI1I=1.0D0/AI1
      AI2I=1.0D0/AI2
      AI3I=1.0D0/AI3
      FAE1=AI3I-AI2I
      FAE2=AI1I-AI3I
      FAE3=AI2I-AI1I
c      
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--- INITIALISATION FOR COMPENSATED SUMMATION
      XCS=0.0D0
      DO I=1,3
      AMCS(I)=0.0D0
      DO J=1,3
      QVARCS(I,J)=0.0D0
      AMVARCS(I,J)=0.0D0
      END DO
      END DO
      DO I=1,4
      QCS(I)=0.0D0
      END DO
c---
c Loop with detection of conjugate points at each step
      ISTEP=0
      DO WHILE ((ISTEP.LT.NSTEP.OR.ISCONJ).AND.(ABS(HB).GT.EPS))
        ISTEP=ISTEP+1
c      DO ISTEP=1,NSTEP
c        X=X+HB
        XCS=XCS+HB
        XOLD=X
        X=X+XCS
        XCS=XCS+(XOLD-X)
        AM1=AM(1)
        AM2=AM(2)
        AM3=AM(3)
        HD=HB/2.0D0
        HC=4.0D0/HB
C ---      COMPUTATION OF THE MODIFIED MOMENTS OF INERTIA
c If there is no external potential, these coefficients 
c can be computed only once (and not at each step).
        IF (ISTEP.EQ.1.OR.ISCONJ) THEN
        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
C Constants for variational eq.
        CSSDERH=HB**2*(S3C2
     &     +(S5C2*ANOR0+2.0D0*S5C3*HAM0)
     &     +(S7C2*ANOR2+2.0D0*S7C3*ANORHAM+3.0D0*S7C4*HAM2)
     &     +(S9C2*ANOR3+2.0D0*S9C3*ANOR2HAM
     &     +3.0D0*S9C4*ANORHAM2+4.0D0*S9C5*HAM3))
        CDDDERH=HB**2*(D3C2
     &     +(D5C2*ANOR0+2.0D0*D5C3*HAM0)
     &     +(D7C2*ANOR2+2.0D0*D7C3*ANORHAM+3.0D0*D7C4*HAM2)
     &     +(D9C2*ANOR3+2.0D0*D9C3*ANOR2HAM
     &     +3.0D0*D9C4*ANORHAM2+4.0D0*D9C5*HAM3))
        CSSDERC=HB**2*(S3C1
     &     +(2.0D0*S5C1*ANOR0+S5C2*HAM0)
     &     +(3.0D0*S7C1*ANOR2+2.0D0*S7C2*ANORHAM+S7C3*HAM2)
     &     +(4.0D0*S9C1*ANOR3+3.0D0*S9C2*ANOR2HAM+2.0D0*S9C3*ANORHAM2
     &     +S9C4*HAM3))
        CDDDERC=HB**2*(D3C1
     &     +(2.0D0*D5C1*ANOR0+D5C2*HAM0)
     &     +(3.0D0*D7C1*ANOR2+2.0D0*D7C2*ANORHAM+D7C3*HAM2)
     &     +(4.0D0*D9C1*ANOR3+3.0D0*D9C2*ANOR2HAM+2.0D0*D9C3*ANORHAM2
     &     +D9C4*HAM3))
        CDSDERH=(CDDDERH*CSS-CDD*CSSDERH)/CSS**2
        CDSDERC=(CDDDERC*CSS-CDD*CSSDERC)/CSS**2
        CSSDER1=AI1I*CSSDERH+CSSDERC
        CSSDER2=AI2I*CSSDERH+CSSDERC
        CSSDER3=AI3I*CSSDERH+CSSDERC
        CDSDER1=AI1I*CDSDERH+CDSDERC
        CDSDER2=AI2I*CDSDERH+CDSDERC
        CDSDER3=AI3I*CDSDERH+CDSDERC
        END IF
C---------------
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
          VAR=1.0D50
          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
            VAROLD=VAR
            VAR=ABS(CM1B-CM1)+ABS(CM2B-CM2)+ABS(CM3B-CM3)
            IF (VAR.EQ.0.0D0.OR.VAR.GT.VAROLD) GOTO 22
          END DO
 22       CONTINUE
C ---      UPDATE Q
          Q0=Q(1)
          Q1=Q(2)
          Q2=Q(3)
          Q3=Q(4)
c Compensated summation for rotation matrix
c          Q(1)=Q0-CM1*Q1-CM2*Q2-CM3*Q3
c          Q(2)=Q1+CM1*Q0+CM3*Q2-CM2*Q3
c          Q(3)=Q2+CM2*Q0+CM1*Q3-CM3*Q1
c          Q(4)=Q3+CM3*Q0+CM2*Q1-CM1*Q2
          QCS(1)=QCS(1)-CM1*Q1-CM2*Q2-CM3*Q3
          QCS(2)=QCS(2)+CM1*Q0+CM3*Q2-CM2*Q3
          QCS(3)=QCS(3)+CM2*Q0+CM1*Q3-CM3*Q1
          QCS(4)=QCS(4)+CM3*Q0+CM2*Q1-CM1*Q2
          Q(1)=Q0+QCS(1)
          Q(2)=Q1+QCS(2)
          Q(3)=Q2+QCS(3)
          Q(4)=Q3+QCS(4)
          QCS(1)=QCS(1)+(Q0-Q(1))
          QCS(2)=QCS(2)+(Q1-Q(2))
          QCS(3)=QCS(3)+(Q2-Q(3))
          QCS(4)=QCS(4)+(Q3-Q(4))
C To avoid large quaternion
C (norm of quaternion is not one to allow using compensated summation)
          IF (ABS(Q(1))+ABS(Q(2))+ABS(Q(3))+ABS(Q(4)).GT.4) THEN
          DO I=1,4
          Q(I)=Q(I)*0.25D0
          QCS(I)=QCS(I)*0.25D0
          END DO
          END IF
C ---       UPDATE MOMENTUM
          CALPHA1=4.0D0/(HB*CALPHA*AI1MODI)
          CALPHA2=4.0D0/(HB*CALPHA*AI2MODI)
          CALPHA3=4.0D0/(HB*CALPHA*AI3MODI)
C COMPENSATED SUMMATION
          AMCS(1)=AMCS(1)+CM2*CM3*(CALPHA2-CALPHA3)
          AMCS(2)=AMCS(2)+CM3*CM1*(CALPHA3-CALPHA1)
          AMCS(3)=AMCS(3)+CM1*CM2*(CALPHA1-CALPHA2)
          AM(1)=AM1+AMCS(1)
          AM(2)=AM2+AMCS(2)
          AM(3)=AM3+AMCS(3)
          AMCS(1)=AMCS(1)+(AM1-AM(1))
          AMCS(2)=AMCS(2)+(AM2-AM(2))
          AMCS(3)=AMCS(3)+(AM3-AM(3))
C --- VARIATIONNAL EQUATION
        AZ1=CSS*AI1MOD*CM1
        AZ2=CSS*AI2MOD*CM2
        AZ3=CSS*AI3MOD*CM3
C SOLVE linear system AMAT dz_n/dy_0 = BMAT and store result in BMAT
c left matrix
        CM1I=HB*AI1MODI*CM1
        CM2I=HB*AI2MODI*CM2
        CM3I=HB*AI3MODI*CM3
        AMAT(1,1)=1.0D0    -CM1I*AM1
        AMAT(1,2)=-FAE1*AZ3-CM2I*AM1
        AMAT(1,3)=-FAE1*AZ2-CM3I*AM1
        AMAT(2,1)=-FAE2*AZ3-CM1I*AM2
        AMAT(2,2)=1.0D0    -CM2I*AM2
        AMAT(2,3)=-FAE2*AZ1-CM3I*AM2
        AMAT(3,1)=-FAE3*AZ2-CM1I*AM3
        AMAT(3,2)=-FAE3*AZ1-CM2I*AM3
        AMAT(3,3)=1.0D0    -CM3I*AM3
c right matrix
        CLAMBDA=CALPHA*HD
        CB=CSS*CLAMBDA
        CMU=HB*(AI1MODI*AZ1**2+AI2MODI*AZ2**2+AI3MODI*AZ3**2)
        CC1=CLAMBDA*CSSDER1+CMU*CDSDER1
        CC2=CLAMBDA*CSSDER2+CMU*CDSDER2
        CC3=CLAMBDA*CSSDER3+CMU*CDSDER3
c compute CMAT(I,J)=AMVAR(I,J)*AMJ I,J=1,2,3
        CMAT(1,1)=AMVAR(1,1)*AM1
        CMAT(2,1)=AMVAR(2,1)*AM2
        CMAT(3,1)=AMVAR(3,1)*AM3
        CMAT(1,2)=AMVAR(1,2)*AM1
        CMAT(2,2)=AMVAR(2,2)*AM2
        CMAT(3,2)=AMVAR(3,2)*AM3
        CMAT(1,3)=AMVAR(1,3)*AM1
        CMAT(2,3)=AMVAR(2,3)*AM2
        CMAT(3,3)=AMVAR(3,3)*AM3
c
        CC0=CC1*CMAT(1,1)+CC2*CMAT(2,1)+CC3*CMAT(3,1)
        BMAT(1,1)=CB*AMVAR(1,1)+AM1*CC0
        BMAT(2,1)=CB*AMVAR(2,1)+AM2*CC0
        BMAT(3,1)=CB*AMVAR(3,1)+AM3*CC0
        CC0=CC1*CMAT(1,2)+CC2*CMAT(2,2)+CC3*CMAT(3,2)
        BMAT(1,2)=CB*AMVAR(1,2)+AM1*CC0
        BMAT(2,2)=CB*AMVAR(2,2)+AM2*CC0
        BMAT(3,2)=CB*AMVAR(3,2)+AM3*CC0
        CC0=CC1*CMAT(1,3)+CC2*CMAT(2,3)+CC3*CMAT(3,3)
        BMAT(1,3)=CB*AMVAR(1,3)+AM1*CC0
        BMAT(2,3)=CB*AMVAR(2,3)+AM2*CC0
        BMAT(3,3)=CB*AMVAR(3,3)+AM3*CC0
C        
        CALL LINSOLVE(AMAT,BMAT)   
C
c computation of d(d/s)/dy_{0,j}
        CDSDERY1=CDSDER1*CMAT(1,1)+CDSDER2*CMAT(2,1)+CDSDER3*CMAT(3,1)
        CDSDERY2=CDSDER1*CMAT(1,2)+CDSDER2*CMAT(2,2)+CDSDER3*CMAT(3,2)
        CDSDERY3=CDSDER1*CMAT(1,3)+CDSDER2*CMAT(2,3)+CDSDER3*CMAT(3,3)
C computation of de_n/dy_0 and store result in AMAT
        CC1=AI1I+CDD/CSS
        CC2=AI2I+CDD/CSS
        CC3=AI3I+CDD/CSS        
        AMAT(1,1)=BMAT(1,1)*CC1+AZ1*CDSDERY1
        AMAT(1,2)=BMAT(1,2)*CC1+AZ1*CDSDERY2
        AMAT(1,3)=BMAT(1,3)*CC1+AZ1*CDSDERY3
        AMAT(2,1)=BMAT(2,1)*CC2+AZ2*CDSDERY1
        AMAT(2,2)=BMAT(2,2)*CC2+AZ2*CDSDERY2
        AMAT(2,3)=BMAT(2,3)*CC2+AZ2*CDSDERY3
        AMAT(3,1)=BMAT(3,1)*CC3+AZ3*CDSDERY1
        AMAT(3,2)=BMAT(3,2)*CC3+AZ3*CDSDERY2
        AMAT(3,3)=BMAT(3,3)*CC3+AZ3*CDSDERY3
c computation of ds/dy_{0,j}
        CSSDERY1=CSSDER1*CMAT(1,1)+CSSDER2*CMAT(2,1)+CSSDER3*CMAT(3,1)
        CSSDERY2=CSSDER1*CMAT(1,2)+CSSDER2*CMAT(2,2)+CSSDER3*CMAT(3,2)
        CSSDERY3=CSSDER1*CMAT(1,3)+CSSDER2*CMAT(2,3)+CSSDER3*CMAT(3,3)
c computation of -h/b^2 db/dy_{0,j}
        SCAL=-HB/CB**2
        CBDER1=SCAL*(CALPHA*CSSDERY1
     &        +2.0D0*CSS*(CM1*AMAT(1,1)+CM2*AMAT(2,1)+CM3*AMAT(3,1)))
        CBDER2=SCAL*(CALPHA*CSSDERY2
     &        +2.0D0*CSS*(CM1*AMAT(1,2)+CM2*AMAT(2,2)+CM3*AMAT(3,2)))
        CBDER3=SCAL*(CALPHA*CSSDERY3
     &        +2.0D0*CSS*(CM1*AMAT(1,3)+CM2*AMAT(2,3)+CM3*AMAT(3,3)))
c computation of dy_n/dy_0
        CBI=2.0D0/CB
        AMVARCS(1,1)=AMVARCS(1,1)+FAE1*(CBDER1*AZ2*AZ3
     &            +CBI*(BMAT(2,1)*AZ3+BMAT(3,1)*AZ2))
        AMVARCS(1,2)=AMVARCS(1,2)+FAE1*(CBDER2*AZ2*AZ3
     &            +CBI*(BMAT(2,2)*AZ3+BMAT(3,2)*AZ2))
        AMVARCS(1,3)=AMVARCS(1,3)+FAE1*(CBDER3*AZ2*AZ3
     &            +CBI*(BMAT(2,3)*AZ3+BMAT(3,3)*AZ2))
        AMVARCS(2,1)=AMVARCS(2,1)+FAE2*(CBDER1*AZ3*AZ1
     &            +CBI*(BMAT(3,1)*AZ1+BMAT(1,1)*AZ3))
        AMVARCS(2,2)=AMVARCS(2,2)+FAE2*(CBDER2*AZ3*AZ1
     &            +CBI*(BMAT(3,2)*AZ1+BMAT(1,2)*AZ3))
        AMVARCS(2,3)=AMVARCS(2,3)+FAE2*(CBDER3*AZ3*AZ1
     &            +CBI*(BMAT(3,3)*AZ1+BMAT(1,3)*AZ3))
        AMVARCS(3,1)=AMVARCS(3,1)+FAE3*(CBDER1*AZ1*AZ2
     &            +CBI*(BMAT(1,1)*AZ2+BMAT(2,1)*AZ1))
        AMVARCS(3,2)=AMVARCS(3,2)+FAE3*(CBDER2*AZ1*AZ2
     &            +CBI*(BMAT(1,2)*AZ2+BMAT(2,2)*AZ1))
        AMVARCS(3,3)=AMVARCS(3,3)+FAE3*(CBDER3*AZ1*AZ2
     &            +CBI*(BMAT(1,3)*AZ2+BMAT(2,3)*AZ1))
C COMPENSATED SUMMATION
        DO I=1,3
        DO J=1,3
        OLD=AMVAR(I,J)
        AMVAR(I,J)=OLD+AMVARCS(I,J)
        AMVARCS(I,J)=AMVARCS(I,J)+(OLD-AMVAR(I,J))
        END DO
        END DO
c UPDATE FOR QVAR     
c matrix -widehat e_n/(Id+widehat e_n) (*alpha)
        BMAT(1,1)=-CM3**2-CM2**2
        BMAT(1,2)=CM3+CM2*CM1
        BMAT(1,3)=CM3*CM1-CM2
        BMAT(2,1)=-CM3+CM2*CM1
        BMAT(2,2)=-CM3**2-CM1**2
        BMAT(2,3)=CM1+CM3*CM2
        BMAT(3,1)=CM3*CM1+CM2
        BMAT(3,2)=-CM1+CM3*CM2
        BMAT(3,3)=-CM2**2-CM1**2
c computation of (I-E) de_n/dy_0 
c matrix Id-widehat e_n
        CMAT(1,1)=1.0D0
        CMAT(1,2)=CM3
        CMAT(1,3)=-CM2
        CMAT(2,1)=-CM3
        CMAT(2,2)=1.0D0
        CMAT(2,3)=CM1
        CMAT(3,1)=CM2
        CMAT(3,2)=-CM1
        CMAT(3,3)=1.0D0
c computation of A_{n+1}
c compute QVARCS:=QVARCS+BMAT*QVAR+CMAT*AMAT
        CALPHAI=2.0D0/CALPHA
        DO I=1,3
        DO J=1,3
        DO K=1,3
        QVARCS(I,J)=QVARCS(I,J)
     &     +CALPHAI*(BMAT(I,K)*QVAR(K,J)+CMAT(I,K)*AMAT(K,J))
        END DO
        END DO
        END DO

c COMPENSATED SUMMATION
        DO I=1,3
        DO J=1,3
        OLD=QVAR(I,J)
        QVAR(I,J)=OLD+QVARCS(I,J)
        QVARCS(I,J)=QVARCS(I,J)+(OLD-QVAR(I,J))
        END DO
        END DO
C
c ------- DETECTION OF CONJUGATE POINTS
      DET=DETMAT(QVAR)
      if (ISTEP.GT.2) THEN
      IF (DET*DETOLD.LE.0.0D0) THEN
      ISCONJ=.TRUE.
c go backward      
      X=XOLD
      XCS=XCSOLD
      DO I=1,4
      Q(I)=QOLD(I)
      QCS(I)=QCSOLD(I)
      END DO
      DO I=1,3
      AM(I)=AMOLD(I) 
      AMCS(I)=AMCSOLD(I)  
      DO J=1,3
      AMVAR(I,J)=AMVAROLD(I,J)
      QVAR(I,J)=QVAROLD(I,J)
      AMVARCS(I,J)=AMVARCSOLD(I,J)
      QVARCS(I,J)=QVARCSOLD(I,J)
      END DO
      END DO
c divide stepsize
      HB=MIN(0.5D0,DETOLD/(DETOLD-DET))*HB
c      HB=0.5D0*HB
      DET=DETOLD
      END IF
      END IF
C save output
      DETOLD=DET
      XCSOLD=XCS
      DO I=1,4
      QOLD(I)=Q(I)
      QCSOLD(I)=QCS(I)
      END DO
      DO I=1,3
      AMOLD(I)=AM(I) 
      AMCSOLD(I)=AMCS(I) 
      DO J=1,3
      AMVAROLD(I,J)=AMVAR(I,J)
      QVAROLD(I,J)=QVAR(I,J)
      AMVARCSOLD(I,J)=AMVARCS(I,J)
      QVARCSOLD(I,J)=QVARCS(I,J)
      END DO
      END DO
c -----------------------------------
      END DO
      
C --- NORMALIZATION OF QUATERNION
C (norm of quaternion is not one to allow using compensated summation)
      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
      RETURN
      END
C
