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