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, 139) 139 FORMAT (' WARNING from Quatfps in MODULE Quaternion:') 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) WRITE (6, *) !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