Patent Publication Number: US-7593835-B2

Title: Reformulated atmospheric band model method for modeling atmospheric propagation at arbitrarily fine spectral resolution and expanded capabilities.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
   This application is a continuation in part of parent application Ser. No. 09/838,801 filed Apr. 20, 2001, now U.S. Pat. No. 7,433,806. Priority is claimed. The disclosure of the parent application is incorporated by reference herein. This application also claims priority of Provisional application Ser. No. 60/668,159 filed on Apr. 4, 2005 

   STATEMENT REGARDING FEDERALLY SPONSORED RESEARCH 
   This invention was made with government support under contracts F19628-02-C-0078 and FA8718-04-C-0073 awarded by the Department of the Air Force. The government has certain rights in this invention. 

   FIELD OF THE INVENTION 
   This invention relates to a radiative transport band model for application to high spectral resolution simulations. 
   BACKGROUND OF THE INVENTION 
   Calculation of in-band radiances via the line-by-line (LBL) approach is computationally intensive. Not only does it require the calculation of spectral molecular absorption coefficients for the series of pressure and temperature pairs encountered in the specified atmosphere, but it also requires the radiation transport equation (RTE) be solved for each spectral point. Because of the fine structure of molecular absorption, especially at higher altitudes (lower pressures), a very finely spaced spectral grid is required. Spectral step size depends upon scenario and accuracy requirements, but values of between 0.01 and 0.00001 cm −1  is common. 
   The MODTRAN™ solution to this problem is to use a narrow statistical band model. In this approach, the RTE is solved for a relatively narrow spectral band all at once. A number of physical and statistical approximations are made. The quality of the MODTRAN™ results depend largely upon the validity of these approximations for any specific scenario. Over the past 30 years of MODTRAN™ development, much effort has been focused on fine tuning the model to improve accuracy for the terrestrial atmosphere. An earlier version of MODTRAN is described in U.S. Pat. No. 5,884,226. 
   The integral RTE defines the line-of-sight (LOS) monochromatic diffuse radiant intensity I v (Ω 0 ;0) at spectral frequency v along a refracted path in the direction specified by solid angle Ω 0 : 
                     I   v     ⁡     (       Ω   0     ;   0     )       =         ⅇ     -       τ   v     ⁡     (   s   )           ⁢       I   v     ⁡     (       Ω   s     ;   s     )         +       ∫   0       τ   v     ⁡     (   s   )         ⁢       ⅇ     -       τ   v     ⁡     (     s   ′     )           ⁢       J   v     ⁡     (       Ω     s   ′       ;     s   ′       )       ⁢           ⁢       ⅆ       τ   v     ⁡     (     s   ′     )         .                   (   1   )               
This equation states that the observed spectral radiant intensity at the sensor (s=0) equals the spectral radiant intensity directed towards the sensor a distance s away along the refracted path attenuated by the foreground path transmittance e −τ     v     (s)  plus the integral of intervening source radiation J v (Ω s ′;s′) attenuated by its foreground transmittance e −τ     v     (s′) . The path spectral transmittances in Eq. (1) are defined in terms of dimensionless optical depths, τ v (s). The optical depth terms themselves are path integrals over extinction coefficients, b v (s′) [units of reciprocal length].
 
                     τ   v     ⁡     (   s   )       ≡       ∫   0   s     ⁢         b   v     ⁡     (     s   ′     )       ⁢           ⁢     ⅆ     s   ′                   (   2   )               
The extinction coefficients consists of absorption b v   abs (s) and scattering b v   sct (s) components, which are computed as sums over atmospheric molecular and particulate species contributions b v,i (s).
 
                           b   v     ⁡     (   s   )       ≡       ⁢         b   v   abs     ⁡     (   s   )       +       b   v   sct     ⁡     (   s   )                     ≡       ⁢       ∑     i   ⁢           ⁢   species       ⁢           ⁢     [         b     v   ,   i     abs     ⁡     (   s   )       +       b     v   ,   i     sct     ⁡     (   s   )         ]                   ≡       ⁢       ∑     i   ⁢           ⁢   species       ⁢           ⁢       [         κ     v   ,   i     abs     ⁡     (   s   )       +       κ     v   ,   i     sct     ⁡     (   s   )         ]     ⁢       ρ   i     ⁡     (   s   )                         (   3   )               
The individual species absorption and scattering coefficients are themselves the product of a cross-section κ v,i (s) [units of area] and a density ρ i (s) [number per unit volume].
 
   The atmospheric source term J v (Ω s ;s) consists of three distinct components: local thermal emission, single scattered directly transmitted solar (and/or lunar) irradiance and diffuse radiation scattered into the LOS 
                           J   v     ⁡     (       Ω   s     ;   s     )       =           b   v   abs         b   v     ⁡     (   s   )         ⁢       B   v     ⁡     (     T   s     )         ⁢           +             local   ⁢           ⁢   thermal   ⁢           ⁢   emission                     b   v   sct         b   v     ⁡     (   s   )         ⁢       p   v     ⁡     (       Ω   s     ,       Ω   s   sun     ;   s       )       ⁢     ⅇ     -       τ   v   sun     ⁡     (   s   )           ⁢     F   v   sun       ⁢           +           single   ⁢           ⁢   scatter   ⁢           ⁢   solar                     b   v   sct     ⁡     (   s   )           b   v     ⁡     (   s   )         ⁢       ∫     4   ⁢   π               ⁢         p   v     ⁡     (       Ω   s     ,       Ω   ′     ;   s       )       ⁢       I   v     ⁡     (       Ω   ′     ;   s     )       ⁢           ⁢     ⅆ     Ω   ′                   multiple   ⁢           ⁢     scatter   .                   (   4   )               
The first term, the thermal emission, is the product of the Planck blackbody function B v (T s ) at local temperature T s  and the local emissivity, which is represented as a ratio of absorption to extinction coefficients. This expression for thermal emission assumes local thermodynamic equilibrium (LTE) conditions—an assumption, valid in the lower terrestrial atmosphere, that molecular collisions occur in sufficient frequency to control upper state population densities.
 
   The latter two terms in Eq. (4) are the path scattered radiance. The single scatter solar radiation is the product of the top-of-atmosphere (TOA) solar irradiance F v   sun , the scattering point to sun transmittance e −τ     v       sun     (s) , the effective scattering phase function p v (Ω s ,Ω s   sun ;s), and the single scattering albedo, equal to the ratio of scattering to extinction coefficients. The multiply scattered radiation is the integral over all directions (4π steradians) of the incoming diffuse radiation scattered into the LOS. The dependence of the multiple scattering on path radiance itself significantly complicates the radiation transport equation, transforming it into a fully coupled problem. MODTRAN™ uses the DISORT model (below), a discrete ordinate algorithm, to solve the coupled equations. The effective phase function that appears in both scattered radiance source terms is a scattering coefficient weighted average of individual species i scattering phase functions p v,i (Ω s ,Ω s ′;s): 
                     p   v     ⁡     (       Ω   s     ,       Ω   s   ′     ;   s       )       ≡         ∑     i   ⁢           ⁢   species       ⁢           ⁢         b     v   ,   i     sct     ⁡     (   s   )       ⁢       p     v   ,   i       ⁡     (       Ω   s     ,       Ω   s   ′     ;   s       )               b   v   sct     ⁡     (   s   )                 (   5   )               
The phase functions of Eq. (4) are all unit normalized, with dimensions of reciprocal steradians (1/sr). It is somewhat more common to define phase functions as dimensionless with a 4π sr normalization, but MODTRAN™ uses the unit normalization to eliminate all constant factors from the source function equation.
 
   Substituting the source term of Eq. (4) into the integral RTE, the spectral radiant intensity of Eq. (1) can be rewritten as absorption optical depth τ v   abs (s) and scattering optical depth τ v   sct (s) integrals: 
                         I   v     ⁡     (       Ω   0     ;   0     )       =         ⅇ     -       τ   v     ⁡     (   s   )           ⁢       I   v     ⁡     (       Ω   s     ;   s     )         +       ∫   0       τ   v   abs     ⁡     (   s   )         ⁢       ⅇ     -       τ   v     ⁡     (     s   ′     )           ⁢       B   v     ⁢           (     T     s   ′       )     ⁢     ⅆ       τ   v   abs     ⁡     (     s   ′     )             +       ∫   0       τ   v   sct     ⁡     (   s   )         ⁢         ⅇ     -       τ   v     ⁡     (     s   ′     )           ⁡     [                 p   v     ⁡     (       Ω     s   ′       ,       Ω     s   ′     sun     ;     s   ′         )       ⁢     ⅇ     -       τ   v   sun     ⁡     (     s   ′     )           ⁢     F   v   sun       +                 ∫     4   ⁢   π               ⁢         I   v     ⁡     (       Ω   ′     ;     s   ′       )       ⁢       p   v     ⁡     (       Ω     s   ′       ,       Ω   ′     ;     s   ′         )       ⁢           ⁢     ⅆ     Ω   ′                 ]       ⁢           ⁢     ⅆ       τ   v   sct     ⁡     (     s   ′     )                 ;     ⁢     
     ⁢         τ   v   abs     ⁡     (   s   )       ≡       ∫   0   s     ⁢         b   v   abs     ⁡     (     s   ′     )       ⁢           ⁢     ⅆ     s   ′       ⁢           ⁢   and   ⁢           ⁢       τ   v   sct     ⁡     (   s   )           ≡       ∫   0   s     ⁢         b   v   sct     ⁡     (     s   ′     )       ⁢           ⁢       ⅆ     s   ′       .                   (   6   )               
For band model theory, the key atmospheric attenuation parameter is the transmittance itself, not the optical depth. For that reason, it is convenient to again rewrite the RTE in terms of transmittance, t v (s)=exp [−τ v (s)], integrals:
 
   
     
       
         
           
             
               
                 
                   
                     I 
                     v 
                   
                   ⁡ 
                   
                     ( 
                     
                       
                         Ω 
                         0 
                       
                       ; 
                       0 
                     
                     ) 
                   
                 
                 = 
                 
                   
                     
                       
                         t 
                         v 
                       
                       ⁡ 
                       
                         ( 
                         s 
                         ) 
                       
                     
                     ⁢ 
                     
                       
                         I 
                         v 
                       
                       ⁡ 
                       
                         ( 
                         
                           
                             Ω 
                             s 
                           
                           ; 
                           s 
                         
                         ) 
                       
                     
                   
                   + 
                   
                     
                       ∫ 
                       
                         
                           t 
                           v 
                           abs 
                         
                         ⁡ 
                         
                           ( 
                           s 
                           ) 
                         
                       
                       1 
                     
                     ⁢ 
                     
                       
                         
                           t 
                           v 
                           sct 
                         
                         ⁡ 
                         
                           ( 
                           
                             s 
                             ′ 
                           
                           ) 
                         
                       
                       ⁢ 
                       
                         
                           B 
                           v 
                         
                         ⁡ 
                         
                           ( 
                           
                             T 
                             
                               s 
                               ′ 
                             
                           
                           ) 
                         
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         ⅆ 
                         
                           
                             t 
                             v 
                             abs 
                           
                           ⁡ 
                           
                             ( 
                             
                               s 
                               ′ 
                             
                             ) 
                           
                         
                       
                     
                   
                   + 
                   
                     
                       ∫ 
                       
                         
                           t 
                           v 
                           sct 
                         
                         ⁡ 
                         
                           ( 
                           s 
                           ) 
                         
                       
                       1 
                     
                     ⁢ 
                     
                       
                         
                           
                             t 
                             v 
                             abs 
                           
                           ⁡ 
                           
                             ( 
                             
                               s 
                               ′ 
                             
                             ) 
                           
                         
                         ⁡ 
                         
                           [ 
                           
                             
                               
                                 
                                   
                                     
                                       
                                         p 
                                         v 
                                       
                                       ⁡ 
                                       
                                         ( 
                                         
                                           
                                             Ω 
                                             
                                               s 
                                               ′ 
                                             
                                           
                                           , 
                                           
                                             
                                               Ω 
                                               
                                                 s 
                                                 ′ 
                                               
                                               sun 
                                             
                                             ; 
                                             
                                               s 
                                               ′ 
                                             
                                           
                                         
                                         ) 
                                       
                                     
                                     ⁢ 
                                     
                                       
                                         t 
                                         v 
                                         sun 
                                       
                                       ⁡ 
                                       
                                         ( 
                                         
                                           s 
                                           ′ 
                                         
                                         ) 
                                       
                                     
                                     ⁢ 
                                     
                                       F 
                                       v 
                                       sun 
                                     
                                   
                                   + 
                                 
                               
                             
                             
                               
                                 
                                   
                                     ∫ 
                                     
                                       4 
                                       ⁢ 
                                       π 
                                     
                                     
                                         
                                     
                                   
                                   ⁢ 
                                   
                                     
                                       
                                         I 
                                         v 
                                       
                                       ⁡ 
                                       
                                         ( 
                                         
                                           
                                             Ω 
                                             ′ 
                                           
                                           ; 
                                           
                                             s 
                                             ′ 
                                           
                                         
                                         ) 
                                       
                                     
                                     ⁢ 
                                     
                                       
                                         p 
                                         v 
                                       
                                       ⁡ 
                                       
                                         ( 
                                         
                                           
                                             Ω 
                                             
                                               s 
                                               ′ 
                                             
                                           
                                           , 
                                           
                                             
                                               Ω 
                                               ′ 
                                             
                                             ; 
                                             
                                               s 
                                               ′ 
                                             
                                           
                                         
                                         ) 
                                       
                                     
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     
                                       ⅆ 
                                       
                                         Ω 
                                         ′ 
                                       
                                     
                                   
                                 
                               
                             
                           
                           ] 
                         
                       
                       ⁢ 
                       
                         
                           ⅆ 
                           
                             
                               t 
                               v 
                               sct 
                             
                             ⁡ 
                             
                               ( 
                               
                                 s 
                                 ′ 
                               
                               ) 
                             
                           
                         
                         . 
                       
                     
                   
                 
               
             
             
               
                 ( 
                 7 
                 ) 
               
             
           
         
       
     
   
   Within MODTRAN™, the boundary radiant intensity term I v (Ω s ;s) for the full LOS path is typically zero unless the path terminates at the spherical Earth surface. There is an option to include grey body emission with no reflected source at the path final altitude. For a down-looking sensor whose LOS intersects the ground, the surface radiance includes both emissive and reflective components: 
                           I   v     ⁡     (       Ω   s     ;   s     )       =           ɛ   v     ⁡     (       Ω   s     ;   s     )       ⁢       B   v     ⁡     (     T   g     )         +             Surface   ⁢           ⁢   emission                   μ   s   sun     ⁢       f   v     ⁡     (       Ω   s     ,       Ω   s   sun     ;   s       )       ⁢       t   v   sun     ⁡     (   s   )       ⁢     F   v   sun       ⁢           +                 Reflected   ⁢           ⁢   directly                 transmitted   ⁢           ⁢   solar     ⁢                               ⁢   irradiance                       ∫     2   ⁢   π               ⁢         I   v     ⁡     (       Ω   ′     ;   s     )       ⁢       f   v     ⁡     (       Ω   s     ,       Ω   ′     ;   s       )       ⁢     μ   ′     ⁢     ⅆ     Ω   ′                         ⁢     Reflected   ⁢           ⁢   downward   ⁢           ⁢   flux                   (   8   )               
In this equation, μ′ is the cosine of the zenith angle of incoming radiation measured from the ground, f v (Ω s ,Ω′;s) is the surface bi-direction reflectance distribution function (BRDF) in units of reciprocal steradians at ground location s, ε v (Ω s ;s) is the dimensionless directional emissivity in the LOS path direction at s, and the surface reflected downward flux is computed as the upper hemisphere integral (2π sr) over incoming diffuse radiation. The BRDF defines the reflected radiance resulting from surface illumination as a function of both incoming and outgoing angles, and it is normalized to the surface albedo α v (s):
 
                     α   v     ⁡     (   s   )       =       1   π     ⁢       ∫     2   ⁢   π               ⁢       ∫     2   ⁢   π               ⁢         f   v     ⁡     (     Ω   ,       Ω   ′     ;   s       )       ⁢     μ   ′     ⁢     ⅆ     Ω   ′       ⁢   μ   ⁢     ⅆ   Ω                     (   9   )               
The directional emissivity, the complement of the directional reflectivity ρ v (Ω s ), dictates the angular distribution of the surface thermal emission, and is computed directly from the BRDF by the following relationship:
 
   
     
       
         
           
             
               
                 
                   
                     
                       ɛ 
                       v 
                     
                     ⁡ 
                     
                       ( 
                       
                         
                           Ω 
                           s 
                         
                         ; 
                         s 
                       
                       ) 
                     
                   
                   ≡ 
                   
                     1 
                     - 
                     
                       
                         ρ 
                         v 
                       
                       ⁡ 
                       
                         ( 
                         
                           
                             Ω 
                             s 
                           
                           ; 
                           s 
                         
                         ) 
                       
                     
                   
                 
                 = 
                 
                   1 
                   - 
                   
                     
                       ∫ 
                       
                         2 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         π 
                       
                     
                     ⁢ 
                     
                       
                         
                           f 
                           v 
                         
                         ⁢ 
                         
                             
                         
                         ( 
                         
                           
                             Ω 
                             s 
                           
                           , 
                           
                             
                               Ω 
                               ′ 
                             
                             ; 
                             s 
                           
                         
                         ) 
                       
                       ⁢ 
                       
                         μ 
                         ′ 
                       
                       ⁢ 
                       
                         ⅆ 
                         
                           Ω 
                           ′ 
                         
                       
                     
                   
                 
               
             
             
               
                 ( 
                 10 
                 ) 
               
             
           
         
       
     
   
   Finally, the Earth viewing radiance is obtained by substituting the surface radiance term, Eq. (8), into the RTE, Eq. (7): 
                         I   v     ⁡     (     Ω   ;   0     )       =           t   v     ⁡     (   s   )       ⁢       I   v     ⁡     (     Ω   ;   s     )         +       ∫       t   v   abs     ⁡     (   s   )       1     ⁢         t   v   sct     ⁡     (     s   ′     )       ⁢       B   v     ⁡     (     T     s   ′       )       ⁢           ⁢     ⅆ       t   v   abs     ⁡     (     s   ′     )             +       ∫       t   v   sct     ⁡     (   s   )       1     ⁢         t   v   abs     ⁡     (     s   ′     )       ⁢     {                 t   v   sun     ⁡     (     s   ′     )       ⁢     F   v   sun     ⁢       p   v     ⁡     (     Ω   ,       Ω   sun     ;     s   ′         )         +                 ∫     4   ⁢   π               ⁢         I   v     ⁡     (       Ω   ′     ;   s     )       ⁢       p   v     ⁡     (     Ω   ,       Ω   ′     ;     s   ′         )       ⁢           ⁢     ⅆ     Ω   ′                 }     ⁢           ⁢     ⅆ       t   v   sct     ⁡     (     s   ′     )                 ;     ⁢     
     ⁢         I   v     ⁡     (     Ω   ;   s     )       =             ɛ   v     ⁡     (     Ω   ;   s     )       ⁢       B   v     ⁡     (     T   g     )         +       ∫     2   ⁢   π               ⁢       [           t   v   sun     ⁡     (   s   )       ⁢     F   v   sun     ⁢     δ   ⁡     (     Ω   -     Ω   sun       )         +       I   v     ⁡     (       Ω   ′     ;   s     )         ]     ⁢       f   v     ⁡     (     Ω   ,       Ω   ′     ;   s       )       ⁢     μ   ′     ⁢     ⅆ     Ω   ′       ⁢           ⁢   if   ⁢           ⁢   s         =       s   g     .                 (   11   )               
In writing this final set of equations, the explicit dependence of the solid angle Ω on position s has been dropped for notational simplicity and the Dirac delta function δ(Ω−Ω sun ) has been introduced to move the surface direct solar illumination term inside the angular integration. There is also an implicit understanding that the boundary term I v (Ω;s) is zero unless the path actually terminates at the ground, s=s g . Eq. (11) serves as the starting point for the derivation of the MODTRAN™ band model.
 
The MODTRAN™ Band Model
 
   The MODTRAN™ band model approach is concerned with calculation of narrow spectral bin radiant intensities &lt;I 0 (Ω)&gt;: 
                   〈       I   0     ⁡     (   Ω   )       〉     ≡       1     Δ   ⁢           ⁢   v       ⁢       ∫     Δ   ⁢           ⁢   v       ⁢         I   v     ⁡     (     Ω   ;   0     )       ⁢           ⁢     ⅆ   v                   (   12   )               
The width Δv of the band model spectral bin should be selected sufficiently narrow to insure no more than a modest variation in all spectral data other than molecular absorption. This definition enables all the spectrally smoothly varying properties to be modeled with an effective band value, computed either as an interval average or central value. A smaller band width typically produces finer accuracy. The reformulated band model incorporated into MODTRAN™, described in the parent application which is incorporated herein by reference, provides 4 band model resolution options, 0.1, 1.0, 5.0 and 15.0 cm −1 . With an over-line used to represent an effective band value, the band model complement to the LBL LOS radiance expression, Eq. (11), has the following form:
 
                           〈       I   0     ⁡     (   Ω   )       〉     =       ⁢       〈       t     0   -&gt;   s       ⁢       I   s     ⁡     (   Ω   )         〉     -       ∫   0   s     ⁢         t     0   -&gt;     s   ′       sct     _     ⁢       B   _     ⁡     (     T     s   ′       )       ⁢       ⅆ     [         t     0   -&gt;     s   ′       part_abs     _     ⁢     〈     t     0   -&gt;     s   ′       mol_abs     〉       ]         ⅆ     s   ′         ⁢           ⁢     ⅆ     s   ′           -                     ⁢         ∫   0   s     ⁢         t     0   -&gt;       s   ′     -&gt;   sum       part_abs     _     ⁢     〈     t     0   -&gt;       s   ′     -&gt;   sum       mol_abs     〉     ⁢       t       s   ′     -&gt;   sum     sct     _     ⁢       F   sun     _     ⁢         p     s   ′       _     ⁡     (     Ω   ,     Ω   sun       )       ⁢     (       ⅆ       t     0   -&gt;     s   ′       sct     _         ⅆ     s   ′         )     ⁢           ⁢     ⅆ     s   ′           -                     ⁢       ∫   0   s     ⁢         t     0   -&gt;     s   ′       part_abs     _     ⁢     〈       t     0   -&gt;     s   ′       mol_abs     ⁢       ∫     4   ⁢   π       ⁢         I     s   ′       ⁡     (     Ω   ′     )       ⁢         p     s   ′       _     ⁡     (     Ω   ,     Ω   ′       )       ⁢           ⁢     ⅆ     Ω   ′             〉     ⁢     (       ⅆ       t     0   -&gt;     s   ′       sct     _         ⅆ     s   ′         )     ⁢           ⁢     ⅆ     s   ′                   ⁢     
     ⁢   BoundaryCondition   ⁢     
     ⁢             〈       t     0   -&gt;   s       ⁢       I   s     ⁡     (   Ω   )         〉     =       ⁢             ɛ   s     _     ⁡     (   Ω   )       ⁢       t     0   -&gt;   s     sct     _     ⁢       B   _     ⁡     (     T   g     )       ⁢       t     0   -&gt;     s   ′       part_abs     _     ⁢     〈     t     0   -&gt;   s     mol_abs     〉       +                     ⁢           t     0   -&gt;     s   -&gt;   sun       sct     _     ⁢       t     0   -&gt;     s   -&gt;   sun       part_abs     _     ⁢     〈     t     0   -&gt;     s   -&gt;   sun       mol_abs     〉     ⁢       F   sun     _     ⁢         f   s     _     ⁡     (     Ω   ,     Ω   sun       )       ⁢     μ   sun       +                     ⁢         t     0   -&gt;   s     sct     _     ⁢       t     0   -&gt;   s     part_abs     _     ⁢       〈       t     0   -&gt;   s     mol_abs     ⁢       ∫     2   ⁢   π       ⁢         I   s     ⁡     (     Ω   ′     )       ⁢         f   s     _     ⁡     (     Ω   ,     Ω   ′       )       ⁢           ⁢     μ   ′     ⁢     ⅆ     Ω   ′             〉     .                       (   13   )               
In these equations, the LOS integrands are expressed as functions of the refractive path length s′. Furthermore, the attenuation resulting from absorption is partitioned into its molecular and particulate contributions. The labeling of the band-averaged absorption as arising solely from particulate sources is somewhat misleading. Detailed line compilation data is unavailable for some molecular species, especially large molecules such as the chiorofluorocarbons (CFC&#39;s). For these species, absorption is modeled via relatively coarse spectral resolution cross-section data. The resulting transmittances are band averaged values and included in the  t part     —     abs    terms.
 
   Eq. (13) is the fundamental narrow band model radiance equation. It explicitly partitions the spectral data for which effective band values are used from the highly spectrally varying molecular absorption contributions. The effective coarse spectral resolution band data falls into three categories.
         1. Ground properties: The surface emissivity  ε s    and the surface BRDF  f s   ;   2. Particulate, Molecular Cross-Section &amp; Rayleigh data: The scattering phase function  p s′    and the atmospheric attenuation resulting from scattering  t sct    and particulate (and molecular cross-section) absorption  t part     —     abs   ; and   3. Source functions: The Planck blackbody function  B  and the TOA solar irradiance  F sun   .
 
The ground is comprised of solids and liquids. Employing lower spectral resolution surface data is justified because condensed matter does not exhibit the fine spectral structure of molecular gas absorption. Similarly, even though scattering data can be highly structured for a collection of particles of a specific size and shape; terrestrial aerosols and clouds consist of particulates of various size, shape and orientation. These distributions smooth out the fine spectral structure. Rayleigh scattering optical properties are also spectrally smoothly varying, and the molecular cross-section data is modeled as spectrally smoothly varying in the absence of accurate fine spectral resolution information. The variation of the Planck blackbody function at atmospheric temperatures (˜180 to 305K) for the MODTRAN™ band model widths is modest for the spectral range over which thermal emission is significant.
       

   The one term for which spectral averaging is somewhat suspect is the TOA solar irradiance. The sun is essentially a ˜6000K blackbody; however, molecular absorption within the solar atmosphere generates Fraunhofer structure. Two facts help support the spectral averaging of the solar data approximation: (1) The spectral correlation between the solar spectrum and the Earth atmospheric absorption is well represented as randomly correlated for most spectral bins; and (2) The current state-of-the-art in measurement and calculation of TOA solar irradiances suggests that spectral structure assignments finer than ˜1 cm −1  are unreliable. 
   Currently, the multiple scattering (MS) algorithms within MODTRAN™ only support specification of surface reflectance via an angularly independent spectral albedo, f v (Ω,Ω′;s)=α v (s)/π. Given this limitation, the surface reflected downward flux is approximated as a product of the directional reflectivity and the hemisphere integrated downward flux, F v   ↓ (s): 
                     ∫     2   ⁢           ⁢   π       ⁢         I     v   ⁢               ⁡     (       Ω   ′     ;   s     )       ⁢       f   v     ⁡     (     Ω   ,       Ω   ′     ;   s       )       ⁢     μ   ′     ⁢     ⅆ     Ω   ′           ≅         ρ     v   ⁢               ⁡     (     Ω   ;   s     )       ⁢       F   v   ↓     ⁡     (   s   )         ≡         ρ   v     ⁡     (     Ω   ;   s     )       ⁢       ∫     2   ⁢           ⁢   π       ⁢         I   v     ⁢           (       Ω   ′     ;   s     )     ⁢     μ   ′     ⁢       ⅆ     Ω   ′       .                   (   14   )               
This expression for the reflected downward flux is quite reasonable if the distribution of downward diffuse radiance is nearly homogeneous or the surface reflectance nearly Lambertian. On the other hand, for clear sky solar illumination of a specular surface such as still water, Eq. (14) is a poor approximation; however, in this case, the observed radiance will be dominated by the direct solar illumination and the error in total observed radiance may still be small. Merging the flux approximation into the narrow band model radiance expression, Eq. (13), the surface reflected diffuse radiance term factors as follows:
 
   
     
       
         
           
             
               
                 
                   〈 
                   
                     
                       t 
                       
                         
                             
                         
                         ⁢ 
                         
                           0 
                           ⁢ 
                           
                               
                           
                           -&gt; 
                           
                               
                           
                           ⁢ 
                           s 
                         
                       
                       
                         
                             
                         
                         ⁢ 
                         mol_abs 
                       
                     
                     ⁢ 
                     
                       
                         ∫ 
                         
                           2 
                           ⁢ 
                           π 
                         
                       
                       ⁢ 
                       
                         
                           
                             I 
                             s 
                           
                           ⁡ 
                           
                             ( 
                             
                               Ω 
                               ′ 
                             
                             ) 
                           
                         
                         ⁢ 
                         
                           
                             
                               f 
                               s 
                             
                             _ 
                           
                           ⁡ 
                           
                             ( 
                             
                               Ω 
                               , 
                               
                                 Ω 
                                 ′ 
                               
                             
                             ) 
                           
                         
                         ⁢ 
                         
                           μ 
                           ′ 
                         
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           ⅆ 
                           
                             Ω 
                             ′ 
                           
                         
                       
                     
                   
                   〉 
                 
                 ≅ 
                 
                   
                     
                       
                         ρ 
                         s 
                       
                       _ 
                     
                     ⁡ 
                     
                       ( 
                       Ω 
                       ) 
                     
                   
                   ⁢ 
                   
                     
                       〈 
                       
                         
                           t 
                           
                             
                                 
                             
                             ⁢ 
                             
                               0 
                               ⁢ 
                               
                                   
                               
                               -&gt; 
                               
                                   
                               
                               ⁢ 
                               s 
                             
                           
                           
                             
                                 
                             
                             ⁢ 
                             mol_abs 
                           
                         
                         ⁢ 
                         
                           
                             F 
                             v 
                             ↓ 
                           
                           ⁡ 
                           
                             ( 
                             s 
                             ) 
                           
                         
                       
                       〉 
                     
                     . 
                   
                 
               
             
             
               
                 ( 
                 15 
                 ) 
               
             
           
         
       
     
   
   When MODTRAN™ utilizes its correlated-k (CK) option, an algorithm that recasts the molecular band model into weighted sum of monochromatic calculations, the path radiant intensity computation follows the prescription of Eq. (13) with the caveat of Eq. (15). However, further approximations are necessary if the band model algorithm is used without CK. The issue here is that the MODTRAN™ atmospherically scattered diffuse radiance and the downward surface flux terms are computed from the MS algorithms within MODTRAN™. These MS algorithms were developed to solve the coupled monochromatic radiation transport problem, and they rely heavily upon Beer&#39;s Law. Beer&#39;s Law states that the spectral transmittance of a multi-segmented path equals the product of the individual segment transmittances:
 
 t   v (0 →s′→s )= t   v (0→ s ′)  t   v ( s′→s ).  (16)
 
Alternatively stated, segment extinction optical depths are additive. This relationship does hold for CK transmittances but not for band model in-band transmittances:
 
&lt; t   0→s′→s   &gt;≠&lt;t   0→s   &gt;&lt;t   s′→s &gt;.  (17)
 
When the CK option is not invoked, MODTRAN™ proceeds with multiple scatter source term and level flux calculations computed from a Beer&#39;s Law model even though the layer optical depths are poorly defined. The inaccuracy of this approach is limited somewhat by combining the MS source terms with foreground transmittances computed for the inhomogeneous path. Thus, the LOS radiance terms containing diffuse radiation are calculated as follows:
 
   
     
       
         
           
             
               
                 
                   
 
                 
                 ⁢ 
                 
                   
                     
                       
                         
                           
                             〈 
                             
                               
                                 t 
                                 
                                   0 
                                   -&gt; 
                                   s 
                                 
                                 mol_abs 
                               
                               ⁢ 
                               
                                 F 
                                 s 
                                 ↓ 
                               
                             
                             〉 
                           
                           ≅ 
                           
                             
                               〈 
                               
                                 t 
                                 
                                   0 
                                   -&gt; 
                                   s 
                                 
                                 mol_abs 
                               
                               〉 
                             
                             ⁢ 
                             
                               〈 
                               
                                 F 
                                 s 
                                 ↓ 
                               
                               〉 
                             
                           
                         
                       
                       
                         
                           band 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           model 
                         
                       
                     
                     
                       
                         
                           
                             〈 
                             
                               
                                 t 
                                 
                                   0 
                                   -&gt; 
                                   
                                     s 
                                     ′ 
                                   
                                 
                                 mol_abs 
                               
                               ⁢ 
                               
                                 
                                   ∫ 
                                   
                                     4 
                                     ⁢ 
                                     
                                         
                                     
                                     ⁢ 
                                     π 
                                   
                                 
                                 ⁢ 
                                 
                                   
                                     
                                       I 
                                       
                                         s 
                                         ′ 
                                       
                                     
                                     ⁢ 
                                     
                                         
                                     
                                     ( 
                                     
                                       Ω 
                                       ′ 
                                     
                                     ) 
                                   
                                   ⁢ 
                                   
                                     
                                       
                                         p 
                                         
                                           s 
                                           ′ 
                                         
                                       
                                       _ 
                                     
                                     ⁡ 
                                     
                                       ( 
                                       
                                         Ω 
                                         , 
                                         
                                           Ω 
                                           ′ 
                                         
                                       
                                       ) 
                                     
                                   
                                   ⁢ 
                                   
                                     ⅆ 
                                     
                                       Ω 
                                       ′ 
                                     
                                   
                                 
                               
                             
                             〉 
                           
                           ≅ 
                           
                             
                               〈 
                               
                                 t 
                                 
                                   0 
                                   -&gt; 
                                   s 
                                 
                                 mol_abs 
                               
                               〉 
                             
                             ⁢ 
                             
                               
                                 ∫ 
                                 
                                   4 
                                   ⁢ 
                                   
                                       
                                   
                                   ⁢ 
                                   π 
                                 
                               
                               ⁢ 
                               
                                 
                                   〈 
                                   
                                     
                                       I 
                                       
                                         s 
                                         ′ 
                                       
                                     
                                     ⁢ 
                                     
                                         
                                     
                                     ( 
                                     
                                       Ω 
                                       ′ 
                                     
                                     ) 
                                   
                                   〉 
                                 
                                 ⁢ 
                                 
                                   
                                     
                                       p 
                                       
                                         s 
                                         ′ 
                                       
                                     
                                     _ 
                                   
                                   ⁡ 
                                   
                                     ( 
                                     
                                       Ω 
                                       , 
                                       
                                         Ω 
                                         ′ 
                                       
                                     
                                     ) 
                                   
                                 
                                 ⁢ 
                                 
                                   ⅆ 
                                   
                                     Ω 
                                     ′ 
                                   
                                 
                               
                             
                           
                         
                       
                       
                         
                           with 
                           ⁢ 
                           
                               
                           
                           ⁢ 
                           no 
                           ⁢ 
                           
                             
                                 
                             
                             ⁢ 
                             
                                 
                             
                           
                           ⁢ 
                           CK 
                         
                       
                     
                   
                   . 
                 
               
             
             
               
                 ( 
                 18 
                 ) 
               
             
           
         
       
     
   
   Together, Eqs. (13), (15) and (18) define the MODTRAN™ path radiance formulation. However, implementation of the approach still requires a description of three key model components: the MODTRAN™ statistical band model molecular transmittance calculation, the multiple scattering algorithms, and the surface water model. Descriptions of these components follow. 
   SUMMARY OF THE INVENTION 
   This invention comprises an improved band model method for modeling atmospheric propagation at arbitrarily fine spectral resolution, sometimes termed herein as “Improvements to MODTRAN4™”. Featured in the invention is a method for predicting atmospheric radiation transport parameters including transmittance and radiance, comprising defining one or more conditions under which the atmospheric transmittance and radiance will be predicted, providing atomic and molecular transition data for a given spectral range and atmospheric conditions, dividing the spectral region into a number of spectral bins that determine the spectral resolution, each bin having a width of less than 1.0 cm −1 , calculating atomic and molecular species line center absorption from at least the equivalent width of the atomic and molecular transitions centered within each spectral bin, calculating line tail absorption within each spectral bin from atomic and molecular transitions not centered within the bin, and determining atomic and molecular species spectral transmittances for each spectral bin, the spectral transmittance having a value which is a function of at least the calculated line center absorptions and the calculated line tail absorptions. 
   The spectral bins may have a width of about 0.1 cm −1 . The step of calculating line center absorption may include calculating, from an exact expansion, the bin Voigt equivalent width of atomic and molecular transitions whose centers lie within each spectral bin. The exact expansion may be an exact modified Bessel functions expansion. The calculating line tail absorption step may include subtracting line-tail absorption as calculated from the column strength, the Lorentz half-width, the Doppler half-width, and the line tail spectral displacement. The calculating line center absorption step may include determining the Voigt line-shape function computed at specific frequencies. 
   The line tail calculation step may include calculating line tail absorption within each bin from atomic and molecular transitions centered outside of the bin using Padé approximant spectral fits to Voigt absorption coefficient curves. The line tail absorption calculation step may include determining a database of temperature and pressure dependent Padé approximant spectral fits to Voigt absorption coefficient curves. There may be five Padé parameters. The Padé parameters may be determined from summed line tail spectral absorption coefficients. One Padé parameter may be determined at the center of the bin, and one at each edge of the bin. One Padé parameter may be the derivative of the absorption coefficient with respect to the normalized spectral variable at the line center. One Padé parameter may be the integral of the spectral absorption coefficient over the spectral band. The Padé parameters database may be generated for a plurality of temperatures. The Padé parameters database may be determined for a plurality of pressures. The line center absorptions may be calculated from atomic and molecular transitions centered no more than half a spectral bin width from the bin, and the line tail absorptions are calculated from atomic and molecular transitions not centered within a half spectral bin from the bin. 
   The conditions may comprise atmospheric species for which transition data is provided. Band model data may be provided for at least some of the species. Absorption cross section data may be provided for at least some of the species. The conditions may comprise an arbitrary quantity of species for which transition data is provided. The conditions may comprise a surface water model. The surface water model may introduce a spectral attenuation layer. The surface water model may predict surface leaving radiances due to surface water. 
   The defining step may comprise allowing the user to define the radiation refractive path through the atmosphere. The refractive path may comprise a circular arc. The conditions may comprise the aerosol spectral optical properties. The properties can preferably be varied as a function of altitude. The properties may comprise spectral extinction, spectral scattering and spectral phase function. The method may further comprise determining the medium embedded diffuse transmittance. 
   Featured in another embodiment of the invention is a method for predicting atmospheric radiation transport parameters including transmittance and radiance, comprising defining one or more conditions under which the atmospheric transmittance and radiance will be predicted, wherein the conditions comprise an arbitrary number of atmospheric species for which transition data is provided, a surface water model, a user-defined radiation refractive path through the atmosphere, and the aerosol spectral optical properties, providing atomic and molecular transition data for a given spectral range and atmospheric conditions, dividing the spectral region into a number of spectral bins that determine the spectral resolution, each bin having a width of less than 1.0 cm −1 , calculating atomic and molecular species line center absorption from at least the equivalent width of the atomic and molecular transitions centered within each spectral bin, calculating line tail absorption within each spectral bin from atomic and molecular transitions not centered within the bin, determining atomic and molecular species spectral transmittances for each spectral bin, the spectral transmittance having a value which is a function of at least the calculated line center absorptions and the calculated line tail absorptions, and determining the medium embedded diffuse transmittance. 
   A primary use of the invention is spectral exploitation. More specifically, this invention greatly enhances the ability to perform the radiation transport necessary to accurately atmospherically correct airborne and satellite imagery of the Earth. Thus, a user can accurately translate the at-sensor observed radiances (i.e. an image captured by an airborne or space-based image sensor) into at-ground spectral reflectances and emissivities. When accurately determined, the spectral exploitation scientists can use the retrieved surface reflectances to classify and identify materials on the Earth. Applications are many including tracking of crop health, monitoring of the littoral zone, mapping of mineral fields, and detecting military targets. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
     Other objects, features and advantages will occur to those skilled in the art from the following description of the preferred embodiments and the accompanying drawings in which: 
       FIG. 1  is a plot of spectral absorbance for the primary sources of atmospheric extinction. These 4 nm resolution curves were generated by an embodiment of the invention for a vertical path from space with a mid-latitude winter model atmosphere. 
       FIG. 2  is an illustration of LOS and solar illumination paths, useful in understanding an aspect of the invention. 
       FIG. 3  is two plots illustrating the auxiliary species aspect of the invention. These are Isolated HCl Nadir Viewing 0.2 cm −1  Spectral Transmittances. The inventive auxiliary species option was exercised with a 1 km thick surface layer of HCl. The left-hand side contains the continuous spectrum from 3.2 to 3.7 μm with 10.0 ppm HCl; the right-hand side illustrates the change in absorption arising from varying HCl concentrations for two of the vibrational rotational bands. The doublet structure arises because of the Cl 35  and Cl 37  isotopes. 
       FIG. 4  is two plots also illustrating the auxiliary species aspect of the invention. These are Spectral Transmittances and Thermal Emission as Predicted by an embodiment of the invention for Nadir-Viewing through a MLS atmosphere with a 1 km High HCl layer at the Surface. On the right ths blackbody ground temperature is 10K below the surface air temperature. 
       FIG. 5  is an illustration of missing photon trajectory, useful in understanding an aspect of the invention. 
       FIG. 6  is a plot of spectral absorption and scattering coefficients of liquid water, useful in understanding aspects of the invention. 
       FIG. 7  is three illustrations of reflections by liquid water, useful in understanding aspects of the invention. 
       FIG. 8  is a plot of liquid water spectral index of refraction, useful in understanding aspects of the invention. 
   

   DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS OF THE INVENTION 
   IMPROVEMENTS TO MODTRAN4™ Molecular Transmittance 
   SERTRAN, the Spectral Enhanced Resolution band model IMPROVEMENT TO MODTRAN4™, is statistical. SERTRAN is described in detail in the parent patent application incorporated herein by reference. MODTRAN™ does not account for the exact location of molecular transitions and the distribution of their line strengths within each spectral bin. Instead, these distributions are modeled statistically. In this spirit, the formulation of the molecular transmittance calculation begins with the stipulation that the spectral bin molecular absorption from distinct species is randomly correlated, so that the combined species molecular transmittance &lt;t mol     —     abs &gt; for a band model spectral bin can be computed as a product of individual band model species transmittances &lt;t k     —     abs &gt;: 
                   〈     t             ⁢   mol_abs       〉     =       ∏     k   ⁢           ⁢   mol       ⁢           ⁢       〈     t   k_abs     〉     .               (   19   )               
This approximation works well for a number of reasons. As illustrated in  FIG. 1 , a single molecular species serves as a dominant absorption source in much of the visible through IR spectral region. Where overlap between molecular species does occur, H 2 O is typically one of the species; since H 2 O is an asymmetry top, its transition frequency and line strength distributions are well represented as random. It is therefore quite reasonable to approximate H 2 O absorption correlation with other molecular species as random. Furthermore, when moderate spectral resolution results are sought, spectral averaging of the relatively fine resolution MODTRAN™ values tends to stochastically cancel statistical variations if there is no overall biasing.
 
   Within MODTRAN™, the individual species spectral bin transmittance is itself partitioned into three distinct components. These components are:
         (1) Line center transmittance: the transmittance arising from molecular transitions centered within the spectral bin;   (2) Line tail transmittance: the transmittance arising from molecular transitions centered outside of the spectral bin but not further than 25 cm −1  from the spectral bin center; and   (3) Continuum transmittance: the transmittance arising from molecular transitions centered more than 25 cm −1  from a spectral bin.
 
Again, the assumption is made that these terms are randomly correlated:
 
&lt; t   k     abs     &gt;=&lt;t   k     —     cen   &gt;&lt;t   k     —     tail   &gt;&lt;t   k     —     cont &gt;.  (20)
 
Since the line tail and continuum contributions tend to be relatively smoothly varying across the spectral interval, the random correlation assumption is quite reasonable. In the following subsections, the modeling approach for each of these components is detailed. Additional sections describe treatment of auxiliary species and of aerosols for which spectral optical properties vary with altitude.
 
Line Center Transmittance
       

   The starting point for the calculation of the MODTRAN™ in-band molecular line center transmittance is the Plass expression. Plass proved that the homogeneous path transmittance resulting from n identical molecular lines randomly distributed within each spectral interval of width Δv for a contiguous collection of spectral bins is given by the following power law expression: 
                     〈     t   k_cen     〉     =       (     1   -       W   sl       Δ   ⁢           ⁢   v         )     n       ,           (   21   )               
where W sl  is the finite-bin equivalent width (EW) of a single Voigt line,
 
                         W   sl     ⁡     (     Su   ,     γ   c     ,     γ   d       )       ≡       ∫     Δ   ⁢           ⁢   v       ⁢       [     1   -     exp   ⁡     (       -   Su     ⁢           ⁢     f   v       )         ]     ⁢           ⁢     ⅆ   v           ;     ⁢     
     ⁢         f             ⁢   v       ⁡     (       γ             ⁢   c       ,     γ             ⁢   d         )       ≡                 ⁢     γ             ⁢   c                   ⁢     π             ⁢     3   /   2             ⁢       ∫     -   ∞               ⁢   ∞       ⁢         exp   ⁢     (     -           ⁢     x             ⁢   2         )                 ⁢       γ             ⁢   c               ⁢   2       ⁢           +           ⁢       (     v   ⁢           -           ⁢       γ             ⁢   d       ⁢           ⁢   x       )     2           ⁢           ⁢       ⅆ   x     .                     (   22   )               
The Voigt spectral line-shape f v (γ c ,γ d ) [cm] at line center displacement frequency v, the convolution of Lorentz and Doppler line-shape functions, characterizes the spectral dependence of molecular line absorption. In these equations, γ c  is a collision (Lorentz) half-width at half maximum [cm −1 ], γ d  is a Doppler half-width at 1/e of the maximum [cm −1 ], S is a line strength [cm −2 /Atm] and u is a path column density [Atm cm].
 
   The Plass expression, Eq. (21), defines the band model curve-of-growth (COG) for the randomly distributed collection of identical molecular lines, i.e., the rate at which the molecular absorption increases with column density. The actual distribution of molecular lines in a specific spectral interval will exhibit a COG unique to the spectral band. The objective of band model theory in general is to find an optimal fit to each COG. In MODTRAN™, the Plass expression is used as the COG functional form. The fitting parameters are the average line strength of each line S, the effective number of molecular lines in the spectral interval n, the spectral bin averaged Lorentz half-width γ c  and the Doppler half-width γ d . 
   Band Model Parameter Determination 
   The traditional band model equations for the temperature dependent line strength, S(T), and effective line number, n(T), parameters are derived by fitting the EW COG in the weak-line (small column density) and Lorentz strong-line (large column density) limits. The Voigt EW varies linearly with column density in the weak-line limit, and it increases with the square root of the column density in the strong line limit. The accuracy of the interpolation between these limits largely dictates the accuracy of the band model. The expressions for the band model parameters resulting from the weak and strong line limit constraints are: 
                     S   ⁡     (   T   )       =       1     n   ⁡     (   T   )         ⁢       ∑     sub   -     bin   ⁢           ⁢   j         ⁢       S   j     ⁡     (   T   )             ;       n   ⁡     (   T   )       =       1       S   ⁡     (   T   )           ⁢       ∑     sub   -     bin   ⁢           ⁢   j         ⁢           S   j     ⁡     (   T   )         .                   (   23   )               
(Note that the actual temperature dependent parameters stored by MODTRAN™ are (S/d)≡nS/Δv and (1/d)≡n/Δv.) Typically, these sums would be over all the transitions centered within the spectral bin for a given molecular species. For MODTRAN™, the spectral bin sums are over spectral sub-intervals of width 0.01 cm −1 , and each S j (T) is itself the sum of temperature-dependent line strengths for lines centered in sub-interval j. This coagulation of molecular transitions within narrow spectral bins has no effect on the calculation the total line strength, n(T) S(T), but it reduces the effective number of lines, n(T). MODTRAN™ temperature dependent band model parameters are generated for 6 atmospheric temperatures, T=180, 205, 230, 255, 280, and 305K. The Lorentz half-width band model parameter, γ c   0 , is computed as a line strength weighted average of molecular transition half-widths at standard temperature (T 0 =273.15K) and pressure (P 0 =1 Atm).
 
                     γ   c   0     =         ∑   p     ⁢       γ     c   ,   p     0     ⁢     S   p             n   ⁡     (     T   0     )       ⁢     S   ⁡     (     T   0     )             ;         γ   c     ⁡     (     P   ,   T     )       =         γ   c   0     ⁡     (     P     P   0       )       ⁢         (       T   0     T     )       3   /   4       .                 (   24   )               
The Doppler half-width at 1/e of maximum is defined using the spectral bin center frequency,  v :
 γ d =√{square root over (2 kT/m )}    v   / c,   (25) 
where k is the Boltzmann constant (1.3807E-16 erg/K), m is the molecular mass (g/molecule), and c is the speed of light in a vacuum (2.9979246E+10 cm/sec).
 
   Sub-intervals are used in the definition of the band model parameters, Eq. (23), for two reasons: (1) the atmospheric absorption from nearly coincident lines (such as those arising from degenerate transitions) is accurately modeled as a single line of combined strength, and (2) the Plass expression was derived for a random distribution of line centers, but near degeneracy occurs at a higher frequency than a random distribution would predict. On a more practical level, it has been shown by extensive comparisons to LBL calculations that the classical expression for the effective number of lines n is too large; too much absorption is predicted. The sub-interval summing decreases the effective line number without changing the total line strength nS, and thereby lowers the band model absorption. 
   To increase accuracy in the IMPROVEMENTS TO MODTRAN4™ band model absorption for the coarser resolution band models (Δv=1.0, 5.0 and 15.0 cm −1 ), the effective line number band model parameter is tuned using the 0.1 cm −1  band model. The two-step approach is relatively straightforward. First, the column amount U is determined for which the 0.1 cm −1  band model predicts a homogeneous path transmittance equal to 1/e when integrated (summed) over a given coarse spectral bin: 
                       1   N     ⁢       ∑     N   ⁢           ⁢   bias       ⁢       〈       t   k_abs     ⁡     (   U   )       〉       0.1   ⁢     cm     -   1               =     1   e       ;     N   =         Δ   ⁢           ⁢   v       0.1   ⁢           ⁢     cm     -   1           .               (   26   )               
Here Δv is the coarse band model bin width and the loop is over the N 0.1 cm −1  spectral sub-bins within the coarse spectral bin. Once the spectrally dependent 1/e transmittance column density has been determined, the coarse resolution effective line number n is perturbed to give the same in-band transmittance:
 &lt; t   k     —     abs ( U;n )&gt; Δv =1 /e.   (27) 
These calculations are performed at each of the MODTRAN™ temperatures and for a specific pressure depending upon the molecule. For H 2 O, which is concentrated near the ground, the reference pressure is set to 1 Atm. The absorption of O 3  is most important at higher altitudes, so its effective line number values are tuned for a pressure of 0.1 Atm. For all other species, the homogeneous path is defined with a reference pressure of 0.5 Atm.
 
Calculation of Equivalent Width
 
   The precise spectral domain of the finite bin EW integral in Eq. (22) was left ambiguous. In traditional band model theory, the finite bin EW is computed for a line located in the center of the spectral bin; thus, the spectral domain of the integral extends from −Δv/2 to +Δv/2. However, MODTRAN™ seeks to determine the EW of lines whose centers are randomly distributed within the spectral bin. For a single pure Lorentzian line in the weak line limit, the exact expression for the offset from bin center D(r c ) that produces a finite bin integrated absorption equal to the random distribution value can be derived as a function of the ratio of the Lorentz half-width to the bin width, r c =γ c /Δv: 
                     (       D   ⁡     (     r   c     )         Δ   ⁢           ⁢   v       )     2     =       1   4     -       (     1   +     r   c   2       )     ⁢         1   -       r   c     ⁢     tan   ⁡     [       r   c     ⁢     ln   ⁡     (     1   +     1     r   c   2         )         ]             2   +       (       1   -     r   c   2         r   c       )     ⁢     tan   ⁡     [       r   c     ⁢     ln   ⁡     (     1   +     1     r   c   2         )         ]             .                 (   28   )               
The function D(r c ) is monotonic for positive r c , decreasing from D(0)=½Δv(the edge of the spectral bin) to D(r c →∞)=Δv/√12. The terrestrial surface Lorentz half-widths are typically ˜0.07 cm −1 . For the MODTRAN™ band model widths of 0.1, 1.0, 5.0 and 15.0 cm −1 , the 0.07 cm −1  Lorentz transition weak-line displacement values are D(0.07/0.1)=0.028 cm −1 , D(0.07/1.0)=0.34 cm −1 , D(0.07/5.0)=1.96 cm −1 , and D(0.07/15.0)=6.2 cm −1 , respectively.
 
   The weak-line Lorentz transition offsets are only a guide for determining the optimal Voigt line-shape transition values; the actual displacements used in MODTRAN™ are D[0.1 cm −1 ]=0.025 cm −1 , D[1.0 cm −1 ]=0.30 cm −1 , D[5.0 cm −1 ]=1.97 cm −1 , and D[15.0 cm −1 ]=6.3 cm −1  with the EW integral having the form: 
                         W   sl     =       ⁢       ∫           -   1     /   2     ⁢   Δ   ⁢           ⁢   v     +   D           1   /   2     ⁢   Δ   ⁢           ⁢   v     +   D       ⁢       (     1   -     ⅇ     -     Suf   v           )     ⁢     ⅆ   v                     =       ⁢         ∫   0         1   /   2     ⁢   Δ   ⁢           ⁢   v     -   D       ⁢       (     1   -     ⅇ     -     Suf   v           )     ⁢     ⅆ   v         +       ∫   0         1   /   2     ⁢   Δ   ⁢           ⁢   v     +   D       ⁢       (     1   -     ⅇ     -     Suf   v           )     ⁢     ⅆ   v                       =       ⁢       2   ⁢       ∫   0   ∞     ⁢       (     1   -     ⅇ     -     Suf   v           )     ⁢     ⅆ   v           -       ∫         1   /   2     ⁢   Δ   ⁢           ⁢   v     -   D     ∞     ⁢       (     1   -     ⅇ     -     Suf   v           )     ⁢     ⅆ   v         -                     ⁢       ∫         1   /   2     ⁢   Δ   ⁢           ⁢   v     +   D     ∞     ⁢       (     1   -     ⅇ     -     Suf   v           )     ⁢     ⅆ   v                     ≡       ⁢       2   ⁢     W   0       -     W         1   /   2     ⁢   Δ   ⁢           ⁢   v     -   D       -     W         1   /   2     ⁢   Δ   ⁢           ⁢   v     +   D                       (   29   )               
The second line above is derived by noting that f v is symmetric in spectral displacement frequency v(f v =f −v ), and the Voigt line tail EW W Δ  is defined in the last line. The methodology employed to compute the Voigt line tail EW is described in the parent SERTRAN patent application, incorporated by reference herein.
 
Spherical Refractive Geometry
 
   From its inception, MODTRAN™ included a spherical refractive geometry package. The modeling of refraction was introduced to model near-horizon sun and/or low-altitude near-horizontal path scenarios. As is well known, the sun is already below the horizon when sunset begins. The sun only appears as if it is above the horizon because of refraction. 
   With the locally spherical Earth approximation used in MODTRAN™, Snell&#39;s Law dictates that the product of three terms, (1) the index of refraction, 1+N   v   (H s ), (2) the Earth center distance (the sum of the Earth radius R e  and the altitude H s ), and (3) the sine of the path zenith angle √{square root over (1−μ s   2 )} at each position s along a line-of-sight (LOS), is a path constant:
 
[1 +N     v   ( H   s )]( R   e   +H   s )√{square root over (1−μ s   2 )}=Constant  (32)
 
Here, μ s  is the cosine of the path zenith at path length s. The MODTRAN™ refractivity, N   v   , vertical profile is computed at the mean spectral frequency  v  of the calculation band pass and is defined as a function of pressure, temperature, H 2 O density, CO 2  density and the central frequency.
 
   For a path specified by an initial altitude H 0  and direction μ 0 &lt;0, MODTRAN™ follows the line-of-sight (LOS) trajectory through the atmosphere continually updating the path zenith angle according to Eq. (1). With normal refraction, the path bends downward from the straight-line path because of the decrease in refractivity with increasing altitude. A limb path is defined as a downward LOS which passes through a minimum tangent altitude, H tan , before beginning an upward trajectory. At the tangent altitude, the path zenith angle is exactly 90°:
 
[1 +N     v   ( H   0 )]( R   e   +H   0 )√{square root over (1−μ 0   2 )}=[1 +N     v   ( H   tan )]( R   e   +H   tan ).  (33)
 
The MODTRAN™ refractive geometry package was specifically designed to model this scenario.
 
   When large temperature and/or H 2 O gradients are present within the atmosphere, non-standard refraction can occur. For optical frequencies, this scenario can arise within the marine boundary layer. If refractivity gradient, dN   v   /d z, is positive at altitude z s  along the optical path, subrefraction, i.e., upward bending of the optical path, occurs. On the other hand, if the decrease in refractivity with altitude is excessive, superrefraction results. Superrefraction is defined as larger than normal bending towards the Earth. If the radius of curvature for the refraction in an upward path exceeds the radius of curvature of the Earth, the path will pass through a tangent maximum and proceed back downward towards the Earth. The MODTRAN™ refractive geometry package was not designed to handle either of these scenarios. In particular, an upward path is always assumed to continue upward; if the superrefraction exceeds the critical refractivity gradient, the sine of the path zenith angle will exceed one at the boundary layer altitude level just above the tangent maximum. IMPROVEMENTS TO MODTRAN4™ detects this failure, and aborts processing. 
   To properly model the subrefraction and superrefraction scenarios, the invention (IMPROVEMENTS TO MODTRAN4™) provides an option to bypass the built-in refractive geometry package and explicitly enter refractive path data. Specification of a user-defined path requires that some stringent rules be followed. Four columns of data are required. As illustrated in Table 1 (This example is kindly been provided by Drs. Denis Dion and Vincent Ross, Defence Research and Development Canada, Valcartier, 2005), these columns contain the cumulative (spherical) Earth surface distance in meters, the altitude in meters, the path zenith angle in degrees, and the incremental path segment length in meters. In Table 1, all lines beginning with an exclamation are ignored and the word “NEXT” is used as the delimiter. Every atmospheric profile altitude level between the beginning of the path (H 1 ) and the end of the path (H 2 ) must be explicitly included in the user-defined path file. If the LOS passes through a minimum tangent point, the altitude of the tangent point and the altitude of the minimum of H 1  and H 2  must be included in the path definition. Similarly, if the LOS passes through a maximum tangent point, the altitude of the tangent point and the altitude of the maximum of H 1  and H 2  must be included in the path definition. In Table 1, a minimum tangent altitude occurs at 0.16075426526368 m, and the final altitude (H 2 =25.821244410585 m) is listed as the third entry in the downward leg. If a MODTRAN™ atmosphere includes a non-zero ground altitude or a cloud, it is best to run MODTRAN™ first without a user-defined path to determine the model relayering. Note that it is not necessary that “mirrored” path segments within a common atmospheric layer have identical path lengths or angles. 
   
     
       
         
             
           
             
               TABLE 1 
             
           
          
             
                 
             
             
               Sample User-Defined Path Parameter File. 
             
          
         
         
             
          
             
               !#Full refraction (Mirage ray) 
             
             
               !#Curved earth representation (R = 6371230.0 m) 
             
             
               ! 
             
          
         
         
             
             
             
             
          
             
               !Earth surface 
               Altitude 
               Zenith Angle 
               Segment Lengths 
             
             
               !distance (m) 
               (m) 
               (deg) 
               (m) 
             
             
                 
             
          
         
         
             
          
             
               ! 
             
          
         
         
             
             
             
             
          
             
               0.0000000000000 
               29.900000000000 
               90.189000000000 
               0.0000000000000 
             
             
               1041.9012954644 
               26.538244292473 
               90.180733953598 
               1041.9113336288 
             
             
               1270.3560625281 
               25.821244410585 
               90.178908112818 
               228.45683147277 
             
             
               2038.7644494551 
               23.463065803670 
               90.172762325063 
               768.41497742361 
             
             
               2964.7678300263 
               20.731295125158 
               90.165288500462 
               926.01062164312 
             
             
               3824.1902900742 
               18.304583613239 
               90.158276863904 
               859.42851892983 
             
             
               4621.1260668439 
               16.148865050039 
               90.151693040154 
               796.94084715816 
             
             
               5359.4938828508 
               14.233877421347 
               90.145503806883 
               738.37205984105 
             
             
               6043.0445474172 
               12.532738097948 
               90.139677023359 
               683.55421721182 
             
             
               6675.3673650761 
               11.021566456806 
               90.134181491957 
               632.32579225008 
             
             
               7259.8955694608 
               9.6791486444739 
               90.128986871664 
               584.53069546048 
             
             
               7799.9109567221 
               8.4866397766488 
               90.124063541044 
               540.01747381047 
             
             
               8298.5481454158 
               7.4272993933662 
               90.119382401136 
               498.63893670373 
             
             
               8758.7986010627 
               6.4862564561342 
               90.114914739278 
               460.25192023864 
             
             
               9183.5148066126 
               5.6503005880400 
               90.110631911152 
               424.71743276353 
             
             
               9575.4155575907 
               4.9076966262520 
               90.106504727038 
               391.90177926706 
             
             
               9937.0915791127 
               4.2480198836023 
               90.102503543286 
               361.67688300182 
             
             
               10271.013265885 
               3.6620098066495 
               90.098596673395 
               333.92240826494 
             
             
               10579.540504055 
               3.1414399758701 
               90.094750408679 
               308.52784207464 
             
             
               10864.936400155 
               2.6790026230402 
               90.090926277237 
               285.39640111885 
             
             
               11129.386251322 
               2.2682060446596 
               90.087080169216 
               264.45027291130 
             
             
               11375.024382517 
               1.9032834713067 
               90.083158113572 
               245.63848268456 
             
             
               11603.982028179 
               1.5791121136302 
               90.079087071713 
               228.95793773649 
             
             
               11818.455868850 
               1.2911412485443 
               90.074773233897 
               214.47408232363 
             
             
               12020.841683321 
               1.0353283360997 
               90.070068860443 
               202.38601311481 
             
             
               12214.001067940 
               .80808227023856 
               90.064744519726 
               193.15954626868 
             
             
               12401.895848521 
               .60621296678608 
               90.058369650734 
               187.89490992932 
             
             
               12591.541884888 
               .42688658099378 
               90.049986342712 
               189.64613663414 
             
             
               12802.014760646 
               .26758572597844 
               90.036744693052 
               210.47294785778 
             
             
               13135.178041776 
               .16075426526368 
               90.000000000000 
               333.16331427414 
             
             
               13468.341309379 
               .26758572597844 
               89.963255306948 
               333.16330074681 
             
             
               13678.814177420 
               .42688658099378 
               89.950013657288 
               210.47294014061 
             
             
               13868.460213787 
               .60621296678608 
               89.941630349266 
               189.64613663414 
             
             
               14056.354994368 
               .80808227023856 
               89.935255480274 
               187.89490992932 
             
             
               14249.514378987 
               1.0353283360997 
               89.929931139557 
               193.15954626868 
             
             
               14451.900193458 
               1.2911412485443 
               89.925226766103 
               202.38601311481 
             
             
               14666.374034128 
               1.5791121136302 
               89.920912928287 
               214.47408232363 
             
             
               14895.331679791 
               1.9032834713067 
               89.916841886428 
               228.95793773649 
             
             
               15140.969810986 
               2.2682060446596 
               89.912919830784 
               245.63848268456 
             
             
               15405.419662153 
               2.6790026230402 
               89.909073722763 
               264.45027291130 
             
             
               15690.815558253 
               3.1414399758701 
               89.905249591321 
               285.39640111885 
             
             
               15999.342796423 
               3.6620098066495 
               89.901403326605 
               308.52784207464 
             
             
               16333.264483195 
               4.2480198836023 
               89.897496456714 
               333.92240826494 
             
             
               16694.940504717 
               4.9076966262520 
               89.893495272962 
               361.67688300182 
             
             
               17086.841255695 
               5.6503005880400 
               89.889368088848 
               391.90177926706 
             
             
               17511.557461466 
               6.4862564561342 
               89.885085260722 
               424.71743298422 
             
             
               17971.807917113 
               7.4272993933662 
               89.880617598864 
               460.25192023864 
             
             
               18470.445105806 
               8.4866397766488 
               89.875936458956 
               498.63893670373 
             
             
               19010.460493068 
               9.6791486444739 
               89.871013128336 
               540.01747381047 
             
             
               19594.988697452 
               11.021566456806 
               89.865818508043 
               584.53069546048 
             
             
               20227.311515111 
               12.532738097948 
               89.860322976641 
               632.32579225008 
             
             
               20910.862179678 
               14.233877421347 
               89.854496193117 
               683.55421721182 
             
             
               21649.229995685 
               16.148865050039 
               89.848306959846 
               738.37205984105 
             
             
               22446.165772454 
               18.304583613239 
               89.841723136096 
               796.94084715816 
             
             
               23305.588232502 
               20.731295125158 
               89.834711499538 
               859.42851892984 
             
             
               24231.591613073 
               23.463065803670 
               89.827237674937 
               926.01062164312 
             
             
               25000.000000000 
               25.821244410585 
               89.821091887182 
               768.41497742361 
             
          
         
         
             
          
             
               NEXT 
             
             
                 
             
          
         
       
     
   
   The full suite of refractive path input requirements do make the use of the user-specified refractive geometry option rather complex, but this is currently necessary for modeling strong index of refraction gradients such as those observed in the marine boundary layer. 
   Inhomogeneous Paths 
   This line center transmittance section began with a description of the Plass expression for the transmittance of a homogeneous path and the use of this expression by MODTRAN™ to estimate line center transmittances as a function of both column density and the temperature and pressure dependent band model parameters (pressure for the Lorentz half-width). However, the molecular transmittances of the band model RTE, Eq. (13), are defined for an inhomogeneous path. More specifically, the molecular spectral cross-sections κ v,i   abs (s) and molecular densities n i (s) are functions of path length s along the line-of-sight, and atmospheric conditions change as the path is traversed. In MODTRAN™, this disconnect is rectified by modeling the Earth atmosphere as stratified, with molecular and particulate densities, total pressure, and ambient temperature defined as a function of altitude (there is an option in MODTRAN™ to alternatively define temperature and density profiles as a function of pressure with altitudes determined from the hydrostatic equation and a user defined ground altitude). The layering is assumed to be relatively finely gridded, so that a homogeneous path can be defined within individual layers. 
   Implementation of the IMPROVEMENTS TO MODTRAN4™ user-specified refractive geometry option required that a completely new approach for calculating MODTRAN™ column densities be incorporated into the model. In MODTRAN4™, column density calculations are performed coincident with the refractive geometry calculation; as the refractive path is incremented in small steps, ds, and the cosine of the path zenith is updated, column amounts are summed by modeling the short incremental segments as a straight lines. With the user-specified refractive geometry option, the full path length within a layer is provided, but not the detailed mapping of path length as a function of altitude within the layer. Like earlier MODTRAN™ versions, IMPROVEMENTS TO MODTRAN4™ assumes an exponential variation of densities and pressure vertically between altitude levels and a linear variation of temperature vertically between altitude levels. However, IMPROVEMENTS TO MODTRAN4™ also models the refracted path between altitude levels as a circular arc, typically with an extremely large radius of curvature—much larger than the Earth radius for normal refraction: 
                     u   σ     =       ρ     l   σ       ⁢       ∫     s     σ   -   1         s   σ       ⁢         (       ρ       l   σ     -   1       /     ρ     l   σ         )         h   σ     ⁡     (   s   )         ⁢     ⅆ   s             ⁢     
     ⁢           P   σ     _     =           ρ     l   σ     air     ⁢     P     l   σ           u   σ   air       ⁢       ∫     s     σ   -   1         s   σ       ⁢         (         ρ     l   σ     air     ⁢     P       l   σ     -   1             ρ     l   σ     air     ⁢     P     l   σ           )         h   σ     ⁡     (   s   )         ⁢     ⅆ   s             ;     ⁢     
     ⁢         h   σ     ⁡     (   s   )       ≡         H   ⁡     (   s   )       -     H     l   σ             H       l   σ     -   1       -     H     l   σ             ⁢     
     ⁢         T   σ     _     =         ρ     l   σ     air       u   σ   air       ⁢       ∫     s     σ   -   1         s   σ       ⁢           (       ρ       l   σ     -   1     air       ρ     l   σ     air       )         h   σ     ⁡     (   s   )         ⁡     [       T     l   σ       +         h   σ     ⁡     (   s   )       ⁢     (       T       l   σ     -   1       -     T     l   σ         )         ]       ⁢     ⅆ   s                     (   34   )               
In these path integrals, s is the path length variable; s σ−1  and s σ  are path lengths to the beginning and end of segment σ, H(s) is the altitude at path length s; l σ  is index of the vertical layer containing segment σ, H l     σ   , ρ l     σ   , P l     σ    and T l     σ    are the altitude, density, pressure and temperature at the bottom of vertical layer l σ ; and H l     σ     −1 , ρ l     σ     −1 , P l     σ     −1  and T l     σ     −1  are the same quantities but at the top of vertical layer l σ .
 
   If Beer&#39;s Law were applicable, the LOS molecular line center transmittance could be computed by simply stringing together individual segment transmittances for each species. However, as detailed above, the in-band absorption does not retain a linear dependence on column density as that density increases. The Curtis-Godson (CG) approximation provides a solution to the band model inhomogeneous path problem. The approach involves approximating the inhomogeneous path by a quasi-equivalent homogeneous path. First, a total path in-band weak-line optical depth (τ cen )   v     weak-line  for the band centered at  v  is computed: 
                           (     τ   cen     )       v   _       weak   -   line       ≡       ⁢       ∑     σ   ⁢           ⁢   segments               ⁢       n   ⁡     (       T   σ     _     )       ⁢     S   ⁡     (       T   σ     _     )       ⁢     u   σ                     =       ⁢       ∑     σ   ⁢           ⁢   segments               ⁢       [       ∑   l             ⁢       S   l     ⁡     (       T   σ     _     )         ]     ⁢       u   σ     .                       (   35   )               
This value is used to normalized a CG path-averaged effective line number,  n   
                   n   _     ≡         ∑     σ   ⁢           ⁢   segments               ⁢       [       n   ⁡     (       T   σ     _     )       ⁢     S   ⁡     (       T   σ     _     )       ⁢     u   σ       ]     ⁢     n   ⁡     (       T   σ     _     )               (     τ   cen     )       v   _       weak   -   line                 (   36   )               
The CG single-line optical depth parameter  Su  is defined as total path in-band weak-line optical depth normalized by the number of lines:
 
                   Su   _     ≡         (     τ   cen     )       v   _       weak   -   line         n   _               (   37   )               
Finally, the CG Lorentz and Doppler line-width parameters are defined as the optical depth weighted sum of the n line widths in each band normalized by the product of the weak-line optical depth and the CG effective line number:
 
                   γ   _     ≡         ∑     σ   ⁢           ⁢   segments               ⁢       [       n   ⁡     (       T   σ     _     )       ⁢     S   ⁡     (       T   σ     _     )       ⁢     u   σ       ]     ⁢     n   ⁡     (       T   σ     _     )       ⁢     γ   ⁡     (         T   σ     _     ,       P   σ     _       )                 (     τ   cen     )       v   _       weak   -   line       ⁢     n   _                 (   38   )               
In MODTRAN™, the CG data is inserted into the Plass expression to compute line center transmittance for the multi-segmented paths.
 
                     〈     t   k_cen     〉       inhomogeneous   ⁢           ⁢   path       =         (     1   -         W   st     ⁡     (       Su   _     ,       γ   _     c     ,       γ   _     d       )         Δ   ⁢           ⁢   v         )     n     .             (   39   )               
The CG approximation works well when the bulk of the absorption occurs within a sub-section of the LOS over which the line shape function does not vary drastically. For most ambient atmospheric species, this condition is satisfied independent of the particular LOS. An important exception for the terrestrial atmosphere is ozone. The absolute O 3  concentration has a relatively constant profile from the Earth&#39;s surface up into the stratosphere. However, the character of Voigt line shape changes from Lorentzian to Doppler over this altitude range. Thus, MODTRAN™ models the O 3  in-band absorption with a higher-order CG correction.
 
   Applying CG path averaging for each of the solar illumination paths arising in the spherical Earth single scatter solar calculation was deemed computationally too expensive when these calculations were introduced into MODTRAN™. As illustrated in  FIG. 2 , CG averages can be computed incrementally for the direct LOS as the calculation traverses outward along the path from the sensor; however, the LOS level intersection to sun paths are each independent because of the Earth&#39;s curvature, even with horizontal homogeneity of the spherical shell layering. To simplify the solar illumination calculations, MODTRAN™ approximates the illumination path at each LOS altitude level intersection (and the tangent point) as a single homogeneous segment using cumulative column densities and species-dependent effective temperatures and pressures: 
                       P   s   k     _     =       ∑           σ   ⁢           ⁢   segments               s   ⁢           ⁢   to   ⁢           ⁢   sum                     ⁢       u   σ   k     ⁢         P   σ     _     /       ∑           σ   ⁢           ⁢   segments               s   ⁢           ⁢   to   ⁢           ⁢   sum                     ⁢     u   σ   k               ⁢     
     ⁢         T   s   k     _     =       ∑           σ   ⁢           ⁢   segments               s   ⁢           ⁢   to   ⁢           ⁢   sum                     ⁢       u   σ   k     ⁢         T   σ     _     /       ∑           σ   ⁢           ⁢   segments               s   ⁢           ⁢   to   ⁢           ⁢   sum                     ⁢     u   σ   k                       (   40   )               
These single illumination path layers are then coupled to the LOS using CG path averaging to calculate the L-shaped path transmittances.
 
Line Tail Transmittance
 
   The MODTRAN™ spectral bin line tail transmittance is the attenuation arising from molecular transitions centered outside of the bin but no more than 25 cm −1  from the bin center. Unlike the line center spectral attenuation, the line tail absorption varies smoothly across each spectral bin with at most a single minimum. Thus, the entire line tail spectral absorption cross-section curve can be accurately fit to temperature and pressure dependent Padé approximants. This fitting procedure is described in detail within the parent patent application. 
   For the IMPROVEMENTS TO MODTRAN4™ invention, the Padé parameter database is generated at pressures of 1.0 and 0.1 Atm and at the 6 evenly spaced MODTRAN™ temperatures between 180 and 305K. The Padé parameter interpolations are linear in temperature and linear in pressure-squared. The P 2  dependence arises because of the γ c   2  term in the denominator of the Voigt line-shape, Eq. (22). 
   An IMPROVEMENTS TO MODTRAN4™ parameter defines the number of monochromatic line tail absorption cross-section calculations per spectral bin so the line tails can be computed at arbitrarily fine spectral resolution. Nominally, the parameter is set to 5. The monochromatic line tail optical depths are added for multiple segment paths since Beer&#39;s Law applies, and these summed optical depths are modeled as varying linearly between spectral points for computing the in-band line tail transmittances. 
   Line Continuum Transmittance 
   The primary atmospheric molecules for which lines tails beyond 25 cm −1  from line center sum to make a significant contribution to atmospheric absorption are H 2 O and CO 2 . MODTRAN™ adopts the continuum algorithms employed in standard LBL models for both of these molecules. 
   The temperature dependent H 2 O continuum spectral absorption cross-section κ v   H     2     O cont (T) is partitioned into a foreign and a temperature-dependent self-broadened component: 
                   κ   v       H   2     ⁢   O   ⁢           ⁢   cont       =     v   ⁢           ⁢     tanh   ⁡     (         c   2     ⁢   v       2   ⁢   T       )       ×     {         (         ρ   total     -     ρ       H   2     ⁢   O           ρ   total       )     ⁢       C   v       H   2     ⁢   O   ⁢           ⁢   foreign       ⁡     (     296   ⁢   K     )         +       (       ρ       H   2     ⁢   O         ρ   total       )     ⁡     [         (       T   -     260   ⁢   K         36   ⁢   K       )     ⁢       C   v       H   2     ⁢   O   ⁢           ⁢   self       ⁡     (     296   ⁢   K     )         +       (         296   ⁢   K     -   T       36   ⁢   K       )     ⁢       C   v       H   2     ⁢   O   ⁢           ⁢   self       ⁡     (     260   ⁢   K     )           ]         }               (   44   )               
The spectral foreign (C v   H     2     O foreign ) and self (C v   H     2     O self ) broadened H 2 O continuum parameters [(cm 2 /mol)/cm −1 ] are tabulated from 0 to 10,000 cm −1  on a 5 cm −1  grid; the self-broadened term endpoint values are used for temperatures below 260K or above 296K rather than extrapolation. The H 2 O continuum data in MODTRAN™ are CKD 2.4, i.e. Version 2.4 of the Clough-Kneizys Data.
 
   The form of the CO 2  continuum spectral absorption cross-section is simpler than that of H 2 O. CO 2  has a zero dipole moment, so partitioning into foreign and self broadened components is unnecessary: 
                     κ   v       CO   2     ⁢           ⁢   cont       ⁡     (   T   )       =     v   ⁢           ⁢     tanh   ⁡     (         c   2     ⁢   v       2   ⁢   T       )       ⁢     C   v     CO   2                 (   45   )               
The CO 2  data is tabulated on a 10 cm −1  grid from 0 to 10,000 cm −1 . Since the CO 2  continuum cross-section values are themselves not functions of mole fraction, they can and are directly incorporated into the CO 2  line tail band model parameters.
 
Auxiliary Species
 
   All versions of MODTRAN™ compute the absorption and emission resulting from a standard set of atmospheric molecules. For MODTRAN4™, this standard set includes 26 molecules. However, the invention (IMPROVEMENTS TO MODTRAN4™) provides a new auxiliary species option which couples an arbitrary number of additional molecular species into MODTRAN&#39;s radiation transport. These molecules can be modeled either as Beer&#39;s Law absorption cross-section species or via the IMPROVEMENTS TO MODTRAN4™ band model. The band model provides greater accuracy if optical depths are substantial, while cross-section data is generally more readily available. The IMPROVEMENTS TO MODTRAN4™ database includes band model parameter data for 24 HITRAN species (OH, HF, HCl, HBr, HI, ClO, OCS, H 2 CO, HOCl, N 2 , HCN, CH 3 Cl, H 2 O 2 , C 2 H 2 , C 2 H 6 , PH 3 , COF 2 , H 2 S, HCOOH, HO 2 , O, NO + , HOBr and C 2 H 4 ) not among the MODTRAN™ set of standard molecules (actually, N 2  is one of the MODTRAN™ standard atmospheric molecules, but it is modeled in MODTRAN™ as a cross-section species, not as a band model molecule). 
   In the IMPROVEMENTS TO MODTRAN4™ invention, auxiliary species molecular densities can be defined via user-specified atmospheric vertical profiles or with ‘default’ model profiles; currently, profiles are provided for 16 of the new band model species (OH through PH 3  in the above list). The molecular species temperature-dependent optical data inputs must conform to MODTRAN™ definitions. Thus, if molecular band model data are user-supplied, this data must be defined in accordance with the MODTRAN™ band model formalism [Eqs. (23), (24) and (41)]; similarly, molecular absorption cross-section data must be defined at the 6 MODTRAN™ temperatures (180, 205, 230, 255, 280 and 305K) and in units of cm 2 /molecule. 
   The mechanism for introducing auxiliary species is similar to that for MODTRAN™ standard atmospheric molecular species when user-defined profiles or radiosonde data are being entered. For each atmospheric altitude profile level, the density of each absorption species is specified. Unlike the standard molecules, the number of auxiliary species is read from the profile table header. This number is a TOTAL number of auxiliary species in the user&#39;s master list. The input to the IMPROVEMENTS TO MODTRAN4™ invention accepts a negative sign “−” prefix for any species name to indicate the species is to be excluded from the current calculation. This input format was instituted to avoid having to completely reconfigure an input stream each time a species is added or removed. A star (*) at the end of species name indicates band model rather than absorption cross-section data is to be used. The optical data for each species is stored in species-dependent auxiliary input files (SF6.BM for a SF 6  absorption cross-section data file, SF6.cBM for a SF 6  line center band model data file, and SF6.tBM for a SF 6  line tail band model data file). 
   To illustrate the auxiliary species option, IMPROVEMENTS TO MODTRAN4™ Top-Of-Atmosphere (TOA) nadir view calculations were run with a 1 km thick surface layer of HCl embedded into the Mid-Latitude Summer (MLS) model atmosphere. The auxiliary species 0.1 cm −1  band model data files, p1_HCl.cBM and p1_HCl.tBM, were read in. Calculations could have been performed with the default HCl profile, but the current objective was to predict the change in spectral signatures with increasing HCl column density. HCl densities of 0.1, 1.0 and 10.0 ppm were input. The transmittance spectra for the HCl component are illustrated in  FIG. 3 . On the left, the envelope of R and P branch transitions from the (0-1) vibrational transition are displayed. At 0.2 cm −1  spectral resolution, the splitting of the HCl 35  and HCl 37  isotopes is resolved. This is illustrated for R2 and P3 transitions on the right in  FIG. 3 . 
   The transmittance spectra for the HCl fundamental R2 and P3 vibrational/rotational lines with the HCl surface layer embedded in the MLS atmosphere are shown in  FIG. 4  (left). Not surprisingly, the spectra are considerably more complex than the isolated spectra because of overlap with other molecular species, most notably CH 4  and H 2 O. 
   The thermal emission spectra are also plotted in  FIG. 5  (right). In this simulation, the grey-body (emissivity=0.95) surface is modeled with a temperature 10K below the surface air temperature. As a result, the increase in thermal emission resulting from higher HCl concentrations more than compensates for the lowering of surface emission due to the HCl attenuation. If the HCl layer temperature were coincident with the ground surface temperature, the thermal emission spectra would vary minimally with changes to the HCl column density. 
   Spectrally Varying Aerosol Profiles 
   MODTRAN™ has always included an option to define up to 4 aerosol types, each specified in a distinct altitude region. Nominally, a boundary layer aerosol density profile is defined for altitude levels up to 2 km, a tropospheric aerosol density profile is defined for altitude levels from 3 to 10 km, a stratospheric aerosol density profile is defined for altitude levels from 11 to 30 km, and a high altitude aerosol density profile is defined above the stratopause. Aerosol spectral optical properties within each of these regions do not vary with altitude; the aerosol concentration profile provide the only variation with altitude. A new Spectral Aerosol Profiles (SAP) option has been developed for the IMPROVEMENTS TO MODTRAN4™ invention. The SAP option allows aerosol spectral extinction coefficients, spectral scattering coefficients, and spectral scattering phase functions all to be entered as a functions of altitude. This enables the IMPROVEMENTS TO MODTRAN4™ invention to more realistically capture the true effect of aerosols on atmospheric radiation transport. As part of this development, the Legendre expansion representation of the SAP aerosol phase function is directly coupled into the DISORT multiple scattering (below). In all the other MODTRAN™ aerosol options, default or user-specified, the scattering phase function is always approximated by a Henyey-Greenstein parameterization for use in multiple scattering computations. This is a critical upgrade for accurate inclusion of aerosols in the scattered radiance calculations. 
   Whenever the SAP option is exercised within IMPROVEMENTS TO MODTRAN4™, spectral aerosol profile data is read in from an external file, a sample of which is included in Table 2(This example is kindly been provided by Drs. Denis Dion and Vincent Ross, Defence Research and Development Canada, Valcartier, 2005). The first line of the SAP file contains 3 integers: the number of spectral grid points, the number of scattering phase function Legendre moments, and the number of scattering phase function angular points. The second line contains the scattering phase function angular grid in degrees. The next group of lines contains the aerosol data for the first MODTRAN altitude level. There is a separate line for each wavelength. The individual lines contain the altitude level, the spectral wavelength, extinction and scattering coefficients, and phase function data. Legendre moments are listed first followed by phase function values at each of the angular grid points. Subsequent groups of lines contain the data for additional altitude levels. 
   
     
       
         
             
           
             
               TABLE 2 
             
           
          
             
                 
             
             
               Sample Spectral Aerosol Profiles (SAP) Input Data File. 
             
          
         
         
             
             
             
          
             
               40 64 181 
                 
               Phase Function vs Angle (°) 
             
          
         
         
             
             
             
             
             
             
             
             
             
             
             
             
          
             
               altitude (m) 
               wavlen 
               Bext(1/km) 
               Bsca(1/km) 
               Pmom1 
               Pmom2 
               . . . 
               Pmom64 
               0.0000 
               1.0000 
               . . . 
               180.0000 
             
             
                 
             
          
         
         
             
             
             
             
             
             
             
             
             
             
             
             
          
             
               .0003653 
               .2000 
               .9173681 
               .9158943 
               1.0000 
               .7966921 
               . . . 
               .1319957 
               208.9782 
               193.6871 
               . . . 
               .7423355 
             
             
               .0003653 
               .3000 
               .8726480 
               .8726359 
               1.0000 
               .8118536 
               . . . 
               .1070675 
               208.8804 
               194.0608 
               . . . 
               .3700137 
             
             
               .0003653 
               .3371 
               .8525812 
               .8525700 
               1.0000 
               .8112088 
               . . . 
               .1012243 
               207.9938 
               193.4231 
               . . . 
               .3866308 
             
             
               . . . 
               . . . 
             
             
               .0003653 
               25.00 
               .2620940 
               .0969520 
               1.0000 
               .7294320 
               . . . 
               .0000000 
               38.06391 
               37.61441 
               . . . 
               .1531952 
             
             
               .0003653 
               30.00 
               .2369326 
               .0895308 
               1.0000 
               .7221835 
               . . . 
               .0000000 
               29.82364 
               29.57460 
               . . . 
               .1475345 
             
             
               .0003653 
               40.00 
               .2156776 
               .0699561 
               1.0000 
               .6948593 
               . . . 
               .0000000 
               21.41463 
               21.30646 
               . . . 
               .1690632 
             
             
               .1260741 
               .2000 
               .3358318 
               .3343440 
               1.0000 
               .7796006 
               . . . 
               .0807715 
               172.1363 
               160.5204 
               . . . 
               .6840323 
             
             
               .1260741 
               .3000 
               .3008995 
               .3008881 
               1.0000 
               .7767678 
               . . . 
               .0675285 
               161.2990 
               150.9579 
               . . . 
               .6125076 
             
             
               .1260741 
               .3371 
               .2897007 
               .2896896 
               1.0000 
               .7733779 
               . . . 
               .0646587 
               157.8538 
               147.8797 
               . . . 
               .6019086 
             
             
               . . . 
               . . . 
             
             
               .1260741 
               25.00 
               .0616450 
               .0247540 
               1.0000 
               .6856267 
               . . . 
               .0000000 
               17.56935 
               17.49937 
               . . . 
               .1463460 
             
             
               .1260741 
               30.00 
               .0562292 
               .0217360 
               1.0000 
               .6611843 
               . . . 
               .0000000 
               14.18904 
               14.14785 
               . . . 
               .1608776 
             
             
               .1260741 
               40.00 
               .0584414 
               .0161276 
               1.0000 
               .5706048 
               . . . 
               .0000000 
               10.02048 
               10.00178 
               . . . 
               .2772588 
             
             
               . . . 
               . . . 
             
             
               1993.880 
               .2000 
               .2613175 
               .2602409 
               1.0000 
               .7697431 
               . . . 
               .0518847 
               152.5157 
               142.7655 
               . . . 
               .6152996 
             
             
               1993.880 
               .3000 
               .2297381 
               .2297099 
               1.0000 
               .7593464 
               . . . 
               .0413473 
               136.5168 
               128.3560 
               . . . 
               .5908967 
             
             
               1993.880 
               .3371 
               .2196234 
               .2196121 
               1.0000 
               .7543240 
               . . . 
               .0390868 
               131.9052 
               124.1712 
               . . . 
               .5874675 
             
             
               . . . 
               . . . 
             
             
               1993.880 
               25.00 
               .0295673 
               .0118886 
               1.0000 
               .6489874 
               . . . 
               .0000000 
               13.93146 
               13.88930 
               . . . 
               .1702973 
             
             
               1993.880 
               30.00 
               .0271156 
               .0099214 
               1.0000 
               .6250932 
               . . . 
               .0000000 
               11.57627 
               11.55026 
               . . . 
               .1904991 
             
             
               1993.880 
               40.00 
               .0317162 
               .0074431 
               1.0000 
               .5096798 
               . . . 
               .0000000 
               8.037910 
               8.026117 
               . . . 
               .3562035 
             
             
                 
             
          
         
       
     
   
   In the RTE, the absorption and scattering coefficients, Eq. (3), and the effective phase function, Eq. (5), are all express as sums of individual species contribution. When the SAP options is exercised in IMPROVEMENTS TO MODTRAN4™, the SAP aerosol is treated as one of these species. The extinction coefficient, scattering coefficient, and required phase function profile data are spectrally interpolated (linearly in wavelength) at the central bin frequency within the band model spectral loop. The extinction and scattering coefficients are modeled as varying exponentially with altitude, while the scattering phase functions and the Legendre coefficient are linearly interpolated in altitude. The result is a first principles treatment of the aerosol contributions. 
   Multiple Scattering via the DISORT Algorithm 
   MODTRAN™ computes LOS spectral radiant intensities by integrating the attenuated spectral source function along the path defined by its spherical refractive geometry algorithm and summing the boundary term if termination occurs at the hard Earth. This is the essence of Eq. (1). However, since the multiple scattering component of the source function [Eq. (4)] and the downward flux component [Eq. (14)] of the boundary term both depend on LOS spectral radiant intensities themselves, there is an internal coupling that complicates the path integration procedure. MODTRAN™ resolves the self-consistency problem by using a plane-parallel multiple scattering algorithm to compute the LOS spectral radiant intensity dependent contributions. Underlying this procedure is the assumption that multiple scattering in the terrestrial atmosphere is predominately a lower atmosphere effect and therefore well characterized with the locally flat Earth approximation. This ansatz is least accurate for low sun conditions. Thus, for solar elevation angles below ˜10°, MODTRAN™ solar radiances carry a relatively high uncertainty. 
   The MODTRAN™ model includes two multiple scattering algorithms. An analytic two-stream algorithm dubbed the Isaacs approach approximates the multiple scattering source terms with minimal computations. These two-stream results are best suited for calculation of spatially integrated radiances (fluxes) and thermal scatter contributions. However, the accuracy of the 2-stream solar scatter radiant intensities is relatively poor, typically ±20%, but factor of 2 errors are not uncommon in the shortwave end of visible spectral region (blue). For higher fidelity results, MODTRAN™ employs the DISORT discrete ordinate algorithm to compute multiple scatter radiances. In this section, the DISORT algorithm and its IMPROVEMENTS TO MODTRAN4™ interface are described. 
   Chandrasekhar&#39;s Solution to the Plane-Parallel Atmosphere RTE 
   DISORT explicitly solves the integral-differential RTE for a horizontally-homogeneous plane-parallel atmosphere with no diffuse radiation impingent on the top of atmosphere (TOA) and with an emitting and reflecting ground surface: 
                     μ   ⁢       ⅆ     I   ⁡     (     Ω   ;     τ   ↓       )           ⅆ     τ   ↓           =       I   ⁡     (     Ω   ;     τ   ↓       )       -     J   ⁡     (     Ω   ;     τ   ↓       )           ⁢     
     ⁢       Source   ⁢           ⁢   Function     _     ⁢     
     ⁢       J   ⁡     (     Ω   ;     τ   ↓       )       ≡         (     1   -     ω     τ   ↓         )     ⁢     B   ⁡     (     T     τ   ↓       )         +       ω     τ   ↓       [           ⁢         p   ⁡     (     Ω   ,       Ω   sun     ;     τ   ↓         )       ⁢     ⅇ       τ   ↓     /     μ   sun         ⁢     F   sun       +       ∫     4   ⁢   π       ⁢       p   ⁡     (     Ω   ,       Ω   ′     ;     τ   ↓         )       ⁢     I   ⁡     (       Ω   ′     ;     τ   ↓       )       ⁢     ⅆ     Ω   ′             ]         ⁢     
     ⁢         ω     τ   ↓       ≡           b   sct     ⁡     (     τ   ↓     )       /     b   ⁡     (     τ   ↓     )         ⁢     
     ⁢       Boundary   ⁢           ⁢   Conditions     _     ⁢     
     ⁢       I   ⁡     (       +   μ     ,     ϕ   ;     τ   g   ↓         )       =         ɛ   ⁡     (       +   μ     ,   ϕ     )       ⁢     B   ⁡     (     T   g     )         +       μ   sun     ⁢     f   ⁡     (       +   μ     ,   ϕ   ,     μ   sun     ,     ϕ   sun       )       ⁢     ⅇ       -     τ   g   ↓       /     μ   sun         ⁢     F   sun               ;     ⁢     
     ⁢       I   ⁡     (       -   μ     ,     ϕ   ;   0       )       =     0   +       ∫   0   1     ⁢       ∫   0     2   ⁢   π       ⁢       I   ⁡     (       -     μ   ′       ,       ϕ   ′     ;     τ   g   ↓         )       ⁢     f   ⁡     (       +   μ     ,   ϕ   ,     -     μ   ′       ,     ϕ   ′       )       ⁢     μ   ′     ⁢     ⅆ     ϕ   ′       ⁢       ⅆ     μ   ′       .                       (   46   )               
The top equation in this set is the path optical depth derivative of the integral RTE, Eq. (1), for a plane-parallel atmosphere. Since refraction and spherical Earth effects are excluded, the path solid angle Ω[≡(μ,φ)] is constant for a given LOS. The reference to position has also been dropped as an argument of the BRDF because the Earth surface is modeled as horizontally homogeneous. For notational simplicity, the spectral frequency subscript has also been dropped from Eq. (46) and all subsequent equations in this section, even though all quantities other than the temperatures and angles are spectrally dependent. Finally, the path optical depth is represented as the quotient of the path nadir optical depth measured from the TOA τ ↓  and the cosine of the path zenith angle with the ground value, τ g   ↓ , equal to the total atmosphere nadir extinction optical depth. Any remaining references to position s have been replaced by references to the nadir optical depth τ ↓  because of the lack of horizontal variability.
 
   In order to facilitate solving the plane parallel atmosphere radiation transport problem, two additional model assumptions are made up front: (1) the scattering phase function is modeled as a function of the angle Θ between the incident and outgoing beams, p(Ω,Ω′;τ ↓ )=p(Θ;τ ↓ ); and (2) the surface BRDF is modeled as a function of the incoming and outgoing zenith angles and of the relative azimuth angle, but not of the incoming and outgoing azimuth angles independently, i.e., f(Ω,Ω′)=f(μ,−μ′,φ−φ′). These assumptions limit DISORT&#39;s ability to precisely model physical phenomena such as preferential orientation of ice particles in cirrus clouds and plowed fields on the Earth&#39;s surface. However, under most terrestrial atmospheric scenarios, it is quite reasonable to model suspended particles as randomly oriented and the Earth surface as locally isotropic. 
   DISORT follows a 4-step approach, originally laid out by Chandrasekhar, to solve the integral-differential RTE. In developing the DISORT model, Stamnes and his coworkers solved the previous unsolved and challenging problem of creating a numerically stable realization of the Chandrasekhar approach. 
   Step 1: Factoring the Azimuth Dependence 
   The azimuth dependence of the RTE is factored from the integral-differential RTE by expanding the radiance and bidirectional reflectivity in Fourier cosine series, and expanding the scattering phase function in Legendre polynomials, P l  (μ), and unit normalized associated Legendre polynomials, Λ l   m (μ): 
                           I   ⁡     (     Ω   ;     τ   ↓       )       =       ∑     m   =   0         2   ⁢   M     -   1       ⁢         I   m     ⁡     (     μ   ;     τ   ↓       )       ⁢   cos   ⁢           ⁢     m   ⁡     (       ϕ   sun     -   ϕ     )             ;     ⁢     
     ⁢     f   ⁡     (     μ   ,     -     μ   ′       ,     ϕ   -     ϕ   ′         )       =       ∑     m   =   0         2   ⁢   M     -   1       ⁢         f   m     ⁡     (     μ   ,     -     μ   ′         )       ⁢   cos   ⁢           ⁢     m   ⁡     (     ϕ   -     ϕ   ′       )             ;     ⁢     
     ⁢   and   ⁢     
     ⁢             p   ⁡     (     Θ   ;     τ   ↓       )       =       ⁢       ∑     ℓ   =   0         2   ⁢   M     -   1       ⁢       (       2   ⁢   ℓ     +   1     )     ⁢       g     τ   ↓     ℓ       4   ⁢   π       ⁢       P   ℓ     ⁡     (     cos   ⁢           ⁢   Θ     )                       =       ⁢       ∑     ℓ   =   0         2   ⁢   M     -   1       ⁢       (       2   ⁢   ℓ     +   1     )     ⁢       g     τ   ↓     ℓ       4   ⁢   π                           ⁢       {           P   ℓ     ⁡     (   μ   )       ⁢       P   ℓ     ⁡     (     μ   ′     )         +     2   ⁢       ∑     m   =   1     ℓ     ⁢         Λ   ℓ   m     ⁡     (   μ   )       ⁢       Λ   ℓ   m     ⁡     (     μ   ′     )       ⁢     cos   ⁡     (     ϕ   -     ϕ   ′       )               }     .                     (   47   )               
In these equations, the I m (μ;τ ↓ ) are the azimuth moment radiant intensities, the f m (μ;−μ′) are the azimuth moment BRDFs, the g τ     ↓     l  are the scattering phase function Legendre moments, and 2M is the number of Legendre moments (l=0 to 2M−1). These expansions yield independent/orthogonal integral differential equations for each m component:
 
                     μ   ⁢       ⅆ       I   m     ⁡     (     μ   ;     τ   ↓       )           ⅆ     τ   ↓           =         I   m     ⁡     (     μ   ;     τ   ↓       )       -         J   m     ⁡     (     μ   ;     τ   ↓       )       ⁢     (       m   =   0     ,   1   ,   …   ⁢           ,       2   ⁢   M     -   1       )           ⁢     
     ⁢       Source   ⁢           ⁢   Function     _     ⁢     
     ⁢         J   m     ⁡     (     μ   ;     τ   ↓       )       =         ω     τ   ↓       ⁢     {         (     2   -     δ     m   ⁢           ⁢   0         )     ⁢       F   sun     ⁡     [       ∑     ℓ   =   m         2   ⁢   M     -     1   ′         ⁢         (     -   1     )       ℓ   +   m       ⁢     (       2   ⁢   ℓ     +   1     )     ⁢       g     τ   ↓     ℓ       4   ⁢   π       ⁢       Λ   ℓ   m     ⁡     (   μ   )       ⁢       Λ   ℓ   m     ⁡     (     μ   sun     )           ]         +       ∫     -   1     1     ⁢           I   m     ⁡     (       μ   ′     ;     τ   ↓       )       ⁡     [       ∑     ℓ   =   m         2   ⁢   M     -     1   ′         ⁢       (       2   ⁢   ℓ     +   1     )     ⁢       g     τ   ↓     ℓ     2     ⁢       Λ   ℓ   m     ⁡     (   μ   )       ⁢       Λ   ℓ   m     ⁡     (     μ   ′     )           ]       ⁢     ⅆ     μ   ′             }       +         δ     m   ⁢           ⁢   0       ⁡     (     1   -     ω     τ   ↓         )       ⁢     B   ⁡     (     T     τ   ↓       )             ⁢     
     ⁢       Boundary   ⁢           ⁢   Conditions     _     ⁢     
     ⁢           I   m     ⁡     (       +   μ     ;     τ   g   ↓       )       =         δ     m   ⁢           ⁢   0       ⁢     ɛ   ⁡     (   μ   )       ⁢     B   ⁡     (     T   g     )         +       μ   sun     ⁢       f   m     ⁡     (     μ   ,     -     μ   sun         )       ⁢     ⅇ       -     τ   g   ↓       /     μ   sun         ⁢     F   sun           ;     ⁢     
     ⁢         I   m     ⁡     (       -   μ     ;   0     )       =     0   +       (     1   +     δ     m   ⁢           ⁢   0         )     ⁢       ∫   0   1     ⁢         I   m     ⁡     (       -     μ   ′       ;     τ   g   ↓       )       ⁢       f   m     ⁡     (       +   μ     ,     -     μ   ′         )       ⁢     μ   ′     ⁢     ⅆ     μ   ′                         (   48   )               
Note that in the absence of solar radiation, only the m=0 problem needs to be solved. This is also true if there is azimuthal symmetry resulting from either a zenith sun or a vertical/nadir LOS.
 
Step 2: Introducing Discrete Ordinates
 
   In the second step of the Chandrasekhar procedure, the independent integral-differential equations are each converted into a system of 2N coupled ordinary differential equations (ODE&#39;s) by applying the discrete ordinate approximation to the cosine zenith integrals. 
                         μ   i     ⁢       ⅆ       I   m     ⁡     (       μ   1     ;     τ   ↓       )           ⅆ     τ   ↓           =         I   m     ⁡     (       μ   i     ;     τ   ↓       )       -         J   m     ⁡     (       μ   i     ;     τ   ↓       )       ⁢     (       i   =     ±   1       ,     ±   2     ,   …   ⁢           ,     ±   N       )           ;     ⁢     
     ⁢       Source   ⁢           ⁢   Function     _     ⁢     
     ⁢           J   m     ⁡     (     μ   ;     τ   ↓       )       =         ω     τ   ↓       ⁢     {         (     2   -     δ     m   ⁢           ⁢   0         )     ⁢       F   sun     ⁡     [       ∑     ℓ   =   m         2   ⁢   M     -     1   ′         ⁢         (     -   1     )       ℓ   +   m       ⁢     (       2   ⁢   ℓ     +   1     )     ⁢       g     τ   ↓     ℓ       4   ⁢   π       ⁢       Λ   ℓ   m     ⁡     (   μ   )       ⁢       Λ   ℓ   m     ⁡     (     μ   i   sun     )           ]         +       ∑       j   =     -   N         j   ≠   0       N     ⁢       w   j     ⁢         I   m     ⁡     (       μ   j     ;     τ   ↓       )       ⁡     [       ∑     ℓ   =   m         2   ⁢   M     -     1   ′         ⁢       (       2   ⁢   ℓ     +   1     )     ⁢       g     τ   ↓     ℓ     2     ⁢       Λ   ℓ   m     ⁡     (   μ   )       ⁢       Λ   ℓ   m     ⁡     (     μ   j     )           ]             }       +         δ     m   ⁢           ⁢   0       ⁡     (     1   -     ω     τ   ↓         )       ⁢     B   ⁡     (     T     τ   ↓       )             ;     ⁢     
     ⁢       Boundary   ⁢           ⁢   Conditions     _     ⁢     
     ⁢           I   m     ⁡     (       +     μ   i       ;     τ   g   ↓       )       =         δ     m   ⁢           ⁢   0       ⁢     ɛ   ⁡     (     μ   i     )       ⁢     B   ⁡     (     T   g     )         +       μ   sun     ⁢       f   m     ⁡     (       μ   i     ,     -     μ   sun         )       ⁢     ⅇ       -     τ   g   ↓       /     μ   sun         ⁢     F   sun           ;     ⁢     
     ⁢         I   m     ⁡     (       -     μ   i       ;   0     )       =     0   +       (     1   +     δ     m   ⁢           ⁢   0         )     ⁢       ∑     j   =   1     N     ⁢       w   j     ⁢     μ   j     ⁢       I   m     ⁡     (       -     μ   j       ;     τ   g   ↓       )       ⁢       f   m     ⁡     (       +   μ     ,     -     μ   j         )                 ⁢     
     ⁢     (       i   =     +   1       ,     +   2     ,   …   ⁢           ,     +   N       )             (   49   )               
DISORT employs a double Gaussian quadrature, with paired cosine angles, μ j =−μ −j , symmetric in each hemisphere about μ=½ and having a common Gaussian weight, w j =w −j .
 
Step 3: Simplifying Optical Depth Dependence via Atmospheric Layering
 
   In the third step of the Chandrasekhar procedure, the optical depth dependence of the atmospheric single scattering albedo, the scattering phase function and the Planck emission is parameterized so that an analytic matrix solution to the system of ODE&#39;s can be derived. This is accomplished by partitioning the atmosphere into l=1 to L layers with l=0 to L level boundaries. Altitudes z 0 =TOA to z L =ground, nadir optical depths τ 0   ↓ =0 to τ L   ↓ =τ g   ↓ , and temperatures T 0  to T L  are specified at each altitude level. Scattering properties are modeled as constants within each layer: 
                       ω   l     ≡       ∫     z     l   -   1         z   l       ⁢         b   sct     ⁡     (   z   )       ⁢       ⅆ   z     /       ∫     z     l   -   1         z   l       ⁢       b   ⁡     (   z   )       ⁢     ⅆ   z                 ;     ⁢     
     ⁢         g   l   ℓ     ≡       ∫     z     l   -   1         z   l       ⁢         g     τ   ↓     ℓ     ⁡     (   z   )       ⁢       b   sct     ⁡     (   z   )       ⁢       ⅆ   z     /       ∫     z     l   -   1         z   l       ⁢         b   sct     ⁡     (   z   )       ⁢     ⅆ   z                 ,             (   50   )               
while Planck emission is characterized as varying linearly with nadir optical depth through each layer:
 
                         B   l     ⁡     (     T     τ   ↓       )       =       B   ⁡     (     T     l   -   1       )       +           τ   ↓     -     τ     l   -   1     ↓           τ   l   ↓     -     τ     l   -   1     ↓         ⁡     [       B   ⁡     (     T   l     )       -     B   ⁡     (     T     l   -   1       )         ]           ;     ⁢     
     ⁢       τ     l   -   1     ↓     ≤     τ   ↓     ≤       τ   l   ↓     .               (   51   )               
A general solution can be derived for the m-dependent RTE within each atmospheric layer given this parameterization of optical depth dependence:
 
                         I   l   m     ⁡     (       μ   i     ;     τ   ↓       )       =         δ     m   ⁢           ⁢   0       ⁡     [         Y     0   ⁢   l       ⁡     (     μ   i     )       +         Y     1   ⁢   l       ⁡     (     μ   i     )       ⁢     τ   ↓         ]       +         Z   l   m     ⁡     (     μ   i     )       ⁢     ⅇ       -     τ   ↓       /     μ   sun           +       ∑       j   =     -   N         j   ≠   0       N     ⁢       C   jl   m     ⁢       G   jl   m     ⁡     (     μ   i     )       ⁢     ⅇ       -     k   jl   m       ⁢     τ   ↓                 ;     ⁢     
     ⁢       τ     l   -   1     ↓     ≤           ⁢     τ   ↓     ≤           ⁢       τ   l   ↓     .               (   52   )               
In Eq. (52), the Y 0l (μ i ) and Y 1l (μ i ) arise as a particular solution to the thermal source problem and are obtained by solving a system 4N linear equations. Similarly, the Z l   m (μ i ) are a particular solution to the solar source problem, obtained by solving a system of 2N linear equations for each value of m. The k jl   m  and G jl   m  are the eigenvalues and eigenvectors, respectively, that arise in solving the hormogeneous plane-parallel RTE problem [F sun =0; B l =0], and the C jl   m  are the constants of integration that must be determined from boundary and layer continuity conditions.
 
Step 4: Applying Boundary and Continuity Equations
 
   For each of the azimuth moments m, the L atmospheric layers each contain 2N constants of integration. These values are determined from the (2N)×L linear boundary and continuity equations: 
                       I   1   m     ⁡     (       -     μ   i       ;       τ   0   ↓     =   0       )       =   0     ⁢     
     ⁢       (       i   =     +   1       ,     +   2     ,   …   ⁢           ,     +   N       )     ;     ⁢     
     ⁢         I     l   -   1     m     ⁡     (       -     μ   i       ;     τ     l   -   1     ↓       )       =       I   l   m     ⁡     (       -     μ   i       ;     τ     l   -   1     ↓       )         ⁢     
     ⁢       (       l   =   2     ,   3   ,   …   ⁢           ,     L   ;     i   =     ±   1         ,     ±   2     ,   …   ⁢           ,     ±   N       )     ;     ⁢     
     ⁢         I   L   m     ⁡     (       +     μ   i       ;     τ   g   ↓       )       =         δ     m   ⁢           ⁢   0       ⁢     ɛ   ⁡     (     μ   i     )       ⁢     B   ⁡     (     T   g     )         +       μ   sun     ⁢       f   m     ⁡     (       μ   i     ,     -     μ   sun         )       ⁢     ⅇ       -     τ   g   ↓       /     μ   sun         ⁢     F   sun       +       (     1   +     δ     m   ⁢           ⁢   0         )     ⁢       ∑     j   =   1     N     ⁢       w   j     ⁢     μ   j     ⁢       I   m     ⁡     (       -     μ   j       ;     τ   g   ↓       )       ⁢       f   m     ⁡     (       +     μ   i       ,     -     μ   j         )       ⁢       (       i   =     +   1       ,     +   2     ,   …   ⁢           ,     +   N       )     .                       (   53   )               
Solving for these constants is the most computationally time intensive step in DISORT.
 
   Once the constants of integration have been calculated, the specific requested spectral radiant intensities and horizontal fluxes can be computed. The interpolated spectral radiant intensity in direction (±μ, φ) and at an arbitrary nadir optical depth τ ↓  within the plane-parallel atmosphere is computed by inserting the integral RTE into the radiant intensity expression of Eq. (47): 
                       I   ⁡     (       +   μ     ,     ϕ   ;     τ   ↓         )       =       ∑     m   =   0         2   ⁢   M     -   1       ⁢       {           I   m     ⁡     (       +   μ     ;     τ   L   ↓       )       ⁢     ⅇ         τ   ↓     -     τ   L   ↓       μ         +       ∫     τ   ↓       τ   L   ↓       ⁢         J   m     ⁡     (       +   μ     ;     τ     ↓   ′         )       ⁢     ⅇ         τ     ↓   ′       -     τ   ↓       μ       ⁢       ⅆ     τ     ↓   ′         μ           }     ⁢   cos   ⁢           ⁢     m   ⁡     (       ϕ   sun     -   ϕ     )             ;     ⁢     
     ⁢     I   ⁡     (       -   μ     ,     ϕ   ;     τ   ↓         )       =       ∑     m   =   0         2   ⁢   M     -   1       ⁢       {       ∫   0     τ   ↓       ⁢         J   m     ⁡     (       -   μ     ;     τ     ↓   ′         )       ⁢     ⅇ         τ   ↓     -     τ     ↓   ′         μ       ⁢       ⅆ     τ     ↓   ′         μ         }     ⁢   cos   ⁢           ⁢       m   ⁡     (       ϕ   sun     -   ϕ     )       .                 (   54   )               
The azimuth moments of the source function J m (±μ;τ ↓ ) are defined by Eq. (49). Substitution of Eqs. (50), (51) and (52) for the nadir optical depth dependent terms yields simple exponential integrals that are evaluated analytically. For the upward and downward flux calculations, only m=0 terms contribute. The fluxes are computed from the Gaussian quadrature integral representation:
 
                         F   +     ⁡     (     τ   ↓     )       ≅     2   ⁢   π   ⁢       ∑     i   =   1     N     ⁢       w   i     ⁢     μ   i     ⁢       I   0     ⁡     (       +   μ     ;     τ   ↓       )               ;     ⁢     
     ⁢           F   -     ⁡     (     τ   ↓     )       ≅         μ   sun     ⁢     F   sun     ⁢     ⅇ       -     τ   ↓       /     μ   sun           +     2   ⁢   π   ⁢       ∑     i   ⁢           =           ⁢   1     N     ⁢       w   i     ⁢     μ   i     ⁢       I   0     ⁡     (       -   μ     ;     τ   ↓       )                 ;     ⁢     
     ⁢       F   -     ⁢       (     τ   g   ↓     )     ≅         μ   sun     ⁢     F   sun     ⁢     ⅇ     -       τ   g   ↓     /     μ   sun             +     F   ⁢               v   ↓     ⁢     (     s   g     )     .                     (   55   )               
In the last line of Eq. (55), the relationship between the total (direct+diffuse) downward flux at the ground and the diffuse downward flux of Eq. (14) is explicitly presented.
 
The MODTRAN™ Interface to DISORT
 
   In principle, MODTRAN™ could be used solely as the optical property and geometry specifier for DISORT, defining the DISORT inputs as a function of wavelength, and DISORT used as the RT engine, the generator of the RT outputs. However, in this mode, the MODTRAN™ ability to incorporate spherical refractive effects would be lost. As described above, DISORT is instead only used to compute the multiple scattering and surface downward flux contributions to Eq. (13). In this sub-section, the precise interface is laid out. The sub-section begins with a description of DISORT convergence issues, the delta-M method and its impact on the MODTRAN™ interface. The IMPROVEMENTS TO MODTRAN4™ method of defining DISORT inputs is described next, and that is followed by a description of special uses of DISORT via the IMPROVEMENTS TO MODTRAN4™. 
   DISORT Convergence, the Delta-M Method and Radiant Intensities 
   Atmospheric aerosol and cloud particulates typically possess highly forward peaked scattering phase functions. Accurate representation of these functions using a Legendre expansion, Eq. (47), requires hundreds of terms. Computation times grow rapidly with the number of discrete ordinate streams, which DISORT couples to the number of azimuth terms. 
   To alleviate the convergence problem, the delta-M method was introduced. In this method, the scattering phase function is represented as a sum of a forward peaked Dirac-delta function δ(1-cos Θ) and a Legendre expansion: 
                         p   ′     ⁡     (     Θ   ;     τ   ↓       )       =       1     4   ⁢   π       ⁡     [       2   ⁢   f   ⁢           ⁢     δ   ⁡     (     1   -     cos   ⁢           ⁢   Θ       )         +       (     1   -   f     )     ⁢       ∑     ℓ   =   0         2   ⁢   M     -   1       ⁢       (       2   ⁢   ℓ     +   1     )     ⁢     g     τ   ↓     ′ℓ     ⁢       P   ℓ     ⁡     (     cos   ⁢           ⁢   Θ     )               ]         ;     ⁢     
     ⁢         g     τ   ↓     ′ℓ     =           g     τ   ↓     ℓ     -   f       1   -   f       ⁢     (       ℓ   =   0     ,   …   ⁢           ,       2   ⁢   M     -   1       )         ;     ⁢     
     ⁢     f   =       g     τ   ↓       2   ⁢   M       .               (   56   )               
Here, the delta-M method expansion coefficients g′ τ     ↓     l  are defined by requiring preservation of the overall Legendre expansion up to order 2M−1, and, as indicated here, the extra parameter f is set by fitting one additional term from the expansion. The delta-M method changes the assumption used for the higher order Legendre coefficients; instead of setting these coefficients to zero, Eq. (47), they are set to g′ τ     ↓     2M . The Legendre expansion coefficients from observed particulate phase functions typically fall off quite slowly, so the constant value that results from the delta-M method is a more realistic representation. From an alternative perspective, the integrals used to compute the traditional scattering phase function Legendre expansion coefficients are dominated by the forward peak, which yields a poor representation for off-peak angles. By factoring out the forward peak, Eq. (56) allows the phase function in the ε&lt;Θ&lt;180° angular region (ε small) to be much better characterized by the expansion.
 
   With the delta-M method, accurate spectral radiant intensities and fluxes can generally be computed with as few as 4, 8 or 16 streams. Substituting of the delta-M expansion for the phase function into the RTE transforms it into an identical equation, but with p, τ ↓ , and ω replaced by p′, τ′ ↓ , and ω′, where the later two quantities are defined by: 
                       d   ⁢           ⁢     τ     ′   ↓         =       (     1   -       ω   l     ⁢   f       )     ⁢   d   ⁢           ⁢     τ   ↓         ;     ⁢     
     ⁢       ω   l   ′     =         1   -   f       1   -       ω   l     ⁢   f         ⁢       ω   l     .                 (   57   )               
DISORT automatically implements these transformations.
 
   The one obvious drawback of the delta-M method is that the transformed radiant intensities are somewhat inaccurate for very near forward scattering scenarios, such as occurs in the study of the solar aureole, because the finite width of the forward peak is not characterized. To improve near forward scattering accuracy, DISORT provides an option to invoke a series of radiant intensity corrections developed by Nakajima and Tanaka. In these corrections, the plane-parallel atmosphere single scatter solar radiant intensity is explicitly subtracted from the diffuse solar radiant intensity and then added back in using a high order Legendre expansion for the scattering phase function. Additional, less dramatic Nakajima and Tanaka corrections improve higher order solar scattering intensity terms. 
   For DISORT runs from within MODTRAN™, it is important that the Nakajima and Tanaka corrections not be invoked. MODTRAN™ itself subtracts the plane-parallel single scatter solar contribution from the DISORT at-source solar radiant intensity for each layer l, ΔI l (μ,φ): 
                     Downward   ⁢           ⁢   Viewing     _     ⁢     
     ⁢         Δ   ⁢           ⁢       I   l   mss     ⁡     (     +     μ   .   ϕ       )         =       Δ   ⁢           ⁢       I   l   sun     ⁡     (       +   μ     ,   ϕ     )         -       (       μ   sun       μ   +     μ   sun         )     ⁢     F   sun     ⁢     ω   l     ⁢         p   l   ′     ⁡     (   Θ   )       ⁡     [       t       l   -   1     -&gt;   sun       -       ⅇ       -     Δτ   l   ↓       /   μ       ⁢     t     l   -&gt;   sun           ]             ;     ⁢     
     ⁢         Δ   ⁢           ⁢       I   l     ⁡     (       +   μ     ,   ϕ     )         ≡       I   ⁡     (       +   μ     ,     ϕ   ;     τ     l   -   1     ↓         )       -       ⅇ       -     Δτ   l   ↓       /   μ       ⁢     I   ⁡     (       +   μ     ,     ϕ   ;     τ   l   ↓         )             ;     ⁢     
     ⁢       Upward   ⁢           ⁢   Viewing     _     ⁢     
     ⁢         Δ   ⁢           ⁢       I   l   mss     ⁡     (       -   μ     ,   ϕ     )         =       Δ   ⁢           ⁢       I   l   sun     ⁡     (       -   μ     ,   ϕ     )         -       (       μ   sun       μ   -     μ   sun         )     ⁢     F   sun     ⁢     ω   l     ⁢         p   l   ′     ⁡     (   Θ   )       ⁡     [       t     l   -&gt;   sun       -       ⅇ       -     Δτ   l   ↓       /   μ       ⁢     t       l   -   1     -&gt;   sun           ]             ;     ⁢     
     ⁢       Δ   ⁢           ⁢       I   l     ⁡     (       -   μ     ,   ϕ     )         ≡       I   ⁡     (       -   μ     ,     ϕ   ;     τ   l   ↓         )       -       ⅇ       -     Δτ   l   ↓       /   μ       ⁢       I   ⁡     (       -   μ     ,     ϕ   ;     τ     l   -   1     ↓         )       .                   (   58   )               
In this equation, ΔI l   mss (±μ,φ) is the at-source multiple (no single) scatter solar radiant intensity from layer l; ω l  and p′ l (Θ) are the single scattering albedo and delta-M scattering phase function within layer l; Δτ l   ↓  is the nadir optical depth through layer l; t l→sun  is the solar path transmittance from the bottom of layer l; and t l−1→sun  is the solar path transmittance from the top of layer l. MODTRAN™ subsequently adds in the spherical refractive geometry single scatter solar radiant intensity. MODTRAN™ computes the species scattering phase functions either from an analytic form (for Rayleigh and Henyey-Greenstein representations) or from a scattering angle dependent table. For user-defined scattering phase functions, both the Legendre coefficients and the tabulated values are required inputs to insure that the Legendre expansion does not have to be relied upon for the single scatter calculation.
 
   MODTRAN incorporates the DISORT multiple scattering layer contributions into its LOS radiant intensity by first defining a DISORT thermal and multiple scattering solar layer source terms, ΔJ l   thm  and ΔJ l   mss : 
                       Δ   ⁢           ⁢     J   l   thm       ≡       Δ   ⁢           ⁢       I   l   thm     ⁡     (       ±   μ     ,   ϕ     )           1   -     ⅇ       -     Δτ   l   ↓       /   μ             ;     ⁢     
     ⁢       Δ   ⁢           ⁢     J   l   mss       ≡         Δ   ⁢           ⁢       I   l   mss     ⁡     (       ±   μ     ,   ϕ     )           1   -     ⅇ       -     Δτ   l   ↓       /   μ           .               (   59   )               
These source terms are then multiplied by refractive path transmittance differences and added to the single scatter solar contributions:
 
                   〈       I   0     ⁡     (   Ω   )       〉     =       〈       t     0   -&gt;   s       ⁢       I   s     ⁡     (   Ω   )         〉     -       ∫   0   s     ⁢       〈     t     0   -&gt;       s   ′     -&gt;   sun       mol_abs     〉     ⁢       t     0   -&gt;       s   ′     -&gt;   sun       part_abs     _     ⁢       t       s   ′     -&gt;   sun     sct     _     ⁢       F   sun     _     ⁢         p     s   ′       _     ⁡     (     Ω   ,     Ω   sun       )       ⁢     (       ⅆ       t     0   -&gt;     s   ′       sct     _         ⅆ     s   ′         )     ⁢     ⅆ     s   ′           +       ∑     σ   ⁢           ⁢   segments       ⁢       {         〈     t     0   -&gt;     s     σ   -   1         mol_abs     〉     ⁢       t     0   -&gt;     s     σ   -   1         part_abs     _     ⁢       t     0   -&gt;     s     σ   -   1         sct     _       -       〈     t     0   -&gt;     s   σ       mol_abs     〉     ⁢       t     0   -&gt;     s   σ       part_abs     _     ⁢       t     0   -&gt;     s   σ       sct     _         }     ⁢     (       Δ   ⁢           ⁢     J     l   σ     thm       +     Δ   ⁢           ⁢     J     l   σ     mss         )                   (   60   )               
This procedure enables the band model molecular transmittance COG to be incorporated into the multiple scatter calculation.
 
DISORT I/O for MODTRAN™
 
   The complete list of inputs to the plane-parallel multiple scattering algorithms in MODTRAN™ is quite modest. As noted above, the atmosphere is stratified with each layer l characterized by a nadir optical depth, Δτ l   ↓ , and layer averaged values for the single scattering albedo, ω l , and the scattering phase function Legendre moments, g l   l . In addition, layer boundaries are characterized by altitude level temperatures, T l . Other than this profile data, the only inputs are spectral range (required for the Planck function evaluation), TOA solar irradiance, solar angle and observer zenith and relative solar azimuth angles. 
   In principle, solving for layer source functions, ΔJ l , is simplest in the optical thin limit. However, this limit is occasionally problematic for DISORT. To simplify the source code, the developers of DISORT made a conscience effort to avoid branching for special cases. As a result, inclusion of layers with a very small nadir optical depth, Δτ l   ↓ , can produce numerical instabilities. This problem is avoided by merging MODTRAN™ layers to insure that Δτ l   ↓  exceeds a 0.0001 threshold. This actually provides the added benefit of speeding up the DISORT calculations, because computation time is proportional to the number of computation layers. The DISORT output of vertical flux profiles are redefined onto the finer gridded MODTRAN™ levels via optical depth weighted interpolations. The coarse layer source function values are assigned without modification to the MODTRAN™ sub-layers with one caveat: the thermal source for thin layers (nadir optical depth of less than 0.001) is set to the layer blackbody emission value (the conservative scattering approximation). 
   The calculation of the layer source functions is numerically more stable if the solar and thermal DISORT scatter calculations are performed independently, especially for azimuth dependent scenarios. One reason for this is that when DISORT solar scatter calculations are performed and the nadir absorption optical depth is large, DISORT only calculates the solar scattering contributions for the upper atmosphere; atmospheric layers for which the nadir absorption optical depth exceeds 10 are excluded. However, the dark layers are not excluded from the solar calculations if thermal scatter contributions are computed simultaneously. Combining of solar and thermal calculations results in numerical instabilities. Thus, whenever MODTRAN performs azimuth dependent scattering DISORT calculations, the solar and thermal contributions are computed from separate calls to DISORT. 
   Additionally, special cases can arise that are problematic for DISORT and either require special handling by the MODTRAN™ interface or complete avoidance. Three of these cases involve a solar geometry. If the solar zenith angle is too close to one of the Gaussian double quadrature angles, the IMPROVEMENTS TO MODTRAN4™ package dithers the solar angle to insure that |μ i −μ sun |&gt;μ sun /10,000. If the view zenith angle is near coincident with the solar zenith angle, the view angle is perturbed as follows: 
                           If   ⁢           ⁢            μ   sun     +     μ   o              &lt;   0.002         ⇒           μ   o   ′     =           ⁢     {             -     (       μ   sun     +   0.002     )               for   ⁢           ⁢     μ   sun       ≤   0.5               -     (       μ   sun     -   0.002     )               for   ⁢           ⁢     μ   sun       &gt;   0.5           ,               (   61   )               
where μ 0  is the original view cosine zenith and μ 0 ′ is the view cosine zenith used in the DISORT algorithm. Finally, if the sun is below the horizon (μ sun ≦0), then solar multiple scattering contributions are set to zero. Spectrally, if the single scattering albedo is too close to unity for any layer, it is lowered to a value that will insure a numerically stable solution; this can lead to inaccurate scattered radiances for a purely Rayleigh atmosphere.
 
MODTRAN™ Atmospheric Correction Data File
 
   The observed solar scatter radiance at nadir optical depth τ ↓  measured from an Earth viewing Visible/Near-Infrared (VNIR) through Shortwave-Infrared (SWIR) hyperspectral imager along the Ω 0  direction can be conveniently partitioned into three components: the scattered only by atmosphere radiance I atm (Ω 0 ;τ ↓ ) the directly transmitted from ground reflected radiance I dir   gnd (Ω 0 ;τ ↓ ), and the diffusely transmitted from ground reflected radiance I dif   gnd (Ω 0 ;τ ↓ ): 
                   I   ⁡     (       Ω   o     ;     τ   ↓       )       =         I   atm     ⁡     (       Ω   o     ;     τ   ↓       )       +       I   dir   gnd     ⁡     (       Ω   o     ;     τ   ↓       )       +         I   dir   gnd     ⁡     (       Ω   o     ;     τ   ↓       )       .               (   62   )               
As before, the spectral index has been omitted for notational simplicity. This partitioning is particularly useful for remote sensing applications because it provides a distribution of ground information content: I atm  contains no information about the ground; I dir   gnd  is directly proportional to the imaged pixel reflectance; and I dif   gnd  provides data about surrounding ground pixels.
 
   For a specified atmosphere, the scattered only by atmosphere solar radiance is computed in DISORT by setting the upper illumination at the lower boundary, I m (+μ;τ g   ↓ ), to zero and dropping the thermal contributions from Eq. (46): 
                       I   atm     ⁡     (       Ω   o     ;     τ   ↓       )       =     I   ⁡     (       Ω   o     ;     τ   ↓       )         ⁢     
     ⁢       Source   ⁢           ⁢   Function     _     ⁢     
     ⁢       J   ⁡     (     Ω   ;     τ   ↓       )       ≡       ω     τ   ↓       ⁡     [         p   ⁡     (     Ω   ,       Ω   sun     ;     τ   ↓         )       ⁢     ⅇ       -     τ   ↓       /     μ   sun         ⁢     F   sun       +       ∫     4   ⁢   π       ⁢       p   ⁡     (     Ω   ,       Ω   ′     ;     τ   ↓         )       ⁢     I   ⁡     (       Ω   ′     ;     τ   ↓       )       ⁢     ⅆ     Ω   ′             ]         ⁢     
     ⁢       Boundary   ⁢           ⁢   Conditions     _     ⁢     
     ⁢         I   ⁡     (       +   μ     ,     ϕ   ;     τ   g   ↓         )       =   0     ;     ⁢     
     ⁢       I   ⁡     (       -   μ     ,     ϕ   ;   0       )       =   0             (   63   )               
In practice, the lower boundary condition is specified by setting a thermal flag to false and setting the surface reflectance to zero.
 
   The directly transmitted from ground reflected radiance, I dir   gnd  can be expressed as a function of ground-dependent only and atmosphere-dependent only terms if the imaged pixel is modeled as a Lambertian reflector with BRDF, f(+μ,−μ′), equal to α/π: 
                     I   dir   gnd     ⁡     (       Ω   o     ;     τ   ↓       )       =           μ   sun     ⁢     F   sun     ⁢     t   total     sun   -&gt;   gnd           1   -       〈   α   〉     ⁢     σ   gnd           ⁢     (     α   π     )     ⁢       t   direct     obs   -&gt;   gnd       .               (   64   )               
In this expression, the product μ sun  F sun  equals the top-of-atmosphere (TOA) solar downward horizontal flux. Multiplying this product by the sun-to-ground total (direct+diffuse) transmittance, t total   sun→gnd  (=t direct   sun→gnd +t diffuse   sun→gnd ) results in the expression for the surface downward solar flux given a black (i.e., non-reflecting) and cold (i.e., non-emitting) ground. If the ground is actually a Lambertian reflector with a positive spatially-averaged surface albedo, i.e., &lt;α&gt;&gt;0, then the product &lt;α&gt;σ gnd  equals the fraction of initial downward illumination at the ground that will return to the bottom of the atmosphere after a single ground reflection, where σ gnd  is the atmospheric spherical albedo from the ground. The atmospheric spherical albedo is defined as the fraction of lower boundary isotropic illumination that will scatter back to the ground. The multiple reflections expression is a geometric sum:
 
                   1   +       〈   α   〉     ⁢     σ   gnd       +       (       〈   α   〉     ⁢     σ   gnd       )     2     +       (       〈   α   〉     ⁢     σ   gnd       )     3     +   ⋯     =       1     1   -       〈   α   〉     ⁢     σ   gnd           .             (   65   )               
Thus, the left-most quotient on the right-hand side of Eq. (64) equals the total downward solar flux at the ground. A final multiplication of this surface flux by the Lambertian reflectance and the direct transmittance to the observer, t direct   obs→gnd , provides the expression for the directly transmitted ground reflected radiance.
 
   MODTRAN™ itself computes direct path transmittances, t direct   sun→gnd  and t direct   obs→gnd . The remaining atmosphere-dependent only terms of Eq. (64) are readily determined using DISORT by appealing to the RT reciprocity principle. As noted by Stamnes, the plane albedo from and diffuse transmissivity through a vertically inhomogeneous planetary atmosphere can be expressed in terms of normalized emergent diffuse intensities from isotropic incident illumination: 
                         t   diffuse     sun   -&gt;   gnd       =         I   0     ⁡     (       +     μ   sun       ;   0     )           I   0     ⁡     (       +   μ     ;     τ   g   ↓       )           ;     ⁢     
     ⁢           σ   gnd     ⁡     (     -     μ   ′       )       =         I   0     ⁡     (       -     μ   ′       ;     τ   g   ↓       )           I   0     ⁡     (       +   μ     ;     τ   g   ↓       )           ;     ⁢     
     ⁢       Source   ⁢           ⁢   Function     _     ⁢     
     ⁢       J   0     ⁡     (     μ   ;     τ   ↓       )       =       ω     τ   ↓       ⁢     {       ∑       j   =     -   N         j   ≠   0       N     ⁢       w   j     ⁢         I   0     ⁡     (       μ   j     ;     τ   ↓       )       [       ∑     ℓ   =   0         2   ⁢   M     -     1   ′         ⁢       (       2   ⁢   ℓ     +   1     )     ⁢       g     τ   ↓     ℓ     2     ⁢       Λ   ℓ   0     ⁡     (   μ   )       ⁢       Λ   ℓ   0     ⁡     (     μ   j     )           ]         }         ⁢     
     ⁢       Boundary   ⁢           ⁢   Conditions     _     ⁢     
     ⁢           I   0     ⁡     (       +   μ     ;     τ   g   ↓       )       =     constant   &gt;   0       ;     ⁢     
     ⁢           I   0     ⁡     (       -     μ   i       ;   0     )       =   0     ,             (   66   )               
where the plane albedo, σ gnd (−μ′), is related to the spherical albedo, σ gnd , by the integral equation
 σ gnd =∫ 2π σ gnd (−μ′)μ′ dΩ′.   (67) 
The beauty of these equations is that there is no azimuth dependence, so only the m=0 solution need be computed, and the particular solution does not have to be calculated because the boundary and the external source and the emitted source terms are all zero. Thus, these calculations are computationally modest.
 
   The final term of Eq. (62) is the diffusely transmitted ground reflected radiance: 
                     I   dif   gnd     ⁡     (       Ω   o     ;     τ   ↓       )       =           μ   sun     ⁢     F   sun     ⁢     t   total     sun   -&gt;   gnd           1   -       〈   α   〉     ⁢     σ   gnd           ⁢     (       〈   α   〉     π     )     ⁢       t   emb_dif     obs   -&gt;   gnd       .               (   68   )               
In this equation, the total downward solar flux, the left-most quotient on the right-hand side, is multiplied by the spatially averaged Lambertian reflectance and a medium embedded diffuse transmittance term, t emb     —     dif   obs→gnd . In remote sensing literature, this transmittance term is typically simply referred to as the observer-to-ground diffuse transmittance. However, the observer-to-ground diffuse transmittance is defined as the fraction of parallel beam illumination along the observer view direction that would transmit through a bounded medium. From the principle of reciprocity, this equals illumination-normalized radiation exiting the medium in the reverse observer view direction when isotropically illuminated from below. If the sensor is embedded within the medium, or more specifically, the Earth&#39;s atmosphere, then Eq. (68) needs to include the back-scattering from atmosphere above the observer, such as that illustrated in the trajectory of  FIG. 5 . The term medium embedded diffuse transmittance is being introduced to capture all the ground leaving radiance that reaches the observer after being scattered by atmosphere without further interacting with the ground surface. This quantity is calculated as follows:
 
                       t   emb_dif     obs   -&gt;   gnd       =         I   0     ⁡     (       +     μ     gnd   -&gt;   obs         ;     τ   ↓       )           I   0     ⁡     (       +   μ     ;     τ   g   ↓       )           ;     ⁢     
     ⁢     0   ≤     τ   ↓     ≤     τ   g   ↓       ⁢     
     ⁢       Source   ⁢           ⁢   Function     _     ⁢     
     ⁢         J   0     ⁡     (     μ   ;     τ   ↓       )       =       ω     τ   ↓       ⁢     {       ∑       j   =     -   N         j   ≠   0       N     ⁢       w   j     ⁢         I   0     ⁡     (       μ   j     ;     τ   ↓       )       ⁡     [       ∑     ℓ   =   0         2   ⁢   M     -     1   ′         ⁢       (       2   ⁢   ℓ     +   1     )     ⁢       g     τ   ↓     ℓ     2     ⁢       Λ   ℓ   0     ⁡     (   μ   )       ⁢       Λ   ℓ   0     ⁡     (     μ   j     )           ]           }         ⁢     
     ⁢       Boundary   ⁢           ⁢   Conditions     _     ⁢     
     ⁢           I   0     ⁡     (       +   μ     ;     τ   g   ↓       )       =     constant   &gt;   0       ;     ⁢     
     ⁢         I   0     ⁡     (       -     μ   i       ;   0     )       =   0.             (   69   )               
As this equation illustrates, the radiative transport problem that must be solved to calculate t emb     —     dif   obs→gnd  is identical to the problem solved for the t diffuse   sun→gnd  calculation, Eq. (66), except that the upwelling radiance, I 0 (+μ gnd→obs ;τ ↓ ) is determined within the medium, not at its upper boundary. [The equations also differ in that Eq. (66) is solved for the solar cosine μ sun , while Eq. (69) is solved for the view cosine μ gnd→obs ].
 
   The calculation of spherical albedo, the solar path diffuse transmittance and the view angle medium-embedded diffuse transmittance by DISORT has been streamlined as part of the IMPROVEMENTS TO MODTRAN4™. DISORT provides an option to calculate spherical albedo and diffuse transmittance in a dedicated run of the model. The IMPROVEMENTS TO MODTRAN4™ use DISORT to calculate both the line-of-sight (LOS) scattered radiance contributions and the atmospheric correction data. These two calculations required the same sets of Legendre polynomials and atmospheric layer Eigenvalues and Eigenvectors; they differ in their constants of integration because of differing boundary conditions. DISORT in IMPROVEMENTS TO MODTRAN4™ has been modified to generate both the LOS and atmospheric correction data (including the medium-embedded diffuse transmittance) in a single call with redundant computations eliminated. With a vector optimized compiler, the extra computational effort required to generate the atmospheric correction data during a LOS scattered radiance calculation only amounts to a few percent of the overall time. 
   Surface Moisture 
   Atmospheric water vapor retrievals from down-looking visible and infrared imaging sensor data are adversely impacted by surface droplets and/or embedded liquid water on otherwise dry backgrounds, as with dew or within chlorophyll structures like the leaves of a forested canopy. The water vapor retrieval algorithms typically assume the background and atmospheric water signatures are uncorrelated; this is untrue if the background contains liquid water. This spectral correlation of the surface signature with atmospheric water also presents a problem for viewing against extended water body backgrounds, as in littoral scenes. The details of the radiation transport (RT) in the background underwater scenario is considerably different from embedded water or water droplet scenario. This section describes approaches recently implemented within IMPROVEMENTS TO MODTRAN4™ for dealing with both situations. 
   Embedded Water Attenuation Model 
   The IMPROVEMENTS TO MODTRAN4™ embedded water or moisture attenuation model incorporates the effects of the liquid water extinction directly into its surface spectral reflectance model. For a dry background, MODTRAN™ parameterizes surface reflectance via a spectral grid λ i (i=0, . . . ,N) of bidirectional reflectance distribution functions (BRDFs), f i (Ω v ,Ω s ) and their angularly integrated products. These reflectance products include (a) ρ i (Ω v ), the surface spectral directional reflectivities in the view solid angle direction, Ω v ≡(μ v ,φ v ); (b) ρ i (Ω s ), the surface spectral directional reflectivities in the source solid angle direction, Ω s ≡(μ s ,φ s ); and (c) α i , the double hemisphere integrated surface spectral albedos: 
                         ρ   i     ⁡     (     Ω   v     )       ≡       ∫     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v     ,     Ω   s       )       ⁢     μ   s     ⁢     ⅆ     Ω   s             ;     ⁢     
     ⁢           ρ   i     ⁡     (     Ω   s     )       ≡       ∫     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v     ,     Ω   s       )       ⁢     μ   v     ⁢     ⅆ     Ω   v             ;     ⁢     
     ⁢         α   i     ≡       1   π     ⁢       ∫     2   ⁢   π       ⁢         ρ   i     ⁡     (     Ω   v     )       ⁢     μ   v     ⁢     ⅆ     Ω   v               =       1   π     ⁢       ∫     2   ⁢   π       ⁢         ρ   i     ⁡     (     Ω   s     )       ⁢     μ   s     ⁢       ⅆ     Ω   s       .                     (   70   )               
As indicated above, the solid angles Ω v  and Ω s  have two components: the cosine of a zenith angle μ and an azimuth angle φ.
 
   The embedded water attenuation model uses a single semi-empirical depth parameter, D, to incorporate liquid water attenuation into the surface spectral reflectance terms: 
                         f   i     ⁡     (       Ω   v     ,       Ω   s     ;   D       )       ≡       ⅇ       -     b   i       ⁢     D   /     μ   v           ⁢       f   i     ⁡     (       Ω   v     ,     Ω   s       )       ⁢     ⅇ       -     b   ⅈ       ⁢     D   /     μ   s               ;     ⁢     
     ⁢               ρ   i     ⁡     (       Ω   v     ;   D     )       ≡       ⁢       ∫     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v     ,       Ω   s     ;   D       )       ⁢     μ   s     ⁢     ⅆ     Ω   s                         =       ⁢       ⅇ       -     b   i       ⁢     D   /     μ   v           ⁡     [       ∫     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v     ,     Ω   s       )       ⁢     ⅇ       -     b   ⅈ       ⁢     D   /     μ   s           ⁢     μ   s     ⁢     ⅆ     Ω   s           ]         ;           ⁢     
     ⁢               ρ   i     ⁡     (       Ω   s     ;   D     )       =       ⁢       ∫     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v     ,       Ω   s     ;   D       )       ⁢     μ   v     ⁢     ⅆ     Ω   v                         ≡       ⁢       [       ∫     2   ⁢   π       ⁢       ⅇ       -     b   ⅈ       ⁢     D   /     μ   v           ⁢       f   i     ⁡     (       Ω   v     ,     Ω   s       )       ⁢     μ   v         ]     ⁢     ⅇ       -     b   ⅈ       ⁢     D   /     μ   s               ;           ⁢     
     ⁢               α   i     ⁡     (   D   )       ≡       ⁢       1   π     ⁢       ∫     2   ⁢   π       ⁢         ρ   i     ⁡     (       Ω   v     ;   D     )       ⁢     μ   v     ⁢     ⅆ     Ω   v                         =       ⁢       1   π     ⁢       ∫     2   ⁢   π       ⁢         ρ   i     ⁡     (       Ω   s     ;   D     )       ⁢     μ   s     ⁢       ⅆ     Ω   s       .                           (   71   )               
In these definitions, b i  is the liquid water spectral extinction coefficient (1/m) at spectral wavelength λ i . The extinction coefficients are computed as the sum of absorption b i   abs  and scattering b i   scat  terms,  FIG. 6 . The absorption coefficients b i   abs  are derived from the Segelstein spectral complex refractive index tables. Smith and Baker scattering coefficient b i   scat  data has been fit to a simple power law expression:
 
                     b   λ   scat     =       1.12926   ×       10     -   4       /   m           [     λ   ⁡     (   µm   )       ]     4.32         ,           (   72   )               
In this expression, λ is the spectral wavelength in microns. As evidenced by  FIG. 7 , the scattering contribution to the overall extinction is negligible for wavelengths above ˜1 μm, and extrapolation of the data is inconsequential.
 
   Since MODTRAN™ only needs the BRDF and directional reflectivities values at the solar and (downward) viewing angles, a generally coarse spectral grid of values are computed up front, i.e., before the primary MODTRAN spectral loop. For the dry surface, the reflectance quantities are then linearly interpolated in wavelength, as needed: 
                           f   λ     ⁡     (       Ω   v     ,     Ω   s       )       =         (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁢       f     i   -   1       ⁡     (       Ω   v     ,     Ω   s       )         +       (       λ   -     λ     i   -   1             λ   i     -     λ     i   -   1           )     ⁢       f   i     ⁡     (       Ω   v     ,     Ω   s       )             ;     ⁢     
     [       λ     i   -   1       ≤   λ   ≤     λ   i       ]     ⁢     
     ⁢           ρ   λ     ⁡     (     Ω   v     )       =         (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁢       ρ     i   -   1       ⁡     (     Ω   v     )         +       (       λ   -     λ     i   -   1             λ   i     -     λ     i   -   1           )     ⁢       ρ   i     ⁡     (     Ω   v     )             ;     ⁢     
     ⁢           ρ   λ     ⁡     (     Ω   s     )       =         (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁢       ρ     i   -   1       ⁡     (     Ω   s     )         +       (       λ   -     λ     i   -   1             λ   i     -     λ     i   -   1           )     ⁢       ρ   i     ⁡     (     Ω   s     )             ;     ⁢     
     ⁢       α   λ     =         (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁢     α     i   -   1         +       (       λ   -     λ     i   -   1             λ   i     -     λ     i   -   1           )     ⁢       α   i     .                   (   73   )               
However, the coarse wavelength grids commonly used for the surface description will not capture the spectral structure of the liquid water absorption,  FIG. 7 .
 
   Incorporating the finer spectral dependence of the water absorption is straightforward for the BRDF term: 
                           f   λ     ⁡     (       Ω   v     ,       Ω   s     ;   D       )       ≡       ⁢       ⅇ       -     b   λ       ⁢     D   /     μ   v           ⁢       f   λ     ⁡     (       Ω   v     ,     Ω   s       )       ⁢     ⅇ       -     b   λ       ⁢     D   /     μ   s                         =       ⁢       ⅇ       -     b   λ       ⁢     D   /     μ   v           [         (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁢       f     i   -   1       ⁡     (       Ω   v     ,     Ω   s       )         +                         ⁢       (       λ   -     λ     i   -   1             λ   i     -     λ     i   -   1           )     ⁢       f   i     ⁡     (       Ω   v     ,     Ω   s       )         ]     ⁢       ⅇ       -     b   λ       ⁢     D   /     μ   s           .                   (   74   )               
Modeling of the directional reflectivity and double hemisphere integrated albedo terms is more difficult because the angular integrations, which are preferentially performed in preprocessing, are spectrally dependent on the water absorption. The following approximation is introduced to compute the spectral directional reflectivities:
 
                             ρ   λ     ⁡     (       Ω   v     ;   D     )       ≡       ⁢       ∫     2   ⁢   π       ⁢         f   λ     ⁡     (       Ω   v     ,       Ω   s     ;   D       )       ⁢     μ   s     ⁢     ⅆ     Ω   s                       =       ⁢       ⅇ       -     b   λ       /     μ   v         [         (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁢       f     i   -   1       ⁡     (       Ω   v     ,       〈     Ω   s     〉       i   -   1         )         +                         ⁢       (       λ   -     λ     i   -   1             λ   i     -     λ     i   -   1           )     ⁢       f   i     ⁡     (       Ω   v     ,       〈     Ω   s     〉     i       )         ]     ⁢       ∫     2   ⁢   π       ⁢       ⅇ       -     b   λ       ⁢     D   /     μ   s           ⁢     μ   s     ⁢     ⅆ     Ω   s                       ≈       ⁢     [         (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁢         ρ     i   -   1       ⁡     (       Ω   v     ;   D     )         2   ⁢       E   3     ⁡     (       b     i   -   1       ⁢   D     )       ⁢     ⅇ       -     b     i   -   1         ⁢     D   /     μ   v                 +                           ⁢       (       λ   -     λ     i   -   1             λ   i     -     λ     i   -   1           )     ⁢         ρ   i     ⁡     (       Ω   v     ;   D     )         2   ⁢       E   3     ⁡     (       b   i     ⁢   D     )       ⁢     ⅇ       -     b   i       ⁢     D   /     μ   v                 ]     ⁢   2   ⁢       E   3     ⁡     (       b   λ     ⁢   D     )       ⁢     ⅇ       -     b   λ       ⁢     D   /     μ   v             ,           ⁢     
     ⁢   where           (   75   )                         f   i     ⁡     (       Ω   v     ,       〈     Ω   s     〉     i       )       ≡       ⁢         ∫     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v     ,     Ω   s       )       ⁢     ⅇ       -     b   λ       ⁢     D   /     μ   s           ⁢     μ   s     ⁢     ⅆ     Ω   s               ∫     2   ⁢   π       ⁢       ⅇ       -     b   λ       ⁢     D   /     μ   s           ⁢     μ   s     ⁢     ⅆ     Ω   s                         ≈       ⁢         ∫     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v     ,     Ω   s       )       ⁢     ⅇ       -     b   ⅈ       ⁢     D   /     μ   s           ⁢     μ   s     ⁢     ⅆ     Ω   s               ∫     2   ⁢   π       ⁢       ⅇ       -     b   ⅈ       ⁢     D   /     μ   s           ⁢     μ   s     ⁢     ⅆ     Ω   s                         =       ⁢         ρ   i     ⁡     (       Ω   s     ;   D     )         2   ⁢   π   ⁢           ⁢       E   3     ⁡     (       b   i     ⁢   D     )       ⁢     ⅇ       -     b   ⅈ       ⁢     D   /     μ   s                             (   76   )               
and E 3 (b λ D) is the exponential integral of order 3:
 
                     2   ⁢   π   ⁢           ⁢       E   3     ⁡     (       b   λ     ⁢   D     )         ≡     2   ⁢   π   ⁢       ∫   1   ∞     ⁢         ⅇ       -     b   λ       ⁢   Dt         t   3       ⁢     ⅆ   t             =       ∫     2   ⁢   π       ⁢       ⅇ       -     b   λ       ⁢     D   /   μ         ⁢   μ   ⁢       ⅆ   Ω     .                 (   77   )               
Here the mean value theorem is being invoked to pull f i (Ω v ,Ω s ) out of the integrand at Ω s =&lt;Ω s &gt; i  and f i (Ω v , &lt;Ω s &gt; i ) is being approximated by a transmittance-weighted average at λ=λ i . Thus, even though f i (Ω v , &lt;Ω s &gt; i ) is actually defined as a spectrally varying quantity because of its dependence on the extinction coefficient b λ , it is approximated by a spectrally independent quantity. This weighted average insures that ρ λ (Ω v ; D) indeed equals the precomputed integral ρ i (Ω v ; D) at the spectral grid point, λ=λ i . Since the effective hemispheric transmittance equals twice the exponential integral [note that 2E 3 (b i D)| D=0  equals unity, as it should], a coefficient of 2 is included with each occurrence of the exponential integral in the final line of Eq. (75). With this expression for the directional reflectivity, the only integral that has to be computed on-the-fly spectrally is the exponential integral of order 3, and well established algorithms are readily available for this evaluation.
 
   The double integral mean value theorem and the corresponding double path transmittance weighted average are used to derive an analogous expression for the spectral albedo: 
                             α   λ     ⁡     (   D   )       ≡       ⁢       1   π     ⁢       ∫     2   ⁢   π       ⁢       [       ∫     2   ⁢   π       ⁢         f   λ     ⁡     (       Ω   v     ,       Ω   s     ;   D       )       ⁢     μ   s     ⁢     ⅆ     Ω   s           ]     ⁢     μ   v     ⁢     ⅆ     Ω   v                         =       ⁢       1   π     [         (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁢       f     i   -   1       ⁡     (       〈     (       Ω   v     ,     Ω   s       )     〉       i   -   1       )         +                         ⁢       (       λ   -     λ     i   -   1             λ   i     -     λ     i   -   1           )     ⁢       f   i     ⁡     (       〈     (       Ω   v     ,     Ω   s       )     〉     i     )         ]     ⁢       (       ∫     2   ⁢   π       ⁢       ⅇ       -     b   λ       ⁢     D   /   μ         ⁢   μ   ⁢     ⅆ   Ω         )     2                 ≈       ⁢     [         (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁢         α     i   -   1       ⁡     (   D   )         4   ⁢         E   3     ⁡     (       b     i   -   1       ⁢   D     )       2           +                           ⁢       (       λ   -     λ     i   -   1             λ   i     -     λ     i   -   1           )     ⁢         α   i     ⁡     (   D   )         4   ⁢         E   3     ⁡     (       b   i     ⁢   D     )       2           ]     ⁢   4   ⁢         E   3     ⁡     (       b   λ     ⁢   D     )       2       ,           ⁢     
     ⁢   where           (   78   )                         1   π     ⁢       f   i     ⁡     (       〈     (       Ω   v     ,     Ω   s       )     〉     λ     )         ≡       ⁢               1   π     ⁢       ∫     2   ⁢   π       ⁢     ⅇ       -     b   λ       ⁢     D   /     μ   v                           [       ∫     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v     ,     Ω   s       )       ⁢     ⅇ       -     b   λ       ⁢     D   /     μ   s           ⁢     μ   s     ⁢     ⅆ     Ω   s           ]     ⁢     μ   v     ⁢     ⅆ     Ω   v                   (       ∫     2   ⁢   π       ⁢       ⅇ       -     b   λ       ⁢     D   /   μ         ⁢   μ   ⁢     ⅆ   Ω         )     2                   ≈       ⁢               1   π     ⁢       ∫     2   ⁢   π       ⁢     ⅇ       -     b   λ       ⁢     D   /     μ   v                           [       ∫     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v     ,     Ω   s       )       ⁢     ⅇ       -     b   i       ⁢     D   /     μ   s           ⁢     μ   s     ⁢     ⅆ     Ω   s           ]     ⁢     μ   v     ⁢     ⅆ     Ω   v                   (       ∫     2   ⁢   π       ⁢       ⅇ       -     b   i       ⁢     D   /   μ         ⁢   μ   ⁢     ⅆ   Ω         )     2                   =       ⁢         α   i     ⁡     (   D   )         4   ⁢         E   3     ⁡     (       b   i     ⁢   D     )       2                       (   79   )               
Again, the mean value theorem allows that BRDF to be pulled out of the integrand, but f i  is evaluated at a pair of solid angles, (Ω s ,Ω v )=&lt;(Ω s ,Ω v )&gt; i  for the surface albedo calculation.
 
   In the preprocessing step, the ρ i (Ω;D) of Eq. (75) and α i (D) of Eq. (78) are computed by numerically integrating the BRDF; these integrals are divided by the spectral grid denominator terms [2 E 3 (b i D) exp(−b i D/μ) and 4 E 2 (b i D) 2 , respectively] and stored for use in the fine spectral resolution calculations. No numeric intergrations over the BRDF are performed during this band model spectral loop—only spectral interpolations and calculations of exponentials and exponential intergrals. 
   Flat Water Layer Model 
   The flat water layer model is quite distinct form the embedded layer model in that the liquid water is modeled as a layer that lies above a surface with specified reflectance properties. This introduces three effects in addition to the attenuation embedded water model:
         Fresnel transmittance and reflectance from the air/water interface;   Snell&#39;s Law refraction across the interface; and   Multiple reflections arising primarily from photons that impinge upon the underside of air/water interface at angles greater than the critical angle, i.e., the total internal reflections contribution.
 
Consider first the directly transmitted contribution,  FIG. 7 , part a, to the BRDF at λ i , f i   direct (Ω v ,Ω s ;D;r i ), defined here as source photons that are
   (1) transmitted down through the water-air interface,   (2) redirected according to Snell&#39;s Law refraction,   (3) attenuated by water layer extinction along the straight downward path,   (4) reflected by the underlying surface,   (5) attenuated by water layer extinction along the straight upward path,   (6) transmitted up through the water-air surface, and again   (7) redirected according to Snell&#39;s Law refraction.
 
The water layer BRDF contribution has the following form
       

                             f   i   direct     ⁡     (       Ω   v     ,       Ω   s     ;   D   ;     r   i         )       ≡       ⁢       r   i   2     ⁢       t   i     ⁡     (       θ   v   ′     ,     θ   v       )       ⁢       f   i     ⁡     (       Ω   v   ′     ,       Ω   s   ′     ;   D       )       ⁢       t   i     ⁡     (       θ   s   ′     ,     θ   s       )       ⁢     r   i     -   2                       =       ⁢         t   i     ⁡     (       θ   v   ′     ,     θ   v       )       ⁢       f   i     ⁡     (       Ω   v   ′     ,       Ω   s   ′     ;   D       )       ⁢       t   i     ⁡     (       θ   s   ′     ,     θ   s       )           ,           ⁢     
     ⁢   where           (   80   )                     t   ⁡     (       θ   ′     ,   θ     )       ≡     1   -       1   2     ⁢     (           sin   2     ⁡     (       θ   ′     -   θ     )           sin   2     ⁡     (       θ   ′     +   θ     )         +         tan   2     ⁡     (       θ   ′     -   θ     )           tan   2     ⁡     (       θ   ′     +   θ     )           )           ;     ⁢     
     ⁢       r   i     ≡       n   ⁡     (     λ   i     )       /       n   ′     ⁡     (     λ   i     )           ⁢     
     ⁢         sin   ⁢           ⁢     θ   ′       =       r   i     ⁢   sin   ⁢           ⁢   θ       ;     ⁢     
     ⁢       μ   ′2     =         cos   2     ⁢     θ   ′       =       1   -       r   i   2     ⁡     (     1   -       cos   2     ⁢   θ       )         =       (     1   -     r   i   2       )     +       r   i   2     ⁢     μ   2               ⁢     
     ⁢         μ   ′     ⁢   d   ⁢           ⁢     Ω   ′       =         μ   ′     ⁢   d   ⁢           ⁢     ϕ   ′     ⁢   d   ⁢           ⁢     μ   ′       =         r   i   2     ⁢   μ   ⁢           ⁢   d   ⁢           ⁢   ϕ   ⁢           ⁢   d   ⁢           ⁢   μ     =       r   i   2     ⁢   μ   ⁢           ⁢   d   ⁢           ⁢     Ω   .                     (   81   )               
In these equations, r i (&lt;1) is the air-to-water index of refraction ratio. In the first line of Eq. (80), the r i   2  and r i   −2  terms equal the radiance decrease from the expanded solid angle and increase from the constricted solid angle, respectively, resulting from Snell&#39;s Law refraction; these terms cancel because the radiation both enters and exits the water. The t(θ′,θ) terms are the spectral grid Fresnel interface transmittances. Primes (′) are used to indicate underwater angles.
 
   Approximate approaches are used to model the initial above water Fresnel reflection and the multiple under water reflections. Since water bodies are seldom well represented as having a strictly flat surface, the Fresnel reflection is averaged over the full hemisphere of incident angles, θ s ; this approximation smoothes out the large contributions that occur for the specular reflection geometry. The reflected underwater contributions are modeled by treating the upwelling radiation from below the water-air interface is nearly isotropic. The fraction of photons whose incident angle will exceed the critical angle is (1−r i   2 ). Additional photons are Fresnel reflected, but this contribution is small because the Fresnel transmittance is close to one for incident angles less than the critical angle. Thus, a fraction equal to r i   2  of the upwelling flux penetrates the air-water interface, and (1−r i   2 ) is reflected. The next question that must be answered is “what fraction of the photons reflected from the underside of the water-air interface will return to that interface?” This water layer bottom spherical albedo term is approximated by the double path attenuated water layer surface albedo α i (D) of Eq. (71). Two somewhat compensating assumptions are made: (1) that the backscatter from within the water layer is relatively minor and (2) that the reflected downward flux is nearly isotropic even though the total internal reflections only occur for larger incident angles. Given these approximations, the flat water layer model BRDF is expressed as 
                             f   i     ⁡     (       Ω   v     ,       Ω   s     ;   D   ;     r   i         )       ≈       ⁢           t   i     ⁡     (       θ   s   ′     ,     θ   s       )       ⁢       r   i   2     ⁡     [   ⋯   ]       ⁢       t   i     ⁡     (       θ   v   ′     ,     θ   v       )         +                     ⁢       1     π   2       ⁢       ∫     2   ⁢   π       ⁢       [     1   -       t   i     ⁡     (       θ   s   ′     ,     θ   s       )         ]     ⁢     μ   s     ⁢     ⅆ     Ω   s                           =       ⁢           r   i   2     ⁢       f   i   direct     ⁡     (       Ω   v     ,       Ω   s     ;   D   ;     r   i         )           1   -       (     1   -     r   i   2       )     ⁢       α   i     ⁡     (   D   )             +       1   π     ⁢     〈     1   -       t   i     ⁡     (       θ   s   ′     ,     θ   s       )         〉           ,           ⁢     
     ⁢   where           (   82   )                         [   ⋯   ]     ≡       ⁢     [               f   i     ⁡     (       Ω   v   ′     ,       Ω   s   ′     ;   D       )       +                   ∫       a   ⁢           ⁢   sin   ⁢           ⁢     r   i       &lt;     θ   ′     &lt;     π   /   2         ⁢       f   i     ⁢     (       Ω   v   ′     ,       Ω   ′     ;   D       )     ⁢       f   i     ⁡     (       Ω   ′     ,       Ω   s   ′     ;   D       )       ⁢     μ   ′     ⁢     ⅆ     Ω   ′           +                 ∫       a   ⁢           ⁢   sin   ⁢           ⁢     r   i       &lt;     θ   ′     &lt;     π   /   2         ⁢         f   i     ⁡     (       Ω   v   ′     ,       Ω   ′     ;   D       )       ⁢     (             ∫       a   ⁢           ⁢   sin   ⁢           ⁢     r   i     ⁢     θ   ″       &lt;     π   /   2         ⁢       f   i     ⁡     (       Ω   ′     ,       Ω   s   ″     ;   D       )                       f   i     ⁡     (       Ω   ″     ,       Ω   s   ′     ;   D       )       ⁢     μ   ″     ⁢     ⅆ     Ω   ″               )                       μ   ′     ⁢     ⅆ     Ω   ′         +   ⋯           ]                 ≈       ⁢         f   i     ⁡     (       Ω   v   ′     ,       Ω   s   ′     ;   D       )       ⁡     [     1   +       (     1   -     r   i   2       )     ⁢       α   i     ⁡     (   D   )         +         (     1   -     r   i   2       )     2     ⁢       α   i   2     ⁡     (   D   )         +   ⋯     ]                   =       ⁢         f   i     ⁡     (       Ω   v   ′     ,       Ω   s   ′     ;   D       )         1   -       (     1   -     r   i   2       )     ⁢       α   i     ⁡     (   D   )                     ⁢     
     ⁢   and           (   83   )                 〈     1   -       t   i     ⁡     (       θ   s   ′     ,     θ   s       )         〉     ≡     1   -       1   π     ⁢       ∫       2   ⁢   π     ⁢               ⁢         t   i     ⁡     (       θ   s   ′     ,     θ   s       )       ⁢     μ   s     ⁢     ⅆ     Ω   s                       (   84   )               
The first two terms in Eq. (83) are illustrated pictorially by  FIGS. 7  part a and  7  part b, respectively;  FIG. 7  part c depicts the remaining terms in the series.
 
   To calculate of the flat water layer BRDF at an arbitrary wavelength, the underwater surface BRDF and the water layer attenuated surface albedo are spectrally interpolated, Eqs. (73) and (76), respectively: 
                           f   λ     ⁡     (       Ω   v     ,       Ω   s     ;   D   ;     r   λ         )       ≡       ⁢           r   λ   2     ⁢       f   λ   direct     ⁡     (       Ω   v     ,       Ω   s     ;   D   ;     r   λ         )           1   -       (     1   -     r   λ   2       )     ⁢       α   λ     ⁡     (   D   )             +                     ⁢       1   π     ⁢     (     1   -       1   π     ⁢       ∫     2   ⁢   π       ⁢         t   λ     ⁡     (       θ   s   ′     ,     θ   s       )       ⁢     μ   s     ⁢     ⅆ     Ω   s               )                   =       ⁢           r   λ   2     ⁢       t   λ     ⁡     (       θ   v   ′     ,     θ   v       )       ⁢       f   λ     ⁡     (       Ω   v   ′     ,       Ω   s   ′     ;   D       )       ⁢       t   λ     ⁡     (       θ   s   ′     ,     θ   s       )           1   -       (     1   -     r   λ   2       )     ⁢       α   λ     ⁡     (   D   )             +                     ⁢       1   π     ⁢       〈     1   -       t   λ     ⁡     (       θ   s   ′     ,     θ   s       )         〉     .                     (   85   )               
The angularly averaged Fresnel reflection, &lt;1−t λ (θ′ s , θ s )&gt;, is solely a function of the air-to-water refractive index ratio, r λ . This ratio falls in a range from 1.09 to 1.56 over the entire far ultra-violet through long-wave infrared spectral region from 0.2 to 32 μm; with this limited range, the angularly averaged Fresnel reflection can be accurately fit to a quadratic in r λ  for efficient computation:
 &lt;1 −t   λ (θ′ s ,θ s )&gt;≈−0.3421+0.5033 r   λ −0.1532 r   λ   2   (86) 
   The spectral grid directional reflectance in the view direction, ρ i (Ω v ;D;r i ), and surface albedo, α i (D;r i ), for the flat water layer model are defined by the single and double hemisphere integrated BRDF, respectively: 
                           ρ   i     ⁡     (       Ω   v     ;   D   ;     r   i       )       ≡       ⁢       ∫     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v     ,       Ω   s     ;   D   ;     r   i         )       ⁢     μ   s     ⁢     ⅆ     Ω   s                       =       ⁢             t   i     ⁡     (       θ   v   ′     ,     θ   v       )       ⁢     r   i   2         1   -       (     1   -     r   i   2       )     ⁢       α   i     ⁡     (   D   )             ⁢       ∫   0   1     ⁢       t   i     ⁡     (       θ   s   ′     ,     θ   s       )                           ⁢         ∫   0     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v   ′     ,       Ω   s   ′     ;   D       )       ⁢     ⅆ     ϕ   s       ⁢     μ   s     ⁢     ⅆ     μ   s           +                     ⁢     〈     1   -       t   i     ⁡     (       θ   s   ′     ,     θ   s       )         〉                   ≡       ⁢           ρ   ~     i     ⁡     (       Ω   v     ;   D   ;     r   i       )       +     〈     1   -       t   i     ⁡     (       θ   s   ′     ,     θ   s       )         〉         ;                 (   87   )                         α   i     ⁡     (     D   ;     r   i       )       ≡       ⁢       1   π     ⁢       ∫     2   ⁢   π       ⁢       [       ∫     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v     ,       Ω   s     ;   D   ;     r   i         )       ⁢     μ   s     ⁢     ⅆ     Ω   s           ]     ⁢     μ   v     ⁢     ⅆ     Ω   v                           =       ⁢       1   π     ⁢       ∫     2   ⁢   π       ⁢         ρ   i     ⁡     (       Ω   v     ;   D   ;     r   i       )       ⁢     μ   v     ⁢     ⅆ     Ω   v               ;                 (   88   )               
where {tilde over (ρ)} is defined as the water leaving radiance contribution to ρ. Since the BRDF is defined for the underwater surface, it is easiest to perform numerical integrations using the underwater source and view angles. A change of variables via Eq. (81) leads to the following expressions:
 
                       ρ   i     ⁡     (       Ω   v     ;   D   ;     r   i       )       =       ⁢       1     r   i   2       ⁢       ∫       1   -     r   i   2         1     ⁢       ∫   0     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v     ,       Ω   s     ;   D   ;     r   i         )       ⁢     ⅆ     ϕ   s   ′       ⁢     μ   s   ′     ⁢     ⅆ     μ   s   ′                 ,     
     ⁢               α   i     ⁡     (     D   ;     r   i       )       =       ⁢       1     π   ⁢           ⁢     r   i   4         ⁢       ∫       1   -     r   i   2         1     ⁢     ∫   0     2   ⁢   π                           ⁢       [       ∫       1   -     r   i   2         1     ⁢       ∫   0     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v     ,       Ω   s     ;   D   ;     r   i         )       ⁢     ⅆ     ϕ   s   ′       ⁢     μ   s   ′     ⁢     ⅆ     μ   s   ′             ]     ⁢     ⅆ     ϕ   v   ′       ⁢     μ   v   ′     ⁢     ⅆ     μ   v   ′                     =       ⁢       1     π   ⁢           ⁢     r   i   2         ⁢       ∫       1   ⁢           -           ⁢     r             ⁢   i               ⁢   2           1     ⁢       ∫   0     2   ⁢   π       ⁢         ρ   i     ⁡     (       Ω   v     ;   D   ;     r   i       )       ⁢     ⅆ     ϕ   v   ′       ⁢     μ   v   ′     ⁢       ⅆ     μ   v   ′       .                             (   89   )               
In these expressions, the angular integral domain itself is spectrally dependent.
 
   For the calculation of directional reflectivity at an arbitrary spectral point, ρ λ (Ω′;D;r λ ), the spectrally slowly varying BRDF assumption used in the embedded water attenuation model is replaced with the ansatz that the product of the BRDF with the Fresnel interface transmittance is spectrally slowly varying: 
                             [   ft   ]     λ     ⁢     (       Ω   v   ′     ,     Ω   s   ′       )       ≡       ⁢         f   λ     ⁡     (       Ω   v   ′     ,     Ω   s   ′       )       ⁢       t   λ     ⁡     (       θ   s   ′     ,     θ   s       )                     ≈       ⁢         (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁢       f     i   -   1       ⁡     (       Ω   v   ′     ,     Ω   s   ′       )       ⁢       t     i   -   1       ⁡     (       θ   s   ′     ,     θ   s       )         +                     ⁢       (       λ   -     λ     i   -   1             λ   i     -     λ     i   -   1           )     ⁢       f   i     ⁡     (       Ω   v   ′     ,     Ω   s   ′       )       ⁢       t   i     ⁡     (       θ   s   ′     ,     θ   s       )                     =       ⁢           (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁡     [   ft   ]         i   -   1       ⁢     (       Ω   v   ′     ,     Ω   s   ′       )                       ⁢           (       λ   ⁢           -           ⁢     λ     i   ⁢           -           ⁢   1             λ   i     ⁢           -           ⁢     λ     i   ⁢           -           ⁢   1           )     ⁡     [   ft   ]       i     ⁢       (       Ω   v   ′     ,           ⁢     Ω   s   ′       )     .                     (   90   )               
As illustrated in  FIG. 8 , the liquid water index of refraction only varies between 1.09 and 1.56 for essentially the entire MODTRAN spectral range; however, there are relatively dramatic variations near the 3.0 μm and 6.1 μm liquid water vibrational fundamental bands and in the heart of the pure rotational bands (&gt;10 μm). Within these spectral sub-regions, it is best if the user specifies the sub-surface BRDF on a relative fine spectral grid to protect against inaccurate interpolations of the Fresnel transmittance.
 
   Once again, the mean value theorem is invoked; here it is used to extract the [ft] λ  product out of the spectral directional reflectivity angular integration: 
                             ρ   ~     λ     ⁡     (       Ω   v     ;   D   ;     r   i       )       ≡       ⁢         1     r   λ   2       ⁢       ∫       1   -     r   λ   2           2   ⁢   π       ⁢       ∫   0     2   ⁢   π       ⁢         f   λ     ⁡     (       Ω   v     ,       Ω   s     ;   D   ;     r   λ         )       ⁢     μ   s   ′     ⁢     ⅆ     ϕ   s   ′       ⁢     ⅆ     μ   s   ′               -                     ⁢     〈     1   -       t   λ     ⁡     (       θ   s   ′     ,     θ   s       )         〉                 =       ⁢           t   λ     ⁡     (       θ   v   ′     ,     θ   v       )         1   -       (     1   -     r   λ   2       )     ⁢       α   λ     ⁡     (   D   )             ⁢       ∫       1   -     r   λ   2         1     ⁢       ∫   0     2   ⁢   π       ⁢       f   λ     ⁡     (       Ω   v   ′     ,       Ω   s   ′     ;   D       )                             ⁢         t   λ     ⁡     (       θ   s   ′     ,     θ   s       )       ⁢     μ   s   ′     ⁢     ⅆ     ϕ   s   ′       ⁢     ⅆ     μ   s   ′                     =       ⁢                 t   λ     ⁡     (       θ   v   ′     ,     θ   v       )       ⁢     ⅇ       -     b   λ       ⁢     D   /     μ   v               1   -       (     1   -     r   λ   2       )     ⁢       α   λ     ⁡     (   D   )             ⁡     [   ft   ]       λ     ⁢     (       Ω   v   ′     ,       〈     Ω   s   ′     〉     λ       )                       ⁢       ∫       1   -     r   λ   2         1     ⁢       ∫   0     2   ⁢   π       ⁢       ⅇ       -     b   λ       ⁢     D   /     μ   s   ′           ⁢     μ   s   ′     ⁢     ⅆ     ϕ   s   ′       ⁢       ⅆ     μ   s   ′       .                         (   91   )               
Like the full hemisphere water layer transmittance, Eq. (77), the conical water layer transmittance of Eq. (91) is a function of elliptical integrals, but with added complexity resulting from the restricted domain:
 
                     ∫       1   -     r   λ   2         1     ⁢       ∫   0     2   ⁢   π       ⁢       ⅇ       -     b   λ       ⁢     D   /     μ   s   ′           ⁢     μ   s   ′     ⁢     ⅆ     ϕ   s   ′       ⁢     ⅆ     μ   s   ′             =     2   ⁢     π   ⁡     [         E   3     ⁡     (       b   λ     ⁢   D     )       -       (     1   -     r   λ   2       )     ⁢       E   3     ⁡     (       b   λ     ⁢     D   /       1   -     r   λ   2             )           ]                 (   92   )               
The spectral grid components of [ft] λ  are approximated by the conical water layer transmittance weighted average evaluated at the spectral grid points:
 
                         ⅇ       -     b   λ       ⁢     S   /     μ   v           ⁡     [   ft   ]       i     ⁢     (       Ω   v   ′     ,       〈     Ω   s   ′     〉     λ       )       ≈         ∫       1   -     r   i   2         1     ⁢       [       ∫   0     2   ⁢   π       ⁢         f   i     ⁡     (       Ω   v   ′     ,       Ω   s   ′     ;   D       )       ⁢     ⅆ     ϕ   s   ′           ]     ⁢       t   i     ⁡     (       θ   s   ′     ,     θ   s       )       ⁢     μ   s   ′     ⁢     ⅆ     μ   s   ′             2   ⁢     π   ⁡     [         E   3     ⁡     (       b   i     ⁢   D     )       -       (     1   -     r   i   2       )     ⁢       E   3     ⁡     (       b   i     ⁢     D   /       1   -     r   i   2             )           ]                   (   93   )               
By first substituting Eqs. (90) and (92) into the expression for the spectral directional reflectivity, Eq. (91); then approximating the [ft] i  terms according to Eq. (93); and finally incorporating the spectral grid directional reflectivity expression, Eq. (89), the spectral directional reflectivity in the view (source) direction has the following form:
 
                             ρ   ~     λ     ⁡     (       Ω   v     ;   D   ;     r   λ       )           C   λ     (   ρ   )       ⁡     (       Ω   v     ;   D   ;     r   λ       )         ≈           ⁢         (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁢           ρ   ~       i   -   1       ⁡     (       Ω   v     ;   D   ;     r     i   -   1         )           C     i   -   1       (   ρ   )       ⁡     (       Ω   v     ;   D   ;     r     i   -   1         )           +       (       λ   -     λ     i   -   1             λ   i     -     λ     i   -   1           )     ⁢           ρ   ~     i     ⁡     (       Ω   v     ;   D   ;     r   i       )           C   i     (   ρ   )       ⁡     (       Ω   v     ;   D   ;     r   i       )               ;     ⁢     
     ⁢               ρ   ~     λ     ⁡     (       Ω   s     ;   D   ;     r   λ       )           C   λ     (   ρ   )       ⁡     (       Ω   s     ;   D   ;     r   λ       )         ≈         (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁢           ρ   ~       i   -   1       ⁡     (       Ω   s     ;   D   ;     r     i   -   1         )           C     i   -   1       (   ρ   )       ⁡     (       Ω   s     ;   D   ;     r     i   -   1         )           +       (       λ   -     λ     i   -   1             λ   i     -     λ     i   -   1           )     ⁢           ρ   ~     i     ⁡     (       Ω   s     ;   D   ;     r   i       )           C   i     (   ρ   )       ⁡     (       Ω   s     ;   D   ;     r   i       )               ,             (   94   )               
where the denominator C λ   (ρ)  is a transmittance enhanced by multiple reflections
 
   
     
       
         
           
             
               
                 
                   
                     C 
                     λ 
                     
                       ( 
                       ρ 
                       ) 
                     
                   
                   ⁡ 
                   
                     ( 
                     
                       
                         Ω 
                         
                           v 
                           , 
                           s 
                         
                       
                       ; 
                       D 
                       ; 
                       
                         r 
                         λ 
                       
                     
                     ) 
                   
                 
                 ≡ 
                 
                   
                     
                       
                         
                           t 
                           λ 
                         
                         ⁡ 
                         
                           ( 
                           
                             
                               θ 
                               
                                 v 
                                 , 
                                 s 
                               
                               ′ 
                             
                             , 
                             
                               θ 
                               
                                 v 
                                 , 
                                 s 
                               
                             
                           
                           ) 
                         
                       
                       ⁡ 
                       
                         [ 
                         
                           
                             
                               E 
                               3 
                             
                             ⁡ 
                             
                               ( 
                               
                                 
                                   b 
                                   λ 
                                 
                                 ⁢ 
                                 D 
                               
                               ) 
                             
                           
                           - 
                           
                             
                               ( 
                               
                                 1 
                                 - 
                                 
                                   r 
                                   λ 
                                   2 
                                 
                               
                               ) 
                             
                             ⁢ 
                             
                               
                                 E 
                                 3 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     b 
                                     λ 
                                   
                                   ⁢ 
                                   
                                     D 
                                     / 
                                     
                                       
                                         1 
                                         - 
                                         
                                           r 
                                           λ 
                                           2 
                                         
                                       
                                     
                                   
                                 
                                 ) 
                               
                             
                           
                         
                         ] 
                       
                     
                     
                       1 
                       - 
                       
                         
                           ( 
                           
                             1 
                             - 
                             
                               r 
                               λ 
                               2 
                             
                           
                           ) 
                         
                         ⁢ 
                         
                           
                             α 
                             λ 
                           
                           ⁡ 
                           
                             ( 
                             D 
                             ) 
                           
                         
                       
                     
                   
                   . 
                 
               
             
             
               
                 ( 
                 95 
                 ) 
               
             
           
         
       
     
   
   A completely analogous procedure is used to define the flat water layer spectral double hemisphere albedo. The water leaving contribution to the spectral grid value has the following fully expanded form: 
                             α   ~     i     ⁡     (     D   ;     r   i       )       ≡       ⁢         α   i     ⁡     (     D   ;     r   i       )       +     〈     1   -       t   λ     ⁡     (       θ   s   ′     ,     θ   s       )         〉                   =       ⁢               ∫       1   ⁢           -           ⁢     r   i   2         1     ⁢               ⁢     ⅇ       -     b   l       ⁢           ⁢     D   /     μ   v   ′             ⁢           ⁢     ∫   0     2   ⁢           ⁢   π                       [       ∫       1   ⁢           -           ⁢     r   i   2         1     ⁢               ⁢     ⅇ       -     b   i       ⁢           ⁢     D   /     μ   s   ′             ⁢           ⁢       t   i     ⁡     (       θ   v   ′     ,           ⁢     θ   v       )       ⁢           ⁢       f   i     ⁡     (       Ω   v   ′     ,           ⁢     Ω   v       )       ⁢           ⁢       t   i     ⁡     (       θ   s   ′     ,           ⁢     θ   s       )       ⁢           ⁢     ⅆ     ϕ   s   ′       ⁢           ⁢             ⁢     μ   s   ′       ⁢           ⁢     ⅆ     μ   s   ′           ]     ⁢           ⁢     ⅆ     ϕ   v   ′       ⁢           ⁢             ⁢     μ   v   ′       ⁢           ⁢     ⅆ     μ   v   ′                 π   ⁢           ⁢       r   i   2     ⁡     [     1   -       (     1   -     r   i   2       )     ⁢       α   i     ⁡     (   D   )           ]                         (   96   )               
Spectral interpolations are performed over a triple product, the BRDF with the Fresnel transmittance in both the view and source directions:
 
                             [   tft   ]     λ     ⁢     (       Ω   v   ′     ,     Ω   s   ′       )       ≡       ⁢         t   λ     ⁡     (       θ   v   ′     ,     θ   v       )       ⁢       f   λ     ⁡     (       Ω   v   ′     ,     Ω   s   ′       )       ⁢       t   λ     ⁡     (       θ   s   ′     ,     θ   s       )                     ≈       ⁢         (         λ   i     -   λ         λ   i     -     λ     i   -   1           )     ⁢       t     i   -   1       ⁡     (       θ   v   ′     ,     θ   v       )       ⁢       f     i   -   1       ⁡     (       Ω   v   ′     ,     Ω   s   ′       )       ⁢       t     i   -   1       ⁡     (       θ   s   ′     ,     θ   s       )         +                     ⁢       (       λ   -     λ     i   -   1             λ   i     -     λ     i   -   1           )     ⁢       t   i     ⁡     (       θ   v   ′     ,     θ   v       )       ⁢       f   i     ⁡     (       Ω   v   ′     ,     Ω   s   ′       )       ⁢         t   i     ⁡     (       θ   s   ′     ,     θ   s       )       .                     (   97   )               
The final expression for the flat water layer spectral albedo contains the double path conical transmittance enhanced by multiple reflections:
 
   
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           α 
                           ~ 
                         
                         λ 
                       
                       ⁡ 
                       
                         ( 
                         
                           D 
                           ; 
                           
                             r 
                             λ 
                           
                         
                         ) 
                       
                     
                     
                       
                         C 
                         λ 
                         
                           ( 
                           α 
                           ) 
                         
                       
                       ⁡ 
                       
                         ( 
                         
                           D 
                           ; 
                           
                             r 
                             λ 
                           
                         
                         ) 
                       
                     
                   
                   ≈ 
                   
                     
                       
                         ( 
                         
                           
                             
                               λ 
                               i 
                             
                             - 
                             λ 
                           
                           
                             
                               λ 
                               i 
                             
                             - 
                             
                               λ 
                               
                                 i 
                                 - 
                                 1 
                               
                             
                           
                         
                         ) 
                       
                       ⁢ 
                       
                         
                           
                             
                               α 
                               ~ 
                             
                             
                               i 
                               - 
                               1 
                             
                           
                           ⁡ 
                           
                             ( 
                             
                               D 
                               ; 
                               
                                 r 
                                 
                                   i 
                                   - 
                                   1 
                                 
                               
                             
                             ) 
                           
                         
                         
                           
                             C 
                             
                               i 
                               - 
                               1 
                             
                             
                               ( 
                               α 
                               ) 
                             
                           
                           ⁡ 
                           
                             ( 
                             
                               D 
                               ; 
                               
                                 r 
                                 
                                   i 
                                   - 
                                   1 
                                 
                               
                             
                             ) 
                           
                         
                       
                     
                     + 
                     
                       
                         ( 
                         
                           
                             λ 
                             - 
                             
                               λ 
                               
                                 i 
                                 - 
                                 1 
                               
                             
                           
                           
                             
                               λ 
                               i 
                             
                             - 
                             
                               λ 
                               
                                 i 
                                 - 
                                 1 
                               
                             
                           
                         
                         ) 
                       
                       ⁢ 
                       
                         
                           
                             
                               α 
                               ~ 
                             
                             i 
                           
                           ⁡ 
                           
                             ( 
                             
                               D 
                               ; 
                               
                                 r 
                                 i 
                               
                             
                             ) 
                           
                         
                         
                           
                             C 
                             i 
                             
                               ( 
                               α 
                               ) 
                             
                           
                           ⁡ 
                           
                             ( 
                             
                               D 
                               ; 
                               
                                 r 
                                 i 
                               
                             
                             ) 
                           
                         
                       
                     
                   
                 
                 , 
                 
                   
 
                 
                 ⁢ 
                 where 
               
             
             
               
                 ( 
                 98 
                 ) 
               
             
           
           
             
               
                 
                   
                     C 
                     λ 
                     
                       ( 
                       α 
                       ) 
                     
                   
                   ⁡ 
                   
                     ( 
                     
                       D 
                       ; 
                       
                         r 
                         λ 
                       
                     
                     ) 
                   
                 
                 ≡ 
                 
                   
                     
                       
                         [ 
                         
                           
                             
                               E 
                               3 
                             
                             ⁡ 
                             
                               ( 
                               
                                 
                                   b 
                                   λ 
                                 
                                 ⁢ 
                                 D 
                               
                               ) 
                             
                           
                           - 
                           
                             
                               ( 
                               
                                 1 
                                 - 
                                 
                                   r 
                                   λ 
                                   2 
                                 
                               
                               ) 
                             
                             ⁢ 
                             
                               
                                 E 
                                 3 
                               
                               ⁡ 
                               
                                 ( 
                                 
                                   
                                     b 
                                     λ 
                                   
                                   ⁢ 
                                   
                                     D 
                                     / 
                                     
                                       
                                         1 
                                         - 
                                         
                                           r 
                                           λ 
                                           2 
                                         
                                       
                                     
                                   
                                 
                                 ) 
                               
                             
                           
                         
                         ] 
                       
                       2 
                     
                     
                       
                         r 
                         λ 
                         2 
                       
                       ⁡ 
                       
                         [ 
                         
                           1 
                           - 
                           
                             
                               ( 
                               
                                 1 
                                 - 
                                 
                                   r 
                                   λ 
                                   2 
                                 
                               
                               ) 
                             
                             ⁢ 
                             
                               
                                 α 
                                 k 
                               
                               ⁡ 
                               
                                 ( 
                                 D 
                                 ) 
                               
                             
                           
                         
                         ] 
                       
                     
                   
                   . 
                 
               
             
             
               
                 ( 
                 99 
                 ) 
               
             
           
         
       
     
   
   During MODTRAN™ pre-processing for the flat water layer model, the BRDF is numerically integrated to compute the spectral grid point directional reflectivities ρ i (Ω;D;r i ) and albedos α i (D;r i ). These values are stored after being divided by transmittances terms C λ   (ρ)  and C λ   (α) , respectively. No numerical integration of the BRDF is required at the finer spectral resolution of the MODTRAN™ band model. The flat water layer model calculations performed at this finer resolution are quick and straightforward; they include calculation of the liquid water and air refractive indices, calculation of the Fresnel transmittance, calculation of the elliptical integrals and calculation of the spectral interpolations. 
   Special Cases 
   It is of interest to examine both the embedded water attenuation and the flat water layer models for the scenario in which the underlying BRDF is Lambertian:
 
 f   i   Lamb (Ω v ,Ω s )≡α i /π.  (100)
 
For the embedded water attenuation model, substitution into Eq. (71) gives the following expressions:
 
                         f   i   Lamb     ⁡     (       Ω   v     ,       Ω   s     ;   D       )       ≡         α   i     π     ⁢     ⅇ       -     b   i       ⁢     D   /     μ   s               ;     ⁢     
     ⁢           ρ   i   Lamb     ⁡     (       Ω   v     ;   D     )       ≡     2   ⁢     α   i     ⁢       E   3     ⁡     (       b   i     ⁢   D     )       ⁢     ⅇ       -     b   i       ⁢     D   /     μ   v               ;     ⁢     
     ⁢           ρ   i   Lamb     ⁡     (       Ω   s     ;   D     )       ≡     2   ⁢     α   i     ⁢       E   3     ⁡     (       b   i     ⁢   D     )       ⁢     ⅇ       -     b   i       ⁢     D   /     μ   s               ;     ⁢     
     ⁢         α   i   Lamb     ⁡     (   D   )       ≡     4   ⁢     α   i     ⁢           E   3     ⁡     (       b   i     ⁢   D     )       2     .                 (   101   )               
Not surprisingly, as the depth D approaches zero, the reflectances each approach the Lambertian limit.
 
   The flat water layer model is more complex because of the Fresnel interface and the multiple underwater reflections. Combining the results of Eq. (101), the definition from Eq. (90), and the BRDF expression of Eq. (82) gives 
                     f   i   Lamb     ⁡     (       Ω   v     ,       Ω   s     ;   D   ;     r   i         )       ≡         α   i     π     ⁢       (         r   i   2     ⁢       t   i     ⁡     (       θ   s   ′     ,     θ   s       )       ⁢       t   i     ⁡     (       θ   v   ′     ,     θ   v       )       ⁢     ⅇ       -     b   i       ⁢     D   /     μ   v   ′           ⁢     ⅇ       -     b   i       ⁢     D   /     μ   s   ′               1   -     4   ⁢       α   i     ⁡     (     1   -     r   i   2       )       ⁢         E   3     ⁡     (       b   i     ⁢   D     )       2           )     .               (   102   )               
To obtain approximate analytic expressions for the integrated reflectance terms, the additional simplification is made that the Fresnel transmittance within the angular integrands is unity for zenith angles below the critical angle, sin θ′&lt;r i :
 
                   t   ⁡     (       θ   ′     ,   θ     )       ≈     {               1   ⁢           ⁢   for   ⁢           ⁢   sin   ⁢           ⁢     θ   ′       &lt;     r   i                   0   ⁢           ⁢   for   ⁢           ⁢   sin   ⁢           ⁢     θ   ′       &gt;     r   i             .               (   103   )               
For the Lambertian limit, Eq. (89) becomes
 
                       ρ   i   Lamb     ⁡     (       Ω   v     ;   D   ;     r   i       )       ≈       2   ⁢     α   i     ⁢     ⅇ       -     b   i       ⁢     D   /     μ   v   ′           ⁢         t   i     ⁡     (       θ   v   ′     ,     θ   v       )       ⁡     [         E   3     ⁡     (       b   i     ⁢   D     )       -       (     1   -     r   i   2       )     ⁢       E   3     ⁡     (       b   i     ⁢     D   /       1   -     r   i   2             )           ]           1   -     4   ⁢       α   i     ⁡     (     1   -     r   i   2       )       ⁢         E   3     ⁡     (       b   i     ⁢   D     )       2             ,     
     ⁢         α   i   Lamb     ⁡     (     D   ;     r   i       )       ≈         4   ⁢         α   i     ⁡     [         E   3     ⁡     (       b   i     ⁢   D     )       -       (     1   -     r   i   2       )     ⁢       E   3     ⁡     (       b   i     ⁢     D   /       1   -     r   i   2             )           ]       2           r   i   2     ⁡     [     1   -     4   ⁢       α   i     ⁡     (     1   -     r   i   2       )       ⁢         E   3     ⁡     (       b   i     ⁢   D     )       2         ]         .               (   104   )               
As the depth of the water goes to zero, these reflectance simplify even further:
 
                       f   i   Lamb     ⁡     (       Ω   v     ,       Ω   s     ;   0   ;     r   i         )       ≡         α   i     π     ⁢     (         r   i   2     ⁢       t   i     ⁡     (       θ   s   ′     ,     θ   s       )       ⁢       t   i     ⁡     (       θ   v   ′     ,     θ   v       )           1   -       α   i     ⁡     (     1   -     r   i   2       )           )         ,     
     ⁢         ρ   i   Lamb     ⁡     (       Ω   v     ;   0   ;     r   i       )       ≈         α   i     ⁢       t   i     ⁡     (       θ   v   ′     ,     θ   v       )       ⁢     r   i   2         1   -       α   i     ⁡     (     1   -     r   i   2       )             ,     
     ⁢         α   i   Lamb     ⁡     (     0   ;     r   i       )       ≈           α   i     ⁢     r   i   2         1   -       α   i     ⁡     (     1   -     r   i   2       )           .               (   105   )               
If one furthers considers the limit of a perfect Lambertian reflector, α i =1, and the Fresnel transmittances are set to one, the Lambertian limit is obtained. The Lambertian limit is also obtained for a black (α i =0) surface. However, Eq. (105) does not reduce to the Lambertian limit for grey bodies. In particular, if the reflectance is 50% (α i =½), then the albedo equals r i   2 /(1+r i   2 ). For water, this ratio is closer to one-third than to one-half.
 
   
     
       
         
             
           
             
                 
             
             
               List of Symbols 
             
             
                 
             
           
          
             
                 
             
          
         
         
             
             
          
             
               B v (T) 
               Planck Blackbody Spectral Emission at Temperature T [W cm −2  sr −1 /cm −1 ]. 
             
             
               b v (s) 
               Spectral Extinction Coefficient at Path Length s [cm −1 ]. 
             
             
               b v,i   
               Spectral Extinction Coefficient for Species i [cm −1 ]. 
             
             
               b v   abs (s) 
               Spectral Absorption Coefficient at Path Length s [cm −1 ]. 
             
             
               b v,i   abs (s) 
               Spectral Absorption Coefficient at Path Length s for Species i [cm −1 ]. 
             
             
               b v   sct (s) 
               Spectral Scattering Coefficient at Path Length s [cm −1 ]. 
             
             
               b v,i   sct (s) 
               Spectral Scattering Coefficient at Path Length s for Species i [cm −1 ]. 
             
             
               C v   CO     2     
               Spectral CO 2  Continuum Parameter [(cm 2 /mol)/cm −1 ]. 
             
             
               C v   H     2     O foreign   
               Spectral Foreign-Broadened H 2 O Continuum Parameter [(cm 2 /mol)/cm −1 ]. 
             
             
               C v   H     2     O self   
               Spectral Self-Broadened H 2 O Continuum Parameter [(cm 2 /mol)/cm −1 ]. 
             
             
               C jp   m   
               Plane-Parallel RTE Homogeneous Problem Constants of Integration. 
             
             
               c 
               Speed of Light in a vacuum (2.9979246E+10 cm/sec). 
             
             
               D(r c ) 
               Offset of Spectral Bin Center [cm −1 ]. 
             
             
               F v   sun   
               Top-Of-Atmosphere (TOA) Spectral Solar Irradiance [W cm −2 /cm −1 ]. 
             
             
               F v   ↓ (s) 
               Downward Spectral Diffuse Flux at Path Length s [W cm −2 /cm −1 ]. 
             
             
               F + (τ ↓ ) 
               Upward Spectral Total Flux at Vertical Optical Depth τ ↓ [W cm −2 /cm −1 ]. 
             
             
               F − (τ ↓ ) 
               Downward Spectral Total Flux at Vertical Optical Depth τ ↓ [W cm −2 /cm −1 ]. 
             
             
               ƒ 
               Delta-M Method Separation Parameter. 
             
             
               ƒ v (γ c ,γ d ) 
               Voigt Spectral Line-Shape Function [cm]. 
             
             
               f v (Ω s ,Ω′ s ;s) 
               Ground BRDF for Outgoing Ω s  and Incoming Ω′ s  Angles at Path Length s [sr −1 ]. 
             
             
               f m (μ;−μ′) 
               m th  Azimuth Moment of the Ground BRDF [sr −1 ]. 
             
             
               G jp   m   
               Plane-Parallel RTE Homogeneous Problem Eigenvectors [W cm −2 sr −1 /cm −1 ]. 
             
             
               g   τ     ↓   λ   
               Scattering Phase Function λ th  Legendre Moment at Vertical Optical Depth τ ↓ . 
             
             
               g′   τ     ↓   λ   
               Delta-M Scattering Phase Function λ th  Legendre Moment at Vertical Optical 
             
             
                 
               Depth τ ↓ . 
             
             
               g p   λ   
               Scattering Phase Function λ th  Legendre Moment Averaged over Layer p. 
             
             
               H s′   
               Above Sea Level Altitude at Path Length s′ [km]. 
             
             
               H (p)   
               Above Sea Level Altitude at Profile Level p [km]. 
             
             
               h s′   p   
               Altitude Interpolation Fraction at Path Length s′ for Atmospheric Layer p. 
             
             
               I n   
               Modified Bessel Function of order n. 
             
             
               I v (Ω s ;s) 
               Spectral Radiant Intensity at Path Length s Along the Ω s  Direction [W cm −2 sr −1 /cm −1 ]. 
             
             
               I th (μ;τ ↓)   
               m th  Azimuth Moment of the Plane-Parallel Atmosphere Spectral Radiant 
             
             
                 
               Intensities in Direction μ at Vertical Optical Depth τ ↓  [W cm −2 sr −1 /cm −1 ]. 
             
             
               I p   m (μ i ;τ ↓ ) 
               Plane-Parallel Spectral Radiant Intensities in Layer p [W cm −2 sr −1 /cm −1 ]. 
             
             
               i 
               Species (Molecular or Particulate) Index. 
             
             
               i 
               Gaussian Moment Index. 
             
             
               j 
               Spectral Sub-Bin Index. 
             
             
               j 
               Gaussian Moment Index of Incoming Scattered or Reflected Radiation. 
             
             
               J v (Ω s ;s) 
               Spectral Source Function at Path Length s Along the Ω s  Direction 
             
             
                 
               [W cm −2 sr −1 /cm −1 ]. 
             
             
               k 
               Boltzmann constant (1.3807E−16 erg/K). 
             
             
               k 
               Molecular Index. 
             
             
               k jp   m   
               Eigenvalues of the Plane-Parallel RTE Homogeneous Problem. 
             
             
               L 
               Total Number of Atmospheric Layers. 
             
             
               l 
               Molecular Transition Index. 
             
             
               λ 
               Legendre Moment Index. 
             
             
               M 
               Spectral Band Index. 
             
             
               
                 M 
               
               Effective Band Value Index. 
             
             
               M 
               One More Than the Total Number of Azimuth Moments. 
             
             
               m 
               Molecular Mass [g/molecule]. 
             
             
               m 
               Azimuth Moment Index. 
             
             
               N 
               Number of Spectral Sub-Interval in a Spectral Band. 
             
             
               N 
               Number of DISORT Streams in each Hemisphere. 
             
             
               n v (h) 
               Index of Refraction at Spectral Frequency v and Altitude h. 
             
             
               n 
               Line Number Band Model Parameter. 
             
             
               
                 n 
               
               Curtis-Godson Path Averaged Line Number Band Model Parameter. 
             
             
               P 
               Pressure [Atm]. 
             
             
               P 0   
               Standard Pressure [1 Atm]. 
             
             
               
                 P Δs   
               
               Path Segment Effective Pressure [Atm]. 
             
             
               
                 P s   k   
               
               Effective Pressure of Solar Illumination Path to s for Molecule k [Atm]. 
             
             
               P n   
               Integrals Used in the Definition of the V n  Fourier Coefficients. 
             
             
               P λ (μ) 
               λ th  Legendre Polynomial at Cosine Zenith μ. 
             
             
               p 
               Molecular Transition Index. 
             
             
               p 
               LOS Segment Index. 
             
             
               p 
               Layer or Level Index for Stratified Atmosphere. 
             
             
               p(Θ;τ ↓ ) 
               Spectral Scattering Phase Function for Scattering Angle Θ at Nadir Optical Depth 
             
             
                 
               τ ↓ . 
             
             
               p′(Θ;τ ↓ ) 
               Delta-M Spectral Scattering Phase Function for Scattering Angle Θ at Nadir 
             
             
                 
               Optical Depth τ ↓ . 
             
             
               p v (Ω s ,Ω′ s ;s) 
               Spectral Scattering Phase Function for Outgoing and Incoming Angles Ω s  and Ω′ s   
             
             
                 
               at Path Length s [sr −1 ]. 
             
             
               p v,i (Ω s ,Ω′ s ;s) 
               Spectral Scattering Phase Function of Species i for Outgoing and Incoming 
             
             
                 
               Angles Ω s  and Ω ′   s  at Path Length s [sr −1 ]. 
             
             
               R e   
               Radius of the Earth [km]. 
             
             
               r c   
               Ratio of Collision (Lorentz) Half-Width to Spectral Bin Width. 
             
             
               S 
               Molecular Line Strength [cm −2 /atm] 
             
             
               s 
               Full Path Length [cm]. 
             
             
               s′ 
               Partial Path Length [cm]. 
             
             
               T 
               Temperature [K]. 
             
             
               T 0   
               Standard Temperature [273.15 K]. 
             
             
               
                 T 
                 Δs 
               
               Path Segment Effective Temperature [K]. 
             
             
               
                 T 
                 s 
                 k 
               
               Effective Temperature of Solar Illumination Path to s for Molecule k [K]. 
             
             
               t v (s) 
               Spectral Extinction Transmittance to Path Length s. 
             
             
               t k   cen   
               Molecular Transmittance from Line Center Absorption. 
             
             
               t k   tail   
               Molecular Transmittance from Line Tail Absorption. 
             
             
               t k   coat   
               Molecular Transmittance from Line Continuum Absorption. 
             
             
               t v   abs (s) 
               Spectral Absorption Transmittance to Path Length s. 
             
             
               t v   sct (s) 
               Spectral Scattering Transmittance to Path Length s. 
             
             
               U 
               Path Column Density Corresponding to a Transmittance of 1/e [Atm cm]. 
             
             
               u 
               Path Column Density [Atm cm]. 
             
             
               V n   
               Fourier Expansion Coefficients Used in the Modified Bessel Function Expansion 
             
             
                 
               of the Voigt Line Tail Equivalent Width. 
             
             
               W sl   
               Finite-Bin Equivalent Width (EW) of a Single Voigt Line [cm −1 ] 
             
             
               W Δ   
               Voigt Line Tail Equivalent Width [cm −1 ]. 
             
             
               w i   
               Discrete Ordinate Gaussian Weight. 
             
             
               x n   
               Denominator n th -Order Spectral Coefficient in Padé Approximant Fit to Line Tail 
             
             
                 
               Absorption Cross-Section [cm 2 ] 
             
             
               Y 0p (μ i ) 
               A Particular Solution to the Constant Thermal Source Plane-Parallel RTE 
             
             
                 
               Problem for Layer p [W cm −2 sr −1 /cm −1 ]. 
             
             
               Y 1p (μ i ) 
               A Particular Solution to the Linear in Nadir Optical Depth Thermal Source Plane- 
             
             
                 
               Parallel RTE Problem for Layer p [W cm −2 sr −1 /cm −1 ]. 
             
             
               Z p   m (μ i ) 
               A Particular Solution to the Solar Source Plane-Parallel RTE Problem for 
             
             
                 
               Azimuth Moment m and Layer p [W cm −2 sr −1 /cm −1 ]. 
             
             
               α v (s) 
               Ground Surface Spectral Albedo intersected a distance s along the LOS. 
             
             
               δ m0   
               Delta function (Equals 1 for m equal to 0; otherwise, equals 0). 
             
             
               δ v   
               Normalized Spectral Bin Domain Variable (−1 to +1). 
             
             
               δ(x) 
               Dirac delta function. 
             
             
               γ c   
               Collision (Lorentz) half-width at half maximum [cm −1 ]. 
             
             
               
                 γ 
                 c 
               
               Curtis-Godson Path Averaged Collision (Lorentz) half-width [cm −1 ]. 
             
             
               γ c   0   
               Standard Pressure and Temperature Collision (Lorentz) half-width at half 
             
             
                 
               maximum [cm −1 ]. 
             
             
               γ d   
               Doppler half-width at 1/e of the maximum [cm −1 ]. 
             
             
               
                 γ 
                 d 
               
               Curtis-Godson Path Averaged Doppler half-width [cm −1 ]. 
             
             
               ε v (Ω s ) 
               Ground Surface Spectral Directional Emissivity in Direction Ω S . 
             
             
               φ 
               Path Azimuth Angle. 
             
             
               κ v,i (s) 
               Species i Spectral Extinction Cross-Section at Path Length s [cm 2 ]. 
             
             
               κ v,i   abs (s) 
               Species i Spectral Absorption Cross-Section at Path Length s [cm 2 ]. 
             
             
               κ v,i   sct (s) 
               Species i Spectral Scattering Cross-Section at Path Length s [cm 2 ]. 
             
             
               κ n   
               Numerator n th -Order Spectral Coefficient in Padé Approximant Fit to Line Tail 
             
             
                 
               Absorption Cross-Section [cm 2 ] 
             
             
               Λ λ   m (μ) 
               Unit Normalized Associated Legendre Polynomial at Cosine Zenith μ. 
             
             
               μ′ 
               Cosine of the Path Zenith Angle at the Ground. 
             
             
               μ i   
               Discrete Ordinate Gaussian Abscissa. 
             
             
               v 
               Spectral Frequency [cm −1 ]. 
             
             
               v m   
               Spectral Bin Center Frequency [cm −1 ]. 
             
             
               
                 v 
               
               Calculation Spectral Domain Center Frequency [cm −1 ]. 
             
             
               ρ i (s) 
               Density of Species i at Path Length s [number/cm 3 ]. 
             
             
               ρ v (Ω s ) 
               Ground Surface Spectral Directional Reflectivity in Direction Ω s.   
             
             
               τ v (s) 
               Spectral Extinction Optical Depth to Path Length s. 
             
             
               τ v   abs (s) 
               Spectral Absorption Optical Depth to Path Length s. 
             
             
               τ v   sct (s) 
               Spectral Scattering Optical Depth to Path Length s. 
             
             
               τ v   vert (s) 
               Vertical Spectral Extinction Optical Depth to Path Length s. 
             
             
               τ ↓   
               Nadir Optical Depth from TOA. 
             
             
               τ p   ↓   
               Nadir Optical Depth from TOA to Level p. 
             
             
               X v   
               CO 2  line shape sub-Lorentzian multiplicative correction factor. 
             
             
               Ω s   
               Solid Angle at Path Length s [sr]. 
             
             
               ω   τ     ↓   
               Single Scattering Albedo at Nadir Optical Depth τ ↓ . 
             
             
               ω p   
               Single Scattering Albedo in Layer p. 
             
             
               ω′ p   
               Delta-M Single Scattering Albedo in Layer p. 
             
             
                 
             
          
         
       
     
   
   Although specific features of the invention are shown in some drawings and not others, this is for convenience only as some feature may be combined with any or all of the other features in accordance with the invention.