|
|
SUBROUTINE CONVERT_GEO_COORD(R_LAT_IN,R_LON_IN,R_HEIGHT_IN,
|
|
|
* R_LAT_OUT,R_LON_OUT,I_FLAG,I_ERROR)
|
|
|
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
|
c $Revision: 1.2 $
|
|
|
c The initial version of this subroutine was written by RADEX, INC.
|
|
|
c for use on VAX/VMS systems.
|
|
|
c
|
|
|
c Subsequent revisions for UNIX systems have been made by KBB at
|
|
|
c The Johns Hopkins Univ. Applied Physics Laboratory. These revisions
|
|
|
c have been managed using the Revision Control System (RCS) and a
|
|
|
c log of the revisions will be found at the end of the comments.
|
|
|
c
|
|
|
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
|
C
|
|
|
C
|
|
|
C Purpose:
|
|
|
C
|
|
|
C This subroutine uses a set of spherical harmonic coefficients to
|
|
|
C perform a coordinate conversion between geographic and corrected
|
|
|
C geomagnetic coordinates (or vice versa).
|
|
|
C
|
|
|
C The spherical harmonics for the geographic to corrected geomagnetic
|
|
|
C coordinate system correspond to a conversion from geographic to
|
|
|
C the corrected geomagnetic coordinates (expressed in terms of
|
|
|
C centered dipole coordinates at the input altitude), and are then
|
|
|
C transformed to ground level.
|
|
|
C
|
|
|
C The spherical harmonic coefficients used for the inverse are
|
|
|
C computed relative to centered dipole coordinates at the input
|
|
|
C altitude. The input CGM coordinates are converted to equivalent
|
|
|
C values in this coordinate sytem using the inverse altitude
|
|
|
C algorithm before the evaluation of the geographic spherical
|
|
|
C harmonic expansion.
|
|
|
C
|
|
|
C Method:
|
|
|
C
|
|
|
C This subroutine uses a five-step process in converting a position
|
|
|
C in one coordinate system to a position in the other.
|
|
|
C The five steps are as follows:
|
|
|
C
|
|
|
C 1. The appropriate spherical harmonic coefficients for the
|
|
|
C coordinate conversion are computed for the input altitude.
|
|
|
C
|
|
|
C 2. The appropriate coordinates for use in the spherical harmonic
|
|
|
C expansion are computed. For the geographic ==> corrected
|
|
|
C geomagnetic coordinate version, these are geographic colatitude
|
|
|
C and longitude. For the inverse coordinate conversion, the
|
|
|
C input coordinates are first converted into equivalent dipole
|
|
|
C coordinates at the input altitude.
|
|
|
C
|
|
|
C 3. The Cartesian coordinates of the unit vector in the desired
|
|
|
C coordinate system are then computed using the appropriate
|
|
|
C spherical harmonic expansion.
|
|
|
C
|
|
|
C 4. For the geographic ==> corrected geomagnetic coordinate
|
|
|
C conversion, the dipole-equivalent coordinates at input
|
|
|
C altitude are converted to their cgm cartesian coordinates.
|
|
|
C
|
|
|
C 5. Standard trigonometric identities are used to compute the
|
|
|
C latitude and longitude in the desired output coordinate system.
|
|
|
C
|
|
|
C Input Arguments:
|
|
|
C
|
|
|
C R_LAT_IN - REAL*4 - The input latitude in degrees.
|
|
|
C This could be either geographic latitude
|
|
|
C or corrected geomagnetic latitude. The
|
|
|
C acceptable range of values is (-90 to +90
|
|
|
C degrees).
|
|
|
C
|
|
|
C R_LON_IN - REAL*4 - The input longitude in degrees east.
|
|
|
C This could be either geographic longitude
|
|
|
C or corrected geomagnetic longitude. The
|
|
|
C acceptable range of values is (0 to 360 degrees).
|
|
|
C
|
|
|
C R_HEIGHT_IN - REAL*4 - The height in kilometers for which the
|
|
|
C coordinate transformation will be accomplished.
|
|
|
C The acceptable range of values is (0 km to
|
|
|
C 7200 km)
|
|
|
C
|
|
|
C I_FLAG - INTEGER - The flag that indicates which way
|
|
|
C the conversion will proceed.
|
|
|
C
|
|
|
C = 1 convert geographic to corrected
|
|
|
C geomagnetic coordinates
|
|
|
C
|
|
|
C = 2 convert corrected geomagnetic
|
|
|
C to geographic coordinates
|
|
|
C
|
|
|
C Output Arguments:
|
|
|
C
|
|
|
C R_LAT_OUT - REAL*4 - The output latitude in degrees.
|
|
|
C This could be either the geographic latitude
|
|
|
C or corrected geomagnetic latitude. This
|
|
|
C value will be between -90 and +90 degrees.
|
|
|
C
|
|
|
C R_LON_OUT - REAL*4 - The output longitude in degrees.
|
|
|
C This could be either geographic longitude or
|
|
|
C corrected geomagnetic longitude. This
|
|
|
C value will be between 0 and 360 degrees.
|
|
|
C
|
|
|
C I_ERROR - INTEGER - The error flag
|
|
|
C
|
|
|
C = 0 normal processing
|
|
|
C = -2 R_HEIGHT_IN is outside the allowable range
|
|
|
C = -4 I_FLAG value is invalid
|
|
|
C = -8 R_LAT_IN is outside the allowable range
|
|
|
C = -16 R_LON_IN is outside the allowable range
|
|
|
C = -32 Magnitude of the "unit vector" of the
|
|
|
C target coordinate system deviates
|
|
|
C significantly (+/- 20% or more) from 1.
|
|
|
C = -64 For altitudes > 0, for the corrected
|
|
|
C geomagnetic to geographic coordinate
|
|
|
C conversion there is a range of latitudes
|
|
|
C for which the transformation is invalid.
|
|
|
C This flag is set when the requested input
|
|
|
C cgm latitude falls within the invalid
|
|
|
C region.
|
|
|
C
|
|
|
C Local Variables:
|
|
|
C
|
|
|
C C_CLASS - Character String, used to identify message
|
|
|
C class in error messages.
|
|
|
C
|
|
|
C C_FLAG - Character String used to store I_FLAG
|
|
|
C as a character string for use in error
|
|
|
C message.
|
|
|
C
|
|
|
C C_HEIGHT_IN - Character String used to store R_HEIGHT_IN
|
|
|
C as a character string for use in error
|
|
|
C message.
|
|
|
C
|
|
|
C C_LAT_IN - Character String used to store R_LAT_IN
|
|
|
C as a character string for use in error
|
|
|
C message.
|
|
|
C
|
|
|
C C_LON_IN - Character String used to store R_LON_IN
|
|
|
C as a character string for use in error
|
|
|
C message.
|
|
|
C
|
|
|
C C_R - Character String used to store D_R
|
|
|
C as a character string for use in error
|
|
|
C message.
|
|
|
C
|
|
|
C C_REP_PARAMS - Character String used to store error
|
|
|
C message.
|
|
|
C
|
|
|
C
|
|
|
C C_ROUTINE_NAME - Character String containing name of
|
|
|
C subroutine generating error message.
|
|
|
C
|
|
|
C C_UNDEFINED - Character String used to store input data
|
|
|
C (R_HEIGHT, R_LAT_IN, R_LON_IN) as a
|
|
|
C character string for use in error message.
|
|
|
C
|
|
|
C D_CINT(,,) - Double Precision Array (3-D) -
|
|
|
C Contains the spherical harmonic coefficients
|
|
|
C interpolated to the input height, R_HEIGHT_IN.
|
|
|
C
|
|
|
C D_COEF(,,,) - Double Precision Array (3-D) -
|
|
|
C Contains the spherical harmonic coefficients
|
|
|
C used to compute the Cartesian coordinates of
|
|
|
C unit vector in the target coordinate system.
|
|
|
C
|
|
|
C First index: sp. harm. coeff. index
|
|
|
C Second x, y, z components of unit vector
|
|
|
C Third altitude indices
|
|
|
C Fourth direction of conversion index
|
|
|
C
|
|
|
C coefficients for a given altitude have the form
|
|
|
C
|
|
|
C a0 + h a1 + h^2 * a2 + h^3 * a3
|
|
|
C where h = alt/7000 [km]
|
|
|
C
|
|
|
C D_COLAT_INPUT - Double Precision - colatitude (radians) in the
|
|
|
C input coordinate system
|
|
|
C
|
|
|
C D_COLAT_OUTPUT - Double Precision - colatitude (radians) in the
|
|
|
C output coordinate system
|
|
|
C
|
|
|
C D_LON_INPUT - Double Precision - longitude (radians) in the
|
|
|
C input input system
|
|
|
C
|
|
|
C D_LON_OUTPUT - Double Precision - longitude (radians) in the
|
|
|
C output coordinate system
|
|
|
C
|
|
|
C D_R - Double Precision - magnitude of the
|
|
|
C unit radius vector in the target coordinate
|
|
|
C system. The target coordinate system is
|
|
|
C the system to which this subroutine is
|
|
|
C converting the latitude and longitude.
|
|
|
C
|
|
|
C D_X - Double Precision - the X-component of the
|
|
|
C unit radius vector in the target coordinate
|
|
|
C system
|
|
|
C
|
|
|
C D_Y - Double Precision - the Y-component of the
|
|
|
C unit radius vector in the target coordinate
|
|
|
C system
|
|
|
C
|
|
|
C D_Z - Double Precision - the Z-component of the
|
|
|
C unit radius vector in the target coordinate
|
|
|
C system
|
|
|
C
|
|
|
C D_YLMVAL - Double Precision - the array of spherical
|
|
|
C harmonic basis functions evaluated at
|
|
|
C a particular colatitude and longitude.
|
|
|
C
|
|
|
C
|
|
|
C I_COND_VALUE - Integer - the condition value returned from the
|
|
|
C call to SFC$$PUT_USER_MSG
|
|
|
C
|
|
|
C I_ERR64 - Logical - error flag for CGM_TO_ALTITUDE
|
|
|
C Routine.
|
|
|
C
|
|
|
C
|
|
|
C I_SYSTEM_ERR - Integer - Set to the value, SFC$$_NORMAL in
|
|
|
C all calls to SFC$$PUT_USER_MSG
|
|
|
C
|
|
|
C K - Integer - first index of the D_CINT
|
|
|
C interpolated spherical harmonic coefficient
|
|
|
C array.
|
|
|
C
|
|
|
C L - Integer - order of each spherical harmonic
|
|
|
C coefficient and spherical harmonic function.
|
|
|
C
|
|
|
C M - Integer - zonal index of the spherical
|
|
|
C harmonic functions.
|
|
|
C
|
|
|
C R_HEIGHT_OLD() - Real*4 Array (1-D) - Variable containing
|
|
|
C previous height (km) used to determine whether
|
|
|
C the interpolation to compute D_CINT() from
|
|
|
C D_COEF() needs to be done.
|
|
|
C
|
|
|
C R_LAT_ADJ - Real*4 Variable used to store result of
|
|
|
C ALTITUDE_TO_CGM or CGM_TO_ALTITUDE
|
|
|
C conversion.
|
|
|
C
|
|
|
C R_LAT_ALT - Real*4 Variable used to store at altitude
|
|
|
C dipole latitude.
|
|
|
C
|
|
|
C R_R - Real*4 - magnitude of unit vector (= D_R)
|
|
|
C (D_X, D_Y, D_Z) in single precision
|
|
|
C
|
|
|
C Constants:
|
|
|
C
|
|
|
C
|
|
|
C I_MAX_LENGTH - Integer - the length of the D_COEF(,,,)
|
|
|
C and D_YLMVAL spherical harmonic function
|
|
|
C arrays, corresponding to the order of the
|
|
|
C spherical harmonic expansion used.
|
|
|
C I_MAX_LENGTH = (I_ORDER + 1) * (I_ORDER + 1)
|
|
|
C
|
|
|
C I_MSGID - Integer - contains a message shell number of
|
|
|
C 4014 which will be used in all calls to
|
|
|
C SFC$$PUT_USER_MSG
|
|
|
C
|
|
|
C I_NUM_AXES - Integer - the number of axes in a Cartesian
|
|
|
C coordinate system (3)
|
|
|
C
|
|
|
C I_NUM_FLAG - Integer - the number of coordinate trans-
|
|
|
C formation flags available (2)
|
|
|
C
|
|
|
C I_NUM_LEVEL - Integer - the number of grid levels
|
|
|
C available. This is the number of
|
|
|
C distinct heights for which there are
|
|
|
C spherical harmonic coefficients available
|
|
|
C
|
|
|
C I_ORDER - Integer - the order of the spherical harmonic
|
|
|
C expansion used.
|
|
|
C
|
|
|
C DEGRAD - Double Precision - the multiplicative
|
|
|
C conversion factor for converting degrees
|
|
|
C to radians
|
|
|
C
|
|
|
C PI - Double Precision - mathematical constant
|
|
|
C
|
|
|
C Subroutines Required:
|
|
|
C
|
|
|
C RYLM - This subroutine returns the spherical harmonic
|
|
|
C function value array for a given colatitude and
|
|
|
C longitude. The spherical harmonics are returned
|
|
|
C in the D_YLMVAL array.
|
|
|
C
|
|
|
C CG_ALT_DIP - Given altitude (km), latitude, and a direction
|
|
|
C flag, returns altitude corrected latitude, and
|
|
|
C an error flag. See comments in subroutine for
|
|
|
C additional information.
|
|
|
C
|
|
|
C
|
|
|
C SFC$$PUT_USER_MSG - This subroutine is used to send error messages
|
|
|
C to the Message Handling System.
|
|
|
C
|
|
|
C Files Used at Compile Time:
|
|
|
C
|
|
|
C SFCLIB$$DEFS$$:SFC$$SFCDEF/LIST
|
|
|
C
|
|
|
C - This include file contains useful system
|
|
|
C constants and definitions. (eg. this is where
|
|
|
C SFC$$_NORMAL is defined)
|
|
|
C
|
|
|
C
|
|
|
C Files Used at Run Time: None
|
|
|
C
|
|
|
C
|
|
|
C Databases Accessed: None
|
|
|
C
|
|
|
C
|
|
|
C Warnings:
|
|
|
C
|
|
|
C 1. For certain values of Corrected Geomagnetic Coordinates, the
|
|
|
C transformation to Geographic Coordinates (Geocentric) is
|
|
|
C undefined. When this situation occurs, an error flag is set
|
|
|
C I_ERROR = -64 and returned to the calling routine.
|
|
|
C
|
|
|
C
|
|
|
C 2. This subroutine contains the non-standard include statement
|
|
|
C
|
|
|
C
|
|
|
C
|
|
|
C
|
|
|
C Revision History:
|
|
|
C
|
|
|
C Original Routine, based upon the Kyle Baker routine was converted
|
|
|
C to a standard library subroutine for SFC$$ by MFS on 14 OCT 1992.
|
|
|
C
|
|
|
C 09/17/93 SPR1993-0203 (SWS) The return condition value needs to
|
|
|
C be initialized before coordinate conversion processing begins
|
|
|
C
|
|
|
C 08/18/93 SPR1993-0175 (SWS) Routine performed a comparison of
|
|
|
C Z with a -1.D0. The variable Z should be D_Z.
|
|
|
C
|
|
|
C The present routine represents the use of an improved algorithm for
|
|
|
C the computation of corrected geomagnetic coordinates and (where
|
|
|
C it exists) the inverse, together with the use of an updated
|
|
|
C magnetic field model (IGRF 1995).
|
|
|
C
|
|
|
C The present routine incorporates the functionality of Kyle Baker
|
|
|
C routine, and, uses the same subroutine calls. Relevent portions of
|
|
|
C the old FORTRAN code have been retained where possible.
|
|
|
C
|
|
|
C The new algorithm and the current code have been developed by
|
|
|
C Radex, Inc, 3 Preston Court, Bedford, MA 01730. A technical
|
|
|
C report describing the improved algorithm, and discussing its
|
|
|
C limitations is in preparation, and will be published in early
|
|
|
C 1997 as a Philips Laboratory technical report entitled:.
|
|
|
C
|
|
|
C An Expanded Algorithm for the Computation of Corrected Geomagnetic
|
|
|
C Coordinates, by C. Hein and K. Bhavnani.
|
|
|
C
|
|
|
C The principal changes implementing the new algorithm and code are
|
|
|
C as follows:
|
|
|
C
|
|
|
C (1) The set of spherical harmonic function values are computed
|
|
|
C recursively using a single call to a new subroutine RYLM
|
|
|
C rather than using repetitive function call. The Kyle Baker
|
|
|
C routine calculated each of the spherical harmonics separately
|
|
|
C as needed.
|
|
|
C
|
|
|
C Computationally the recursive caluclation is much faster.
|
|
|
C
|
|
|
C (2) The spherical harmonic coefficients are based upon
|
|
|
C the IGRF 95 magnetic field model. The order of the
|
|
|
C spherical harmonic expansion used here is 10 (previously
|
|
|
C was 4).
|
|
|
C
|
|
|
C (3) The coefficients represent a cubic fit to the spherical
|
|
|
C harmonic coefficients over a range of 24 altitudes from 0 to
|
|
|
C 7200 km. The previous (7/17/95) version used a quadratic fit
|
|
|
C at 0, 300 and 1200 km altitude, and had an allowable
|
|
|
C range of altitudes is 0 - 2000 km. The original Kyle Baker
|
|
|
C routine used polynomial interpolations of the coefficients
|
|
|
C at 0, 150, 300 and 450 km.
|
|
|
C
|
|
|
C (4) The spherical harmonic coefficients are computed in an
|
|
|
C auxilliary coordinate system, which is aligned with the
|
|
|
C desired target coordinate system. This results in a
|
|
|
C considerable simplification of the code, eliminating the
|
|
|
C need to perform multiple coordinate system rotations
|
|
|
C as was the case in the Kyle Baker routine.
|
|
|
C
|
|
|
C 07/17/95
|
|
|
C
|
|
|
C MODIFIED FOR THE IGRF 95 MAGNETIC FIELD MODEL.
|
|
|
C
|
|
|
C CHANGES: BLOCKDATA SECTION HAS BEEN REPLACED WITH COEFFICIENTS
|
|
|
C APPROPRIATE FOR THE IGRF 1995 MODEL
|
|
|
C
|
|
|
C REFERENCES TO IGRF 90 IN THE COMMENT LINES HAVE BEEN
|
|
|
C CHANGED TO IGRF 95.
|
|
|
C
|
|
|
C 07/01/96
|
|
|
C
|
|
|
C The BLOCKDATA section has been expanded to accomodate spherical
|
|
|
C harmonic coefficients for a cubic interpolation, in order to
|
|
|
C accomodate the expanded 0 - 7200 km altitude range. The
|
|
|
C spherical harmonic coefficients were caluclated on the basis
|
|
|
C of the IGRF 95 magnetic field model.
|
|
|
C
|
|
|
C **********************************************************************
|
|
|
C
|
|
|
C The version provided here does not include the BLOCK DATA source
|
|
|
C code. Instead individual BLOCK DATA files are provided as follows
|
|
|
C for the following Magnetic Field Models:
|
|
|
C
|
|
|
C Field Model File Name
|
|
|
c
|
|
|
C DGRF 1975 BLKDAT75.FOR
|
|
|
C DGRF 1980 BLKDAT80.FOR
|
|
|
C DGRF 1985 BLKDAT85.FOR
|
|
|
C DGRF 1990 BLKDAT90.FOR
|
|
|
C IGRF 1995 BLKDAT95.FOR
|
|
|
C
|
|
|
C Line 3 of the BLKDATxx.FOR files contains an identification field
|
|
|
C for the magnetic field model used to generate the coefficients.
|
|
|
C
|
|
|
C **********************************************************************
|
|
|
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
|
c
|
|
|
c RCS revision log for POSIX systems
|
|
|
c
|
|
|
c $Log: not supported by cvs2svn $
|
|
|
c Revision 1.1 2004/07/08 18:01:56 brideout
|
|
|
c files added to support AACGM conversion
|
|
|
c
|
|
|
c Revision 1.8 1997/11/14 19:27:49 baker
|
|
|
c implementation of modifications to extend the coordinate
|
|
|
c system to 7200 km in altitude.
|
|
|
c
|
|
|
c Revision 1.7 1997/10/28 21:04:24 baker
|
|
|
c There was an error in the value of PI that has been fixed.
|
|
|
c
|
|
|
c Revision 1.6 1996/03/20 22:53:03 baker
|
|
|
c Explicitly declared the sizes of the arguments to the
|
|
|
c subroutine. This should help avoid compatibility
|
|
|
c problems when mixing C and Fortran with different
|
|
|
c machines and different compilers.
|
|
|
c
|
|
|
c Revision 1.5 1996/03/12 18:59:18 baker
|
|
|
c Added code to force a recomputation of the
|
|
|
c conversion coefficients at a given ghheight if the
|
|
|
c coordinates model has changed from the last call.
|
|
|
c
|
|
|
c Revision 1.4 1996/03/11 19:25:36 baker
|
|
|
c Modifications for the 1995 version of AACGM.
|
|
|
c
|
|
|
c
|
|
|
c Revision 1.3 94/10/17 12:35:32 12:35:32 baker (Kile Baker S1G)
|
|
|
cadded error code -64 to indicate invalid magnetic coordinates specified
|
|
|
c as input. This also requires a change in the call to cg_alt_dip
|
|
|
c
|
|
|
c Revision 1.2 94/10/14 10:53:36 10:53:36 baker (Kile Baker S1G)
|
|
|
c Added the SAVE instruction to make variables static
|
|
|
c
|
|
|
c Revision 1.1 94/10/12 15:28:38 15:28:38 baker (Kile Baker S1G)
|
|
|
c Initial revision
|
|
|
c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
|
|
|
C
|
|
|
C
|
|
|
C
|
|
|
C
|
|
|
C Set the order used in this algorithm to 10.
|
|
|
C
|
|
|
C
|
|
|
C Parameterize the maximum values of the indices of the
|
|
|
C spherical harmonic coefficient array D_COEF( , , , ) and the
|
|
|
C interpolated spherical harmonic coefficient array, D_CINT( , , ).
|
|
|
C
|
|
|
C
|
|
|
C
|
|
|
C
|
|
|
C
|
|
|
C
|
|
|
c
|
|
|
c The APL version of this code needs to keep track of
|
|
|
c the previous set of coefficients defining the conversion
|
|
|
c model (date). This is done by looking at the first coefficient
|
|
|
c in the model.
|
|
|
c
|
|
|
C
|
|
|
C CHARACTER*300 C_REP_PARAMS
|
|
|
C CHARACTER*8 C_LAT_IN
|
|
|
C CHARACTER*8 C_LON_IN
|
|
|
C CHARACTER*8 C_HEIGHT_IN
|
|
|
C CHARACTER*2 C_FLAG
|
|
|
C CHARACTER*6 C_R
|
|
|
C CHARACTER*70 C_UNDEFINED
|
|
|
C
|
|
|
C .. Parameters ..
|
|
|
DOUBLE PRECISION PI
|
|
|
PARAMETER (PI=3.141592653589793D0)
|
|
|
DOUBLE PRECISION DEGRAD
|
|
|
PARAMETER (DEGRAD=1.745329251994330D-2)
|
|
|
INTEGER I_ORDER
|
|
|
PARAMETER (I_ORDER=10)
|
|
|
INTEGER I_NUM_TERMS
|
|
|
PARAMETER (I_NUM_TERMS=121)
|
|
|
INTEGER I_NUM_AXES
|
|
|
PARAMETER (I_NUM_AXES=3)
|
|
|
INTEGER I_NUM_LEVEL
|
|
|
PARAMETER (I_NUM_LEVEL=5)
|
|
|
INTEGER I_NUM_FLAG
|
|
|
PARAMETER (I_NUM_FLAG=2)
|
|
|
C ..
|
|
|
C .. Scalar Arguments ..
|
|
|
DOUBLE PRECISION R_HEIGHT_IN,R_LAT_IN,R_LAT_OUT,R_LON_IN,R_LON_OUT
|
|
|
INTEGER I_ERROR,I_FLAG
|
|
|
C ..
|
|
|
C .. Local Scalars ..
|
|
|
DOUBLE PRECISION ALT_VAR,ALT_VAR_CU,ALT_VAR_QU,ALT_VAR_SQ,
|
|
|
* D_COLAT_INPUT,D_COLAT_OUTPUT,D_COLAT_TEMP,
|
|
|
* D_LON_INPUT,D_LON_OUTPUT,D_LON_TEMP,D_R,D_X,D_Y,
|
|
|
* D_Z,FIRST_COEFF_OLD,R_LAT_ADJ,R_LAT_ALT,R_R
|
|
|
INTEGER I,J,K,L,M
|
|
|
LOGICAL I_ERR64
|
|
|
C ..
|
|
|
C .. Local Arrays ..
|
|
|
DOUBLE PRECISION D_CINT(I_NUM_TERMS,I_NUM_AXES,I_NUM_FLAG),
|
|
|
* D_YLMVAL(I_NUM_TERMS),R_HEIGHT_OLD(2)
|
|
|
C ..
|
|
|
C .. External Subroutines ..
|
|
|
EXTERNAL ALTITUDE_TO_CGM,CGM_TO_ALTITUDE,RYLM
|
|
|
C ..
|
|
|
C .. Intrinsic Functions ..
|
|
|
INTRINSIC ABS,DACOS,DATAN2,DBLE,SQRT
|
|
|
C ..
|
|
|
C .. Common blocks ..
|
|
|
COMMON /SPH_HARM_MODEL/D_COEF
|
|
|
DOUBLE PRECISION D_COEF(I_NUM_TERMS,I_NUM_AXES,I_NUM_LEVEL,
|
|
|
* I_NUM_FLAG)
|
|
|
C ..
|
|
|
C .. Save statement ..
|
|
|
SAVE D_CINT,R_HEIGHT_OLD,FIRST_COEFF_OLD
|
|
|
C ..
|
|
|
C .. Data statements ..
|
|
|
C
|
|
|
C Initialize R_HEIGHT_OLD() to impossible values.
|
|
|
C
|
|
|
C
|
|
|
C
|
|
|
C
|
|
|
DATA R_HEIGHT_OLD(1)/-1.D0/
|
|
|
DATA R_HEIGHT_OLD(2)/-1.D0/
|
|
|
C ..
|
|
|
C **** DATA I_SYSTEM_ERR /SFC$$_NORMAL/
|
|
|
C
|
|
|
C
|
|
|
C
|
|
|
C The COMMON block / SPH_HARM_MODEL / array D_COEF contains
|
|
|
C the spherical harmonic coefficients used to generate the
|
|
|
C Cartesian coordinates of the unit vector in the target coordinate
|
|
|
C system. The target coordinate system is the system to which
|
|
|
C the algorithm is converting and is the coordinate system
|
|
|
C corresponding to the output variables, R_LAT_OUT and R_LON_OUT.
|
|
|
C
|
|
|
C The coefficient set is stored in a 4-dimensional array with
|
|
|
C indices defined as follows:
|
|
|
C
|
|
|
C First Index - Represents the number of terms used in the
|
|
|
C spherical harmonic expansion. Equal to
|
|
|
C (I_ORDER + 1) * (I_ORDER + 1)
|
|
|
C
|
|
|
C Second Index - Represents the Cartesian coordinate a particular
|
|
|
C coefficient will be used to generate. Indices are
|
|
|
C defined as follows:
|
|
|
C
|
|
|
C 1 - X-coordinate (the X-axis points from the
|
|
|
C center of the earth to where the prime
|
|
|
C meridian crosses the equator.
|
|
|
C
|
|
|
C 2 - Y-coordinate (the Y-axis points from the
|
|
|
C center of the earth to 90 degrees east
|
|
|
C of the X-axis)
|
|
|
C
|
|
|
C 3 - Z-coordinate (the Z-axis points from the
|
|
|
C center of the earth to the north pole)
|
|
|
C
|
|
|
C Third Index - Represents the terms of a quadratic fit to
|
|
|
C altitude (independent variable h =
|
|
|
C altitude [km]/1200) for a given spharical
|
|
|
C harmonic coefficient. The indices are defined
|
|
|
C as follows:
|
|
|
C
|
|
|
C 1 - Constant term: corresponds to the fit at
|
|
|
C 0 km altitude
|
|
|
C 2 - linear term
|
|
|
C 3 - quadratic term
|
|
|
C 4 - cubic term
|
|
|
C
|
|
|
C Fourth Index - Represents the direction of the coordinate
|
|
|
C transformation. Indices are defined as follows:
|
|
|
C
|
|
|
C 1 - Conversion of geographic coordinates to
|
|
|
C corrected geomagnetic coordinates.
|
|
|
C
|
|
|
C 2 - Conversion of corrected geomagnetic coordinates
|
|
|
C to geographic coordinates.
|
|
|
C
|
|
|
C
|
|
|
C The Data for D_COEF is provided in a BLOCK DATA file (see below)
|
|
|
C
|
|
|
C
|
|
|
C Initialize the error return condition indicator
|
|
|
C
|
|
|
I_ERROR = 0
|
|
|
C
|
|
|
C The APL version of this program allows the magnetic field
|
|
|
c model to change during the execution of the program.
|
|
|
c We therefore have to check to make sure that if the model
|
|
|
c has been changed that we recalculate the coordinate
|
|
|
c conversion coefficients.
|
|
|
C
|
|
|
C This IF statement checks to see if the magnetic
|
|
|
C coordinates model has been changed. If so, we
|
|
|
C have to force the recalculation of the conversion
|
|
|
C coefficients, by resetting the values of the old height
|
|
|
C
|
|
|
IF (FIRST_COEFF_OLD.NE.D_COEF(1,1,1,1)) THEN
|
|
|
R_HEIGHT_OLD(1) = -1
|
|
|
R_HEIGHT_OLD(2) = -1
|
|
|
END IF
|
|
|
FIRST_COEFF_OLD = D_COEF(1,1,1,1)
|
|
|
C
|
|
|
C
|
|
|
C The following IF-THEN-ELSE block checks the input arguments to
|
|
|
C ensure they are within allowable ranges. If any argument is
|
|
|
C outside the acceptable range, the error flag, I_ERROR will be
|
|
|
C set to some non-zero value, and control will be returned to the
|
|
|
C calling program.
|
|
|
C
|
|
|
IF ((R_HEIGHT_IN.LT.0.D0) .OR. (R_HEIGHT_IN.GT.7200.D0)) THEN
|
|
|
C
|
|
|
C The height in kilometers is outside the allowable range.
|
|
|
C
|
|
|
C Convert the height value R_HEIGHT_IN to a character string
|
|
|
C and send an error message to the user.
|
|
|
C
|
|
|
C WRITE(UNIT = C_HEIGHT_IN,
|
|
|
C $ FMT = '(F8.2)',
|
|
|
C $ IOSTAT = I_STATUS) R_HEIGHT_IN
|
|
|
C
|
|
|
C C_REP_PARAMS = C_ROUTINE_NAME//'\\HEIGHT\\'//C_HEIGHT_IN
|
|
|
C PRINT *, C_REP_PARAMS
|
|
|
C CALL PUT_USER_MSG(I_MSGID, C_CLASS, C_REP_PARAMS,
|
|
|
C $ I_SYSTEM_ERR, I_COND_VALUE)
|
|
|
C
|
|
|
C Set the error flag to the appropriate value and return
|
|
|
C
|
|
|
I_ERROR = -2
|
|
|
RETURN
|
|
|
C
|
|
|
ELSE IF ((I_FLAG.LT.1) .OR. (I_FLAG.GT.2)) THEN
|
|
|
C
|
|
|
C The conversion flag is neither 1 nor 2.
|
|
|
C
|
|
|
C Convert I_FLAG to a character string C_FLAG and send an error
|
|
|
C message to the user
|
|
|
C
|
|
|
C WRITE(UNIT = C_FLAG,
|
|
|
C $ FMT = '(I2)',
|
|
|
C $ IOSTAT = I_STATUS) I_FLAG
|
|
|
C
|
|
|
C C_REP_PARAMS = C_ROUTINE_NAME//'\\CONVERSION FLAG\\'//C_FLAG
|
|
|
C PRINT*,C_REP_PARAMS
|
|
|
C CALL PUT_USER_MSG(I_MSGID, C_CLASS, C_REP_PARAMS,
|
|
|
C $ I_SYSTEM_ERR, I_COND_VALUE)
|
|
|
C
|
|
|
C Set the error flag to the appropriate value and return
|
|
|
C
|
|
|
I_ERROR = -4
|
|
|
RETURN
|
|
|
C
|
|
|
ELSE IF (ABS(R_LAT_IN).GT.90.D0) THEN
|
|
|
C
|
|
|
C The latitude is outside the allowable range.
|
|
|
C
|
|
|
C Convert R_LAT_IN to a character string C_LAT_IN and send
|
|
|
C an error message to the user.
|
|
|
C
|
|
|
C WRITE(UNIT = C_LAT_IN,
|
|
|
C $ FMT = '(F8.2)',
|
|
|
C $ IOSTAT = I_STATUS) R_LAT_IN
|
|
|
C
|
|
|
C C_REP_PARAMS = C_ROUTINE_NAME//'\\LATITUDE\\'//C_LAT_IN
|
|
|
C PRINT*,C_REP_PARAMS
|
|
|
C CALL PUT_USER_MSG(I_MSGID, C_CLASS, C_REP_PARAMS,
|
|
|
C $ I_SYSTEM_ERR, I_COND_VALUE)
|
|
|
C
|
|
|
C Set the error flag to the appropriate value and return
|
|
|
C
|
|
|
I_ERROR = -8
|
|
|
RETURN
|
|
|
C
|
|
|
ELSE IF ((R_LON_IN.LT.0.D0) .OR. (R_LON_IN.GT.360.D0)) THEN
|
|
|
C
|
|
|
C The longitude is outside the allowable range.
|
|
|
C
|
|
|
C Convert R_LON_IN into a character string C_LON_IN and send
|
|
|
C an error message to the user.
|
|
|
C
|
|
|
C WRITE(UNIT = C_LON_IN,
|
|
|
C $ FMT = '(F8.2)',
|
|
|
C $ IOSTAT = I_STATUS) R_LON_IN
|
|
|
C
|
|
|
C C_REP_PARAMS = C_ROUTINE_NAME//'\\LONGITUDE\\'//C_LON_IN
|
|
|
C PRINT*,C_REP_PARAMS
|
|
|
C CALL PUT_USER_MSG(I_MSGID, C_CLASS, C_REP_PARAMS,
|
|
|
C $ I_SYSTEM_ERR, I_COND_VALUE)
|
|
|
C
|
|
|
C Set the error flag to the appropriate value and return
|
|
|
C
|
|
|
I_ERROR = -16
|
|
|
RETURN
|
|
|
C
|
|
|
END IF
|
|
|
C
|
|
|
C All input arguments are within allowable ranges.
|
|
|
C
|
|
|
C Compute Spherical Harmonic Coefficients for current
|
|
|
C altitude if required
|
|
|
C
|
|
|
IF (R_HEIGHT_IN.NE.R_HEIGHT_OLD(I_FLAG)) THEN
|
|
|
C
|
|
|
ALT_VAR = R_HEIGHT_IN/7200.0D0
|
|
|
ALT_VAR_SQ = ALT_VAR*ALT_VAR
|
|
|
ALT_VAR_CU = ALT_VAR*ALT_VAR_SQ
|
|
|
ALT_VAR_QU = ALT_VAR*ALT_VAR_CU
|
|
|
C
|
|
|
DO 20 I = 1,3
|
|
|
C
|
|
|
DO 10 J = 1,121
|
|
|
C
|
|
|
D_CINT(J,I,I_FLAG) = D_COEF(J,I,1,I_FLAG) +
|
|
|
* ALT_VAR*D_COEF(J,I,2,I_FLAG) +
|
|
|
* ALT_VAR_SQ*D_COEF(J,I,3,I_FLAG) +
|
|
|
* ALT_VAR_CU*D_COEF(J,I,4,I_FLAG) +
|
|
|
* ALT_VAR_QU*D_COEF(J,I,5,I_FLAG)
|
|
|
C
|
|
|
10 CONTINUE
|
|
|
C
|
|
|
20 CONTINUE
|
|
|
C
|
|
|
R_HEIGHT_OLD(I_FLAG) = R_HEIGHT_IN
|
|
|
C
|
|
|
END IF
|
|
|
C
|
|
|
C Zero Sums for Spherical Harmonic Expansion Computation
|
|
|
C
|
|
|
D_X = 0.0D0
|
|
|
D_Y = 0.0D0
|
|
|
D_Z = 0.0D0
|
|
|
C
|
|
|
C
|
|
|
C Prepare for Spherical Harmonic Expansion Computation
|
|
|
C
|
|
|
C
|
|
|
D_LON_INPUT = DBLE(R_LON_IN)*DEGRAD
|
|
|
C
|
|
|
IF (I_FLAG.EQ.1) THEN
|
|
|
C
|
|
|
C Computing CGM from Geographic Coordinates. No altitude
|
|
|
C correction required
|
|
|
C
|
|
|
D_COLAT_INPUT = (90.0D0-DBLE(R_LAT_IN))*DEGRAD
|
|
|
C
|
|
|
ELSE
|
|
|
C
|
|
|
C Computing Geographic Coordinates from CGM Coordinates.
|
|
|
C
|
|
|
C Convert CGM Latitude to Dipole Latitude at Altitude R_HEIGHT
|
|
|
C
|
|
|
CALL CGM_TO_ALTITUDE(R_HEIGHT_IN,R_LAT_IN,R_LAT_ADJ,I_ERR64)
|
|
|
C
|
|
|
IF (I_ERR64 .EQV. .TRUE.) THEN
|
|
|
C
|
|
|
C WRITE(UNIT = C_UNDEFINED,
|
|
|
C $ FMT = 90,
|
|
|
C $ IOSTAT = I_STATUS) R_LAT_IN, R_LON_IN, R_HEIGHT_IN
|
|
|
C
|
|
|
C C_REP_PARAMS = C_ROUTINE_NAME//C_UNDEFINED
|
|
|
C PRINT*,C_REP_PARAMS
|
|
|
C CALL PUT_USER_MSG(I_MSGID, C_CLASS, C_REP_PARAMS,
|
|
|
C $ I_SYSTEM_ERR, I_COND_VALUE)
|
|
|
C
|
|
|
I_ERROR = -64
|
|
|
RETURN
|
|
|
END IF
|
|
|
C
|
|
|
D_COLAT_INPUT = (90.0D0-DBLE(R_LAT_ADJ))*DEGRAD
|
|
|
C
|
|
|
END IF
|
|
|
C
|
|
|
C Generate the spherical harmonics at the coordinate point
|
|
|
C
|
|
|
CALL RYLM(D_COLAT_INPUT,D_LON_INPUT,I_ORDER,D_YLMVAL)
|
|
|
C
|
|
|
C Calculate the Cartesian coordinates of the unit vector
|
|
|
C in the target coordinate system.
|
|
|
C
|
|
|
DO 40 L = 0,I_ORDER
|
|
|
DO 30 M = -L,L
|
|
|
C
|
|
|
K = L*(L+1) + M + 1
|
|
|
C
|
|
|
C Add the contribution of each spherical harmonic to the
|
|
|
C Cartesian components of the unit vector in the
|
|
|
C appropriate coordinate system.
|
|
|
C
|
|
|
C
|
|
|
D_X = D_X + D_CINT(K,1,I_FLAG)*D_YLMVAL(K)
|
|
|
D_Y = D_Y + D_CINT(K,2,I_FLAG)*D_YLMVAL(K)
|
|
|
D_Z = D_Z + D_CINT(K,3,I_FLAG)*D_YLMVAL(K)
|
|
|
C
|
|
|
30 CONTINUE
|
|
|
40 CONTINUE
|
|
|
C
|
|
|
C Compute the magnitude of the Cartesian unit vector of the
|
|
|
C magnetic dipole coordinate system.
|
|
|
C
|
|
|
D_R = SQRT(D_X*D_X+D_Y*D_Y+D_Z*D_Z)
|
|
|
C
|
|
|
C If the magnitude of the unit vector differs significantly
|
|
|
C from 1, set the error flag and continue processing.
|
|
|
C
|
|
|
IF ((D_R.LT.0.9D0) .OR. (D_R.GT.1.1D0)) THEN
|
|
|
C
|
|
|
C The magnitude of the target unit vector differs significantly
|
|
|
C from 1. Convert the value D_R into the character string C_R,
|
|
|
C send a routine message to the user, and continue processing.
|
|
|
C
|
|
|
R_R = D_R
|
|
|
C
|
|
|
C WRITE(UNIT = C_R,
|
|
|
C $ FMT = '(F6.2)',
|
|
|
C $ IOSTAT = I_STATUS) R_R
|
|
|
C
|
|
|
C C_REP_PARAMS = C_ROUTINE_NAME//'\\UNIT VECTOR\\'//C_R
|
|
|
C PRINT*,C_REP_PARAMS
|
|
|
C CALL PUT_USER_MSG(I_MSGID, C_CLASS, C_REP_PARAMS,
|
|
|
C $ I_SYSTEM_ERR, I_COND_VALUE)
|
|
|
C
|
|
|
C Set the error flag to the appropriate value and return.
|
|
|
C
|
|
|
I_ERROR = -32
|
|
|
RETURN
|
|
|
END IF
|
|
|
C
|
|
|
C Adjust the components so they do represent the components of
|
|
|
C a unit vector. If D_R is equal to 1.0D0, this step will not
|
|
|
C change the values of D_X, D_Y, or D_Z.
|
|
|
C
|
|
|
D_Z = D_Z/D_R
|
|
|
D_X = D_X/D_R
|
|
|
D_Y = D_Y/D_R
|
|
|
C
|
|
|
C Obtain output co_latitude and longitude from the unit
|
|
|
C vector using standard formulas
|
|
|
C
|
|
|
IF (D_Z.GT.1.0D0) THEN
|
|
|
D_COLAT_TEMP = 0.0D0
|
|
|
ELSE IF (D_Z.LT.-1.0D0) THEN
|
|
|
D_COLAT_TEMP = PI
|
|
|
ELSE
|
|
|
D_COLAT_TEMP = DACOS(D_Z)
|
|
|
END IF
|
|
|
C
|
|
|
IF ((ABS(D_X).LT.1.0D-8) .AND. (ABS(D_Y).LT.1.0D-8)) THEN
|
|
|
D_LON_TEMP = 0.0D0
|
|
|
ELSE
|
|
|
D_LON_TEMP = DATAN2(D_Y,D_X)
|
|
|
END IF
|
|
|
C
|
|
|
C Prepare Latitude Data for Output
|
|
|
C
|
|
|
D_LON_OUTPUT = D_LON_TEMP
|
|
|
C
|
|
|
IF (I_FLAG.EQ.1) THEN
|
|
|
C
|
|
|
R_LAT_ALT = 90.0D0 - D_COLAT_TEMP/DEGRAD
|
|
|
C
|
|
|
CALL ALTITUDE_TO_CGM(R_HEIGHT_IN,R_LAT_ALT,R_LAT_ADJ)
|
|
|
C
|
|
|
D_COLAT_OUTPUT = (90.0D0-DBLE(R_LAT_ADJ))*DEGRAD
|
|
|
C
|
|
|
ELSE
|
|
|
C
|
|
|
D_COLAT_OUTPUT = D_COLAT_TEMP
|
|
|
C
|
|
|
END IF
|
|
|
C
|
|
|
C Convert colatitude into latitude and convert both coordinates
|
|
|
C from DOUBLE PRECISION radians to REAL degrees.
|
|
|
C
|
|
|
R_LAT_OUT = 90.0D0 - D_COLAT_OUTPUT/DEGRAD
|
|
|
R_LON_OUT = D_LON_OUTPUT/DEGRAD
|
|
|
C
|
|
|
IF (R_LON_OUT.LT.0.0D0) R_LON_OUT = R_LON_OUT + 360.D0
|
|
|
C
|
|
|
C
|
|
|
RETURN
|
|
|
9000 FORMAT (' INVERSE UNDEFINED AT CGM LAT ',F6.2,' CGM LON ',F7.2,
|
|
|
* ' ALTITUDE ',F8.2)
|
|
|
END
|
|
|
|
