*********************************************************************
*                                                                   *
*     Atmospheric Absorption for Rayleigh, Mie and Ozone.           * 
*                                                                   *
*     Created: May-98                                               *
*     Author: Aitor Ibarra Ibaibarriaga                             *
*             Jose Carlos Gonzalez                                  *
*                                                                   *
*     Modified December-2002, January-2003                          *
*     Author: Abelardo Moralejo (moralejo@pd.infn.it)               *
*                                                                   *
*     Now this accounts only for Rayleigh scattering. Mie and       *
*     Ozone absorption are now treated separatedly.                 *
*                                                                   *
*********************************************************************
*
* December 2002, Abelardo Moralejo: checked algorithms, removed 
* old/unnecessary code, fixed some bugs, added comments.
*
* Fixed bugs (of small influence) in Mie absorption implementation: there were
* errors in the optical depths table, as well as a confusion: height a.s.l.
* was used as if it was height above the telescope level. The latter error was
* also present in the Ozone aborption implementation.
*
* On the other hand, now we have the tables AE_ABI and OZ_ABI with optical 
* depths between sea level and a height h (before it was between 2km a.s.l 
* and a height h). So that we can simulate also in the future a different 
* observation level. 
*
* AM: WARNING: IF VERY LARGE zenith angle simulations are to be done (say 
* above 85 degrees, for neutrino primaries or any other purpose) this code 
* will have to be adapted accordingly and checked, since in principle it has 
* been written and tested to simulate the absorption of Cherenkov photons 
* arriving at the telescope from above.
*
* AM: WARNING 2: not to be used for wavelengths outside the range 250-700 nm
*
* January 2003, Abelardo Moralejo: found error in Ozone absorption treatment.
* At large zenith angles, the air mass used was the one calculated for 
* Rayleigh scattering, but since the Ozone distribution is rather different
* from the global distribution of air molecules, this is not a good 
* approximation. Now I have left in this code only the Rayleigh scattering,
* and moved to atm.c the Mie scattering and Ozone absorption calculated in
* a better way.
*
* AM, Jan 2003: added obslev as an argument of the function. Changed the 
* meaning of the argument height: now it is the true height above sea level 
* at which a photon has been emitted, before it was the height given by
* Corsika, measured in the vertical of the observer and not in the vertical 
* of the emitting particle.
*
* @endtext
c @code

      SUBROUTINE attenu(wavelength, height, obslev, theta, tr_atmos)

      REAL wavelength, height, obslev, theta, tr_atmos
      COMMON /ATMOS/   BATM,CATM,LAHG,LONG
      DOUBLE PRECISION BATM(4),CATM(4),LAHG(5),Lm(12)
      DOUBLE PRECISION H,RT, m
      DOUBLE PRECISION XR, T_Ray, RHO_TOT, RHO_FI
      
      DOUBLE PRECISION Rcos2, Rsin2, pi, hscale
      DOUBLE PRECISION x1, x2, x3, x4
      DOUBLE PRECISION e1, e2, e3, e4

      REAL LONG(12)
      INTEGER I

c--   BATM, CATM, LAHG:
c--   some parameters of the US standard atmosphere (see Corsika manual, 
c--   appendix C). LAHG contains the limits of the 4 exponential layers, 
c--   in which the mass overburden is given by T = AATM + BATM * exp(-h/CATM)
c--   The parameters AATM are not included in this code because they are not 
c--   needed. The last layer of the US standard atmosphere (in which T varies
c--   linearly with h) is above 100 km and has not been included here for the
c--   same reason.
c--
      DATA BATM / 1222.6562D0,1144.9069D0,1305.5948D0,540.1778D0 /
      DATA CATM / 994186.38D0,878153.55D0,636143.04D0,772170.16D0 /
      DATA LAHG / 0.,4.0D5,1.0D6,4.0D6,1.0D7 /

      DATA Lm /3.70D5,3.85D5,4.0D5,4.17D5,4.17D5,4.35D5,5.0D5,5.56D5,
     +         5.99D5,6.33D5, 6.67D5,7.04D5 /
      DATA LONG /
     +    280.0, 
     +    300.0, 
     +    320.0,
     +    340.0,
     +    360.0,
     +    380.0,
     +    400.0,
     +    450.0,
     +    500.0,
     +    550.0,
     +    600.0,
     +    650.0 /
      DATA PI / 3.141592654D0 /

c--   Take same Earth radius as in CORSIKA
      parameter (rt = 6371315.D2)

c--   Scale-height for Exponential density profile
      parameter (hscale = 7.4d5)

c--   Mean free path for scattering Rayleigh XR (g/cm^2)
      parameter (XR=2.970D3)

c-- Set low limit of first atmospheric layer to the observation level
c-- so that the traversed atmospheric depth in the Rayleigh scattering 
c-- will be calculated correctly.

      LAHG(1) = obslev

c-- For the case of simulating a telescope higher than 4 km...

      if (obslev .gt. LAHG(2)) then
         LAHG(2) = obslev
      end if

      T_Ray = 1.0


c--   AM, Jan 2002: now the argument height is directly the height above 
c--   sea level, calculated in atm.c:

      h = height

***********************************************************************
*
*     LARGE ZENITH ANGLE FACTOR (AIR MASS FACTOR):
      
c--   Air mass factor "m" calculated using a one-exponential density
c--   profile for the atmosphere, rho = rho_0 exp(-h/hscale) with 
c--   hscale = 7.4 km. The air mass factor is defined as I(theta)/I(0), 
c--   the ratio of the optical paths I (in g/cm2) traversed between two 
c--   given heights, at theta and at 0 deg z.a. 

      Rcos2 = rt * cos(theta)**2
      Rsin2 = rt * sin(theta)**2

*     
*     Obsolete lines:
*     x1 = sqrt((2. * obslev + Rcos2) / (2. * hscale))
*     x2 = sqrt((2. * h      + Rcos2) / (2. * hscale))
*     x3 = sqrt((2. * obslev + rt)    / (2. * hscale))
*     x4 = sqrt((2. * h      + rt)    / (2. * hscale))
*     

c--   AM, Dec 2002: slightly changed the calculation of the air mass factor,
c--   for what I think is a better approximation (numerically the results are 
c--   almost exactly the same, this change is irrelevant!):

      x1 = sqrt(Rcos2 / (2*hscale))
      x2 = sqrt((2*(h-obslev) + Rcos2) / (2*hscale))
      x3 = sqrt(rt / (2*hscale))
      x4 = sqrt((2*(h-obslev) + rt) / (2*hscale))

c--   AM Dec 2001, to avoid crash! Sometime a few photons seem to be 
c--   "corrupted" (have absurd values) in cer files... Next lines avoid the
c--   crash. They will also return a -1 transmittance in the case the photon 
c--   comes exactly from the observation level, because in that case the 
c--   "air mass factor" would become infinity and the calculation is not valid!

      if (abs(x3-x4) .lt. 1.e-10) then
         tr_atmos = -1.
         RETURN
      endif

      e1 = derfc(x1)
      e2 = derfc(x2)
      e3 = derfc(x3)
      e4 = derfc(x4)

      m = exp(-Rsin2 / (2. * hscale)) * ((e1 - e2) / (e3 - e4))

**********************************************************************   
*     RAYLEIGH SCATTERING (by the molecules in the atmosphere)
*     SCATTERING M.F.P. FOR RAYLEIGH: XR = 2.970D3
**********************************************************************

c--   Calculate the traversed "vertical thickness" of air using the
c--   US Standard atmosphere:

      RHO_TOT = 0.0

      DO 91 I=2,5
        IF ( H .LT. LAHG(I) ) THEN
          RHO_TOT = RHO_TOT +
     +        BATM(I-1)*(EXP(-LAHG(I-1)/CATM(I-1) ) -
     +        EXP(-H/CATM(I-1)))
          GOTO 92
        ELSE
          RHO_TOT = RHO_TOT +
     +        BATM(I-1)*(EXP(-LAHG(I-1)/CATM(I-1) ) -
     +        EXP(-LAHG(I)/CATM(I-1)))
        ENDIF
 91   CONTINUE
 
c--   We now convert from "vertical thickness" to "slanted thickness"
c--   traversed by the photon on its way to the telescope, simply
c--   multiplying by the air mass factor m:

 92   RHO_FI = m*RHO_TOT

c--   Finally we calculate the transmission coefficient for the Rayleigh
c--   scattering:
c--   AM Dec 2002, introduced ABS below to account (in the future) for 
c--   possible photons coming from below the observation level.

      T_Ray = EXP(-ABS(RHO_FI/XR)*(400.D0/wavelength)**4)
   
      tr_atmos = T_Ray

      RETURN

      END
