!*********************************************************************************************************************************** ! ! D E N S I T Y 2 4 ! ! Program: DENSITY24 ! ! Programmer: David G. Simpson ! NASA Goddard Space Flight Center ! Greenbelt, Maryland 20771 ! ! Date: August 7, 2003 ! ! Language: Fortran-90 ! ! Version: 1.00a ! ! Description: ISN density model. ! ! Files: Source files: ! ! density24.f90 Main program ! ! Notes: Re-worked to follow notes in Pad 12. ! Changed finding TPP from Golden Search to celestial mechanics method. ! Marked all points _ISN or _CEL to indicate reference frame. ! 3D added. ! Cleaned up. ! ISN, ECL, EQ coordinates. ! Replaced calculation of parallelepiped volume with Pappus-Guldinus theorem for torus volume. Fixed bug in ZTP_ISN. ! !*********************************************************************************************************************************** !*********************************************************************************************************************************** ! Global variables !*********************************************************************************************************************************** MODULE DENSITY_GLOBAL DOUBLE PRECISION, PARAMETER :: PI = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803D0 DOUBLE PRECISION, PARAMETER :: TWOPI = 6.28318530717958647692528676655900576839433879875021164194988918461563281257241799726D0 DOUBLE PRECISION, PARAMETER :: HALFPI = 1.5707963267948966192313216916397514420985846996875529104874722961539082031431044993D0 DOUBLE PRECISION, PARAMETER :: AU = 1.49597870D11 ! astronomical unit (m) DOUBLE PRECISION, PARAMETER :: GM = 1.32712438D20 ! GM(sun) (m**3 s**-2) DOUBLE PRECISION, PARAMETER :: BDEL = 1.0D-4 ! AU DOUBLE PRECISION, PARAMETER :: INFINITY = 100.0D0 ! AU DOUBLE PRECISION :: XPP_ISN, YPP_ISN DOUBLE PRECISION :: N, T0, ECC1, ECC2, A, XPP_CEL, YPP_CEL, MATA(2,2), V0, MU DOUBLE PRECISION :: VELOCITY_ISN(3), VELOCITY_ISN_MAG DOUBLE PRECISION :: DBG_SEP, DBG_ACCEL, DBG_DELTA_T DOUBLE PRECISION :: XPP_ISN_TRUE, YPP_ISN_TRUE, ZPP_ISN_TRUE INTEGER :: DBG_CODE LOGICAL :: DEBUG END MODULE DENSITY_GLOBAL !*********************************************************************************************************************************** ! CALC_DENSITY ! ! HYP = 1: top hyperbola ! 2: bottom hyperbola ! !*********************************************************************************************************************************** FUNCTION CALC_DENSITY (X, Y, Z, V0_LCL, MU_LCL, BETA, HYP, CODE) RESULT (DENSITY) USE DENSITY_GLOBAL IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: X, Y, Z DOUBLE PRECISION, INTENT(IN) :: V0_LCL DOUBLE PRECISION, INTENT(IN) :: MU_LCL DOUBLE PRECISION, INTENT(IN) :: BETA INTEGER, INTENT(IN) :: HYP INTEGER, INTENT(IN) :: CODE DOUBLE PRECISION :: DENSITY DOUBLE PRECISION :: B, VOLUME_FACTOR DOUBLE PRECISION :: DENSITY1, DENSITY2 DOUBLE PRECISION :: C, R, THETA ! ! Store global variables ! XPP_ISN = X YPP_ISN = SQRT(Y**2+Z**2) ! reduce 3-D to 2-D V0 = V0_LCL MU = MU_LCL C = (1.0D0-MU)*(GM/(V0*V0)) ! m A = C/AU ! AU N = SQRT((1.0D0-MU)*GM/C**3) ! sec**-1 XPP_ISN_TRUE = X YPP_ISN_TRUE = Y ZPP_ISN_TRUE = Z DBG_CODE = CODE ! ! Find polar coordinates ! R = SQRT(XPP_ISN**2+YPP_ISN**2)*AU ! m THETA=ATAN2(YPP_ISN,XPP_ISN) IF (HYP == 1) THEN ! direct orbit B = 0.5D0*(R*SIN(THETA)+SQRT((R*SIN(THETA))**2 + & 4.0D0*R*C*(1.0D0-COS(THETA))))/AU ! AU ELSE ! retrograde orbit B = 0.5D0*(R*SIN(THETA)-SQRT((R*SIN(THETA))**2 + & 4.0D0*R*C*(1.0D0-COS(THETA))))/AU ! AU END IF CALL CALC_VOLUME_FACTOR (B, VOLUME_FACTOR) DENSITY = 1.0D0/VOLUME_FACTOR RETURN END FUNCTION CALC_DENSITY !*********************************************************************************************************************************** ! CALC_VOLUME_FACTOR !*********************************************************************************************************************************** SUBROUTINE CALC_VOLUME_FACTOR (B, VOLUME_FACTOR) USE DENSITY_GLOBAL IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: B DOUBLE PRECISION, INTENT(OUT) :: VOLUME_FACTOR DOUBLE PRECISION :: ALPHA, DIST_MIN, DPJUNK, B2 DOUBLE PRECISION :: VEC2D(2), XP_CEL, YP_CEL, TRUE_ANOMP, RP, FP, MEAN_ANOMP, TP DOUBLE PRECISION :: DEL_TRUE_ANOM, TRUE_ANOMQ, RQ, XQ_CEL, YQ_CEL, FQ DOUBLE PRECISION :: MEAN_ANOMQ, TQ, F_MIN, MA_MIN, T_MIN, F_MAX, MA_MAX, T_MAX DOUBLE PRECISION :: CLOSEST, MEAN_ANOMPP, FPP, TRUE_ANOMPP, RPP, TPP DOUBLE PRECISION :: TQP, MEAN_ANOMQP, FQP, TRUE_ANOMQP, RQP, XQP_CEL, YQP_CEL DOUBLE PRECISION :: XR_CEL, YR_CEL, TR, TRUE_ANOMR, RR, FR, MEAN_ANOMR DOUBLE PRECISION :: TRUE_ANOMS, RS, XS_CEL, YS_CEL, FS DOUBLE PRECISION :: MEAN_ANOMS, TS DOUBLE PRECISION :: MEAN_ANOMRP, FRP, TRUE_ANOMRP, RRP, XRP_CEL, YRP_CEL, TRP DOUBLE PRECISION :: TSP, MEAN_ANOMSP, FSP, TRUE_ANOMSP, RSP, XSP_CEL, YSP_CEL DOUBLE PRECISION :: VOLUME_V, VOLUME_VP DOUBLE PRECISION :: XTEMP, YTEMP, RTEMP, THETATEMP DOUBLE PRECISION :: XP_ISN, YP_ISN, XQ_ISN, YQ_ISN, XR_ISN, YR_ISN, XS_ISN, YS_ISN DOUBLE PRECISION :: XQP_ISN, YQP_ISN, XRP_ISN, YRP_ISN, XSP_ISN, YSP_ISN DOUBLE PRECISION :: XT_ISN, YT_ISN, XTP_ISN, YTP_ISN DOUBLE PRECISION :: ZP_ISN, ZQ_ISN, ZR_ISN, ZT_ISN, ZPP_ISN, ZQP_ISN, ZRP_ISN, ZTP_ISN DOUBLE PRECISION :: T_HEIGHT, A_BASE, Y_CENTROID, Z_CENTROID DOUBLE PRECISION :: DELTA_T DOUBLE PRECISION :: XQP_ISN_TRUE, YQP_ISN_TRUE, ZQP_ISN_TRUE DOUBLE PRECISION :: UPSILON DOUBLE PRECISION :: ACOSH, QUADAREA, TOROID_SECTOR_VOLUME, CALC_ISN_Y INTEGER :: SIGNUM ! ! Find eccentricity of H. ! ECC1 = SQRT(A**2+B**2)/A IF (ECC1 <= 1.0D0) THEN WRITE (UNIT=*, FMT='(A)') ' Error: not a hyperbolic orbit.' WRITE (UNIT=*, FMT='(A,F10.6,A)') ' Impact parameter = ', B, ' AU' WRITE (UNIT=*, FMT='(A,F10.6)') ' Eccentricity = ', ECC1 STOP END IF ! ! Define points P and Q (ISN frame). ! XP_ISN = INFINITY YP_ISN = CALC_ISN_Y (XP_ISN, A, B, SIGNUM(B)) XQ_ISN = INFINITY + BDEL YQ_ISN = CALC_ISN_Y (XQ_ISN, A, B, SIGNUM(B)) WRITE (UNIT=34, FMT='(/A,ES25.12)') ' YP_ISN = ', YP_ISN WRITE (UNIT=34, FMT='(A,ES25.12)') ' YQ_ISN = ', YQ_ISN WRITE (UNIT=34, FMT='(A,ES25.12)') ' B = ', B WRITE (UNIT=34, FMT='(A,I5)') ' CODE = ', DBG_CODE ! ! Convert P', P, and Q from ISN frame to CEL frame. ! CALL CEL_TO_ISN ((/XPP_ISN, YPP_ISN/), A, B, -1, VEC2D) XPP_CEL = VEC2D(1) YPP_CEL = VEC2D(2) CALL CEL_TO_ISN ((/XP_ISN, YP_ISN/), A, B, -1, VEC2D) XP_CEL = VEC2D(1) YP_CEL = VEC2D(2) CALL CEL_TO_ISN ((/XQ_ISN, YQ_ISN/), A, B, -1, VEC2D) XQ_CEL = VEC2D(1) YQ_CEL = VEC2D(2) ! ! Set the epoch time T0. ! T0 = 0.0D0 ! ! Find anomalies and times for P', P, and Q. ! TRUE_ANOMPP = ATAN2(YPP_CEL,XPP_CEL) RPP = SQRT(XPP_CEL**2+YPP_CEL**2) FPP = ACOSH((RPP/A+1.0D0)/ECC1) FPP = SIGN(FPP,TRUE_ANOMPP) MEAN_ANOMPP = ECC1*SINH(FPP)-FPP TPP = MEAN_ANOMPP/N + T0 TRUE_ANOMP = ATAN2(YP_CEL,XP_CEL) RP = SQRT(XP_CEL**2+YP_CEL**2) FP = ACOSH((RP/A+1.0D0)/ECC1) FP = SIGN(FP,TRUE_ANOMP) MEAN_ANOMP = ECC1*SINH(FP)-FP TP = MEAN_ANOMP/N + T0 TRUE_ANOMQ = ATAN2(YQ_CEL,XQ_CEL) RQ = SQRT(XQ_CEL**2+YQ_CEL**2) FQ = ACOSH((RQ/A+1.0D0)/ECC1) FQ = SIGN(FQ,TRUE_ANOMQ) MEAN_ANOMQ = ECC1*SINH(FQ)-FQ TQ = MEAN_ANOMQ/N + T0 ! ! Find time for Q'. ! TQP = TPP - TP + TQ ! ! Find anomalies for Q'. ! CALL CALC_ANOMALIES (TQP, ECC1, TRUE_ANOMQP, RQP) ! ! Find CEL coordinates of Q'. ! XQP_CEL = RQP*COS(TRUE_ANOMQP) YQP_CEL = RQP*SIN(TRUE_ANOMQP) ! ! Rotate Q' from CEL frame back to ISN frame. ! CALL CEL_TO_ISN ((/XQP_CEL, YQP_CEL/), A, B, +1, VEC2D) XQP_ISN = VEC2D(1) YQP_ISN = VEC2D(2) ! ! Now find impact parameter B2 and eccentricity ECC2 of second hyperbola, H2. ! B2 = B - BDEL ECC2 = SQRT(A**2+B2**2)/A IF (ECC2 <= 1.0D0) THEN WRITE (UNIT=*, FMT='(A)') ' Error: not a hyperbolic orbit.' WRITE (UNIT=*, FMT='(A,F10.6,A)') ' Impact parameter = ', B2, ' AU' WRITE (UNIT=*, FMT='(A,F10.6)') ' Eccentricity = ', ECC2 STOP END IF ! ! Define points R and S (ISN frame). ! XR_ISN = INFINITY YR_ISN = CALC_ISN_Y (XR_ISN, A, B2, SIGNUM(B2)) XS_ISN = INFINITY + BDEL YS_ISN = CALC_ISN_Y (XS_ISN, A, B2, SIGNUM(B2)) ! ! Rotate R and S from ISN frame to CEL frame. ! CALL CEL_TO_ISN ((/XR_ISN, YR_ISN/), A, B2, -1, VEC2D) ! <<< B or B2? XR_CEL = VEC2D(1) YR_CEL = VEC2D(2) CALL CEL_TO_ISN ((/XS_ISN, YS_ISN/), A, B2, -1, VEC2D) ! <<< B or B2? XS_CEL = VEC2D(1) YS_CEL = VEC2D(2) ! ! Find anomalies and times for R and S. ! TRUE_ANOMR = ATAN2(YR_CEL,XR_CEL) RR = SQRT(XR_CEL**2+YR_CEL**2) FR = ACOSH((RR/A+1.0D0)/ECC2) FR = SIGN(FR,TRUE_ANOMR) MEAN_ANOMR = ECC2*SINH(FR)-FR TRUE_ANOMS = ATAN2(YS_CEL,XS_CEL) RS = SQRT(XS_CEL**2+YS_CEL**2) FS = ACOSH((RS/A+1.0D0)/ECC2) FS = SIGN(FS,TRUE_ANOMS) MEAN_ANOMS = ECC2*SINH(FS)-FS ! ! Find time for R (on H2's timescale). ! TR = MEAN_ANOMR/N ! time from H2's perihelion to R ! ! Set TDEL, the time from H's perihelion to H2's perihelion. This is required ! so that times on both hyperbolae are measured with respect to the same time ! scale (H's timescale; i.e. time elapsed from perihelion on H). ! T0 = TR - TP ! time from H's perihelion to H2's perihelion ! (time of H2's perihelion, in H's timescale) ! ! Now convert time TR to H's timescale. ! TR = TR + T0 ! time from H's perihelion to R ! ! Find time for S (H's timescale). ! TS = MEAN_ANOMS/N + T0 ! time from H's perihelion to S ! ! Find times for R' and S'. ! TRP = TPP - TP + TR ! time from H's perihelion to R' TSP = TPP - TP + TS ! time from H's perihelion to S' ! ! Find anomalies for R' and S'. ! CALL CALC_ANOMALIES (TRP, ECC2, TRUE_ANOMRP, RRP) CALL CALC_ANOMALIES (TSP, ECC2, TRUE_ANOMSP, RSP) ! ! Find CEL coordinates of R' and S'. ! XRP_CEL = RRP*COS(TRUE_ANOMRP) YRP_CEL = RRP*SIN(TRUE_ANOMRP) XSP_CEL = RSP*COS(TRUE_ANOMSP) YSP_CEL = RSP*SIN(TRUE_ANOMSP) ! ! Rotate R' and S' from CEL frame back to ISN frame. ! CALL CEL_TO_ISN ((/XRP_CEL, YRP_CEL/), A, B2, +1, VEC2D) ! <<< B or B2? XRP_ISN = VEC2D(1) YRP_ISN = VEC2D(2) CALL CEL_TO_ISN ((/XSP_CEL, YSP_CEL/), A, B2, +1, VEC2D) ! <<< B or B2? XSP_ISN = VEC2D(1) YSP_ISN = VEC2D(2) ! ! Store Z coordinates of P, Q, and R. (P, Q, R, and S all lie in the Z=0 plane.) ! ZP_ISN = 0.0D0 ZQ_ISN = 0.0D0 ZR_ISN = 0.0D0 ZPP_ISN = 0.0D0 ZQP_ISN = 0.0D0 ZRP_ISN = 0.0D0 ! ! Define points T and T' for treatment of 3-D. ! T_HEIGHT = 1.0D-1 * BDEL XT_ISN = XP_ISN YT_ISN = YP_ISN / COS(ATAN(T_HEIGHT/YP_ISN)) ZT_ISN = T_HEIGHT XTP_ISN = XPP_ISN YTP_ISN = YPP_ISN ZTP_ISN = ZT_ISN * (YTP_ISN / YT_ISN) ! ! Find volumes. ! A_BASE = QUADAREA (XP_ISN, YP_ISN, XQ_ISN, YQ_ISN, XS_ISN, YS_ISN, XR_ISN, YR_ISN) Y_CENTROID = 0.25D0 * (YP_ISN + YQ_ISN + YR_ISN + YS_ISN) Z_CENTROID = ZT_ISN * (Y_CENTROID / YT_ISN) VOLUME_V = TOROID_SECTOR_VOLUME (A_BASE, Y_CENTROID, Z_CENTROID) A_BASE = QUADAREA (XPP_ISN, YPP_ISN, XQP_ISN, YQP_ISN, XSP_ISN, YSP_ISN, XRP_ISN, YRP_ISN) Y_CENTROID = 0.25D0 * (YPP_ISN + YQP_ISN + YRP_ISN + YSP_ISN) Z_CENTROID = ZTP_ISN * (Y_CENTROID / YTP_ISN) VOLUME_VP = TOROID_SECTOR_VOLUME (A_BASE, Y_CENTROID, Z_CENTROID) ! ! Compute volume factor (ratio of volumes). ! VOLUME_FACTOR = VOLUME_VP / VOLUME_V ! ! Find the true ISN coordinates of Q'. ! UPSILON = ATAN2(ZPP_ISN_TRUE,YPP_ISN_TRUE) XQP_ISN_TRUE = XQP_ISN YQP_ISN_TRUE = YQP_ISN * COS(UPSILON) ZQP_ISN_TRUE = YQP_ISN * SIN(UPSILON) ! ! Compute the velocity vector. ! DELTA_T = TPP - TQP DBG_DELTA_T = DELTA_T VELOCITY_ISN(1) = (XPP_ISN_TRUE-XQP_ISN_TRUE)*AU/DELTA_T VELOCITY_ISN(2) = (YPP_ISN_TRUE-YQP_ISN_TRUE)*AU/DELTA_T VELOCITY_ISN(3) = (ZPP_ISN_TRUE-ZQP_ISN_TRUE)*AU/DELTA_T VELOCITY_ISN_MAG = SQRT(DOT_PRODUCT(VELOCITY_ISN,VELOCITY_ISN)) WRITE (UNIT=33, FMT='(/A,ES23.15,A,I3)') ' CALCULATED V = ', SQRT(DOT_PRODUCT(VELOCITY_ISN,VELOCITY_ISN)), & ' CODE = ', DBG_CODE WRITE (UNIT=33, FMT='(A,ES23.15,A,I3)') ' EXPECTED V = ', SQRT(GM*(2.0D0/(RPP*AU) + 1.0D0/(A*AU))), & ' CODE = ', DBG_CODE ! ! Save separation and acceleration distances (secondary information) ! DBG_SEP = SQRT((XQP_ISN-XPP_ISN)**2+(YQP_ISN-YPP_ISN)**2)/SQRT((XQ_ISN-XP_ISN)**2+(YQ_ISN-YP_ISN)**2) DBG_ACCEL = SQRT((XRP_ISN-XPP_ISN)**2+(YRP_ISN-YPP_ISN)**2)/SQRT((XR_ISN-XP_ISN)**2+(YR_ISN-YP_ISN)**2) RETURN END SUBROUTINE CALC_VOLUME_FACTOR !*********************************************************************************************************************************** ! HYP_KEPLER ! ! Solve Kepler's equation for the hyperbolic case (e > 1) (Newton-Raphson method): ! ! M = e sinh F - F ! !*********************************************************************************************************************************** FUNCTION HYP_KEPLER (MA, ECC) RESULT (FN) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: MA ! mean anomaly (rad) DOUBLE PRECISION, INTENT(IN) :: ECC ! eccentricity DOUBLE PRECISION :: FN, FN_PREV ! eccentric anomaly (rad) DOUBLE PRECISION, PARAMETER :: EPS = 1.0D-15 DOUBLE PRECISION :: M INTEGER :: J FN = MA FN_PREV = MA DO J = 1, 100 FN = FN - (MA-ECC*SINH(FN)+FN)/(-ECC*COSH(FN)+1.0D0) ! one iteration of Newton-Raphson IF (ABS(FN-FN_PREV) < EPS) EXIT ! exit loop if converged FN_PREV = FN END DO RETURN END FUNCTION HYP_KEPLER !*********************************************************************************************************************************** ! LOG1P ! ! Return LOG(1+X) (natural logarithm). !*********************************************************************************************************************************** FUNCTION LOG1P (X) RESULT (Y) DOUBLE PRECISION, INTENT(IN) :: X DOUBLE PRECISION :: Y DOUBLE PRECISION :: Z Z = 1.0D0 + X Y = LOG(Z) - ((Z-1.0D0)-X)/Z ! cancels errors with IEEE arithmetic RETURN END FUNCTION LOG1P !*********************************************************************************************************************************** ! ACOSH ! ! Inverse hyperbolic cosine. !*********************************************************************************************************************************** FUNCTION ACOSH (X) RESULT (Y) DOUBLE PRECISION, INTENT(IN) :: X DOUBLE PRECISION :: Y DOUBLE PRECISION :: T DOUBLE PRECISION, PARAMETER :: LN2 = 0.6931471805599453094172321214581765680755001343602552541206800094933936219696947156059D0 DOUBLE PRECISION :: LOG1P IF (X > (1.0D0 / SQRT(EPSILON(1.0D0)))) THEN Y = LOG (X) + LN2 ELSE IF (X > 2.0D0) THEN Y = LOG (2.0D0 * X - 1.0D0 / (SQRT (X * X - 1.0D0) + X)) ELSE IF (X > 1.0D0) THEN T = X - 1.0D0 Y = LOG1P (T + SQRT (2.0D0 * T + T * T)) ELSE IF (X == 1.0D0) THEN Y = 0.0D0 ELSE Y = 0.0D0 WRITE (UNIT=*, FMT='(A)') ' ACOSH Error.' END IF RETURN END FUNCTION ACOSH !*********************************************************************************************************************************** ! CEL_TO_ISN ! ! Convert "celestial mechanics" (CEL) coordinates to or from ISN coordinates. ! ! Set DIRECTION=+1 for "forward" conversion (CEL->ISN), or DIRECTION=-1 for "reverse" conversion (ISN->CEL). !*********************************************************************************************************************************** SUBROUTINE CEL_TO_ISN (VEC_IN, A_LCL, B_LCL, DIRECTION, VEC_OUT) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: VEC_IN(2) DOUBLE PRECISION, INTENT(IN) :: A_LCL DOUBLE PRECISION, INTENT(IN) :: B_LCL INTEGER, INTENT(IN) :: DIRECTION DOUBLE PRECISION, INTENT(OUT) :: VEC_OUT(2) DOUBLE PRECISION, PARAMETER :: PI = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803D0 DOUBLE PRECISION :: ECC DOUBLE PRECISION :: ALPHA DOUBLE PRECISION :: MATA(2,2) ECC = SQRT(A_LCL**2+B_LCL**2)/A_LCL ALPHA = ACOS(1.0D0/ECC) - PI IF (B_LCL < 0.0D0) ALPHA = -ALPHA MATA = RESHAPE ( (/ COS(ALPHA), SIN(ALPHA), & ! construct rotation matrix -SIN(ALPHA), COS(ALPHA) /), & SHAPE = (/ 2, 2 /), ORDER = (/ 2, 1 /) ) IF (DIRECTION > 0) THEN VEC_OUT = MATMUL(MATA, VEC_IN) ELSE VEC_OUT = MATMUL(TRANSPOSE(MATA), VEC_IN) END IF RETURN END SUBROUTINE CEL_TO_ISN !*********************************************************************************************************************************** ! QUADAREA ! ! Compute the area of a quadrilateral whose vertices have coordinates (X1,Y1), (X2,Y2), (X3,Y3), and (X4,Y4). !*********************************************************************************************************************************** FUNCTION QUADAREA (X1, Y1, X2, Y2, X3, Y3, X4, Y4) RESULT (AREA) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: X1, Y1, X2, Y2, X3, Y3, X4, Y4 DOUBLE PRECISION :: AREA AREA = ABS(0.5D0*(X1*Y2 + X2*Y3 + X3*Y4 + X4*Y1 - Y1*X2 - Y2*X3 - Y3*X4 - Y4*X1)) RETURN END FUNCTION QUADAREA !*********************************************************************************************************************************** ! TOROID_SECTOR_VOLUME ! ! Compute the volume of a sector of a toroid, using the Pappus-Guldinus Theorem. !*********************************************************************************************************************************** FUNCTION TOROID_SECTOR_VOLUME (A_BASE, Y_CENTROID, Z_CENTROID) RESULT (VOLUME) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: A_BASE DOUBLE PRECISION, INTENT(IN) :: Y_CENTROID, Z_CENTROID DOUBLE PRECISION :: VOLUME DOUBLE PRECISION :: THETA THETA = ATAN(Z_CENTROID/Y_CENTROID) VOLUME = A_BASE * ABS(Y_CENTROID * THETA) RETURN END FUNCTION TOROID_SECTOR_VOLUME !*********************************************************************************************************************************** ! CALC_ISN_Y ! ! Return the Y coordinate of an ISN orbit (ISN frame) given X. ! SGN=+1 for 'top', -1 for 'bottom'. !*********************************************************************************************************************************** FUNCTION CALC_ISN_Y (X, A, B, SGN) RESULT (Y) IMPLICIT NONE DOUBLE PRECISION, INTENT (IN) :: X DOUBLE PRECISION, INTENT (IN) :: A DOUBLE PRECISION, INTENT (IN) :: B INTEGER, INTENT (IN) :: SGN DOUBLE PRECISION :: Y IF (SGN > 0) THEN Y = +B + A*B * (+(X+A)-SQRT((X+A)**2+B**2-A**2))/(B**2-A**2) ELSE Y = -B + A*B * (-(X+A)+SQRT((X+A)**2+B**2-A**2))/(B**2-A**2) Y = -Y ! same as B = -B END IF RETURN END FUNCTION CALC_ISN_Y !*********************************************************************************************************************************** ! SIGNUM !*********************************************************************************************************************************** FUNCTION SIGNUM (X) RESULT (Y) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: X INTEGER :: Y IF (X < 0.0D0) THEN Y = -1 ELSE Y = +1 END IF RETURN END FUNCTION SIGNUM !*********************************************************************************************************************************** ! ECL_TO_ISN ! ! Convert ecliptic (ECL) coordinates to or from ISN coordinates. ! ! Set DIRECTION=+1 for "forward" conversion (ECL->ISN), or DIRECTION=-1 for "reverse" conversion (ISN->ECL). !*********************************************************************************************************************************** SUBROUTINE ECL_TO_ISN (VEC_IN, DIRECTION, VEC_OUT) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: VEC_IN(3) ! input 3-vector INTEGER, INTENT(IN) :: DIRECTION ! conversion direction DOUBLE PRECISION, INTENT(OUT) :: VEC_OUT(3) ! output 3-vector DOUBLE PRECISION, PARAMETER :: PI = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803D0 DOUBLE PRECISION, PARAMETER :: APEX_LAMBDA = 252.0D0*PI/180.0D0 ! LISM solar apex right asc(rad) DOUBLE PRECISION, PARAMETER :: APEX_BETA = +7.0D0*PI/180.0D0 ! LISM solar apex declin (rad) DOUBLE PRECISION :: ALPHA, BETA DOUBLE PRECISION, DIMENSION(3,3) :: MATA, MATB ! ! Start of code ! ALPHA = APEX_LAMBDA BETA = -APEX_BETA MATA = RESHAPE ( (/ COS(ALPHA), SIN(ALPHA), 0.0D0, & ! calculate rotation matrix about z -SIN(ALPHA), COS(ALPHA), 0.0D0, & 0.0D0, 0.0D0, 1.0D0 /), & SHAPE=(/3,3/), ORDER=(/2,1/) ) MATB = RESHAPE ( (/ COS(BETA), 0.0D0, -SIN(BETA), & ! calculate rotation matrix about y 0.0D0, 1.0D0, 0.0D0, & SIN(BETA) , 0.0D0, COS(BETA) /), & SHAPE=(/3,3/), ORDER=(/2,1/) ) IF (DIRECTION > 0) THEN VEC_OUT = MATMUL(MATA,VEC_IN) VEC_OUT = MATMUL(MATB,VEC_OUT) ELSE VEC_OUT = MATMUL(TRANSPOSE(MATB),VEC_IN) VEC_OUT = MATMUL(TRANSPOSE(MATA),VEC_OUT) END IF RETURN END SUBROUTINE ECL_TO_ISN !*********************************************************************************************************************************** ! EQ_TO_ECL ! ! Convert input vector from Earth equatorial (EQ) coordinates to ecliptic (ECL) coordinates. !*********************************************************************************************************************************** SUBROUTINE EQ_TO_ECL (VEC_IN, DIRECTION, VEC_OUT) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: VEC_IN(3) ! input 3-vector INTEGER, INTENT(IN) :: DIRECTION ! conversion direction DOUBLE PRECISION, INTENT(OUT) :: VEC_OUT(3) ! output 3-vector DOUBLE PRECISION, PARAMETER :: PI = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899862803D0 ! ! Variable declarations ! DOUBLE PRECISION, PARAMETER :: EPS = 23.4392911D0 * PI/180.0D0 ! obliquity of ecliptic (J2000) DOUBLE PRECISION :: MATA(3,3) ! rotation matrix DOUBLE PRECISION :: ALPHA ! rotation angle between EQ and ECL frames ! ! Start of code ! ALPHA = EPS ! find frame rotation angle IF (DIRECTION < 0) ALPHA = -ALPHA ! change sign if reverse conversion MATA = RESHAPE ( (/ 1.0D0, 0.0D0, 0.0D0, & ! calculate rotation matrix 0.0D0, COS(ALPHA), SIN(ALPHA), & 0.0D0, -SIN(ALPHA), COS(ALPHA) /), & SHAPE=(/3,3/), ORDER=(/2,1/) ) VEC_OUT = MATMUL(MATA, VEC_IN) ! calculate output 3-vector RETURN END SUBROUTINE EQ_TO_ECL !*********************************************************************************************************************************** ! CALC_ANOMALIES !*********************************************************************************************************************************** SUBROUTINE CALC_ANOMALIES (T1, ECC, TRUE_ANOM, R) USE DENSITY_GLOBAL IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: T1 DOUBLE PRECISION, INTENT(IN) :: ECC DOUBLE PRECISION, INTENT(OUT) :: TRUE_ANOM DOUBLE PRECISION, INTENT(OUT) :: R DOUBLE PRECISION, PARAMETER :: HYPERBOLIC_THRESHOLD = 1.0D-6 DOUBLE PRECISION, PARAMETER :: R1 = 180.0D0/PI DOUBLE PRECISION, PARAMETER :: K = 0.01720209895D0 DOUBLE PRECISION, PARAMETER :: D1 = 1.0D4 DOUBLE PRECISION, PARAMETER :: C = 0.333333333333333333333333333333333333333333333333333333333333333333333333333333333333333D0 DOUBLE PRECISION, PARAMETER :: D = 1.0D-9 DOUBLE PRECISION, PARAMETER :: BETA(0:7) = (/ & 1.0D0, & 0.0D0, & -3.0D0/175.0D0, & -2.0D0/525.0D0, & -3.72D2/3.36875D5, & -8.044D3/2.1896875D7, & -1.81646D5/1.379503125D9, & -3.229596D6/6.5143203125D10 /) DOUBLE PRECISION, PARAMETER :: GAMMA(0:7) = (/ & 1.0D0, & 0.4D0, & 38.0D0/175.0D0, & 328.0D0/2625.0D0, & 2.4726D4/3.36875D5, & 1.14344D5/2.627625D6, & 5.9601184D7/2.299171875D9, & 4.32281228D8/2.7918515625D10 /) DOUBLE PRECISION :: MEAN_ANOM, F, Q DOUBLE PRECISION :: A1, B1, C1, AA, BB, CC, TANW2, TANV2 DOUBLE PRECISION :: COTS, COTS2, COTU, COT2U, AAJ INTEGER :: I, J DOUBLE PRECISION :: HYP_KEPLER, CUBEROOT IF ((ECC-1.0D0) > HYPERBOLIC_THRESHOLD) THEN ! hyperbolic orbit, eccentricity > threshold MEAN_ANOM = N*T1 F = HYP_KEPLER(MEAN_ANOM, ECC) TRUE_ANOM = 2.0D0 * ATAN(SQRT((ECC+1.0D0)/(ECC-1.0D0))*TANH(0.5D0*F)) R = A*(ECC*COSH(F)-1.0D0) ELSE IF (ECC > 1.0D0) THEN ! hyperbolic orbit, but near-parabolic Q = A*(ECC-1.0D0) A1 = SQRT(0.1D0*(1.0D0+9.0D0*ECC)) B1 = 5.0D0*(1.0D0-ECC)/(1.0D0+9.0D0*ECC) C1 = SQRT(5.0D0*(1.0D0+ECC)/(1.0D0+9.0D0*ECC)) BB = 1.0D0 DO I = 1, 25 COTS = 3.0D0*BB*A1*K*T1/SQRT(8.0D0*Q**3) COTS2 = SIGN(SQRT(1.0D0+COTS**2),COTS) + COTS COTU = CUBEROOT(COTS2) COT2U = (COTU**2-1.0D0)/(2.0D0*COTU) TANW2 = 2.0D0*COT2U AA = B1*TANW2**2 BB = 0.0D0 AAJ = 1.0D0 DO J = 0, 7 BB = BB + BETA(J) * AAJ AAJ = AAJ * AA END DO END DO CC = 0.0D0 AAJ = 1.0D0 DO J = 0, 7 CC = CC + GAMMA(J) * AAJ AAJ = AAJ * AA END DO TANV2 = C1*CC*TANW2 TRUE_ANOM = 2.0D0*ATAN(TANV2) R = (2.0D0*Q)/(1.0D0+COS(TRUE_ANOM)) ELSE ! non-hyperbolic orbit WRITE (UNIT=*, FMT='(A)') ' Error: non-hyperbolic orbit found in CALC_ANOMALIES' STOP END IF RETURN END SUBROUTINE CALC_ANOMALIES !*********************************************************************************************************************************** ! BARKER ! ! Solves Barker's Equation (direct method). ! The code below uses GM for the Sun. If the orbit is not heliocentric, GM should be replaced by the appropriate value; e.g. ! ! GM(SUN) = 1.32712438D20 m**3 s**-2 ! GM(EARTH) = 3.986005D14 m**3 s**-2 !*********************************************************************************************************************************** FUNCTION BARKER (TP, Q) RESULT (V) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: TP DOUBLE PRECISION, INTENT(IN) :: Q DOUBLE PRECISION :: V DOUBLE PRECISION, PARAMETER :: GM = 1.32712438D20 ! GM(sun) (m**3 s**-2) DOUBLE PRECISION :: COTS, COTS2, COTW, COT2W, TANV2 DOUBLE PRECISION :: CUBEROOT COTS = 3.0D0*SQRT(GM)*TP/SQRT(8.0D0*Q**3) COTS2 = SIGN(SQRT(1.0D0+COTS**2),COTS) + COTS COTW = CUBEROOT(COTS2) COT2W = (COTW**2-1.0D0)/(2.0D0*COTW) TANV2 = 2.0D0*COT2W V = 2.0D0*ATAN(TANV2) RETURN END FUNCTION BARKER !*********************************************************************************************************************************** ! CUBEROOT ! ! Computes the cube root. !*********************************************************************************************************************************** FUNCTION CUBEROOT (X) RESULT (Y) IMPLICIT NONE DOUBLE PRECISION, INTENT(IN) :: X DOUBLE PRECISION :: Y DOUBLE PRECISION, PARAMETER :: THIRD = 0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333D0 Y = SIGN((ABS(X))**THIRD,X) RETURN END FUNCTION CUBEROOT