MODULE Quaternion

!     see Kagan, Y. Y., 1991.
!     3-D rotation of double-couple earthquake sources,
!     Geophys. J. Int., 106, 709-716.

!     Peter Bird converted DCROT.FOR into Quaternion.f90,
!     in order to ease incorporation into Fortran90 programs,
!     and to provide a function that returns the minimum
!     rotation angle without activating any WRITE statements.
!     The original MAIN program was renamed Demonstration(),
!     and new DOUBLE PRECISION FUNCTION DCROT was created which
!     returns the minimum rotation angle in (positive) degrees.
!     There were no other (intentional) changes to the logic.
!     2003.02.17

!     Suggested uses of this module:
!
!     CALL Demonstration()
!
!           [or]
!
!     INTEGER*2 eq1(4), eq2(4)
!     REAL*8 minimum_angle
!
!     eq1(1) =  24 ! plunge  of T axis for earthquake #1
!     eq1(2) = 120 ! azimuth of T axis for earthquake #1
!     eq1(3) =  41 ! plunge  of P axis for earthquake #1
!     eq1(4) = 232 ! azimuth of P axis for earthquake #1

!     eq2(1) =  55 ! plunge  of T axis for earthquake #2
!     eq2(2) = 295 ! azimuth of T axis for earthquake #2
!     eq2(3) =  17 ! plunge  of P axis for earthquake #2
!     eq2(4) =  51 ! azimuth of P axis for earthquake #2
!
!     minimum_angle = DCROT ( eq1, eq2 ) ! (should be 102.8 degrees)
!===============================================================

      IMPLICIT REAL * 8 (a - h, o - z)
      LOGICAL :: verbose = .FALSE.
      COMMON / mom / rad, perp
      COMMON / pb2003 / verbose

CONTAINS

       SUBROUTINE Demonstration()

!      BEGIN OF THE DCROT PROGRAM

!      ROTATION ANGLE AND AXIS DIRECTION CALCULATIONS
!   FOR ALL FOUR POSSIBLE ROTATIONS OF DOUBLE-COUPLE SOURCE

!           DRIVER PROGRAM -- 6/12/1990 (February 15, 1991)
!    Corrected Aug. 21, 1994 (according to records Nov 18, 93)


      IMPLICIT REAL * 8 (a - h, o - z)

!   1985  2   16   16   28   14.90    39.690  142.690   39      66 264
!   1987  4    7    0   40   49.50    37.300  141.750   31      61 296
      INTEGER * 2 eqa1(4) / 66, 264, 22, 109 / , eqa2(4) / 61, 296, 29, 114 /

!   1985  7    2   12   34   56.10    40.630  143.900   24      57 300
!   1987  1    9    6   14   50.00    39.800  141.380   60      53  89
      INTEGER * 2 eqb1(4) / 57, 300, 33, 113 / , eqb2(4) / 53, 89, 35, 289 /

!   1977  1    6    6   11   50.50    -2.990  144.790   11      24 120
!   1980  9   26   15   20   42.50    -3.000  142.430   25      55 295
      INTEGER * 2 eqc1(4) / 24, 120, 41, 232 / , eqc2(4) / 55, 295, 17, 51 /
!   1986  6    5    9    1   20.30   -36.290  -97.120   15       0 226
!   1986  6   24   23   53   32.30   -36.550 -100.450   15       0 233
      INTEGER * 2 eqd1(4) / 0, 226, 0, 136 / , eqd2(4) / 0, 233, 0, 143 /

      INTEGER * 2 eqe1(4) / 0, 0, 0, 90 / , eqe2(4) / 90, 0, 0, 0 /

      LOGICAL verbose
      COMMON / pb2003 / verbose
      COMMON / mom / rad, perp

      verbose = .TRUE.

      hpi = DACOS(0.0D0)
      rad = 90.0D0 / hpi
      WRITE (6, 15)
      CALL Fps4r (eqa1, eqa2)
      WRITE (6, 15)
      CALL Fps4r (eqb1, eqb2)
      WRITE (6, 15)
      CALL Fps4r (eqc1, eqc2)
      WRITE (6, 15)
      CALL Fps4r (eqd1, eqd2)
      WRITE (6, 15)
      CALL Fps4r (eqe1, eqe2)
      WRITE (6, 15)
   15 FORMAT ('1')

      !STOP
      END SUBROUTINE Demonstration



      REAL*8 FUNCTION DCROT (eqh1, eqh2)

     !added by Peter Bird, 2003, as a silent function
     !based on the code of SUBROUTINE Fps4r.

      IMPLICIT REAL * 8 (a - h, o - z)
      REAL * 8 quat(4), quat1(4), quat2(4), quatr(4), quatr1(4), &
     & quatr2(4), quatt1(4), quatt2(4), quatc1(4)
      REAL * 8 q1a(4), q2a(4), q3a(4), q1b(4), q2b(4), q3b(4), q4(4)
      INTEGER * 2 eqh1(4), eqh2(4)
      LOGICAL matched_Ps, matched_Ts, verbose
      COMMON / pb2003 / verbose
      COMMON / mom / rad, perp

     !Initialize rad in COMMON mom:
      hpi = DACOS(0.0D0)
      rad = 90.0D0 / hpi

     !Initialize verbose in COMMON pb2003:
      verbose = .FALSE.

     !Check for case of two identical earthquakes
     !(which can cause numerical problems)
     !and return angle of 0:
      IF (eqh1(1) == 0) THEN ! horizontal T axis for eqh1
          matched_Ts = (eqh2(1) == 0).AND. &
                     & (MOD(ABS(eqh1(2) - eqh2(2)), 180) == 0)
      ELSE ! plunging T axis for eqh1
          matched_Ts = (eqh2(1) == eqh1(1)).AND. &
                     & (MOD(ABS(eqh1(2) - eqh2(2)), 360) == 0)
      END IF
      IF (eqh1(3) == 0) THEN ! horizontal P axis for eqh1
          matched_Ps = (eqh2(3) == 0).AND. &
                     & (MOD(ABS(eqh1(4) - eqh2(4)), 180) == 0)
      ELSE ! plunging P axis for eqh1
          matched_Ps = (eqh2(3) == eqh1(3)).AND. &
                     & (MOD(ABS(eqh1(4) - eqh2(4)), 360) == 0)
      END IF
      IF (matched_Ts.AND.matched_Ps) THEN
          DCROT = 0.0D0
          RETURN
      END IF

      CALL Quatfps (quat1, eqh1, 0)
      CALL Sphcoor (quat1, angl, theta, phi)
      icode = 0
      CALL Boxtest (quat1, quatr1, qm, icode)
      CALL Sphcoor (quatr1, angl, theta, phi)

      CALL Quatfps (quat2, eqh2, 0)
      CALL Sphcoor (quat2, angl, theta, phi)
      icode = 0
      CALL Boxtest (quat2, quatr2, qm, icode)
      CALL Sphcoor (quatr2, angl, theta, phi)

      angle1 = F4r1_pb (quat1, quat2, q1a, 1)

      angle2 = F4r1_pb (quat1, quat2, q2a, 2)

      angle3 = F4r1_pb (quat1, quat2, q3a, 3)

      angle4 = F4r1_pb (quat1, quat2, q4, 4)

      DCROT = MIN( angle1, angle2, angle3, angle4 )

      END FUNCTION DCROT



      SUBROUTINE Fps4r (eqh1, eqh2)

      IMPLICIT REAL * 8 (a - h, o - z)
      REAL * 8 quat(4), quat1(4), quat2(4), quatr(4), quatr1(4), &
     & quatr2(4), quatt1(4), quatt2(4), quatc1(4)
      REAL * 8 q1a(4), q2a(4), q3a(4), q1b(4), q2b(4), q3b(4), q4(4)
      INTEGER * 2 eqh1(4), eqh2(4)
      LOGICAL verbose
      COMMON / pb2003 / verbose
      COMMON / mom / rad, perp

      IF (verbose) WRITE (6, 20) eqh1, eqh2
      IF (verbose) WRITE (6, 10)

      CALL Quatfps (quat1, eqh1, 0)
      CALL Sphcoor (quat1, angl, theta, phi)
      IF (verbose) WRITE (6, 30) angl, theta, phi, perp
      icode = 0
      CALL Boxtest (quat1, quatr1, qm, icode)
      IF (verbose) WRITE (6, 40) quat1, quatr1
      IF (verbose) WRITE (6, 50) icode, qm
      CALL Sphcoor (quatr1, angl, theta, phi)
      IF (verbose) WRITE (6, 30) angl, theta, phi
      IF (verbose) WRITE (6, 10)

      CALL Quatfps (quat2, eqh2, 0)
      CALL Sphcoor (quat2, angl, theta, phi)
      IF (verbose) WRITE (6, 30) angl, theta, phi, perp
      icode = 0
      CALL Boxtest (quat2, quatr2, qm, icode)
      IF (verbose) WRITE (6, 40) quat2, quatr2
      IF (verbose) WRITE (6, 50) icode, qm
      CALL Sphcoor (quatr2, angl, theta, phi)
      IF (verbose) WRITE (6, 30) angl, theta, phi

      IF (verbose) WRITE (6, 10)
      CALL F4r1 (quat1, quat2, q1a, 1)
      CALL F4r2 (quat1, quat2, q1b, 1)
      IF (verbose) WRITE (6, 10)
      IF (verbose) WRITE (6, 40) q1a, q1b

      IF (verbose) WRITE (6, 10)
      CALL F4r1 (quat1, quat2, q2a, 2)
      CALL F4r2 (quat1, quat2, q2b, 2)
      IF (verbose) WRITE (6, 10)
      IF (verbose) WRITE (6, 40) q2a, q2b

      IF (verbose) WRITE (6, 10)
      CALL F4r1 (quat1, quat2, q3a, 3)
      CALL F4r2 (quat1, quat2, q3b, 3)
      IF (verbose) WRITE (6, 10)
      IF (verbose) WRITE (6, 40) q3a, q3b

      IF (verbose) WRITE (6, 10)
      CALL F4r1 (quat1, quat2, q4, 4)
      CALL F4r2 (quat1, quat2, q4, 4)

      IF (verbose) WRITE (6, 10)
      IF (verbose) WRITE (6, 40) q4
      IF (verbose) WRITE (6, 10)

   10 FORMAT (' ')
   20 FORMAT ('0 EQH1, EQH2 = ', 4I5, 4x, 4I5)
   30 FORMAT ('  ANGL, THETA, PHI = ', 7F14.7)
   40 FORMAT (' ', 9G14.6)
   50 FORMAT (' ICODE, QM = ', I5, 3F14.7)
      RETURN
      END SUBROUTINE Fps4r


      DOUBLE PRECISION FUNCTION F4r1_pb (q1, q2, q, icode)

!     Q = Q2*(Q1*(I,J,K,1))**(-1)

!     created by Peter Bird, 2003, to provide a silent
!     alternative to SUBROUTINE F4r1.
!     Instead of WRITEing the result ANGL,
!     it returns it as the value of F4r1_pb.

      IMPLICIT REAL * 8 (a - h, o - z)
      REAL * 8 q(4), q1(4), q2(4), qr1(4), &
     & qt1(4), qt2(4), qc1(4)

      CALL Boxtest (q1, qr1, qm, icode)

      CALL Quatd (qr1, q2, q)

      CALL Sphcoor (q, angl, theta, phi)

      F4r1_pb = angl

      END FUNCTION F4r1_pb



      SUBROUTINE F4r1 (q1, q2, q, icode)

!     Q = Q2*(Q1*(I,J,K,1))**(-1)

!     Since F4R1 and F4R2 yield the same results, only one subroutine
!     is needed; both programs are kept here for testing purposes.

      IMPLICIT REAL * 8 (a - h, o - z)
      REAL * 8 q(4), q1(4), q2(4), qr1(4), &
     & qt1(4), qt2(4), qc1(4)
      LOGICAL verbose
      COMMON / pb2003 / verbose

      CALL Boxtest (q1, qr1, qm, icode)
!      IF (verbose) WRITE (6, 20) QR1

      CALL Quatd (qr1, q2, q)
!      IF (verbose) WRITE (6, 20) Q1, Q2
!      CALL QUATP (Q1, Q, QT2)
!      IF (verbose) WRITE (6, 20) QT2, Q

      CALL Sphcoor (q, angl, theta, phi)

      IF (verbose) WRITE (6, 10) angl, theta, phi
!      IF (verbose) WRITE (6, 20) Q, QR, QM
   10 FORMAT ('  ANGL, THETA, PHI = ', 3F14.7)
   20 FORMAT (' ', 9G14.6)
      RETURN
      END SUBROUTINE F4r1


      SUBROUTINE F4r2 (q1, q2, q, icode)

!     Q = (Q2*(I,J,K,1))*Q1**(-1)

      IMPLICIT REAL * 8 (a - h, o - z)
      REAL * 8 q(4), q1(4), q2(4), qr2(4), &
     & qt1(4), qt2(4), qc1(4)
      LOGICAL verbose
      COMMON / pb2003 / verbose

      CALL Boxtest (q2, qr2, qm, icode)
!      IF (verbose) WRITE (6, 20) QR2

      CALL Quatd (q1, qr2, q)
!      IF (verbose) WRITE (6, 20) Q1, Q2
!      CALL QUATP (Q1, Q, QT2)
!      IF (verbose) WRITE (6, 20) QT2, Q

      CALL Sphcoor (q, angl, theta, phi)

      IF (verbose) WRITE (6, 10) angl, theta, phi
!      IF (verbose) WRITE (6, 20) Q, QR, QM
   10 FORMAT ('  ANGL, THETA, PHI = ', 3F14.7)
   20 FORMAT (' ', 9G14.6)
      RETURN
      END SUBROUTINE F4r2


      SUBROUTINE Boxtest (q1, q2, qm, icode)

!     for ICODE=0 finds minimal rotation quaternion
!     for ICODE=N finds rotation quaternion Q2 = Q1*(i,j,k,1),
!     for N=1,2,3,4

      IMPLICIT REAL * 8 (a - h, o - z)
      REAL * 8 q1(4), q2(4), quatt(4)
      REAL * 8 quat(4, 3) / 1.0D0, 0.0D0, 0.0D0, 0.0D0, &
     &                   0.0D0, 1.0D0, 0.0D0, 0.0D0, &
     &                   0.0D0, 0.0D0, 1.0D0, 0.0D0 /

      IF (icode /= 0) GO TO 15
      icode = 1
      qm = DABS(q1(1))
      DO 10 ixc = 2, 4
      IF (qm >= DABS(q1(ixc))) GO TO 10
      qm = DABS(q1(ixc))
      icode = ixc
   10 CONTINUE
   15 CONTINUE

      DO 20 ixc = 1, 4
      q2(ixc) = q1(ixc)
   20 CONTINUE

      IF (icode == 4) GO TO 40
      DO 30 ixc = 1, 4
      quatt(ixc) = quat(ixc, icode)
   30 CONTINUE
      CALL Quatp (quatt, q1, q2)
   40 CONTINUE

      IF (q2(4) > 0.0D0) GO TO 60
      DO 50 ixc = 1, 4
      q2(ixc) = - q2(ixc)
   50 CONTINUE
   60 CONTINUE
      qm = q2(4)

      RETURN
      END SUBROUTINE Boxtest


      SUBROUTINE Sphcoor (quat, angl, theta, phi)

!     for the rotation quaternion QUAT the subroutine finds the
!     rotation angle (ANGL) of a counterclockwise rotation and
!     spherical coordinates (colatitude THETA, and azimuth PHI) of the
!     rotation pole (intersection of the axis with reference sphere);
!     THETA=0 corresponds to the vector pointing down.

      IMPLICIT REAL * 8 (a - h, o - z)
      REAL * 8 quat(4)

      IF (quat(4) < 0.0D0) THEN
      DO 10 ism = 1, 4
      quat(ism) = - quat(ism)
   10 CONTINUE
      END IF

      q4n = DSQRT(1.0D0 - quat(4)**2)
      costh = 1.0
      IF (DABS(q4n) > 1.0D-10) costh = quat(3) / q4n
      IF (DABS(costh) > 1.0) costh = jidint(costh)
      theta = dacosd(costh)
      angl = 2.0D0 * dacosd(quat(4))
      phi = 0.0D0
      IF (DABS(quat(1)) > 1.0D-10.OR.DABS(quat(2)) > 1.0D-10) &
     & phi = datan2d(quat(2), quat(1))
      IF (phi < 0.0D0) phi = phi + 360.0D0

      RETURN
      END SUBROUTINE Sphcoor


      SUBROUTINE Quatfps (quat, eqh, icode)
      IMPLICIT REAL * 8 (a - h, o - z)
      REAL * 8 quat(4)
      INTEGER * 2 eqh(4)
      LOGICAL verbose
      COMMON / mom / rad, perp
      COMMON / pb2003 / verbose

!     THIS ROUTINE CALCULATES ROTATION QUATERNION CORRESPONDING TO
!                   EARTHQUAKE FOCAL MECHANISM

!     icode=0 -- four input data:  plunge and azimuth of T-axis
!                                  plunge and azimuth of P-axis
!     Since plunge and azimuth of 2 axes are redundant for calculation,
!     (four degrees of freedom vs three degrees that are necessary)
!     and have low accuracy (integer angular degrees), we calculate
!     plane_normal (V) and slip_vector (S) axes, in order that all axes
!     be orthogonal.

!     icode=1 -- three input data: slip angle (SA), dip angle (DA),
!                                  dip direction (DD)

!     PERP variable checks orthogonality
!     of T- and P-axes, it should be small (@<0.01 or so).

      ERR = 1.0D-15
      ic = 1
      IF (icode == 1) GO TO 200
      plg_t_ax = eqh(1)
      azm_t_ax = eqh(2)
      plg_p_ax = eqh(3)
      azm_p_ax = eqh(4)

      t1 = DCOS(azm_t_ax / rad) * DCOS(plg_t_ax / rad)
      t2 = DSIN(azm_t_ax / rad) * DCOS(plg_t_ax / rad)
      t3 = DSIN(plg_t_ax / rad)
      p1 = DCOS(azm_p_ax / rad) * DCOS(plg_p_ax / rad)
      p2 = DSIN(azm_p_ax / rad) * DCOS(plg_p_ax / rad)
      p3 = DSIN(plg_p_ax / rad)
      perp = t1 * p1 + t2 * p2 + t3 * p3
      IF (perp > 2.0D-02) THEN
      WRITE (6, 140) eqh, t1, t2, t3, p1, p2, p3, perp
  140 FORMAT (' **** T- AND P-AXES ARE NOT ORTHOGONAL: &
     & EQH, T1, T2, T3, P1, P2, P3, PERP = ', /, 4I4, 7G14.5)
      STOP 35
      END IF

      v1 = t1 + p1
      v2 = t2 + p2
      v3 = t3 + p3
      s1 = t1 - p1
      s2 = t2 - p2
      s3 = t3 - p3
      anormv = DSQRT(v1 * v1 + v2 * v2 + v3 * v3)
      v1 = v1 / anormv
      v2 = v2 / anormv
      v3 = v3 / anormv
      anorms = DSQRT(s1 * s1 + s2 * s2 + s3 * s3)
      s1 = s1 / anorms
      s2 = s2 / anorms
      s3 = s3 / anorms

      GO TO 250
  200 CONTINUE

      dd = eqh(1)
      da = eqh(2)
      sa = eqh(3)
      dd = dd / rad
      da = da / rad
      sa = sa / rad
      cdd = DCOS(dd)
      sdd = DSIN(dd)
      cda = DCOS(da)
      sda = DSIN(da)
      csa = DCOS(sa)
      ssa = DSIN(sa)
      s1 = csa * cdd + cda * sdd * ssa
      s2 = csa * sdd - ssa * cda * cdd
      s3 = - ssa * sda
      v1 = - sda * sdd
      v2 = sda * cdd
      v3 = - cda

!      S1 = CSA*SDD - SSA*CDA*CDD
!      S2 = - CSA*CDD - SSA*CDA*SDD
!      S3 = - SSA*SDA
!      V1 = SDA*CDD
!      V2 = SDA*SDD
!      V3 = - CDA

  250 CONTINUE
      an1 = s2 * v3 - v2 * s3
      an2 = v1 * s3 - s1 * v3
      an3 = s1 * v2 - v1 * s2
!      SINV3 = S1*V2*AN3 + S2*V3*AN1 + V1*AN2*S3 -
!     1 S3*V2*AN1 - S1*AN2*V3 - AN3*V1*S2
      d2 = 1.0D0 / DSQRT(2.0D0)
      t1 = (v1 + s1) * d2
      t2 = (v2 + s2) * d2
      t3 = (v3 + s3) * d2
      p1 = (v1 - s1) * d2
      p2 = (v2 - s2) * d2
      p3 = (v3 - s3) * d2
      IF (verbose) WRITE (6, 100) t1, t2, t3, p1, p2, p3, an1, an2, an3
  100 FORMAT (' T1, T2, T3, P1, P2, P3, AN1, AN2, AN3 = ', /, 9G13.4)

      u0 = (t1 + p2 + an3 + 1.0D0) / 4.0D0
      u1 = (t1 - p2 - an3 + 1.0D0) / 4.0D0
      u2 = (-t1 + p2 - an3 + 1.0D0) / 4.0D0
      u3 = (-t1 - p2 + an3 + 1.0D0) / 4.0D0
      um = DMAX1(u0, u1, u2, u3)
      IF (um == u0) GO TO 10
      IF (um == u1) GO TO 20
      IF (um == u2) GO TO 30
      IF (um == u3) GO TO 40
      IF (verbose) WRITE (6, 150)

   10 CONTINUE
      icod = 1 * ic
      u0 = DSQRT(u0)
      u3 = (t2 - p1) / (4.0D0 * u0)
      u2 = (an1 - t3) / (4.0D0 * u0)
      u1 = (p3 - an2) / (4.0D0 * u0)
      GO TO 50
   20 CONTINUE
      icod = 2 * ic
      u1 = DSQRT(u1)
      u2 = (t2 + p1) / (4.0D0 * u1)
      u3 = (t3 + an1) / (4.0D0 * u1)
      u0 = (p3 - an2) / (4.0D0 * u1)
      GO TO 50
   30 CONTINUE
      icod = 3 * ic
      u2 = DSQRT(u2)
      u1 = (t2 + p1) / (4.0D0 * u2)
      u0 = (an1 - t3) / (4.0D0 * u2)
      u3 = (p3 + an2) / (4.0D0 * u2)
      GO TO 50
   40 CONTINUE
      icod = 4 * ic
      u3 = DSQRT(u3)
      u0 = (t2 - p1) / (4.0D0 * u3)
      u1 = (t3 + an1) / (4.0D0 * u3)
      u2 = (p3 + an2) / (4.0D0 * u3)
   50 CONTINUE
      temp = u0 * u0 + u1 * u1 + u2 * u2 + u3 * u3

      IF (DABS(temp - 1.0D0) > ERR) THEN
      WRITE (6, 150)
  150 FORMAT (' **** ERROR *****')
      WRITE (6, 90) t1, t2, t3, p1, p2, p3
   90 FORMAT (' T1, T2, T3, P1, P2, P3 = ', /, 6G18.9)
      WRITE (6, 80) an1, an2, an3
   80 FORMAT (' AN1, AN2, AN3 = ', 3G18.9)
      WRITE (6, 120) temp, u1, u2, u3, u0
  120 FORMAT (' TEMP, U1, U2, U3, U0 = ', 5G18.9)
      END IF
      quat(1) = u1
      quat(2) = u2
      quat(3) = u3
      quat(4) = u0
      IF (verbose) WRITE (6, 130) quat, icod
  130 FORMAT (' QUAT = ', 4G18.9, ' ICOD =', I5)

      RETURN
      END SUBROUTINE Quatfps


      SUBROUTINE Quatp (q1, q2, q3)

!  Calculates product of two quaternions Q3 = Q2*Q1,
!  see F. Klein v.1 p.61, or Altmann, 1986, p.156,
!  or Biedenharn and Louck, 1981, p. 185.

!  Quaternion is taken here  --  q1*i + q2*j + q3*k + q4

      IMPLICIT REAL * 8 (a - h, o - z)

      REAL * 8 q1(4), q2(4), q3(4)

      q3(1) =  q1(4) * q2(1) + q1(3) * q2(2) - q1(2) * q2(3) + q1(1) * q2(4)
      q3(2) = -q1(3) * q2(1) + q1(4) * q2(2) + q1(1) * q2(3) + q1(2) * q2(4)
      q3(3) =  q1(2) * q2(1) - q1(1) * q2(2) + q1(4) * q2(3) + q1(3) * q2(4)
      q3(4) = -q1(1) * q2(1) - q1(2) * q2(2) - q1(3) * q2(3) + q1(4) * q2(4)
      RETURN
      END SUBROUTINE Quatp




      SUBROUTINE Quatd (q1, q2, q3)
      IMPLICIT REAL * 8 (a - h, o - z)

!     Quaternion division Q3 = Q2*(Q1)**(-1), or Q2 = Q3*Q1

      REAL * 8 q1(4), qc1(4), q2(4), q3(4)

      DO 10 i = 1, 3
      qc1(i) = - q1(i)
   10 CONTINUE
      qc1(4) = q1(4)
      CALL Quatp (qc1, q2, q3)

      RETURN
      END SUBROUTINE Quatd

!      END OF THE DCROT PROGRAM

END MODULE Quaternion