********************************************************************* * * * 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