| 1 | *********************************************************************
|
|---|
| 2 | * *
|
|---|
| 3 | * Atmospheric Absorption for Rayleigh, Mie and Ozone. *
|
|---|
| 4 | * *
|
|---|
| 5 | * Created: May-98 *
|
|---|
| 6 | * Author: Aitor Ibarra Ibaibarriaga *
|
|---|
| 7 | * Jose Carlos Gonzalez *
|
|---|
| 8 | * *
|
|---|
| 9 | * Modified December-2002, January-2003 *
|
|---|
| 10 | * Author: Abelardo Moralejo (moralejo@pd.infn.it) *
|
|---|
| 11 | * *
|
|---|
| 12 | * Now this accounts only for Rayleigh scattering. Mie and *
|
|---|
| 13 | * Ozone absorption are now treated separatedly. *
|
|---|
| 14 | * *
|
|---|
| 15 | *********************************************************************
|
|---|
| 16 | *
|
|---|
| 17 | * December 2002, Abelardo Moralejo: checked algorithms, removed
|
|---|
| 18 | * old/unnecessary code, fixed some bugs, added comments.
|
|---|
| 19 | *
|
|---|
| 20 | * Fixed bugs (of small influence) in Mie absorption implementation: there were
|
|---|
| 21 | * errors in the optical depths table, as well as a confusion: height a.s.l.
|
|---|
| 22 | * was used as if it was height above the telescope level. The latter error was
|
|---|
| 23 | * also present in the Ozone aborption implementation.
|
|---|
| 24 | *
|
|---|
| 25 | * On the other hand, now we have the tables AE_ABI and OZ_ABI with optical
|
|---|
| 26 | * depths between sea level and a height h (before it was between 2km a.s.l
|
|---|
| 27 | * and a height h). So that we can simulate also in the future a different
|
|---|
| 28 | * observation level.
|
|---|
| 29 | *
|
|---|
| 30 | * AM: WARNING: IF VERY LARGE zenith angle simulations are to be done (say
|
|---|
| 31 | * above 85 degrees, for neutrino primaries or any other purpose) this code
|
|---|
| 32 | * will have to be adapted accordingly and checked, since in principle it has
|
|---|
| 33 | * been written and tested to simulate the absorption of Cherenkov photons
|
|---|
| 34 | * arriving at the telescope from above.
|
|---|
| 35 | *
|
|---|
| 36 | * AM: WARNING 2: not to be used for wavelengths outside the range 250-700 nm
|
|---|
| 37 | *
|
|---|
| 38 | * January 2003, Abelardo Moralejo: found error in Ozone absorption treatment.
|
|---|
| 39 | * At large zenith angles, the air mass used was the one calculated for
|
|---|
| 40 | * Rayleigh scattering, but since the Ozone distribution is rather different
|
|---|
| 41 | * from the global distribution of air molecules, this is not a good
|
|---|
| 42 | * approximation. Now I have left in this code only the Rayleigh scattering,
|
|---|
| 43 | * and moved to atm.c the Mie scattering and Ozone absorption calculated in
|
|---|
| 44 | * a better way.
|
|---|
| 45 | *
|
|---|
| 46 | * AM, Jan 2003: added obslev as an argument of the function. Changed the
|
|---|
| 47 | * meaning of the argument height: now it is the true height above sea level
|
|---|
| 48 | * at which a photon has been emitted, before it was the height given by
|
|---|
| 49 | * Corsika, measured in the vertical of the observer and not in the vertical
|
|---|
| 50 | * of the emitting particle.
|
|---|
| 51 | *
|
|---|
| 52 | * @endtext
|
|---|
| 53 | c @code
|
|---|
| 54 |
|
|---|
| 55 | SUBROUTINE attenu(wavelength, height, obslev, theta, tr_atmos)
|
|---|
| 56 |
|
|---|
| 57 | REAL wavelength, height, obslev, theta, tr_atmos
|
|---|
| 58 | COMMON /ATMOS/ BATM,CATM,LAHG,LONG
|
|---|
| 59 | DOUBLE PRECISION BATM(4),CATM(4),LAHG(5),Lm(12)
|
|---|
| 60 | DOUBLE PRECISION H,RT, m
|
|---|
| 61 | DOUBLE PRECISION XR, T_Ray, RHO_TOT, RHO_FI
|
|---|
| 62 |
|
|---|
| 63 | DOUBLE PRECISION Rcos2, Rsin2, pi, hscale
|
|---|
| 64 | DOUBLE PRECISION x1, x2, x3, x4
|
|---|
| 65 | DOUBLE PRECISION e1, e2, e3, e4
|
|---|
| 66 |
|
|---|
| 67 | REAL LONG(12)
|
|---|
| 68 | INTEGER I
|
|---|
| 69 |
|
|---|
| 70 | c-- BATM, CATM, LAHG:
|
|---|
| 71 | c-- some parameters of the US standard atmosphere (see Corsika manual,
|
|---|
| 72 | c-- appendix C). LAHG contains the limits of the 4 exponential layers,
|
|---|
| 73 | c-- in which the mass overburden is given by T = AATM + BATM * exp(-h/CATM)
|
|---|
| 74 | c-- The parameters AATM are not included in this code because they are not
|
|---|
| 75 | c-- needed. The last layer of the US standard atmosphere (in which T varies
|
|---|
| 76 | c-- linearly with h) is above 100 km and has not been included here for the
|
|---|
| 77 | c-- same reason.
|
|---|
| 78 | c--
|
|---|
| 79 | DATA BATM / 1222.6562D0,1144.9069D0,1305.5948D0,540.1778D0 /
|
|---|
| 80 | DATA CATM / 994186.38D0,878153.55D0,636143.04D0,772170.16D0 /
|
|---|
| 81 | DATA LAHG / 0.,4.0D5,1.0D6,4.0D6,1.0D7 /
|
|---|
| 82 |
|
|---|
| 83 | DATA Lm /3.70D5,3.85D5,4.0D5,4.17D5,4.17D5,4.35D5,5.0D5,5.56D5,
|
|---|
| 84 | + 5.99D5,6.33D5, 6.67D5,7.04D5 /
|
|---|
| 85 | DATA LONG /
|
|---|
| 86 | + 280.0,
|
|---|
| 87 | + 300.0,
|
|---|
| 88 | + 320.0,
|
|---|
| 89 | + 340.0,
|
|---|
| 90 | + 360.0,
|
|---|
| 91 | + 380.0,
|
|---|
| 92 | + 400.0,
|
|---|
| 93 | + 450.0,
|
|---|
| 94 | + 500.0,
|
|---|
| 95 | + 550.0,
|
|---|
| 96 | + 600.0,
|
|---|
| 97 | + 650.0 /
|
|---|
| 98 | DATA PI / 3.141592654D0 /
|
|---|
| 99 |
|
|---|
| 100 | c-- Take same Earth radius as in CORSIKA
|
|---|
| 101 | parameter (rt = 6371315.D2)
|
|---|
| 102 |
|
|---|
| 103 | c-- Scale-height for Exponential density profile
|
|---|
| 104 | parameter (hscale = 7.4d5)
|
|---|
| 105 |
|
|---|
| 106 | c-- Mean free path for scattering Rayleigh XR (g/cm^2)
|
|---|
| 107 | parameter (XR=2.970D3)
|
|---|
| 108 |
|
|---|
| 109 | c-- Set low limit of first atmospheric layer to the observation level
|
|---|
| 110 | c-- so that the traversed atmospheric depth in the Rayleigh scattering
|
|---|
| 111 | c-- will be calculated correctly.
|
|---|
| 112 |
|
|---|
| 113 | LAHG(1) = obslev
|
|---|
| 114 |
|
|---|
| 115 | c-- For the case of simulating a telescope higher than 4 km...
|
|---|
| 116 |
|
|---|
| 117 | if (obslev .gt. LAHG(2)) then
|
|---|
| 118 | LAHG(2) = obslev
|
|---|
| 119 | end if
|
|---|
| 120 |
|
|---|
| 121 | T_Ray = 1.0
|
|---|
| 122 |
|
|---|
| 123 |
|
|---|
| 124 | c-- AM, Jan 2002: now the argument height is directly the height above
|
|---|
| 125 | c-- sea level, calculated in atm.c:
|
|---|
| 126 |
|
|---|
| 127 | h = height
|
|---|
| 128 |
|
|---|
| 129 | ***********************************************************************
|
|---|
| 130 | *
|
|---|
| 131 | * LARGE ZENITH ANGLE FACTOR (AIR MASS FACTOR):
|
|---|
| 132 |
|
|---|
| 133 | c-- Air mass factor "m" calculated using a one-exponential density
|
|---|
| 134 | c-- profile for the atmosphere, rho = rho_0 exp(-h/hscale) with
|
|---|
| 135 | c-- hscale = 7.4 km. The air mass factor is defined as I(theta)/I(0),
|
|---|
| 136 | c-- the ratio of the optical paths I (in g/cm2) traversed between two
|
|---|
| 137 | c-- given heights, at theta and at 0 deg z.a.
|
|---|
| 138 |
|
|---|
| 139 | Rcos2 = rt * cos(theta)**2
|
|---|
| 140 | Rsin2 = rt * sin(theta)**2
|
|---|
| 141 |
|
|---|
| 142 | *
|
|---|
| 143 | * Obsolete lines:
|
|---|
| 144 | * x1 = sqrt((2. * obslev + Rcos2) / (2. * hscale))
|
|---|
| 145 | * x2 = sqrt((2. * h + Rcos2) / (2. * hscale))
|
|---|
| 146 | * x3 = sqrt((2. * obslev + rt) / (2. * hscale))
|
|---|
| 147 | * x4 = sqrt((2. * h + rt) / (2. * hscale))
|
|---|
| 148 | *
|
|---|
| 149 |
|
|---|
| 150 | c-- AM, Dec 2002: slightly changed the calculation of the air mass factor,
|
|---|
| 151 | c-- for what I think is a better approximation (numerically the results are
|
|---|
| 152 | c-- almost exactly the same, this change is irrelevant!):
|
|---|
| 153 |
|
|---|
| 154 | x1 = sqrt(Rcos2 / (2*hscale))
|
|---|
| 155 | x2 = sqrt((2*(h-obslev) + Rcos2) / (2*hscale))
|
|---|
| 156 | x3 = sqrt(rt / (2*hscale))
|
|---|
| 157 | x4 = sqrt((2*(h-obslev) + rt) / (2*hscale))
|
|---|
| 158 |
|
|---|
| 159 | c-- AM Dec 2001, to avoid crash! Sometime a few photons seem to be
|
|---|
| 160 | c-- "corrupted" (have absurd values) in cer files... Next lines avoid the
|
|---|
| 161 | c-- crash. They will also return a -1 transmittance in the case the photon
|
|---|
| 162 | c-- comes exactly from the observation level, because in that case the
|
|---|
| 163 | c-- "air mass factor" would become infinity and the calculation is not valid!
|
|---|
| 164 |
|
|---|
| 165 | if (abs(x3-x4) .lt. 1.e-10) then
|
|---|
| 166 | tr_atmos = -1.
|
|---|
| 167 | RETURN
|
|---|
| 168 | endif
|
|---|
| 169 |
|
|---|
| 170 | e1 = derfc(x1)
|
|---|
| 171 | e2 = derfc(x2)
|
|---|
| 172 | e3 = derfc(x3)
|
|---|
| 173 | e4 = derfc(x4)
|
|---|
| 174 |
|
|---|
| 175 | m = exp(-Rsin2 / (2. * hscale)) * ((e1 - e2) / (e3 - e4))
|
|---|
| 176 |
|
|---|
| 177 | **********************************************************************
|
|---|
| 178 | * RAYLEIGH SCATTERING (by the molecules in the atmosphere)
|
|---|
| 179 | * SCATTERING M.F.P. FOR RAYLEIGH: XR = 2.970D3
|
|---|
| 180 | **********************************************************************
|
|---|
| 181 |
|
|---|
| 182 | c-- Calculate the traversed "vertical thickness" of air using the
|
|---|
| 183 | c-- US Standard atmosphere:
|
|---|
| 184 |
|
|---|
| 185 | RHO_TOT = 0.0
|
|---|
| 186 |
|
|---|
| 187 | DO 91 I=2,5
|
|---|
| 188 | IF ( H .LT. LAHG(I) ) THEN
|
|---|
| 189 | RHO_TOT = RHO_TOT +
|
|---|
| 190 | + BATM(I-1)*(EXP(-LAHG(I-1)/CATM(I-1) ) -
|
|---|
| 191 | + EXP(-H/CATM(I-1)))
|
|---|
| 192 | GOTO 92
|
|---|
| 193 | ELSE
|
|---|
| 194 | RHO_TOT = RHO_TOT +
|
|---|
| 195 | + BATM(I-1)*(EXP(-LAHG(I-1)/CATM(I-1) ) -
|
|---|
| 196 | + EXP(-LAHG(I)/CATM(I-1)))
|
|---|
| 197 | ENDIF
|
|---|
| 198 | 91 CONTINUE
|
|---|
| 199 |
|
|---|
| 200 | c-- We now convert from "vertical thickness" to "slanted thickness"
|
|---|
| 201 | c-- traversed by the photon on its way to the telescope, simply
|
|---|
| 202 | c-- multiplying by the air mass factor m:
|
|---|
| 203 |
|
|---|
| 204 | 92 RHO_FI = m*RHO_TOT
|
|---|
| 205 |
|
|---|
| 206 | c-- Finally we calculate the transmission coefficient for the Rayleigh
|
|---|
| 207 | c-- scattering:
|
|---|
| 208 | c-- AM Dec 2002, introduced ABS below to account (in the future) for
|
|---|
| 209 | c-- possible photons coming from below the observation level.
|
|---|
| 210 |
|
|---|
| 211 | T_Ray = EXP(-ABS(RHO_FI/XR)*(400.D0/wavelength)**4)
|
|---|
| 212 |
|
|---|
| 213 | tr_atmos = T_Ray
|
|---|
| 214 |
|
|---|
| 215 | RETURN
|
|---|
| 216 |
|
|---|
| 217 | END
|
|---|
