!*************************************************************************************************** ! Copyright 2007,2009 S. Bourdarie ! ! This file is part of IRBEM-LIB. ! ! IRBEM-LIB is free software: you can redistribute it and/or modify ! it under the terms of the GNU Lesser General Public License as published by ! the Free Software Foundation, either version 3 of the License, or ! (at your option) any later version. ! ! IRBEM-LIB is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU Lesser General Public License for more details. ! ! You should have received a copy of the GNU Lesser General Public License ! along with IRBEM-LIB. If not, see . ! !*************************************************************************************************** ! CREATION: S. Bourdarie - ONERA-DESP ! ! FILE CONTENT: ! SUBROUTINE coord_trans1: Generic coordinate transformation from one Earth or Heliospheric coordinate ! system to another one ! SUBROUTINE coord_trans_vec1: Generic coordinate transformation from one Earth or Heliospheric coordinate ! system to another one (handle up to ! ntime_max positions) ! !*************************************************************************************************** C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE PNM80 C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009. Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C C Form the matrix of precession/nutation for a given date, IAU 1976 C precession model, IAU 1980 nutation model. C C The matrix operates in the sense V(date) = RMATPN * V(J2000), C where the p-vector V(date) is with respect to the true C equatorial triad of date DATE1+DATE2 and the p-vector C V(J2000) is with respect to the mean equatorial triad of C epoch J2000. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE PNM80 ( DATE1, DATE2, RMATPN ) IMPLICIT NONE REAL*8 DATE1, DATE2, RMATPN(3,3) REAL*8 RMATP(3,3), RMATN(3,3) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Precession matrix, J2000 to date. CALL PMAT76 ( DATE1, DATE2, RMATP ) * Nutation matrix. CALL NUTM80 ( DATE1, DATE2, RMATN ) * Combine the matrices: PN = N x P. CALL RXR ( RMATN, RMATP, RMATPN ) END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE PMAT76 C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C C Precession matrix from J2000 to a specified date, IAU 1976 model. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE PMAT76 ( DATE1, DATE2, RMATP ) IMPLICIT NONE DOUBLE PRECISION DATE1, DATE2, RMATP(3,3) * Reference epoch (J2000), JD REAL*8 DJ00 PARAMETER ( DJ00 = 2451545D0 ) REAL*8 ZETA, Z, THETA, WMAT(3,3) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Precession Euler angles, J2000 to specified date. CALL PREC76 ( DJ00, 0D0, DATE1, DATE2, ZETA, Z, THETA ) * Form the rotation matrix. CALL IR ( WMAT ) CALL RZ ( -ZETA, WMAT ) CALL RY ( THETA, WMAT ) CALL RZ ( -Z, WMAT ) CALL CR ( WMAT, RMATP ) END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE IR C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 and C INTEGER by INTEGER*4 by S. Bourdarie C C Initialize an r-matrix to the identity matrix. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE IR ( R ) IMPLICIT NONE REAL*8 R(3,3) INTEGER*4 I * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - CALL ZR ( R ) DO 1 I=1,3 R(I,I) = 1D0 1 CONTINUE END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE ZR C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 and C INTEGER by INTEGER*4 by S. Bourdarie C C Initialize an r-matrix to the null matrix. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE ZR ( R ) IMPLICIT NONE REAL*8 R(3,3) INTEGER*4 I, J * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - DO 2 J=1,3 DO 1 I=1,3 R(I,J) = 0.D0 1 CONTINUE 2 CONTINUE END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE RZ C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C CRotate an r-matrix about the z-axis. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE RZ ( PSI, R ) IMPLICIT NONE REAL*8 PSI, R(3,3) REAL*8 S, C, A(3,3), W(3,3) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Matrix representing new rotation. S = SIN(PSI) C = COS(PSI) CALL IR ( A ) A(1,1) = C A(2,1) = -S A(1,2) = S A(2,2) = C * Rotate. CALL RXR ( A, R, W ) * Return result. CALL CR ( W, R ) end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE RY C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C C Rotate an r-matrix about the y-axis. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE RY ( THETA, R ) IMPLICIT NONE REAL*8 THETA, R(3,3) REAL*8 S, C, A(3,3), W(3,3) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Matrix representing new rotation. S = SIN(THETA) C = COS(THETA) CALL IR ( A ) A(1,1) = C A(3,1) = S A(1,3) = -S A(3,3) = C * Rotate. CALL RXR ( A, R, W ) * Return result. CALL CR ( W, R ) end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE RY C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C C Rotate an r-matrix about the x-axis. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE RX ( PHI, R ) IMPLICIT NONE REAL*8 PHI, R(3,3) REAL*8 S, C, A(3,3), W(3,3) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Matrix representing new rotation. S = SIN(PHI) C = COS(PHI) CALL IR ( A ) A(2,2) = C A(3,2) = -S A(2,3) = S A(3,3) = C * Rotate. CALL RXR ( A, R, W ) * Return result. CALL CR ( W, R ) end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE RXR C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C C Multiply two r-matrices. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE RXR ( A, B, ATB ) IMPLICIT NONE REAL*8 A(3,3), B(3,3), ATB(3,3) INTEGER*4 I, J, K REAL*8 W, WM(3,3) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - DO 3 I=1,3 DO 2 J=1,3 W = 0D0 DO 1 K=1,3 W = W + A(I,K)*B(K,J) 1 CONTINUE WM(I,J) = W 2 CONTINUE 3 CONTINUE CALL CR ( WM, ATB ) end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE CR C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The changes to original routine are: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C Subroutine CR and CP have been combined together C C Copy an r-matrix. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE CR ( R, C ) IMPLICIT NONE REAL*8 R(3,3), C(3,3) INTEGER*4 I,J * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - DO 1 I=1,3 DO J=1,3 C(I,J)=R(I,J) ENDDO 1 CONTINUE end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE PREC76 C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE PREC76 ( EP01, EP02, EP11, EP12, ZETA, Z, THETA ) *+ *----------------------------------------------------------------------- IMPLICIT NONE REAL*8 EP01, EP02, EP11, EP12, ZETA, Z, THETA * Arcseconds to radians REAL*8 DAS2R PARAMETER ( DAS2R = 4.848136811095359935899141D-6 ) * Reference epoch (J2000), JD REAL*8 DJ00 PARAMETER ( DJ00 = 2451545D0 ) * Days per Julian century REAL*8 DJC PARAMETER ( DJC = 36525D0 ) REAL*8 T0, T, TAS2R, W * Interval between fundamental epoch J2000.0 and beginning epoch (JC). T0 = ( ( EP01-DJ00 ) + EP02 ) / DJC * Interval over which precession required (JC). T = ( ( EP11-EP01 ) + ( EP12-EP02 ) ) / DJC * Euler angles. TAS2R = T * DAS2R W = 2306.2181D0 + ( : 1.39656D0 : - 0.000139D0 * T0 ) * T0 ZETA = ( W + ( ( 0.30188D0 : - 0.000344D0 * T0 ) : + 0.017998D0 * T ) * T ) * TAS2R Z = ( W + ( ( 1.09468D0 : + 0.000066D0 * T0 ) : + 0.018203D0 * T ) * T ) * TAS2R THETA = ( ( 2004.3109D0 + ( : - 0.85330D0 : - 0.000217D0 * T0 ) * T0 ) + ( ( : - 0.42665D0 : - 0.000217D0 * T0 ) : - 0.041833D0 * T ) * T ) * TAS2R END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE NUTM80 C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C C Form the matrix of nutation for a given date, IAU 1980 model. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE NUTM80 ( DATE1, DATE2, RMATN ) IMPLICIT NONE REAL*8 DATE1, DATE2, RMATN(3,3) REAL*8 DPSI, DEPS, EPSA REAL*8 OBL80 * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Nutation components and mean obliquity. CALL NUT80 ( DATE1, DATE2, DPSI, DEPS ) EPSA = OBL80 ( DATE1, DATE2 ) * Build the rotation matrix. CALL NUMAT ( EPSA, DPSI, DEPS, RMATN ) end c C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE NUTM80 C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C C Nutation, IAU 1980 model. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE NUT80 ( DATE1, DATE2, DPSI, DEPS ) IMPLICIT NONE real*8 DATE1, DATE2, DPSI, DEPS * Arcseconds to radians real*8 DAS2R PARAMETER ( DAS2R = 4.848136811095359935899141D-6 ) * 2Pi real*8 D2PI PARAMETER ( D2PI = 6.283185307179586476925287D0 ) * Units of 0.1 milliarcsecond to radians real*8 U2R PARAMETER ( U2R = DAS2R/1D4 ) * Reference epoch (J2000), JD real*8 DJ00 PARAMETER ( DJ00 = 2451545D0 ) * Days per Julian century real*8 DJC PARAMETER ( DJC = 36525D0 ) real*8 T, EL, ELP, F, D, OM, DP, DE, ARG, S, C INTEGER*4 I, J real*8 ANPM * ------------------------------------------------ * Table of multiples of arguments and coefficients * ------------------------------------------------ * * The coefficient values are in 0.1 mas units and the rates of change * are in mas per Julian millennium. REAL X(9,106) * Multiple of Longitude Obliquity * L L' F D Omega coeff. of sin coeff. of cos * 1 t 1 t DATA ((X(I,J),I=1,9),J=1,10) / : 0., 0., 0., 0., 1., -171996., -1742., 92025., 89., : 0., 0., 0., 0., 2., 2062., 2., -895., 5., : -2., 0., 2., 0., 1., 46., 0., -24., 0., : 2., 0., -2., 0., 0., 11., 0., 0., 0., : -2., 0., 2., 0., 2., -3., 0., 1., 0., : 1., -1., 0., -1., 0., -3., 0., 0., 0., : 0., -2., 2., -2., 1., -2., 0., 1., 0., : 2., 0., -2., 0., 1., 1., 0., 0., 0., : 0., 0., 2., -2., 2., -13187., -16., 5736., -31., : 0., 1., 0., 0., 0., 1426., -34., 54., -1. / DATA ((X(I,J),I=1,9),J=11,20) / : 0., 1., 2., -2., 2., -517., 12., 224., -6., : 0., -1., 2., -2., 2., 217., -5., -95., 3., : 0., 0., 2., -2., 1., 129., 1., -70., 0., : 2., 0., 0., -2., 0., 48., 0., 1., 0., : 0., 0., 2., -2., 0., -22., 0., 0., 0., : 0., 2., 0., 0., 0., 17., -1., 0., 0., : 0., 1., 0., 0., 1., -15., 0., 9., 0., : 0., 2., 2., -2., 2., -16., 1., 7., 0., : 0., -1., 0., 0., 1., -12., 0., 6., 0., : -2., 0., 0., 2., 1., -6., 0., 3., 0. / DATA ((X(I,J),I=1,9),J=21,30) / : 0., -1., 2., -2., 1., -5., 0., 3., 0., : 2., 0., 0., -2., 1., 4., 0., -2., 0., : 0., 1., 2., -2., 1., 4., 0., -2., 0., : 1., 0., 0., -1., 0., -4., 0., 0., 0., : 2., 1., 0., -2., 0., 1., 0., 0., 0., : 0., 0., -2., 2., 1., 1., 0., 0., 0., : 0., 1., -2., 2., 0., -1., 0., 0., 0., : 0., 1., 0., 0., 2., 1., 0., 0., 0., : -1., 0., 0., 1., 1., 1., 0., 0., 0., : 0., 1., 2., -2., 0., -1., 0., 0., 0. / DATA ((X(I,J),I=1,9),J=31,40) / : 0., 0., 2., 0., 2., -2274., -2., 977., -5., : 1., 0., 0., 0., 0., 712., 1., -7., 0., : 0., 0., 2., 0., 1., -386., -4., 200., 0., : 1., 0., 2., 0., 2., -301., 0., 129., -1., : 1., 0., 0., -2., 0., -158., 0., -1., 0., : -1., 0., 2., 0., 2., 123., 0., -53., 0., : 0., 0., 0., 2., 0., 63., 0., -2., 0., : 1., 0., 0., 0., 1., 63., 1., -33., 0., : -1., 0., 0., 0., 1., -58., -1., 32., 0., : -1., 0., 2., 2., 2., -59., 0., 26., 0./ DATA ((X(I,J),I=1,9),J=41,50) / : 1., 0., 2., 0., 1., -51., 0., 27., 0., : 0., 0., 2., 2., 2., -38., 0., 16., 0., : 2., 0., 0., 0., 0., 29., 0., -1., 0., : 1., 0., 2., -2., 2., 29., 0., -12., 0., : 2., 0., 2., 0., 2., -31., 0., 13., 0., : 0., 0., 2., 0., 0., 26., 0., -1., 0., : -1., 0., 2., 0., 1., 21., 0., -10., 0., : -1., 0., 0., 2., 1., 16., 0., -8., 0., : 1., 0., 0., -2., 1., -13., 0., 7., 0., : -1., 0., 2., 2., 1., -10., 0., 5., 0. / DATA ((X(I,J),I=1,9),J=51,60) / : 1., 1., 0., -2., 0., -7., 0., 0., 0., : 0., 1., 2., 0., 2., 7., 0., -3., 0., : 0., -1., 2., 0., 2., -7., 0., 3., 0., : 1., 0., 2., 2., 2., -8., 0., 3., 0., : 1., 0., 0., 2., 0., 6., 0., 0., 0., : 2., 0., 2., -2., 2., 6., 0., -3., 0., : 0., 0., 0., 2., 1., -6., 0., 3., 0., : 0., 0., 2., 2., 1., -7., 0., 3., 0., : 1., 0., 2., -2., 1., 6., 0., -3., 0., : 0., 0., 0., -2., 1., -5., 0., 3., 0. / DATA ((X(I,J),I=1,9),J=61,70) / : 1., -1., 0., 0., 0., 5., 0., 0., 0., : 2., 0., 2., 0., 1., -5., 0., 3., 0., : 0., 1., 0., -2., 0., -4., 0., 0., 0., : 1., 0., -2., 0., 0., 4., 0., 0., 0., : 0., 0., 0., 1., 0., -4., 0., 0., 0., : 1., 1., 0., 0., 0., -3., 0., 0., 0., : 1., 0., 2., 0., 0., 3., 0., 0., 0., : 1., -1., 2., 0., 2., -3., 0., 1., 0., : -1., -1., 2., 2., 2., -3., 0., 1., 0., : -2., 0., 0., 0., 1., -2., 0., 1., 0. / DATA ((X(I,J),I=1,9),J=71,80) / : 3., 0., 2., 0., 2., -3., 0., 1., 0., : 0., -1., 2., 2., 2., -3., 0., 1., 0., : 1., 1., 2., 0., 2., 2., 0., -1., 0., : -1., 0., 2., -2., 1., -2., 0., 1., 0., : 2., 0., 0., 0., 1., 2., 0., -1., 0., : 1., 0., 0., 0., 2., -2., 0., 1., 0., : 3., 0., 0., 0., 0., 2., 0., 0., 0., : 0., 0., 2., 1., 2., 2., 0., -1., 0., : -1., 0., 0., 0., 2., 1., 0., -1., 0., : 1., 0., 0., -4., 0., -1., 0., 0., 0. / DATA ((X(I,J),I=1,9),J=81,90) / : -2., 0., 2., 2., 2., 1., 0., -1., 0., : -1., 0., 2., 4., 2., -2., 0., 1., 0., : 2., 0., 0., -4., 0., -1., 0., 0., 0., : 1., 1., 2., -2., 2., 1., 0., -1., 0., : 1., 0., 2., 2., 1., -1., 0., 1., 0., : -2., 0., 2., 4., 2., -1., 0., 1., 0., : -1., 0., 4., 0., 2., 1., 0., 0., 0., : 1., -1., 0., -2., 0., 1., 0., 0., 0., : 2., 0., 2., -2., 1., 1., 0., -1., 0., : 2., 0., 2., 2., 2., -1., 0., 0., 0. / DATA ((X(I,J),I=1,9),J=91,100) / : 1., 0., 0., 2., 1., -1., 0., 0., 0., : 0., 0., 4., -2., 2., 1., 0., 0., 0., : 3., 0., 2., -2., 2., 1., 0., 0., 0., : 1., 0., 2., -2., 0., -1., 0., 0., 0., : 0., 1., 2., 0., 1., 1., 0., 0., 0., : -1., -1., 0., 2., 1., 1., 0., 0., 0., : 0., 0., -2., 0., 1., -1., 0., 0., 0., : 0., 0., 2., -1., 2., -1., 0., 0., 0., : 0., 1., 0., 2., 0., -1., 0., 0., 0., : 1., 0., -2., -2., 0., -1., 0., 0., 0. / DATA ((X(I,J),I=1,9),J=101,106) / : 0., -1., 2., 0., 1., -1., 0., 0., 0., : 1., 1., 0., -2., 1., -1., 0., 0., 0., : 1., 0., -2., 2., 0., -1., 0., 0., 0., : 2., 0., 0., 2., 0., 1., 0., 0., 0., : 0., 0., 2., 4., 2., -1., 0., 0., 0., : 0., 1., 0., 1., 0., 1., 0., 0., 0. / * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Interval between fundamental epoch J2000.0 and given date (JC). T = ( ( DATE1-DJ00 ) + DATE2 ) / DJC * * FUNDAMENTAL ARGUMENTS in the FK5 reference system * * Mean longitude of the Moon minus mean longitude of the Moon's * perigee. EL = ANPM ( ( 485866.733D0 + ( 715922.633D0 + : ( 31.310D0 + 0.064D0 * T ) * T ) * T ) * DAS2R : + MOD(1325D0*T, 1D0) * D2PI ) * Mean longitude of the Sun minus mean longitude of the Sun's perigee. ELP = ANPM ( ( 1287099.804D0 + ( 1292581.224D0 + : ( -0.577D0 -0.012D0 * T ) * T ) * T ) * DAS2R : + MOD(99D0*T, 1D0) * D2PI ) * Mean longitude of the Moon minus mean longitude of the Moon's node. F = ANPM ( ( 335778.877D0 + ( 295263.137D0 + : ( -13.257D0 + 0.011D0 * T ) * T ) * T ) * DAS2R : + MOD(1342D0*T, 1D0) * D2PI ) * Mean elongation of the Moon from the Sun. D = ANPM ( ( 1072261.307D0 + ( 1105601.328D0 + : ( -6.891D0 + 0.019D0 * T ) * T ) * T ) * DAS2R : + MOD(1236D0*T, 1D0) * D2PI ) * Longitude of the mean ascending node of the lunar orbit on the * ecliptic, measured from the mean equinox of date. OM = ANPM( ( 450160.280D0 + ( -482890.539D0 + : ( 7.455D0 + 0.008D0 * T ) * T ) * T ) * DAS2R : + MOD( -5D0*T, 1D0) * D2PI ) * --------------- * Nutation series * --------------- * Change time argument from centuries to millennia. T = T / 10D0 * Initialize nutation components. DP = 0D0 DE = 0D0 * Sum the nutation terms, ending with the biggest. DO 1 J=106,1,-1 * Form argument for current term. ARG = DBLE(X(1,J)) * EL : + DBLE(X(2,J)) * ELP : + DBLE(X(3,J)) * F : + DBLE(X(4,J)) * D : + DBLE(X(5,J)) * OM * Accumulate current nutation term. S = DBLE(X(6,J)) + DBLE(X(7,J)) * T C = DBLE(X(8,J)) + DBLE(X(9,J)) * T IF ( S .NE. 0D0 ) DP = DP + S * SIN(ARG) IF ( C .NE. 0D0 ) DE = DE + C * COS(ARG) * Next term. 1 CONTINUE * Convert results from 0.1 mas units to radians. DPSI = DP * U2R DEPS = DE * U2R end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Function ANPM C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C C Normalize angle into the range -pi <= A < +pi. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ REAL*8 FUNCTION ANPM ( A ) IMPLICIT NONE REAL*8 A * Pi REAL*8 DPI PARAMETER ( DPI = 3.141592653589793238462643D0 ) * 2Pi REAL*8 D2PI PARAMETER ( D2PI = 6.283185307179586476925287D0 ) REAL*8 W * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - W = MOD(A,D2PI) IF ( ABS(W) .GE. DPI ) W = W - SIGN(D2PI,A) ANPM = W END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Function OBL80 C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C C Mean obliquity of the ecliptic, IAU 1980 model. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ REAL*8 FUNCTION OBL80 ( DATE1, DATE2 ) IMPLICIT NONE DOUBLE PRECISION DATE1, DATE2 * Arcseconds to radians REAL*8 DAS2R PARAMETER ( DAS2R = 4.848136811095359935899141D-6 ) * Reference epoch (J2000), JD REAL*8 DJ00 PARAMETER ( DJ00 = 2451545D0 ) * Days per Julian century REAL*8 DJC PARAMETER ( DJC = 36525D0 ) REAL*8 T * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Interval between fundamental epoch J2000.0 and given date (JC). T = ( ( DATE1-DJ00 ) + DATE2 ) / DJC * Mean obliquity of date. OBL80 = DAS2R * ( 84381.448D0 + : ( -46.8150D0 + : ( -0.00059D0 + : 0.001813D0 * T ) * T ) * T ) end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Subroutine NUMAT C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The only change to original routine is: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C C Form the matrix of nutation. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE NUMAT ( EPSA, DPSI, DEPS, RMATN ) IMPLICIT NONE REAL*8 EPSA, DPSI, DEPS, RMATN(3,3) * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Build the rotation matrix. CALL IR ( RMATN ) CALL RX ( EPSA, RMATN ) CALL RZ ( -DPSI, RMATN ) CALL RX ( -(EPSA+DEPS), RMATN ) end C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Function EQEQ94 C This routine has been extracted from the SOFA software C (Standards of Fundamental Astronomy, http://www.iau-sofa.rl.ac.uk/ ). C by S. Bourdarie in July 2009.Therefore it does not itself * constitute software provided by and/or endorsed by SOFA. C C The changes to original routine are: C DOUBLE PRECISION changed to REAL*8 by S. Bourdarie C Routine ANP has been included in this one by S. Bourdarie C C Equation of the equinoxes, IAU 1994 model. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ REAL*8 FUNCTION EQEQ94 ( DATE1, DATE2 ) IMPLICIT NONE REAL*8 DATE1, DATE2 * Arcseconds to radians REAL*8 DAS2R PARAMETER ( DAS2R = 4.848136811095359935899141D-6 ) * 2Pi REAL*8 D2PI PARAMETER (D2PI = 6.283185307179586476925287D0) * Reference epoch (J2000), JD REAL*8 DJ00 PARAMETER ( DJ00 = 2451545D0 ) * Days per Julian century REAL*8 DJC PARAMETER ( DJC = 36525D0 ) REAL*8 T, OM, DPSI, DEPS, EPS0 REAL*8 ANPM, OBL80 REAL*8 W * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * Interval between fundamental epoch J2000.0 and given date (JC). T = ( ( DATE1-DJ00 ) + DATE2 ) / DJC * Longitude of the mean ascending node of the lunar orbit on the * ecliptic, measured from the mean equinox of date. OM = ANPM( ( 450160.280D0 + ( -482890.539D0 + : ( 7.455D0 + 0.008D0 * T ) * T ) * T ) * DAS2R : + MOD(-5D0*T,1D0) * D2PI ) * Nutation components and mean obliquity. CALL NUT80 ( DATE1, DATE2, DPSI, DEPS ) EPS0 = OBL80 ( DATE1, DATE2 ) * Equation of the equinoxes. EQEQ94 = DPSI * COS(EPS0) + DAS2R * ( 0.00264D0 * SIN(OM) + : 0.000063D0 * SIN(OM+OM)) * Normalize angle into the range 0 <= A < 2pi. W = MOD(EQEQ94,D2PI) IF ( W .LT. 0.D0 ) W = W + D2PI EQEQ94 = W end !--------------------------------------------------------------------------------------------------- ! Introduced in version 4.... ! ! CREATION: S. Bourdarie - July 2009 ! MODIFICATION: None ! ! DESCRIPTION: Compute the dates following the J2000 method. ! ! INPUT: iyr -> year (long integer) ! idoy -> day of year (long integer) ! secs -> UT in seconds (double) ! ! OUTPUT: datE1, date2 -> julday ! ! CALLING SEQUENCE: call date2dateJ2000(iyr,idoy,UT,date1,date2) !--------------------------------------------------------------------------------------------------- SUBROUTINE date2dateJ2000(iyr,idoy,UT,date1,date2) IMPLICIT NONE c c input variables INTEGER*4 iyr,idoy REAL*8 UT c c internal variables REAL*8 DJ00 PARAMETER ( DJ00 = 2451545D0 ) INTEGER*4 Month, Day,hour, minute, second,julday,jd REAL*8 eps,eps2 c c Output variables REAL*8 DATE1,DATE2 c ------------------------------------------------------------------- c Compute the dates following J2000 method DATE1=DJ00 ! julday for 1st january 2000 at 12 UT (J2000 method) call DOY_AND_UT2DATE_AND_TIME(iyr,idoy,UT, Month, Day, & hour, minute, second) jd=JULDAY(iyr,month,day) eps=2.2204460d-16 eps2 = eps*ABS(jd) if (eps2 .gt. eps) eps=eps2 DATE2 = dble(jd) + ( (Hour/24.d0 - 0.5d0) + & (UT-hour*3600.d0)/86400.d0 + eps )-DATE1 !(J2000 method) c end !--------------------------------------------------------------------------------------------------- ! Introduced in version 4.... ! ! CREATION: S. Bourdarie - July 2009 ! MODIFICATION: None ! ! DESCRIPTION: Rotates coordinates from ECI J2000 to ECI true of date ! ! INPUT: iyr -> year (long integer) ! idoy -> day of year (long integer) ! secs -> UT in seconds (double) ! xIN -> position in ECI J2000 system (double array(3)) ! ! OUTPUT: xOUT -> position in ECI True of Date (TOD) system (double array(3)) ! ! CALLING SEQUENCE: call ECIJ2000_2_ECITOD(iyr,idoy,UT,xIN,xOUT) !--------------------------------------------------------------------------------------------------- SUBROUTINE ECIJ2000_2_ECITOD(iyr,idoy,UT,xIN,xOUT) c IMPLICIT NONE c c input variables INTEGER*4 iyr,idoy REAL*8 UT,xIN(3) c c internal variables REAL*8 RMATPN(3,3) REAL*8 DATE1,DATE2 c Output variables REAL*8 xOUT(3) c Compute the dates following J2000 method call date2dateJ2000(iyr,idoy,UT,date1,date2) c Form the matrix of precession/nutation for a given date, IAU 1976 c precession model, IAU 1980 nutation model. call PNM80 ( DATE1, DATE2, RMATPN) c Rotate input coordinates to outputs xOUT(1)=RMATPN(1,1)*xIN(1)+RMATPN(1,2)*xIN(2)+RMATPN(1,3)*xIN(3) xOUT(2)=RMATPN(2,1)*xIN(1)+RMATPN(2,2)*xIN(2)+RMATPN(2,3)*xIN(3) xOUT(3)=RMATPN(3,1)*xIN(1)+RMATPN(3,2)*xIN(2)+RMATPN(3,3)*xIN(3) end !--------------------------------------------------------------------------------------------------- ! Introduced in version 4.... ! ! CREATION: S. Bourdarie - July 2009 ! MODIFICATION: None ! ! DESCRIPTION: Rotates coordinates from ECI true of date to ECI J2000 ! ! INPUT: iyr -> year (long integer) ! idoy -> day of year (long integer) ! secs -> UT in seconds (double) ! xIN -> position in ECI True of Date (TOD) system (double array(3)) ! ! OUTPUT: xOUT -> position in ECI J2000 system (double array(3)) ! ! CALLING SEQUENCE: call ECITOD_2_ECIJ2000(iyr,idoy,UT,xIN,xOUT) !--------------------------------------------------------------------------------------------------- SUBROUTINE ECITOD_2_ECIJ2000(iyr,idoy,UT,xIN,xOUT) c IMPLICIT NONE c c input variables INTEGER*4 iyr,idoy REAL*8 UT,xIN(3) c c internal variables REAL*8 RMATPN(3,3),invRMATPN(3,3),det REAL*8 DATE1,DATE2 c c Output variables REAL*8 xOUT(3) c Compute the dates following J2000 method call date2dateJ2000(iyr,idoy,UT,date1,date2) c Form the matrix of precession/nutation for a given date, IAU 1976 c precession model, IAU 1980 nutation model. call PNM80 ( DATE1, DATE2, RMATPN) c c invert RMATPN matrix det=RMATPN(1,1)*(RMATPN(2,2)*RMATPN(3,3)-RMATPN(2,3)*RMATPN(3,2)) & +RMATPN(1,2)*(RMATPN(2,1)*RMATPN(3,3)-RMATPN(2,3)*RMATPN(3,1)) & +RMATPN(1,3)*(RMATPN(2,1)*RMATPN(3,2)-RMATPN(2,2)*RMATPN(3,1)) invRMATPN(1,1)=(RMATPN(2,2)*RMATPN(3,3)-RMATPN(2,3)*RMATPN(3,2)) & /det invRMATPN(1,2)=(RMATPN(1,3)*RMATPN(3,2)-RMATPN(1,2)*RMATPN(3,3)) & /det invRMATPN(1,3)=(RMATPN(1,2)*RMATPN(2,3)-RMATPN(1,3)*RMATPN(2,2)) & /det invRMATPN(2,1)=(RMATPN(2,3)*RMATPN(3,1)-RMATPN(2,1)*RMATPN(3,3)) & /det invRMATPN(2,2)=(RMATPN(1,1)*RMATPN(3,3)-RMATPN(1,3)*RMATPN(3,1)) & /det invRMATPN(2,3)=(RMATPN(1,3)*RMATPN(2,1)-RMATPN(1,1)*RMATPN(2,3)) & /det invRMATPN(3,1)=(RMATPN(2,1)*RMATPN(3,2)-RMATPN(2,2)*RMATPN(3,1)) & /det invRMATPN(3,2)=(RMATPN(1,2)*RMATPN(3,1)-RMATPN(1,1)*RMATPN(3,2)) & /det invRMATPN(3,3)=(RMATPN(1,1)*RMATPN(2,2)-RMATPN(1,2)*RMATPN(2,1)) & /det c c Rotate input coordinates to outputs xOUT(1)=invRMATPN(1,1)*xIN(1)+invRMATPN(1,2)*xIN(2)+ & invRMATPN(1,3)*xIN(3) xOUT(2)=invRMATPN(2,1)*xIN(1)+invRMATPN(2,2)*xIN(2)+ & invRMATPN(2,3)*xIN(3) xOUT(3)=invRMATPN(3,1)*xIN(1)+invRMATPN(3,2)*xIN(2)+ & invRMATPN(3,3)*xIN(3) end !--------------------------------------------------------------------------------------------------- ! Introduced in version 4.... ! ! CREATION: S. Bourdarie - July 2009 ! MODIFICATION: None ! ! DESCRIPTION: Rotates coordinates from TEME to ECI true of date ! TEME=True Equator Mean Equinox ! ! INPUT: iyr -> year (long integer) ! idoy -> day of year (long integer) ! secs -> UT in seconds (double) ! xIN -> position in TEME system (double array(3)) ! ! OUTPUT: xOUT -> position in ECI TOD system (double array(3)) ! ! CALLING SEQUENCE: call TEME2ECITOD(iyr,idoy,UT,xIN,xOUT) !--------------------------------------------------------------------------------------------------- SUBROUTINE TEME2ECITOD(iyr,idoy,UT,xIN,xOUT) c IMPLICIT NONE c c input variables INTEGER*4 iyr,idoy REAL*8 UT,xIN(3) c internal variables REAL*8 RMATPN(3,3),invRMATPN(3,3),det REAL*8 DATE1,DATE2 REAL*8 EQEQ94,psi,S,C,A(3,3) c c Output variables REAL*8 xOUT(3) c c Compute the dates following J2000 method call date2dateJ2000(iyr,idoy,UT,date1,date2) !c Form the rotation matrix around z axis by angle -Eqe82 psi=-EQEQ94 ( DATE1, DATE2 ) * Matrix representing new rotation. S = SIN(PSI) C = COS(PSI) CALL IR ( A ) A(1,1) = C A(2,1) = -S A(1,2) = S A(2,2) = C xOUT(1)=A(1,1)*xIN(1)+A(1,2)*xIN(2)+A(1,3)*xIN(3) xOUT(2)=A(2,1)*xIN(1)+A(2,2)*xIN(2)+A(2,3)*xIN(3) xOUT(3)=A(3,1)*xIN(1)+A(3,2)*xIN(2)+A(3,3)*xIN(3) end c !--------------------------------------------------------------------------------------------------- ! Introduced in version 4.... ! ! CREATION: S. Bourdarie - July 2009 ! MODIFICATION: None ! ! DESCRIPTION: Rotates coordinates from ECI true of date to TEME ! TEME=True Equator Mean Equinox ! ! INPUT: iyr -> year (long integer) ! idoy -> day of year (long integer) ! secs -> UT in seconds (double) ! xIN -> position in ECI True of Date (TOD) system (double array(3)) ! ! OUTPUT: xOUT -> position in TEME system (double array(3)) ! ! CALLING SEQUENCE: call ECITOD2TEME(iyr,idoy,UT,xIN,xOUT) !--------------------------------------------------------------------------------------------------- SUBROUTINE ECITOD2TEME(iyr,idoy,UT,xIN,xOUT) c IMPLICIT NONE c c input variables INTEGER*4 iyr,idoy REAL*8 UT,xIN(3) c internal variables REAL*8 RMATPN(3,3),invRMATPN(3,3),det REAL*8 DATE1,DATE2 REAL*8 EQEQ94,psi,S,C,A(3,3) c c Output variables REAL*8 xOUT(3) c c Compute the dates following J2000 method call date2dateJ2000(iyr,idoy,UT,date1,date2) !c Form the rotation matrix around z axis by angle Eqe82 psi=EQEQ94 ( DATE1, DATE2 ) * Matrix representing new rotation. S = SIN(PSI) C = COS(PSI) CALL IR ( A ) A(1,1) = C A(2,1) = -S A(1,2) = S A(2,2) = C xOUT(1)=A(1,1)*xIN(1)+A(1,2)*xIN(2)+A(1,3)*xIN(3) xOUT(2)=A(2,1)*xIN(1)+A(2,2)*xIN(2)+A(2,3)*xIN(3) xOUT(3)=A(3,1)*xIN(1)+A(3,2)*xIN(2)+A(3,3)*xIN(3) end c c--------------------------------------------------------------------------------------------------- ! Introduced in version 4.0 ! ! CREATION: S. Bourdarie - January 2007 ! MODIFICATION: None ! ! DESCRIPTION: Coordinate transformation from one Earth or Heliospheric coordinate ! system to another one ! ! INPUT: sysaxesIN -> designed input coordinate system (long integer) ! sysaxesOUT> designed output coordinate system (long integer) ! iyr -> year (long integer) ! idoy -> day of year (long integer) ! secs -> UT in seconds (double) ! xIN -> position in input coordinate system (double array(3)) ! ! OUTPUT: xOUT -> position in output coordinate system (double array(3)) ! ! CALLING SEQUENCE: call coord_trans1(sysaxesIN,sysaxesOUT,iyr,idoy,secs,xIN,xOUT) !--------------------------------------------------------------------------------------------------- c SUBROUTINE coord_trans1(sysaxesIN,sysaxesOUT,iyr,idoy, & secs,xIN,xOUT) c IMPLICIT NONE c INTEGER*4 sysaxesIN,sysaxesOUT,iyr,idoy INTEGER*4 i REAL*8 secs,psi REAL*8 xIN(3),xOUT(3),xTMP(3),alti call initize ! sets rad, pi used by various routines if (sysaxesIN.EQ.sysaxesOUT) then c write(6,*)'sysaxesIN = sysaxesOUT!' do i=1,3 xOUT(i)=xIN(i) enddo return endif if (sysaxesIN.LT.0) then write(6,*)'sysaxesIN out of range !' do i=1,3 xOUT(i)=-1.D31 enddo return endif if (sysaxesIN.GT.13) then write(6,*)'sysaxesIN out of range !' do i=1,3 xOUT(i)=-1.D31 enddo return endif if (sysaxesOUT.LT.0) then write(6,*)'sysaxesOUT out of range !' do i=1,3 xOUT(i)=-1.D31 enddo return endif if (sysaxesOUT.GT.13) then write(6,*)'sysaxesOUT out of range !' do i=1,3 xOUT(i)=-1.D31 enddo return endif c c input=GDZ if (sysaxesIN.EQ.0) then if (sysaxesOUT.EQ.1) then !GEO call gdz_geo(xIN(2),xIN(3),xIN(1),xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.2) then !GSM call gdz_geo(xIN(2),xIN(3),xIN(1),xTMP(1),xTMP(2),xTMP(3)) call geo2gsm1(iyr,idoy,secs,psi,xTMP,xOUT) endif if (sysaxesOUT.EQ.3) then !GSE call gdz_geo(xIN(2),xIN(3),xIN(1),xTMP(1),xTMP(2),xTMP(3)) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.4) then !SM call gdz_geo(xIN(2),xIN(3),xIN(1),xTMP(1),xTMP(2),xTMP(3)) call geo2sm1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.5) then !GEI call gdz_geo(xIN(2),xIN(3),xIN(1),xTMP(1),xTMP(2),xTMP(3)) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.6) then !MAG call gdz_geo(xIN(2),xIN(3),xIN(1),xTMP(1),xTMP(2),xTMP(3)) call geo2mag1(iyr,xTMP,xOUT) endif if (sysaxesOUT.EQ.7) then !SPH call gdz_geo(xIN(2),xIN(3),xIN(1),xTMP(1),xTMP(2),xTMP(3)) call CAR_SPH(xTMP,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.8) then !RLL call gdz_geo(xIN(2),xIN(3),xIN(1),xTMP(1),xTMP(2),xTMP(3)) xOUT(1)=SQRT(xTMP(1)*xTMP(1)+xTMP(2)*xTMP(2)+ & xTMP(3)*xTMP(3)) xOUT(2)=xIN(2) xOUT(3)=xIN(3) endif if (sysaxesOUT.EQ.9) then !HEE call gdz_geo(xIN(2),xIN(3),xIN(1),xTMP(1),xTMP(2),xTMP(3)) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.10) then !HAE call gdz_geo(xIN(2),xIN(3),xIN(1),xTMP(1),xTMP(2),xTMP(3)) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.11) then !HEEQ call gdz_geo(xIN(2),xIN(3),xIN(1),xTMP(1),xTMP(2),xTMP(3)) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hae2heeq1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI ToD call gdz_geo(xIN(2),xIN(3),xIN(1),xTMP(1),xTMP(2),xTMP(3)) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.13) then !J2000 call gdz_geo(xIN(2),xIN(3),xIN(1),xTMP(1),xTMP(2),xTMP(3)) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo CALL ECITOD_2_ECIJ2000(iyr,idoy,secs,xTMP,xOUT) endif endif c c input=GEO if (sysaxesIN.EQ.1) then if (sysaxesOUT.EQ.0) then !GDZ call geo_gdz(xIN(1),xIN(2),xIN(3),xOUT(2),xOUT(3),xOUT(1)) endif if (sysaxesOUT.EQ.2) then !GSM call geo2gsm1(iyr,idoy,secs,psi,xIN,xOUT) endif if (sysaxesOUT.EQ.3) then !GSE call geo2gse1(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.4) then !SM call geo2sm1(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.5) then !GEI call geo2gei1(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.6) then !MAG call geo2mag1(iyr,xIN,xOUT) endif if (sysaxesOUT.EQ.7) then !SPH call CAR_SPH(xIN,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.8) then !RLL call geo_gdz(xIN(1),xIN(2),xIN(3),xTMP(1),xTMP(2),xTMP(3)) xOUT(1)=SQRT(xIN(1)*xIN(1)+xIN(2)*xIN(2)+xIN(3)*xIN(3)) xOUT(2)=xTMP(1) xOUT(3)=xTMP(2) endif if (sysaxesOUT.EQ.9) then !HEE call geo2gse1(iyr,idoy,secs,xIN,xTMP) call gse2hee1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.10) then !HAE call geo2gse1(iyr,idoy,secs,xIN,xTMP) call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.11) then !HEEQ call geo2gse1(iyr,idoy,secs,xIN,xTMP) call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hae2heeq1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI ToD call geo2gei1(iyr,idoy,secs,xIN,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.13) then !J2000 call geo2gei1(iyr,idoy,secs,xIN,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo CALL ECITOD_2_ECIJ2000(iyr,idoy,secs,xTMP,xOUT) endif endif c c input=GSM if (sysaxesIN.EQ.2) then if (sysaxesOUT.EQ.0) then !GDZ call gsm2geo1(iyr,idoy,secs,psi,xIN,xTMP) call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(2),xOUT(3),xOUT(1)) endif if (sysaxesOUT.EQ.1) then !GEO call gsm2geo1(iyr,idoy,secs,psi,xIN,xOUT) endif if (sysaxesOUT.EQ.3) then !GSE call gsm2geo1(iyr,idoy,secs,psi,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.4) then !SM call gsm2sm1(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.5) then !GEI call gsm2geo1(iyr,idoy,secs,psi,xIN,xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.6) then !MAG call gsm2geo1(iyr,idoy,secs,psi,xIN,xTMP) call geo2mag1(iyr,xTMP,xOUT) endif if (sysaxesOUT.EQ.7) then !SPH call gsm2geo1(iyr,idoy,secs,psi,xIN,xTMP) call CAR_SPH(xTMP,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.8) then !RLL call gsm2geo1(iyr,idoy,secs,psi,xIN,xTMP) call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(1),xOUT(2),xOUT(3)) do i=1,3 xTMP(i)=xOUT(i) enddo xOUT(1)=SQRT(xIN(1)*xIN(1)+xIN(2)*xIN(2)+xIN(3)*xIN(3)) xOUT(2)=xTMP(1) xOUT(3)=xTMP(2) endif if (sysaxesOUT.EQ.9) then !HEE call gsm2geo1(iyr,idoy,secs,psi,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.10) then !HAE call gsm2geo1(iyr,idoy,secs,psi,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.11) then !HEEQ call gsm2geo1(iyr,idoy,secs,psi,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hae2heeq1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI ToD call gsm2geo1(iyr,idoy,secs,psi,xIN,xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.13) then !J2000 call gsm2geo1(iyr,idoy,secs,psi,xIN,xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo CALL ECITOD_2_ECIJ2000(iyr,idoy,secs,xTMP,xOUT) endif endif c c input=GSE if (sysaxesIN.EQ.3) then if (sysaxesOUT.EQ.0) then !GDZ call gse2geo1(iyr,idoy,secs,xIN,xTMP) call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(2),xOUT(3),xOUT(1)) endif if (sysaxesOUT.EQ.1) then !GEO call gse2geo1(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.2) then !GSM call gse2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gsm1(iyr,idoy,secs,psi,xTMP,xOUT) endif if (sysaxesOUT.EQ.4) then !SM call gse2geo1(iyr,idoy,secs,xIN,xTMP) call geo2sm1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.5) then !GEI call gse2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.6) then !MAG call gse2geo1(iyr,idoy,secs,xIN,xTMP) call geo2mag1(iyr,xTMP,xOUT) endif if (sysaxesOUT.EQ.7) then !SPH call gse2geo1(iyr,idoy,secs,xIN,xTMP) call CAR_SPH(xTMP,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.8) then !RLL call gse2geo1(iyr,idoy,secs,xIN,xTMP) call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(1),xOUT(2),xOUT(3)) do i=1,3 xTMP(i)=xOUT(i) enddo xOUT(1)=SQRT(xIN(1)*xIN(1)+xIN(2)*xIN(2)+xIN(3)*xIN(3)) xOUT(2)=xTMP(1) xOUT(3)=xTMP(2) endif if (sysaxesOUT.EQ.9) then !HEE call gse2hee1(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.10) then !HAE call gse2hee1(iyr,idoy,secs,xIN,xTMP) call hee2hae1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.11) then !HEEQ call gse2hee1(iyr,idoy,secs,xIN,xTMP) call hee2hae1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hae2heeq1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI ToD call gse2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.13) then !J2000 call gse2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo CALL ECITOD_2_ECIJ2000(iyr,idoy,secs,xTMP,xOUT) endif endif c c input=SM if (sysaxesIN.EQ.4) then if (sysaxesOUT.EQ.0) then !GDZ call sm2geo1(iyr,idoy,secs,xIN,xTMP) call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(2),xOUT(3),xOUT(1)) endif if (sysaxesOUT.EQ.1) then !GEO call sm2geo1(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.2) then !GSM call sm2gsm1(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.3) then !GSE call sm2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.5) then !GEI call sm2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.6) then !MAG call sm2geo1(iyr,idoy,secs,xIN,xTMP) call geo2mag1(iyr,xTMP,xOUT) endif if (sysaxesOUT.EQ.7) then !SPH call sm2geo1(iyr,idoy,secs,xIN,xTMP) call CAR_SPH(xTMP,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.8) then !RLL call sm2geo1(iyr,idoy,secs,xIN,xTMP) call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(1),xOUT(2),xOUT(3)) do i=1,3 xTMP(i)=xOUT(i) enddo xOUT(1)=SQRT(xIN(1)*xIN(1)+xIN(2)*xIN(2)+xIN(3)*xIN(3)) xOUT(2)=xTMP(1) xOUT(3)=xTMP(2) endif if (sysaxesOUT.EQ.9) then !HEE call sm2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.10) then !HAE call sm2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.11) then !HEEQ call sm2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hae2heeq1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI ToD call sm2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.13) then !J2000 call sm2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo CALL ECITOD_2_ECIJ2000(iyr,idoy,secs,xTMP,xOUT) endif endif c input=GEI if (sysaxesIN.EQ.5) then if (sysaxesOUT.EQ.0) then !GDZ call gei2geo1(iyr,idoy,secs,xIN,xTMP) call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(2),xOUT(3),xOUT(1)) endif if (sysaxesOUT.EQ.1) then !GEO call gei2geo1(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.2) then !GSM call gei2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gsm1(iyr,idoy,secs,psi,xTMP,xOUT) endif if (sysaxesOUT.EQ.3) then !GSE call gei2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.4) then !SM call gei2geo1(iyr,idoy,secs,xIN,xTMP) call geo2sm1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.6) then !MAG call gei2geo1(iyr,idoy,secs,xIN,xTMP) call geo2mag1(iyr,xTMP,xOUT) endif if (sysaxesOUT.EQ.7) then !SPH call gei2geo1(iyr,idoy,secs,xIN,xTMP) call CAR_SPH(xTMP,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.8) then !RLL call gei2geo1(iyr,idoy,secs,xIN,xTMP) call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(1),xOUT(2),xOUT(3)) do i=1,3 xTMP(i)=xOUT(i) enddo xOUT(1)=SQRT(xIN(1)*xIN(1)+xIN(2)*xIN(2)+xIN(3)*xIN(3)) xOUT(2)=xTMP(1) xOUT(3)=xTMP(2) endif if (sysaxesOUT.EQ.9) then !HEE call gei2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.10) then !HAE call gei2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.11) then !HEEQ call gei2geo1(iyr,idoy,secs,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hae2heeq1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI ToD call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.13) then !J2000 call teme2ecitod(iyr,idoy,secs,xIN,xTMP) CALL ECITOD_2_ECIJ2000(iyr,idoy,secs,xTMP,xOUT) endif endif c c c input=MAG if (sysaxesIN.EQ.6) then if (sysaxesOUT.EQ.0) then !GDZ call mag2geo1(iyr,xIN,xTMP) call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(2),xOUT(3),xOUT(1)) endif if (sysaxesOUT.EQ.1) then !GEO call mag2geo1(iyr,xIN,xOUT) endif if (sysaxesOUT.EQ.2) then !GSM call mag2geo1(iyr,xIN,xTMP) call geo2gsm1(iyr,idoy,secs,psi,xTMP,xOUT) endif if (sysaxesOUT.EQ.3) then !GSE call mag2geo1(iyr,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.4) then !SM call mag2geo1(iyr,xIN,xTMP) call geo2sm1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.5) then !GEI call mag2geo1(iyr,xIN,xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.7) then !SPH call mag2geo1(iyr,xIN,xTMP) call CAR_SPH(xTMP,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.8) then !RLL call mag2geo1(iyr,xIN,xTMP) call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(1),xOUT(2),xOUT(3)) do i=1,3 xTMP(i)=xOUT(i) enddo xOUT(1)=SQRT(xIN(1)*xIN(1)+xIN(2)*xIN(2)+xIN(3)*xIN(3)) xOUT(2)=xTMP(1) xOUT(3)=xTMP(2) endif if (sysaxesOUT.EQ.9) then !HEE call mag2geo1(iyr,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.10) then !HAE call mag2geo1(iyr,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.11) then !HEEQ call mag2geo1(iyr,xIN,xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hae2heeq1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI ToD call mag2geo1(iyr,xIN,xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.13) then !J2000 call mag2geo1(iyr,xIN,xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo CALL ECITOD_2_ECIJ2000(iyr,idoy,secs,xTMP,xOUT) endif endif c c input=SPH if (sysaxesIN.EQ.7) then if (sysaxesOUT.EQ.0) then !GDZ call SPH_CAR(xIN(1),xIN(2),xIN(3),xTMP) call geo_gdz(xTMP(1),xTMP(2),xTMP(3) & ,xOUT(2),xOUT(3),xOUT(1)) endif if (sysaxesOUT.EQ.1) then !GEO call SPH_CAR(xIN(1),xIN(2),xIN(3),xOUT) endif if (sysaxesOUT.EQ.2) then !GSM call SPH_CAR(xIN(1),xIN(2),xIN(3),xTMP) call geo2gsm1(iyr,idoy,secs,psi,xTMP,xOUT) endif if (sysaxesOUT.EQ.3) then !GSE call SPH_CAR(xIN(1),xIN(2),xIN(3),xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.4) then !SM call SPH_CAR(xIN(1),xIN(2),xIN(3),xTMP) call geo2sm1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.5) then !GEI call SPH_CAR(xIN(1),xIN(2),xIN(3),xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.6) then !MAG call SPH_CAR(xIN(1),xIN(2),xIN(3),xTMP) call geo2mag1(iyr,xTMP,xOUT) endif if (sysaxesOUT.EQ.8) then !RLL call SPH_CAR(xIN(1),xIN(2),xIN(3),xTMP) call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(1),xOUT(2),xOUT(3)) do i=1,3 xTMP(i)=xOUT(i) enddo xOUT(1)=xIN(1) xOUT(2)=xTMP(1) xOUT(3)=xTMP(2) endif if (sysaxesOUT.EQ.9) then !HEE call SPH_CAR(xIN(1),xIN(2),xIN(3),xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.10) then !HAE call SPH_CAR(xIN(1),xIN(2),xIN(3),xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.11) then !HEEQ call SPH_CAR(xIN(1),xIN(2),xIN(3),xTMP) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hae2heeq1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI ToD call SPH_CAR(xIN(1),xIN(2),xIN(3),xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.13) then !J2000 call SPH_CAR(xIN(1),xIN(2),xIN(3),xTMP) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo CALL ECITOD_2_ECIJ2000(iyr,idoy,secs,xTMP,xOUT) endif endif c c input=RLL if (sysaxesIN.EQ.8) then if (sysaxesOUT.EQ.0) then !GDZ call RLL_GDZ(xIN(1),xIN(2),xIN(3),alti) call gdz_geo(xIN(2),xIN(3),alti,xTMP(1),xTMP(2),xTMP(3)) call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(2),xOUT(3),xOUT(1)) endif if (sysaxesOUT.EQ.1) then !GEO call RLL_GDZ(xIN(1),xIN(2),xIN(3),alti) call gdz_geo(xIN(2),xIN(3),alti,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.2) then !GSM call RLL_GDZ(xIN(1),xIN(2),xIN(3),alti) call gdz_geo(xIN(2),xIN(3),alti,xTMP(1),xTMP(2),xTMP(3)) call geo2gsm1(iyr,idoy,secs,psi,xTMP,xOUT) endif if (sysaxesOUT.EQ.3) then !GSE call RLL_GDZ(xIN(1),xIN(2),xIN(3),alti) call gdz_geo(xIN(2),xIN(3),alti,xTMP(1),xTMP(2),xTMP(3)) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.4) then !SM call RLL_GDZ(xIN(1),xIN(2),xIN(3),alti) call gdz_geo(xIN(2),xIN(3),alti,xTMP(1),xTMP(2),xTMP(3)) call geo2sm1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.5) then !GEI call RLL_GDZ(xIN(1),xIN(2),xIN(3),alti) call gdz_geo(xIN(2),xIN(3),alti,xTMP(1),xTMP(2),xTMP(3)) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.6) then !MAG call RLL_GDZ(xIN(1),xIN(2),xIN(3),alti) call gdz_geo(xIN(2),xIN(3),alti,xTMP(1),xTMP(2),xTMP(3)) call geo2mag1(iyr,xTMP,xOUT) endif if (sysaxesOUT.EQ.7) then !SPH call RLL_GDZ(xIN(1),xIN(2),xIN(3),alti) call gdz_geo(xIN(2),xIN(3),alti,xTMP(1),xTMP(2),xTMP(3)) call CAR_SPH(xTMP,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.9) then !HEE call RLL_GDZ(xIN(1),xIN(2),xIN(3),alti) call gdz_geo(xIN(2),xIN(3),alti,xTMP(1),xTMP(2),xTMP(3)) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.10) then !HAE call RLL_GDZ(xIN(1),xIN(2),xIN(3),alti) call gdz_geo(xIN(2),xIN(3),alti,xTMP(1),xTMP(2),xTMP(3)) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.11) then !HEEQ call RLL_GDZ(xIN(1),xIN(2),xIN(3),alti) call gdz_geo(xIN(2),xIN(3),alti,xTMP(1),xTMP(2),xTMP(3)) call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hae2heeq1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI ToD call RLL_GDZ(xIN(1),xIN(2),xIN(3),alti) call gdz_geo(xIN(2),xIN(3),alti,xTMP(1),xTMP(2),xTMP(3)) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.13) then !J2000 call RLL_GDZ(xIN(1),xIN(2),xIN(3),alti) call gdz_geo(xIN(2),xIN(3),alti,xTMP(1),xTMP(2),xTMP(3)) call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo CALL ECITOD_2_ECIJ2000(iyr,idoy,secs,xTMP,xOUT) endif endif c c input=HEE if (sysaxesIN.EQ.9) then if (sysaxesOUT.EQ.0) then !GDZ call hee2gse1(iyr,idoy,secs,xIN,xTMP) call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(2),xOUT(3),xOUT(1)) endif if (sysaxesOUT.EQ.1) then !GEO call hee2gse1(iyr,idoy,secs,xIN,xTMP) call gse2geo1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.2) then !GSM call hee2gse1(iyr,idoy,secs,xIN,xTMP) call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gsm1(iyr,idoy,secs,psi,xTMP,xOUT) endif if (sysaxesOUT.EQ.3) then !GSE call hee2gse1(iyr,idoy,secs,xIN,xTMP) call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gse1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.4) then !SM call hee2gse1(iyr,idoy,secs,xIN,xTMP) call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2sm1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.5) then !GEI call hee2gse1(iyr,idoy,secs,xIN,xTMP) call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gei1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.6) then !MAG call hee2gse1(iyr,idoy,secs,xIN,xTMP) call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2mag1(iyr,xTMP,xOUT) endif if (sysaxesOUT.EQ.7) then !SPH call hee2gse1(iyr,idoy,secs,xIN,xTMP) call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call CAR_SPH(xTMP,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.8) then !RLL call hee2gse1(iyr,idoy,secs,xIN,xTMP) call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(1),xOUT(2),xOUT(3)) xOUT(3)=xOUT(2) xOUT(2)=xOUT(1) xOUT(1)=SQRT(xTMP(1)*xTMP(1)+xTMP(2)*xTMP(2) & +xTMP(3)*xTMP(3)) endif if (sysaxesOUT.EQ.10) then !HAE call hee2hae1(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.11) then !HEEQ call hee2hae1(iyr,idoy,secs,xIN,xTMP) call hae2heeq1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI ToD call hee2gse1(iyr,idoy,secs,xIN,xTMP) call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.13) then !J2000 call hee2gse1(iyr,idoy,secs,xIN,xTMP) call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo CALL ECITOD_2_ECIJ2000(iyr,idoy,secs,xTMP,xOUT) endif endif c c input=HAE if (sysaxesIN.EQ.10) then if (sysaxesOUT.EQ.0) then !GDZ call hae2hee1(iyr,idoy,secs,xIN,xTMP) call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(2),xOUT(3),xOUT(1)) endif if (sysaxesOUT.EQ.1) then !GEO call hae2hee1(iyr,idoy,secs,xIN,xTMP) call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.2) then !GSM call hae2hee1(iyr,idoy,secs,xIN,xTMP) call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gsm1(iyr,idoy,secs,psi,xTMP,xOUT) endif if (sysaxesOUT.EQ.3) then !GSE call hae2hee1(iyr,idoy,secs,xIN,xTMP) call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gse1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.4) then !SM call hae2hee1(iyr,idoy,secs,xIN,xTMP) call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2sm1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.5) then !GEI call hae2hee1(iyr,idoy,secs,xIN,xTMP) call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gei1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.6) then !MAG call hae2hee1(iyr,idoy,secs,xIN,xTMP) call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2mag1(iyr,xTMP,xOUT) endif if (sysaxesOUT.EQ.7) then !SPH call hae2hee1(iyr,idoy,secs,xIN,xTMP) call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call CAR_SPH(xTMP,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.8) then !RLL call hae2hee1(iyr,idoy,secs,xIN,xTMP) call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(1),xOUT(2),xOUT(3)) xOUT(3)=xOUT(2) xOUT(2)=xOUT(1) xOUT(1)=SQRT(xTMP(1)*xTMP(1)+xTMP(2)*xTMP(2) & +xTMP(3)*xTMP(3)) endif if (sysaxesOUT.EQ.9) then !HEE call hae2hee1(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.11) then !HEEQ call hae2heeq1(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI ToD call hae2hee1(iyr,idoy,secs,xIN,xTMP) call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.13) then !J2000 call hae2hee1(iyr,idoy,secs,xIN,xTMP) call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo CALL ECITOD_2_ECIJ2000(iyr,idoy,secs,xTMP,xOUT) endif endif c c input=HEEQ if (sysaxesIN.EQ.11) then if (sysaxesOUT.EQ.0) then !GDZ call heeq2hae1(iyr,idoy,secs,xIN,xTMP) call hae2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(2),xOUT(3),xOUT(1)) endif if (sysaxesOUT.EQ.1) then !GEO call heeq2hae1(iyr,idoy,secs,xIN,xTMP) call hae2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.2) then !GSM call heeq2hae1(iyr,idoy,secs,xIN,xTMP) call hae2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gsm1(iyr,idoy,secs,psi,xTMP,xOUT) endif if (sysaxesOUT.EQ.3) then !GSE call heeq2hae1(iyr,idoy,secs,xIN,xTMP) call hae2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gse1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.4) then !SM call heeq2hae1(iyr,idoy,secs,xIN,xTMP) call hae2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2sm1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.5) then !GEI call heeq2hae1(iyr,idoy,secs,xIN,xTMP) call hae2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gei1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.6) then !MAG call heeq2hae1(iyr,idoy,secs,xIN,xTMP) call hae2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2mag1(iyr,xTMP,xOUT) endif if (sysaxesOUT.EQ.7) then !SPH call heeq2hae1(iyr,idoy,secs,xIN,xTMP) call hae2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call CAR_SPH(xTMP,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.8) then !RLL call heeq2hae1(iyr,idoy,secs,xIN,xTMP) call hae2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(1),xOUT(2),xOUT(3)) xOUT(3)=xOUT(2) xOUT(2)=xOUT(1) xOUT(1)=SQRT(xTMP(1)*xTMP(1)+xTMP(2)*xTMP(2) & +xTMP(3)*xTMP(3)) endif if (sysaxesOUT.EQ.9) then !HEE call heeq2hae1(iyr,idoy,secs,xIN,xTMP) call hae2hee1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.10) then !HAE call heeq2hae1(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI ToD call heeq2hae1(iyr,idoy,secs,xIN,xTMP) call hae2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.13) then !J2000 call heeq2hae1(iyr,idoy,secs,xIN,xTMP) call hae2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gei1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call teme2ecitod(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo CALL ECITOD_2_ECIJ2000(iyr,idoy,secs,xTMP,xOUT) endif endif c c input=ECI ToD (True of Date) if (sysaxesIN.EQ.12) then if (sysaxesOUT.EQ.0) then !GDZ call ECITOD2TEME(iyr,idoy,secs,xIN,xTMP) call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(2),xOUT(3),xOUT(1)) endif if (sysaxesOUT.EQ.1) then !GEO call ECITOD2TEME(iyr,idoy,secs,xIN,xTMP) call gei2geo1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.2) then !GSM call ECITOD2TEME(iyr,idoy,secs,xIN,xTMP) call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gsm1(iyr,idoy,secs,psi,xTMP,xOUT) endif if (sysaxesOUT.EQ.3) then !GSE call ECITOD2TEME(iyr,idoy,secs,xIN,xTMP) call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gse1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.4) then !SM call ECITOD2TEME(iyr,idoy,secs,xIN,xTMP) call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2sm1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.5) then !GEI (J2000) call ECITOD2TEME(iyr,idoy,secs,xIN,xOUT) endif if (sysaxesOUT.EQ.6) then !MAG call ECITOD2TEME(iyr,idoy,secs,xIN,xTMP) call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2mag1(iyr,xTMP,xOUT) endif if (sysaxesOUT.EQ.7) then !SPH call ECITOD2TEME(iyr,idoy,secs,xIN,xTMP) call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call CAR_SPH(xTMP,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.8) then !RLL call ECITOD2TEME(iyr,idoy,secs,xIN,xTMP) call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(1),xOUT(2),xOUT(3)) do i=1,3 xTMP(i)=xOUT(i) enddo xOUT(1)=SQRT(xIN(1)*xIN(1)+xIN(2)*xIN(2)+xIN(3)*xIN(3)) xOUT(2)=xTMP(1) xOUT(3)=xTMP(2) endif if (sysaxesOUT.EQ.9) then !HEE call ECITOD2TEME(iyr,idoy,secs,xIN,xTMP) call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.10) then !HAE call ECITOD2TEME(iyr,idoy,secs,xIN,xTMP) call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.11) then !HEEQ call ECITOD2TEME(iyr,idoy,secs,xIN,xTMP) call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hae2heeq1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI TOD CALL ECITOD2TEME(iyr,idoy,secs,xIN,xOUT) ENDIF endif c input=J2000 if (sysaxesIN.EQ.13) then if (sysaxesOUT.EQ.0) then !GDZ call ECIJ2000_2_ECITOD(iyr,idoy,secs,xIN,xTMP) call ECITOD2TEME(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(2),xOUT(3),xOUT(1)) endif if (sysaxesOUT.EQ.1) then !GEO call ECIJ2000_2_ECITOD(iyr,idoy,secs,xIN,xTMP) call ECITOD2TEME(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gei2geo1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.2) then !GSM call ECIJ2000_2_ECITOD(iyr,idoy,secs,xIN,xTMP) call ECITOD2TEME(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gsm1(iyr,idoy,secs,psi,xTMP,xOUT) endif if (sysaxesOUT.EQ.3) then !GSE call ECIJ2000_2_ECITOD(iyr,idoy,secs,xIN,xTMP) call ECITOD2TEME(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gse1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.4) then !SM call ECIJ2000_2_ECITOD(iyr,idoy,secs,xIN,xTMP) call ECITOD2TEME(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2sm1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.5) then !GEI (TEME) call ECIJ2000_2_ECITOD(iyr,idoy,secs,xIN,xTMP) call ECITOD2TEME(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.6) then !MAG call ECIJ2000_2_ECITOD(iyr,idoy,secs,xIN,xTMP) call ECITOD2TEME(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2mag1(iyr,xTMP,xOUT) endif if (sysaxesOUT.EQ.7) then !SPH call ECIJ2000_2_ECITOD(iyr,idoy,secs,xIN,xTMP) call ECITOD2TEME(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call CAR_SPH(xTMP,xOUT(1),xOUT(2),xOUT(3)) endif if (sysaxesOUT.EQ.8) then !RLL call ECIJ2000_2_ECITOD(iyr,idoy,secs,xIN,xTMP) call ECITOD2TEME(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo_gdz(xTMP(1),xTMP(2),xTMP(3), & xOUT(1),xOUT(2),xOUT(3)) do i=1,3 xTMP(i)=xOUT(i) enddo xOUT(1)=SQRT(xIN(1)*xIN(1)+xIN(2)*xIN(2)+xIN(3)*xIN(3)) xOUT(2)=xTMP(1) xOUT(3)=xTMP(2) endif if (sysaxesOUT.EQ.9) then !HEE call ECIJ2000_2_ECITOD(iyr,idoy,secs,xIN,xTMP) call ECITOD2TEME(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.10) then !HAE call ECIJ2000_2_ECITOD(iyr,idoy,secs,xIN,xTMP) call ECITOD2TEME(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.11) then !HEEQ call ECIJ2000_2_ECITOD(iyr,idoy,secs,xIN,xTMP) call ECITOD2TEME(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gei2geo1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call geo2gse1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call gse2hee1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hee2hae1(iyr,idoy,secs,xTMP,xOUT) do i=1,3 xTMP(i)=xOUT(i) enddo call hae2heeq1(iyr,idoy,secs,xTMP,xOUT) endif if (sysaxesOUT.EQ.12) then !ECI TOD call ECIJ2000_2_ECITOD(iyr,idoy,secs,xIN,xTMP) endif endif c c end C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE coord_trans_vec C Routine to vectorize coord_trans1, the single-element C coordinate transformation utility. C Calling sequence: C C call coord_trans_vec(ntime,sysaxesIN,sysaxesOUT,iyear,idoy,secs, C + xIN,xOUT) C C with the following types: C INTEGER*4 sysaxesIN, sysaxesOUT C INTEGER*4 iyear(nmax),idoy(nmax) C REAL*8 secs(nmax) C REAL*8 xINV(3,nmax), xOUTV(3,nmax) C INTEGER*4 numpoints C C As with all (most?) onera library calls, the maximum array size C is limited to ntime_max elements C Contributed by Timothy Guild, 2.2.07 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE coord_trans_vec1(ntime,sysaxesIN,sysaxesOUT, & iyear,idoy,secs,xINV,xOUTV) IMPLICIT NONE INCLUDE 'ntime_max.inc' INCLUDE 'variables.inc' INTEGER*4 nmax,i,ntime, sysaxesIN, sysaxesOUT INTEGER*4 iyear(ntime_max),idoy(ntime_max),y,d REAL*8 secs(ntime_max),xINV(3,ntime_max),xOUTV(3,ntime_max) ! local vars REAL*8 xIN(3),xOUT(3),s ! Loop over the number of points specified, calling ! coord_trans1 each iteration DO i = 1,ntime y = iyear(i) d = idoy(i) s = secs(i) if (xINV(1,i) .eq. baddata .and. xINV(2,i) .eq. baddata & .and. xINV(3,i) .eq.baddata) then xOUTV(1,i) = baddata xOUTV(2,i) = baddata xOUTV(3,i) = baddata goto 10 endif xIN(1) = xINV(1,i) ! copy each array element into 3x1 array xIN(2) = xINV(2,i) xIN(3) = xINV(3,i) call coord_trans1( sysaxesIN,sysaxesOUT,y,d,s,xIN,xOUT) xOUTV(1,i) = xOUT(1) ! copy back to 3xN array to pass back xOUTV(2,i) = xOUT(2) xOUTV(3,i) = xOUT(3) 10 continue ENDDO END