Euler Rotation - Fortran Code

Fortran Module

      MODULE C_MODULE_EULER_ROTATION

!-----------------------------------------------------------------------
! euler rotation of model geometry
!-----------------------------------------------------------------------
!
! C_ALPHA : rotation about z-axis
! C_BETA : rotation about new y-axis
! C_GAMMA : rotation about new z-axis
!
! C_FORWARD : rotation matricies to geographic coordinates
! C_REVERSE : rotation matricies to model coordinates
!
!-----------------------------------------------------------------------

      REAL (TYPE_REAL_8) ::
     &      C_ALPHA = 000.0

      REAL (TYPE_REAL_8) ::
     &      C_BETA   = 000.0

      REAL (TYPE_REAL_8) ::
     &      C_GAMMA = 000.0

      REAL (TYPE_REAL_8),
     &      DIMENSION (1:3, 1:3) ::
     &      C_FORWARD = 0.0

      REAL (TYPE_REAL_8),
     &      DIMENSION (1:3, 1:3) ::
     &      C_REVERSE = 0.0

!-----------------------------------------------------------------------

      END MODULE C_MODULE_EULER_ROTATION

Euler Matrix

SUBROUTINE C_EULER_ROTATION_MATRIX

!-----------------------------------------------------------------------
! purpose: compute euler rotation matrix
!-----------------------------------------------------------------------
!
! C_ALPHA   : rotation about z-axis
! C_BETA    : rotation about new y-axis
! C_GAMMA   : rotation about new z-axis
! C_REVERSE : rotation to model coordinates
! C_FORWARD : rotation to geographic coordinates
!
!-----------------------------------------------------------------------

      !local vars
      REAL (TYPE_REAL_8) ::
     &
     &      COSA,
     &      SINA,
     &      COSB,
     &      SINB,
     &      COSG,
     &      SING

!-----------------------------------------------------------------------

      !compute sin and cos of all angles in degrees
      COSA = C_COSD (C_ALPHA )
      SINA = C_SIND (C_ALPHA )
      COSB = C_COSD (C_BETA )
      SINB = C_SIND (C_BETA )
      COSG = C_COSD (C_GAMMA )
      SING = C_SIND (C_GAMMA )

      !compute rotation matrix into model coordinates
      C_REVERSE (1,1) =   COSA * COSB * COSG - SINA * SING
      C_REVERSE (1,2) =   SINA * COSB * COSG + COSA * SING
      C_REVERSE (1,3) =        - SINB * COSG
      C_REVERSE (2,1) = - COSA * COSB * SING - SINA * COSG
      C_REVERSE (2,2) = - SINA * COSB * SING + COSA * COSG
      C_REVERSE (2,3) =          SINB * SING
      C_REVERSE (3,1) =   COSA * SINB
      C_REVERSE (3,2) =   SINA * SINB
      C_REVERSE (3,3) =          COSB

      !compute rotation matrix into geographical coordinates
      C_FORWARD (1,1) =   COSG * COSB * COSA - SING * SINA
      C_FORWARD (1,2) = - SING * COSB * COSA - COSG * SINA
      C_FORWARD (1,3) =          SINB * COSA
      C_FORWARD (2,1) =   COSG * COSB * SINA + SING * COSA
      C_FORWARD (2,2) = - SING * COSB * SINA + COSG * COSA
      C_FORWARD (2,3) =          SINB * SINA
      C_FORWARD (3,1) = - COSG * SINB
      C_FORWARD (3,2) =   SING * SINB
      C_FORWARD (3,3) =          COSB

!-----------------------------------------------------------------------

      END SUBROUTINE C_EULER_ROTATION_MATRIX

Native-Grid -> AOMIP-Grid

     SUBROUTINE C_GEO_TO_MODEL (X_IN_ANGLE,
     &                          Y_IN_ANGLE,
     &                          X_OUT_ANGLE,
     &                          Y_OUT_ANGLE)

!--------------------------------------------------------------------
! purpose: a coordinate transformation from geographic coordinates
!          to model coordinates
!--------------------------------------------------------------------
!
! variables
!
! X_IN_ANGLE : X-COORDINATE IN GEOGRAPHICAL SYSTEM.
! Y_IN_ANGLE : Y-COORDINATE IN GEOGRAPHICAL SYSTEM.
! X_OUT_ANGLE : X-COORDINATE IN MODEL SYSTEM.
! Y_OUT_ANGLE : Y-COORDINATE IN MODEL SYSTEM.
!
!---------------------------------------------------------------------

      !passed arguments
      REAL (TYPE_REAL_8),
     &      INTENT (IN) ::
     &
     &      X_IN_ANGLE,
     &      Y_IN_ANGLE

      !passed arguments
      REAL (TYPE_REAL_8),
     &      INTENT (OUT) ::
     &
     &      X_OUT_ANGLE,
     &      Y_OUT_ANGLE

      !local vars
      REAL (TYPE_REAL_8) ::
     &
     &      X_TMP_ANGLE,
     &      Y_TMP_ANGLE,
     &      XX,
     &      YY,
     &      ZZ,
     &      XNEW,
     &      YNEW,
     &      ZNEW,
     &      COSTN,
     &      PHI_NEW,
     &      THETA_NEW

!---------------------------------------------------------------------

      !initial angles of rotation
      X_TMP_ANGLE = X_IN_ANGLE
      Y_TMP_ANGLE = Y_IN_ANGLE

      !avoid trouble of an exactly zero angle by adding offset
      X_TMP_ANGLE = X_TMP_ANGLE + C_EPSILON
      Y_TMP_ANGLE = Y_TMP_ANGLE + C_EPSILON

      !spherical to cartesian
      XX = C_COSD (Y_TMP_ANGLE) * C_COSD (X_TMP_ANGLE)
      YY = C_COSD (Y_TMP_ANGLE) * C_SIND (X_TMP_ANGLE)
      ZZ = C_SIND (Y_TMP_ANGLE)

      !new cartesian coordinates are given by
      XNEW = C_REVERSE (1,1) * XX
     &     + C_REVERSE (1,2) * YY
     &     + C_REVERSE (1,3) * ZZ
      YNEW = C_REVERSE (2,1) * XX
     &     + C_REVERSE (2,2) * YY
     &     + C_REVERSE (2,3) * ZZ
      ZNEW = C_REVERSE (3,1) * XX
     &     + C_REVERSE (3,2) * YY
     &     + C_REVERSE (3,3) * ZZ

      THETA_NEW = C_INV_SIND (ZNEW)
      COSTN = SQRT (1.0 - ZNEW ** 2)

      IF ((XNEW > 0.0).AND.(YNEW > 0.0)) THEN

         IF (XNEW < YNEW) THEN
            PHI_NEW= C_INV_COSD (XNEW/COSTN)
         ELSE
            PHI_NEW= C_INV_SIND (YNEW/COSTN)
         ENDIF

      ELSEIF ((XNEW < 0.0).AND.(YNEW > 0.0)) THEN

         IF (ABS(XNEW) < YNEW) THEN
            PHI_NEW = 180.0 - C_INV_COSD (ABS(XNEW)/COSTN)
         ELSE
            PHI_NEW = 180.0 - C_INV_SIND (YNEW/COSTN)
         ENDIF

      ELSEIF ((XNEW < 0.0).AND.(YNEW < 0.0)) THEN

         IF (ABS(XNEW) < ABS(YNEW)) THEN
            PHI_NEW =-180.0 + C_INV_COSD (ABS(XNEW)/COSTN)
         ELSE
            PHI_NEW =-180.0 + C_INV_SIND (ABS(YNEW)/COSTN)
         ENDIF

      ELSEIF ((XNEW > 0.0) .AND. (YNEW < 0.0)) THEN

         IF (    XNEW < ABS(YNEW)) THEN
            PHI_NEW =-C_INV_COSD (ABS(XNEW)/COSTN)
         ELSE
            PHI_NEW =-C_INV_SIND (ABS(YNEW)/COSTN)
         ENDIF

      ENDIF

      !new spherical coordinates are
      Y_OUT_ANGLE = THETA_NEW
      X_OUT_ANGLE = PHI_NEW

      IF (X_OUT_ANGLE < 0.0) X_OUT_ANGLE = X_OUT_ANGLE + 360.0

     !avoid trouble of an exactly zero angle by subtracting offset
      X_OUT_ANGLE = X_OUT_ANGLE - C_EPSILON
      Y_OUT_ANGLE = Y_OUT_ANGLE - C_EPSILON

!---------------------------------------------------------------------

      END SUBROUTINE C_GEO_TO_MODEL

AOMIP-Grid -> Native_Grid

     SUBROUTINE C_MODEL_TO_GEO (X_IN_ANGLE,
     &                          Y_IN_ANGLE,
     &                          X_OUT_ANGLE,
     &                          Y_OUT_ANGLE)

!---------------------------------------------------------------------
! purpose: a coordinate trnansformation from model coordinates
!          to geographic coordinates
!--------------------------------------------------------------------
!
! variables
!
! X_IN_ANGLE : X-COORDINATE IN MODEL SYSTEM.
! Y_IN_ANGLE : Y-COORDINATE IN MODEL SYSTEM.
! X_OUT_ANGLE : X-COORDINATE IN GEOGRAPHICAL SYSTEM.
! Y_OUT_ANGLE : Y-COORDINATE IN GEOGRAPHICAL SYSTEM.
!
! THE ROTATION FROM THE COORDINATE SYSTEM,
! X''',Y''',Z''' BACK TO THE ORIGINAL COORDINATE SYTEM
! CAN BE PERFORMED BY DOING THE ROTATION IN
! REVERSED ORDER AND WITH OPPOSITE SIGN ON THE ROTATED ANGLES.
!
!---------------------------------------------------------------------

      !passed arguments
      REAL (TYPE_REAL_8),
     &      INTENT (IN) ::
     &
     &      X_IN_ANGLE,
     &      Y_IN_ANGLE

      !passed arguments
      REAL (TYPE_REAL_8),
     &      INTENT (OUT) ::
     &
     &      X_OUT_ANGLE,
     &      Y_OUT_ANGLE

      !local vars
      REAL (TYPE_REAL_8) ::
     &
     &      X_TMP_ANGLE,
     &      Y_TMP_ANGLE,
     &      XX,
     &      YY,
     &      ZZ,
     &      XNEW,
     &      YNEW,
     &      ZNEW,
     &      COSTN,
     &      PHI_NEW,
     &      THETA_NEW

!---------------------------------------------------------------------

      !initial angles of rotation
      X_TMP_ANGLE = X_IN_ANGLE
      Y_TMP_ANGLE = Y_IN_ANGLE

      !avoid trouble of an exactly zero angle by adding offset
      X_TMP_ANGLE = X_TMP_ANGLE + C_EPSILON
      Y_TMP_ANGLE = Y_TMP_ANGLE + C_EPSILON

      !spherical coordiantes to cartesian coordinates
      XX = C_COSD (Y_TMP_ANGLE) * C_COSD (X_TMP_ANGLE)
      YY = C_COSD (Y_TMP_ANGLE) * C_SIND (X_TMP_ANGLE)
      ZZ = C_SIND (Y_TMP_ANGLE)

      !new cartesian coordinates are given by
      XNEW = C_FORWARD (1,1) * XX
     &     + C_FORWARD (1,2) * YY
     &     + C_FORWARD (1,3) * ZZ
      YNEW = C_FORWARD (2,1) * XX
     &     + C_FORWARD (2,2) * YY
     &     + C_FORWARD (2,3) * ZZ
      ZNEW = C_FORWARD (3,1) * XX
     &     + C_FORWARD (3,2) * YY
     &     + C_FORWARD (3,3) * ZZ

      !obtain new angles THETA_NEW,COSTN,PHI_NEW
      THETA_NEW = C_INV_SIND (ZNEW)
      COSTN = SQRT (1.0 - ZNEW**2)

      IF    ((XNEW > 0.0) .AND. (YNEW > 0.0)) THEN

         IF (XNEW < YNEW) THEN
            PHI_NEW = C_INV_COSD (XNEW/COSTN)
         ELSE
            PHI_NEW = C_INV_SIND (YNEW/COSTN)
         ENDIF

      ELSEIF ((XNEW < 0.0) .AND. (YNEW > 0.0)) THEN

         IF (ABS (XNEW) < YNEW) THEN
            PHI_NEW = 180.0 - C_INV_COSD (ABS (XNEW)/COSTN)
         ELSE
            PHI_NEW = 180.0 - C_INV_SIND (YNEW/COSTN)
         ENDIF

      ELSEIF ((XNEW < 0.0) .AND. (YNEW < 0.0)) THEN

         IF (ABS (XNEW) < ABS (YNEW)) THEN
            PHI_NEW =-180.0 + C_INV_COSD (ABS (XNEW)/COSTN)
         ELSE
            PHI_NEW =-180.0 + C_INV_SIND (ABS (YNEW)/COSTN)
         ENDIF

      ELSEIF ((XNEW > 0.0) .AND. (YNEW < 0.0)) THEN

        IF(    XNEW < ABS (YNEW)) THEN
           PHI_NEW =    - C_INV_COSD (ABS (XNEW)/COSTN)
        ELSE
           PHI_NEW =    - C_INV_SIND (ABS (YNEW)/COSTN)
        ENDIF

     ENDIF

     !new spherical coordinates
     X_OUT_ANGLE = PHI_NEW
     Y_OUT_ANGLE = THETA_NEW

     IF (X_OUT_ANGLE < 0.0) X_OUT_ANGLE = X_OUT_ANGLE + 360.0

     !avoid trouble of an exactly zero angle by subtracting offset
      X_OUT_ANGLE = X_OUT_ANGLE - C_EPSILON
      Y_OUT_ANGLE = Y_OUT_ANGLE - C_EPSILON

!-----------------------------------------------------------------------

      END SUBROUTINE C_MODEL_TO_GEO

Last updated: March 5, 2010