
  subroutine cross_section_Rayleigh(   &
    nsmax,       wavelength_ns,  relN2,    &
    relO2,       relAr,          eps_N2,   &      
    eps_O2,      eps_Ar,         csray_N2, &
    csray_O2,    csray_Ar)

!-----------------------------------------------------------
! Subroutine:
!    csray
!  
! Purpose:
!    Calculate Rayleigh cross section for 
!    N2, O2, Ar, and air
!------------------------------------------------------------------

  implicit none

  integer,parameter::dp=selected_real_kind(8) !double precision
 
! nr of wavelengths
  integer,intent(in):: &
     nsmax
! wavelengths around incident light (ns_elas) [nm]
  real(kind=dp):: &
    wavelength_ns(0:nsmax)    

  real(kind=dp),intent(in):: & 
     relN2, & !N2 volume mixing ratio
     relO2, & !O2 volume mixing ratio
     relAr    !Ar volume mixing ratio

  real(kind=dp),dimension(0:nsmax):: &
    rn_N2,      & !refraction index of N2
    rn_O2,      & !refraction index of O2
    rn_Ar         !refraction index of Ar

  real(kind=dp):: &
    Wcm,        & !wavelengthlength [cm]
    WNmu,       & !1/wavelengthlength [um-1]
    frac_rn_N2, &
    frac_rn_O2, &
    frac_rn_Ar, &
    rn_N2m1,    & 
    rn_O2m1,    &
    rn_Arm1,    &
    fk_N2,      & !King factor
    fk_O2,      &
    fk_Ar,      &
    dens

! parametrization of refraction index of N2 and O2 ( Bates 1984 )
  real(kind=dp):: &
     a_N2(3),   &
     b_N2(3),   & 
     a_O2(4),   &
     b_O2(4)
  data &
     a_N2 / 6855.200d0,  5989.242d0,  6998.749d0             /, &
     b_N2 / 3243157.0d0, 3363266.3d0, 3233582.0d0            /, &
     a_O2 / 21351.1d0,   20564.8d0,   22120.4d0,  23796.7d0  /, &
     b_O2 / 218567.0d0,  248089.9d0,  203187.6d0, 168988.4d0 / 
 
  real(kind=dp),parameter::    & 
     pi       = 3.141592653589793115997963468544d0, &
     temp_std = 273.15d0,    & !standard temperature [K]
     Dair     = 2.54743D+19    !molar nrdensity air [part. cm-3]

  integer::ni

  real(kind=dp),dimension(0:nsmax):: &
     csray_air    !Rayleigh cross sections for air

  real(kind=dp),dimension(0:nsmax),intent(out):: &
     eps_N2,   &  !epsilon for N2 
     eps_O2,   &  !epsilon for O2
     eps_AR,   &  !epsilon for Ar = 0
     csray_O2, &  !Rayleigh cross sections for O2 molecules
     csray_N2, &  !Rayleigh cross sections for N2 molecules
     csray_Ar     !Rayleigh cross sections for Ar atoms

  integer::ns

  do ns=0,nsmax

     WNmu = (1.d+07/wavelength_ns(ns)) * 1.d-4

!    King correction factor
     fk_N2 = 1.034d0 + 3.17d-4 * WNmu**2
     fk_O2 = 1.096d0 + 1.385d-03*WNmu**2 + 1.448d-4*WNmu**4
     fk_Ar = 1.d0

!    anisotropic polarizability factor eps
     eps_N2(ns) = 9.d0/2.d0 * (fk_N2-1.d0)
     eps_O2(ns) = 9.d0/2.d0 * (fk_O2-1.d0)
     eps_Ar(ns) = 0.d0

!    refraction index of N2  (n(N2)-1) [Bates1984]
     ni = 1 
     if(wavelength_ns(ns).lT.468.d0) ni = 2 
     if(wavelength_ns(ns).lT.254.d0) ni = 3 
     rn_N2m1  = (a_N2(ni) + b_N2(ni)/(144.d0-WNmu**2))*1.d-08
     rn_N2(ns) = 1.d0 + rn_N2m1
         
!    refraction index of O2  (n(O2)-1) [Bates1984]
     ni = 1
     if(wavelength_ns(ns).lT.546.d0)  ni = 2 
     if(wavelength_ns(ns).lT.288.d0)  ni = 3 
     if(wavelength_ns(ns).lT.221.d0)  ni = 4  
     rn_O2m1  = (a_O2(ni) + b_O2(ni)/(40.90d0-WNmu**2))*1.d-08
     rn_O2(ns) = 1.d0 + rn_O2m1
     
!    refraction index of Ar (n(Ar)-1)  [Peck and Fisher1964]
     rn_Arm1  = (678.6711d0 + 301829.43d0/(144.d0-WNmu**2))*1.d-07
     rn_Ar(ns) = 1.d0 + rn_Arm1

!    (n**2-1)/(n**2+2)
     frac_rn_N2 = (rn_N2(ns)**2-1.d0)**2/(rn_N2(ns)**2+2.d0)**2
     frac_rn_O2 = (rn_O2(ns)**2-1.d0)**2/(rn_O2(ns)**2+2.d0)**2
     frac_rn_Ar = (rn_Ar(ns)**2-1.d0)**2/(rn_Ar(ns)**2+2.d0)**2
      
!    Dair is given for 15 degree Celius
!    scale to standard temp T = 273.15 K because of refractive index

     dens = (288.16d0/temp_std) * Dair

     Wcm = wavelength_ns(ns) * 1.d-7 !change units nm -> cm

     csray_N2(ns) = 24.d0*pi**3/(Wcm**4*dens**2)*frac_rn_N2*fk_N2
     csray_O2(ns) = 24.d0*pi**3/(Wcm**4*dens**2)*frac_rn_O2*fk_O2
     csray_Ar(ns) = 24.d0*pi**3/(Wcm**4*dens**2)*frac_rn_Ar*fk_Ar
         
     csray_air(ns) = relN2 * csray_N2(ns) + &
                     relO2 * csray_O2(ns) + &
                     relAr * csray_Ar(ns)     

  enddo

  return

end subroutine cross_section_Rayleigh
