C$Procedure CRM C SUBROUTINE CRM C C C******************************************************************************* C C Copyright (C) 1993, California Institute of Technology. U.S. C Government Sponsorhip under NASA Contract NAS7-918 is C acknowledged. C C******************************************************************************* C C$ Log C C Date Name Description C ----------------------------------------------------------------------------- C 12-JUL-1990 Eric Cannell creation of CRM C 7-MAY-1991 Bruce Shapiro invoke MEAN orbital density C 28-OCT-1991 B.E.S include error bars on solar/geomag data C C$ Purpose C C CRM is the main subroutine point for entry points CRMINI and CRMGET. C The user should NEVER call CRM, only the entry points CRMINI and CRMGET. C C CRMINI computed and saved daily ratios of: C C Cd * density C ------------ C mass C C where density is based upon the information in the FLUX_DATA file. This C can be done right after the FLUX_DATA file is read. I am computing this C ratio, rather than just density, in order to avoid having to pass CD and C MASS all over Hell. Note that the daily ratios values are stored in a C SAVEd array. C C CRMGET determines the ratio value for a particular epoch. Note that C daily values are somewhat of a misnomer, since density is constantly C changing. Also, the input DATE will likely be at some arbitrary time C during a particular day. CRMGET ignores the hour of the day and returns C the ratio for that day, regardless of the time of day. C C$ Input_Arguments C C Entry Point: CRMINI C Name Type Dim Units Description C ----------------------------------------------------------------------------- C DAYS I 1 - number of daily values C DAYONE C*(*) 1 - epoch ('dd-mmm-yyyy') of the first day C FLX DP DAYS --> daily 10.7 solar flux values in C 10**-22 watts/m**2/cycle/sec C FLXBAR DP DAYS --> daily 81 day moving average of FLX C in 10**-22 watts/m**2/cycle/sec C KP DP DAYS - geomagnetic Kp planetary index C CD DP 1 - constant of atmospheric drag on spacecraft C MASS DP 1 kg mass of spacecraft C DFLX DP DAYS - error bar on FLX C DFLXBAR DP DAYS - error bar on FLX C DKP DP DAYS - error bar on FLX C C C Entry Point: CRMGET C Name Type Dim Units Description C ----------------------------------------------------------------------------- C DATE C*(*) 1 - day of desired density, 'dd-mmm-yyyy' C C$ Output_Arguments C C Entry Point: CRMGET C Name Type Dim Units Description C ----------------------------------------------------------------------------- C CRMVAL DP 1 km**-3 Cd * density / MASS for DATE C CRMVALHI DP 1 km**-3 high value C CRMVALLO DP 1 km**-3 lo value C C$ Restrictions C C 1] If CRMGET is called with a DATE that is not within the time period C covered by the FLUX_DATA file, then an error message is printed and C GTARG is terminated. C C 2] Both CD and MASS must be > 0. C C$ References C C 1] Cannell, P.E., "JRSMPL- A Simple Fortran Function to Approximate C Atmospheric Density", JPL-IOM 314.9/90-481, 13 Jun 1990. C C$ Library_Links C C Entry Point Name Location C ----------------------------------------------------------------------------- C CAL2JD TIMETRANS C CH2JD TIMETRANS C jrsmpl2 jrsmpl2 C C$ Files C C File Name Unit Number Description C ----------------------------------------------------------------------------- C * * standard I/O C OFILE 8 text output file C C$ Parameters C INTEGER MXFLUX PARAMETER ( MXFLUX = 1000 ) C C$ Declarations_of_Input_and_Output_Arguments C C Input for CRMINI. C DOUBLE PRECISION CD INTEGER DAYS CHARACTER*(*) DAYONE DOUBLE PRECISION DFLX ( DAYS ) DOUBLE PRECISION DFLXBAR ( DAYS ) DOUBLE PRECISION DKP ( DAYS ) DOUBLE PRECISION FLX ( DAYS ) DOUBLE PRECISION FLXBAR ( DAYS ) DOUBLE PRECISION KP ( DAYS ) DOUBLE PRECISION MASS C Input for CRMGET. CHARACTER*(*) DATE DOUBLE PRECISION CRMVAL DOUBLE PRECISION CRMVALHI DOUBLE PRECISION CRMVALLO C C$ Declarations_of_Local_Variables C CHARACTER*11 DMY DOUBLE PRECISION LOFLUX, LOFBAR, LOKP DOUBLE PRECISION HIFLUX, HIFBAR, HIKP double precision DCRM_flux,DCRM_fbar,dCRM_kp LOGICAL ERROR INTEGER IDAY integer iday_now integer ihour_now integer iminute_now integer imonth_now integer isecond_now INTEGER iyear_now real*8 jd_jan0 real*8 jd_now real*8 jd_one real*8 frac_now real*8 fractional_year C C$ Save_Statements C INTEGER OFFSET SAVE OFFSET INTEGER MAXDAY SAVE MAXDAY DOUBLE PRECISION CRMDAT ( MXFLUX ) SAVE CRMDAT DOUBLE PRECISION CRMDATHI ( MXFLUX ) SAVE CRMDATHI DOUBLE PRECISION CRMDATLO ( MXFLUX ) SAVE CRMDATLO C C$ External_Statements C DOUBLE PRECISION CAL2JD EXTERNAL CAL2JD DOUBLE PRECISION CH2JD EXTERNAL CH2JD DOUBLE PRECISION JRSMPL2 EXTERNAL JRSMPL2 C C$ Method C-& C******************************************************************************* C******************************************************************************* C*** *** C*** Main Entry Point CRM *** C*** *** C******************************************************************************* C******************************************************************************* C1 The user must never call CRM, only the entry points CRMINI and CRMGET. WRITE(*,301) WRITE(8,301) 301 FORMAT(//,1X,'GTARG: Call CRMINI or CRMGET, but never call CRM.') STOP C******************************************************************************* C******************************************************************************* C*** *** C*** Entry Point CRMINI *** C*** *** C******************************************************************************* C******************************************************************************* C1 Entry point CRMINI: ENTRY CRMINI( DAYS , DAYONE , FLX , FLXBAR , KP , CD , MASS, & DFLX, DFLXBAR, DKP ) C2 Check that CD is > 0. IF ( CD .LE. 0D0 ) THEN WRITE(*,302) CD WRITE(8,302) CD 302 FORMAT(//,1X,'GTARG: CD (',D24.18,') must be > 0.') STOP END IF C2 Check that MASS is > 0. IF ( MASS .LE. 0D0 ) THEN WRITE(*,303) MASS WRITE(8,303) MASS 303 FORMAT(//,1X,'GTARG: MASS (',D24.18,') must be > 0.') STOP END IF C2===> 5/7/91 Mod B.E.S. C C Determine julian day at the start of the calculation C DMY = DAYONE jd_one = CH2JD(DMY,.TRUE.,ERROR) C2 Compute and save ratio of CD * rho / MASS. DO 101 IDAY = 1 , DAYS C C Determine the fraction of the current year C jd_now = jd_one + dfloat ( iday - 1 ) call jd2cal ( jd_now, & iyear_now, imonth_now, iday_now, ihour_now, & iminute_now, isecond_now, frac_now ) jd_jan0 = cal2jd (iyear_now, 1, 1, 0, 0, 0, 0.0) - 1.0 if ( mod (iyear_now, 4) .eq. 0 ) then fractional_year = (jd_now - jd_jan0)/366.0 else fractional_year = (jd_now - jd_jan0)/365.0 end if CD write (8,99) iday, iyear_now, imonth_now, iday_now, CD & jd_now-jd_jan0, fractional_year 99 format (1x, 4i10, 2g15.8) CRMDAT( IDAY ) = CD & * JRSMPL2( FLX ( IDAY ) , & FLXBAR( IDAY ) , & KP ( IDAY ) , & fractional_year & ) & / MASS C C calculate the CRM errors due to the sigmas in each of C the components F, FBAR, KP C HIFLUX = FLX ( IDAY )+DFLX(IDAY) HIFBAR = FLXBAR( IDAY )+DFLXBAR(IDAY) HIKP = KP (IDAY) + DKP(IDAY) DCRM_FLUX = CD * JRSMPL2( HIFLUX, FLXBAR(IDAY), KP(IDAY), & fractional_year) / MASS - CRMDAT(IDAY) DCRM_FBAR = CD * JRSMPL2( FLX(IDAY), HIFBAR, KP(IDAY) , & fractional_year) / MASS - CRMDAT(IDAY) DCRM_KP = CD * JRSMPL2( FLX(IDAY), FLXBAR(IDAY), HIKP, & fractional_year) / MASS - CRMDAT(IDAY) C C RSS the three density variations to get the total C density sigma C CRMDATHI( IDAY ) = CRMDAT(IDAY) + & SQRT ( DCRM_FLUX**2 + DCRM_FBAR**2 + DCRM_KP**2 ) C CRMDATHI( IDAY ) = CD C & * JRSMPL2( FLX ( IDAY )+DFLX(IDAY) , C & FLXBAR( IDAY )+DFLXBAR(IDAY) , C & KP ( IDAY )+DKP(IDAY) , C & fractional_year C & ) C & / MASS C C calculate the CRM errors due to the sigmas in each of C the components F, FBAR, KP in the low density C direction C LOFLUX = DMIN1( DMAX1( FLX(IDAY)-DFLX(IDAY), 70.0d0), & FLX(IDAY) ) LOFBAR = DMIN1( DMAX1( FLXBAR(IDAY)-DFLXBAR(IDAY), 70.0d0 ), & FLXBAR(IDAY)) LOKP = DMIN1( DMAX1( KP(IDAY)-DKP(IDAY), 1.0d0 ), KP(IDAY)) C CRMDATLO( IDAY ) = CD C & * JRSMPL2( LOFLUX, C & LOFBAR, LOKP, C & fractional_year C & ) C & / MASS DCRM_FLUX = CD * JRSMPL2( LOFLUX, FLXBAR(IDAY), KP(IDAY), & fractional_year) / MASS - CRMDAT(IDAY) DCRM_FBAR = CD * JRSMPL2( FLX(IDAY), LOFBAR, KP(IDAY) , & fractional_year) / MASS - CRMDAT(IDAY) DCRM_KP = CD * JRSMPL2( FLX(IDAY), FLXBAR(IDAY), LOKP, & fractional_year) / MASS - CRMDAT(IDAY) C C RSS the three density variations to get the total C density sigma C CRMDATLO( IDAY ) = CRMDAT(IDAY) - & SQRT ( DCRM_FLUX**2 + DCRM_FBAR**2 + DCRM_KP**2 ) 101 CONTINUE C2=====> end mod 5/7/91 C2 Determine and save the Julian date offset into the data array C2 so that for any other date: index = date - offset. DMY = DAYONE OFFSET = IDNINT( CH2JD(DMY,.TRUE.,ERROR) - .5D0 ) C2 Save the number of days in the date array. MAXDAY = DAYS C2 Return from CRMINI. RETURN C******************************************************************************* C******************************************************************************* C*** *** C*** Entry Point CRMGET *** C*** *** C******************************************************************************* C******************************************************************************* C1 Entry point CRMGET: ENTRY CRMGET( DATE , CRMVAL, CRMVALHI, CRMVALLO ) C2 Determine the index into the ratio array. DMY = DATE IDAY = IDNINT( CH2JD(DMY,.TRUE.,ERROR) + .5D0 ) - OFFSET C2 Check that DATE is within the time period covered by the data. IF ( IDAY .LT. 1 .OR. MAXDAY .LT. IDAY ) THEN WRITE(8,304) DMY , IDAY WRITE(*,304) DMY , IDAY 304 FORMAT(//,1X,'GTARG: FLUX_DATA does not cover "',A11, & '" (index ',I5,').') STOP END IF C2 Access ratio from saved data array. CRMVAL = CRMDAT( IDAY ) CRMVALHI = CRMDATHI( IDAY ) CRMVALLO = CRMDATLO( IDAY ) C2 Return from CRMGET. RETURN C1 End of CRM. END