hdiff output

r29980/nnutils.f90 2016-03-16 18:32:36.498474042 +0000 r29979/nnutils.f90 2016-03-16 18:32:36.710476222 +0000
947:         ! Lists of parallel and perpendicular quaternions for body J1 in all the images.947:         ! Lists of parallel and perpendicular quaternions for body J1 in all the images.
948:         ALLOCATE(QTN(NIMAGE+2,4))948:         ALLOCATE(QTN(NIMAGE+2,4))
949:         ALLOCATE(PTN(NIMAGE+2,4))949:         ALLOCATE(PTN(NIMAGE+2,4))
950: 950: 
951: !       CONVERT FROM AA -> QUATERNION FOR THE INITIAL AND FINAL FRAMES951: !       CONVERT FROM AA -> QUATERNION FOR THE INITIAL AND FINAL FRAMES
952:         DO K = 1, (NIMAGE+2), (NIMAGE+1)952:         DO K = 1, (NIMAGE+2), (NIMAGE+1)
953:            ! P is the angle-axis vector for rigid body J1 in the initial or final frame953:            ! P is the angle-axis vector for rigid body J1 in the initial or final frame
954:            J = (K-1)*NOPT + RBOFFSET + 3*(J1-1)954:            J = (K-1)*NOPT + RBOFFSET + 3*(J1-1)
955:            P(:)   = XYZ(J+1:J+3)955:            P(:)   = XYZ(J+1:J+3)
956:            THETA  = DSQRT(DOT_PRODUCT(P(:),P(:)))956:            THETA  = DSQRT(DOT_PRODUCT(P(:),P(:)))
957:            IF (ABS(THETA) .LT. 1.0D-8) THEN957:            IF (THETA == 0.D0) THEN
958:               QTN(K,1)   = 1.D0958:               QTN(K,1)   = 1.D0
959:               QTN(K,2:4) = 0.D0959:               QTN(K,2:4) = 0.D0
960:            ELSE960:            ELSE
961:               THETAH     = 0.5D0*THETA961:               THETAH     = 0.5D0*THETA
962:               ST         = SIN(THETAH)962:               ST         = SIN(THETAH)
963:               QTN(K,1)   = COS(THETAH)963:               QTN(K,1)   = COS(THETAH)
964:               QTN(K,2:4) = P(:)*ST/THETA964:               QTN(K,2:4) = P(:)*ST/THETA
965:            ENDIF965:            ENDIF
966:         ENDDO966:         ENDDO
967: 967: 
969:       ! A routine to obtain a series of quaternions which interpolate between two endpoint969:       ! A routine to obtain a series of quaternions which interpolate between two endpoint
970:       ! quaternions at regular intervals.970:       ! quaternions at regular intervals.
971:       ! Throughout, we use theta as a generic name for an angle. The geometric interpretation971:       ! Throughout, we use theta as a generic name for an angle. The geometric interpretation
972:       ! of this angle at any point is specified by \alpha, \beta etc.972:       ! of this angle at any point is specified by \alpha, \beta etc.
973: 973: 
974: !     NOW THETA = \ALPHA (WHICH IS THE ANGLE BETWEEN THE TWO END-QUATERNIONS)974: !     NOW THETA = \ALPHA (WHICH IS THE ANGLE BETWEEN THE TWO END-QUATERNIONS)
975:          ! cos(\alpha) is the dot product between the angle-axis vectors corresponding975:          ! cos(\alpha) is the dot product between the angle-axis vectors corresponding
976:          ! to rigid body J1 in the initial and final structures976:          ! to rigid body J1 in the initial and final structures
977:          ! (because the AA vectors are normalised)977:          ! (because the AA vectors are normalised)
978:          CT       = DOT_PRODUCT(QTN(1,:),QTN(NIMAGE+2,:))978:          CT       = DOT_PRODUCT(QTN(1,:),QTN(NIMAGE+2,:))
979:  
980:          IF (CT < 0.D0) THEN979:          IF (CT < 0.D0) THEN
981:            CT =-CT980:            CT =-CT
982:            QTN(NIMAGE+2,:) =-QTN(NIMAGE+2,:)981:            QTN(NIMAGE+2,:) =-QTN(NIMAGE+2,:)
983:          ENDIF982:          ENDIF
984: 983: 
985:         ! We may sometimes wish to interpolate around the longer of the two angles, for example to984:         ! We may sometimes wish to interpolate around the longer of the two angles, for example to
986:         ! avoid clash of two rigid bodies. Inverting the final quarternion causes iSLERP to985:         ! avoid clash of two rigid bodies. Inverting the final quarternion causes iSLERP to
987:         ! interpolate in the other direction, so this is what we do here.986:         ! interpolate in the other direction, so this is what we do here.
988:          IF (LONG_WAY) THEN987:          IF (LONG_WAY) THEN
989:            QTN(NIMAGE+2,:) =-QTN(NIMAGE+2,:)988:            QTN(NIMAGE+2,:) =-QTN(NIMAGE+2,:)
990:            CT =-CT  ! Inverting the quarternion obviously inverts the inner product989:            CT =-CT  ! Inverting the quarternion obviously inverts the inner product
991:          ENDIF990:          ENDIF
992:  991:  
993:          THETA    = ACOS(CT)992:         THETA    = ACOS(CT)
994:          ST       = SIN(THETA)993:          ST       = SIN(THETA)
995:  
996:          IF(ABS(THETA).LT.1.0D-10) THEN ! This happens with ALIGNRBST, for example 
997:             DO I = 2, NIMAGE+1 
998:                 J = NOPT*(I-1) + RBOFFSET + 3*(J1-1) 
999:                 XYZ(J+1:J+3) = P(:) 
1000:             ENDDO 
1001:             RETURN 
1002:          ENDIF 
1003: !     INCREMENTAL APPROACH: TANGENT QUATERNION994: !     INCREMENTAL APPROACH: TANGENT QUATERNION
1004:          ! The tangent quaternion PTN is a unit quaternion perpendicular to the995:          ! The tangent quaternion PTN is a unit quaternion perpendicular to the
1005:          ! interpolated quaternion996:          ! interpolated quaternion
1006:          PTN(1,:) = (QTN(NIMAGE+2,:) - CT*QTN(1,:))/ST997:          PTN(1,:) = (QTN(NIMAGE+2,:) - CT*QTN(1,:))/ST
1007:  
1008: !     NOW THETA = \BETA = \ALPHA/(NIMAGE+1)998: !     NOW THETA = \BETA = \ALPHA/(NIMAGE+1)
1009:          ! We redefine theta at this point.999:          ! We redefine theta at this point.
1010:          ! \beta is the separation in angle between the evenly-spaced consecutive images.1000:          ! \beta is the separation in angle between the evenly-spaced consecutive images.
1011:          THETA    = THETA/(NIMAGE+1)1001:          THETA    = THETA/(NIMAGE+1)
1012:          ST       = SIN(THETA)1002:          ST       = SIN(THETA)
1013:          CT       = COS(THETA)1003:          CT       = COS(THETA)
1014: 1004: 
1015:          DO I = 2, NIMAGE+11005:          DO I = 2, NIMAGE+1
1016:              ! For each image in turn, we construct the parallel and perpendicular1006:              ! For each image in turn, we construct the parallel and perpendicular
1017:              ! quaternions for this rigid body using the following (recursive) formula.1007:              ! quaternions for this rigid body using the following (recursive) formula.
1018:              QTN(I,:) = CT*QTN(I-1,:) + ST*PTN(I-1,:)1008:              QTN(I,:) = CT*QTN(I-1,:) + ST*PTN(I-1,:)
1019:              PTN(I,:) = CT*PTN(I-1,:) - ST*QTN(I-1,:)1009:              PTN(I,:) = CT*PTN(I-1,:) - ST*QTN(I-1,:)
1020: 1010: 
1021: !     CONVERT FROM QUATERNION -> AA1011: !     CONVERT FROM QUATERNION -> AA
1022:              THETA  = 2.D0*ACOS(QTN(I,1))1012:              THETA  = 2.D0*ACOS(QTN(I,1))
1023:  
1024:              J = NOPT*(I-1) + RBOFFSET + 3*(J1-1)1013:              J = NOPT*(I-1) + RBOFFSET + 3*(J1-1)
1025:              IF (ABS(THETA) .LT. 1.0D-8) THEN1014:              IF (THETA == 0.D0) THEN
1026:                  XYZ(J+1:J+3) = 0.D01015:                  XYZ(J+1:J+3) = 0.D0
1027:              ELSE1016:              ELSE
1028:                  FCT    = DSQRT(DOT_PRODUCT(QTN(I,2:4),QTN(I,2:4)))1017:                  FCT    = DSQRT(DOT_PRODUCT(QTN(I,2:4),QTN(I,2:4)))
1029:                  XYZ(J+1:J+3) = THETA*QTN(I,2:4)/FCT1018:                  XYZ(J+1:J+3) = THETA*QTN(I,2:4)/FCT
1030:              ENDIF1019:              ENDIF
1031: 1020: 
1032:          ENDDO1021:          ENDDO
1033: 1022: 
1034:          DEALLOCATE(QTN)1023:          DEALLOCATE(QTN)
1035:          DEALLOCATE(PTN)1024:          DEALLOCATE(PTN)


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0