Please note: You are viewing the unstyled version of this website. Either your browser does not support CSS (cascading style sheets) or it has been disabled. Skip navigation.

# Euler Rotation - Fortran Code

 Print PDF

## 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