C =============================================================================

C EXOSPHERE DENSITY SUBROUTINE *** COMBINED *** (14 Sept 93)
C THIS CODE COMPUTES TOTAL DENSITY USING EITHER THE ORIGINAL
C CHAMBERLAIN [1963] ANALYTIC EXOSPHERE FORMULATION OR THE
C GEOCORONAL ANALYTIC FORMULATION OF BISHOP [1991]:
C       SPHERICALLY UNIFORM EXOBASE TEMPERATURE & DENSITY
C       EXOPAUSE HANDLING OF SATELLITE COMPONENT
C       NO ROTATION OR WINDS
C INCLUDES BOLTZMANN (OR "BAROMETRIC") FACTOR.
C SELECTION FLAG IS IGEO:
C     0 . . . ORIGINAL CHAMBERLAIN MODEL WITH SATELLITE CRITICAL
C             RADIUS, PASSED IN COMMON BLOCK AS FTSAT IN UNITS OF
C             PLANETARY RADIUS
C     1 . . . ANALYTIC MODEL INCLUDING EXOPAUSE, "EVAPORATIVE" SATELLITE
C             POPULATION PARAMETERS, AND SOLAR IONIZATION DECAY.
C
C [FEBRUARY 2002]
C EVALUATION OF KINETIC TEMPERATURE ADDED.  SOME RESTRUCTURING DONE.

      SUBROUTINE CORONA(RR, DKIN, TKIN)

      IMPLICIT REAL*8 (A-H,O-Z)

      PARAMETER (GG   = 0.8862269255D0)           !!  $\Gamma(3/2)$
      PARAMETER (FF   = 0.1329340388D1)           !!  $\Gamma(5/2)$

      COMMON /EXOS/ GM,CONL,PLANETR,RBASE,RC,RUPR,RP,TEXO,DEXO,VELT,
     &              RADPF,CENTER,ABSCSX,FTSAT,FDSAT,IGEO
      COMMON /NMBR/ PI,RTPI,PID2,OFFSET

      DATA TSOL   / 1.10D6 /                      !!  SOLAR IONIZATION LIFETIME
 
C -----------------------------------------------------------------------------
C SOME BASIC STUFF & NOTES
C -----------------------------------------------------------------------------
C notes:
C   density:      in case of analytic [1991] model, the direct escape component
C                 is included in evaluation of net densities but not the
C                 "extended escape" component
C   temperature:  escape component quantity is calculated but not used.
C                 TKIN determined in terms of bound population only.
C   temperature:  no real kinetic distinction between satellite and ballistic
C                 (or escape) components in Chamberlain 1963 model;  however,
C                 temperature profile evaluated in terms of mixed populations.
C
C ___ densities
C      GX = gamma(3/2,X)
C      ZETA_1 . . . ballistic partition component
C      ZETA_2 . . . satellite partition component
C      ZETA_B . . . bound partition component
C      ZETA_E . . . escape partition component
C
C ___ kinetic temperature factors
C      FX = gamma(5/2,X) = (3/2)*gamma(3/2,X) - X^(3/2) * EXP(-X)
C      ETA_1 . . . ballistic partition component
C      ETA_2 . . . satellite partition component
C      ETA_B . . . bound partition component
C      ETA_E . . . escape partition component
C
C TEMPERATURE OF MIXED POPULATIONS (W/O BULK FLOW SPEEDS)
C      T_mix = f1^2*T1 + f2^2*T2 + f1*f2*(T1 + T2)
C where f1, f2 are mixing ratios

C BASIC RADII PARAMETERS

      AC = CONL/(RC*TEXO)
      AA = CONL/(RR*TEXO)

C -----------------------------------------------------------------------------
C CHECK FOR ERRORS IN CALL
C -----------------------------------------------------------------------------

      IF ( (IGEO.NE.0) .AND. (IGEO.NE.1) ) THEN
        DKIN = 0.0D0
        TKIN = -1.0D0
C        PRINT, 'UNRECOGNIZED EXOSPHERE MODEL IDENTIFIER: ',
C     &         ' IGEO MUST BE 0 OR 1 (WARNING VALUES RETURNED)'
        WRITE(06,*), 'UNRECOGNIZED EXOSPHERE MODEL IDENTIFIER: ',
     &         ' IGEO MUST BE 0 OR 1 (WARNING VALUES RETURNED)'

        GO TO 909

      ELSE IF (AA .GE. (AC*(1.0D0 - OFFSET))) THEN
        DKIN = DEXO * EXP(AA-AC)
        TKIN = TEXO
C        PRINT, 'NOT IN EXOSPHERE -- BELOW EXOBASE ', 
C     &         ' (SIMPLE DEFAULT VALUES RETURNED): ',
C     &         RR, DKIN, TKIN
        WRITE(06,*), 'NOT IN EXOSPHERE -- BELOW EXOBASE ', 
     &         ' (SIMPLE DEFAULT VALUES RETURNED): ',
     &         RR, DKIN, TKIN

        GO TO 909

      ELSE IF (AA .LE. 0.0D0) THEN
        DKIN = 0.0D0
        TKIN = -1.0D0
C        PRINT, 'NOT IN EXOSPHERE -- OUTSIDE EXOPAUSE ',
C     &         ' (WARNING VALUES RETURNED)'
        WRITE(06,*), 'NOT IN EXOSPHERE -- OUTSIDE EXOPAUSE ',
     &         ' (WARNING VALUES RETURNED)'

        GO TO 909

      END IF

C -----------------------------------------------------------------------------
C CHAMBERLAIN 1963
C -----------------------------------------------------------------------------

      IF (IGEO.EQ.0) THEN
 
        RS = FTSAT
        AS = CONL/(RS*TEXO)

        Y1 = (AA**2) / (AC+AA)
        A1 = AA - Y1
        GA = GAMMA(AA)
        G1 = GAMMA(A1)

        CCOEF = SQRT(AC**2 - AA**2) * EXP(-Y1) / AC

C ___ densities

        ZETA_1 = (GA - CCOEF*G1) * 2.0D0/RTPI
        ZETA_E = (1.0D0 - CCOEF) * GG/RTPI - ZETA_1/2.0D0

        IF (AA.LE.AS) THEN
          YS     = (AA**2) / (AS+AA)
          SY     = AA - YS
          GS     = GAMMA(SY)
          SCOEF  = SQRT(AS**2 - AA**2) * EXP(-YS) / AS
          ZETA_B = (GA - SCOEF*GS) * 2.0D0/RTPI
        ELSE
          ZETA_B = 2.0D0*GA/RTPI
        END IF

        ZETA_2 = ZETA_B - ZETA_1

        PART_1 = ZETA_1 * EXP(AA-AC)                 ! BALLISTICS
        PART_2 = ZETA_2 * EXP(AA-AC)                 ! SATELLITES
        PART_E = ZETA_E * EXP(AA-AC)                 ! ESCAPING

        DKIN = DEXO * (PART_1 + PART_2 + PART_E)

C ___ kinetic temperatures

        FA = 1.5D0*GA - SQRT(AA**3) * EXP(-AA)
        F1 = 1.5D0*G1 - SQRT(A1**3) * EXP(-A1)

        ETA_1 = (FA - CCOEF * (Y1*G1 + F1)) * 2.0D0/RTPI
        ETA_E = (1.0D0-CCOEF)*FF/RTPI - CCOEF*Y1*GG/RTPI - ETA_1/2.0D0

        IF (AA.LE.AS) THEN
          ETA_B = (FA - SCOEF * (YS*G1 + F1)) * 2.0D0/RTPI
        ELSE
          ETA_B = 2.0D0*FA/RTPI
        END IF

        ETA_2 = ETA_B - ETA_1

        T_1 = 0.66667D0 * TEXO * ETA_1 / ZETA_1
        T_2 = 0.66667D0 * TEXO * ETA_2 / ZETA_2

        FRAC_1 = PART_1 / (PART_1 + PART_2)
        FRAC_2 = PART_2 / (PART_1 + PART_2)

        TKIN =   (FRAC_1 ** 2) * T_1
     &         + (FRAC_2 ** 2) * T_2
     &         + FRAC_1 * FRAC_2 * (T_1 + T_2)

C -----------------------------------------------------------------------------
C BISHOP 1991
C -----------------------------------------------------------------------------

      ELSE IF (IGEO.EQ.1) THEN

        AR = CONL/(RP*TEXO)

        Y1 = (AA**2) / (AC+AA)
        YA = AA - AR
        YB = AA - AC*AR/(AC+AR)
        YE = (AA**2) / (AA+AR)

        BS = (YB - Y1) / FTSAT
        ES = (YE - YB) / FTSAT

        GA = GAMMA(YA)
        GB = GAMMA(YB-Y1)
        GS = GAMMA(BS)

        ZA = ZAMMA(YE-YA)
        ZB = ZAMMA(YE-YB)
        ZS = ZAMMA(ES)

        CCOEF = SQRT(AC**2 - AA**2) / AC
        ACOEF = SQRT(AA**2 - AR**2) / AR

C ___ densities

        ZETA_1 = GA
     &           - CCOEF * EXP(-Y1) * GB
     &           + ACOEF * EXP(-YE) * (ZA-ZB)
        ZETA_1 = ZETA_1 * 2.0D0 / RTPI

        ZETA_E = GG - GA
     &           - CCOEF * EXP(-Y1) * (GG-GB)
     &           - ACOEF * EXP(-YE) * (ZA-ZB)
        ZETA_E = ZETA_E / RTPI

        ZETA_2 =   CCOEF * EXP(-Y1/FTSAT) * GS
     &           + ACOEF * EXP(-YE/FTSAT) * ZS
        ZETA_2 = ZETA_2 * 2.0D0 / RTPI

        DK    = 2.97D6 * ASIN(1.0D0 - RC/RR) / SQRT(RR) / RADPF
        ZETA_2 = ZETA_2 * EXP(-DK/TSOL)

        PART_1 = ZETA_1 * EXP(AA-AC)
        PART_E = ZETA_E * EXP(AA-AC)
        PART_2 = FDSAT * ZETA_2 * EXP((AA-AC)/FTSAT)

        DKIN  = DEXO * (PART_1 + PART_2 + PART_E)

C ___ kinetic temperature factors

        FA = 1.50D0*GA - SQRT(YA**3)      * EXP(-YA)
        FB = 1.50D0*GB - SQRT((YB-Y1)**3) * EXP(-(YB-Y1))
        FS = 1.50D0*GS - SQRT(BS**3)      * EXP(-BS)

        VA = SQRT((YE-YA)**3) * EXP(YE-YA) - 1.50D0*ZA
        VB = SQRT((YE-YB)**3) * EXP(YE-YB) - 1.50D0*ZB
        VS = SQRT(ES**3)      * EXP(ES)    - 1.50D0*ZS

        ETA_1 = FA
     &          - CCOEF * EXP(-Y1) * (Y1*GB + FB)
     &          + ACOEF * EXP(-YE) * (YE*(ZA-ZB) - (VA-VB))
        ETA_1 = ETA_1 * 2.0D0 / RTPI

        ETA_E = FF - FA
     &          - CCOEF * EXP(-Y1) * (FF + Y1*(GG-GB) - FB)
     &          - ACOEF * EXP(-YE) * (YE*(ZA-ZB) - (VA-VB))
        ETA_E = ETA_E / RTPI

        ETA_2 =   CCOEF * EXP(-Y1/FTSAT) * (Y1*GS + FTSAT*FS)
     &          + ACOEF * EXP(-YE/FTSAT) * (YE*ZS - FTSAT*VS)
        ETA_2 = ETA_2 * 2.0D0 / RTPI

        T_1 = 0.66667D0 * TEXO * ETA_1 / ZETA_1
        T_2 = 0.66667D0 * TEXO * ETA_2 / ZETA_2

        FRAC_1 = PART_1 / (PART_1 + PART_2)
        FRAC_2 = PART_2 / (PART_1 + PART_2)

        TKIN =   (FRAC_1 ** 2) * T_1
     &         + (FRAC_2 ** 2) * T_2
     &         + FRAC_1 * FRAC_2 * (T_1 + T_2)

      END IF

C -----------------------------------------------------------------------------

 909  CONTINUE

      RETURN
      END

C =============================================================================
C =============================================================================
 
C  EVALUATION OF THE INCOMPLETE GAMMA FUNCTION

      FUNCTION GAMMA(AAA)

      IMPLICIT REAL*8 (A-H,O-Z)

      COMMON / GAUS08 / WW08(08),XX08(08)

C -----------------------------------------------------------------------------
C INTEGRATION ROUTINE
C -----------------------------------------------------------------------------

      ROOT = SQRT(AAA)

      ERF  = 0.0D0
      DO JJ = 1,8
        XPNT= ((1.0D0+XX08(JJ))*ROOT/2.0D0)**2
        ERF = ERF + WW08(JJ)*EXP(-XPNT)*ROOT/2.0D0
      END DO

      GAMMA = ERF - ROOT*EXP(-AAA)

C -----------------------------------------------------------------------------

      RETURN
      END

C =============================================================================
C =============================================================================
 
C  EVALUATION OF THE ZAMMA FUNCTION

      FUNCTION ZAMMA(AAA)

      IMPLICIT REAL*8 (A-H,O-Z)

      COMMON / GAUS08 / WW08(08),XX08(08)

C -----------------------------------------------------------------------------
C INTEGRATION ROUTINE
C -----------------------------------------------------------------------------

      ZAMMA = 0.0D0
      DO JJ = 1,8
        XXX   = (1.0D0+XX08(JJ))*AAA/2.0D0
        ZAMMA = ZAMMA + WW08(JJ)*EXP(XXX)*SQRT(XXX)*AAA/2.0D0
      END DO

C -----------------------------------------------------------------------------

      RETURN
      END

C =============================================================================