Ignore:
Timestamp:
01/21/03 15:02:35 (22 years ago)
Author:
moralejo
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Simulation/Detector/ReflectorII/attenu.f

    r1431 r1722  
    77*             Jose Carlos Gonzalez                                  *
    88*                                                                   *
     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*                                                                   *
    915*********************************************************************
    10 
    11 c @T \newpage
    12 
    13 
    14 c @section Source code of {\tt attenu.f}
    15 
    16 * @text
    17 * Copyright $\copyright$ 1998, Aitor Ibarra Ibaibarriaga
     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*
    1852* @endtext
    1953c @code
    2054
    21       SUBROUTINE attenu(wavelength, height, theta, tr_atmos)
    22 
    23 C----------------------------------------------------------------------C
    24 C  RHO (DENSITY) F(UNCTION)                                            C
    25 C                                                                      C
    26 C  CALCULATES DENSITY (G/CM**3) OF ATMOSPHERE DEPENDING ON HEIGHT (CM) C
    27 C  (US STANDARD ATMOSPHERE)                                            C
    28 C  THIS FUNCTION IS CALLED FROM ININKG, UPDATE, CERENE, CERENH         C
    29 C  ARGUMENT:                                                           C
    30 C   ARG    = HEIGHT IN CM                                              C
    31 C----------------------------------------------------------------------C
    32 
    33       COMMON /ATMOS/   AATM,BATM,CATM,LAHG,RHOF,LONG,OZ_ABI,AE_ABI
    34       DOUBLE PRECISION AATM(5),BATM(5),CATM(5),LAHG(5),RHOF(5),Lm(12)
    35       DOUBLE PRECISION H,RT, m,OZ_ABI(48,12),AE_ABI(28,12)
    36       DOUBLE PRECISION XR, TOT_OZ, TOT_AE, T_Ray,T_Mie, T_Oz,
    37      + RHO_TOT, RHO_FI, RHOFP
     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
    3862     
    39       DOUBLE PRECISION Rcos2, Rsin2
     63      DOUBLE PRECISION Rcos2, Rsin2, pi, hscale
    4064      DOUBLE PRECISION x1, x2, x3, x4
    4165      DOUBLE PRECISION e1, e2, e3, e4
    4266
    43       REAL wavelength, height, theta, tr_atmos
    44       real trr, trm, tro
    4567      REAL LONG(12)
    46 c fs: define obervation level
    47       double precision obslev     
    48       INTEGER I,CON_OZ,CON_MI J, ROW
    49 
    50       DATA OZ_ABI /
    51      + 0.2880000D+00,0.5405000D+00,0.7775000D+00,0.1009000D+01,
    52      + 0.1241500D+01,0.1480500D+01,0.1750500D+01,0.2085000D+01,
    53      + 0.2514500D+01,0.3087500D+01,0.3864500D+01,0.4817500D+01,
    54      + 0.5847500D+01,0.6917500D+01,0.8052500D+01,0.9287499D+01,
    55      + 0.1068750D+02,0.1231250D+02,0.1415750D+02,0.1617750D+02,
    56      + 0.1827250D+02,0.2034750D+02,0.2232750D+02,0.2414750D+02,
    57      + 0.2575750D+02,0.2715250D+02,0.2836750D+02,0.2941100D+02,
    58      + 0.3031000D+02,0.3109250D+02,0.3176300D+02,0.3232750D+02,
    59      + 0.3281200D+02,0.3323200D+02,0.3358350D+02,0.3387750D+02,
    60      + 0.3412650D+02,0.3434000D+02,0.3451900D+02,0.3466250D+02,
    61      + 0.3477480D+02,0.3486355D+02,0.3493355D+02,0.3498775D+02,
    62      + 0.3503010D+02,0.3506360D+02,0.3509020D+02,0.3511185D+02,
    63      + 0.2740000D-01,0.5140000D-01,0.7395000D-01,0.9600000D-01,
    64      + 0.1181500D+00,0.1409000D+00,0.1666000D+00,0.1984500D+00,
    65      + 0.2393500D+00,0.2939500D+00,0.3679500D+00,0.4589500D+00,
    66      + 0.5573000D+00,0.6593000D+00,0.7673000D+00,0.8848000D+00,
    67      + 0.1017800D+01,0.1172300D+01,0.1348300D+01,0.1540800D+01,
    68      + 0.1740300D+01,0.1937800D+01,0.2126300D+01,0.2299800D+01,
    69      + 0.2453300D+01,0.2586300D+01,0.2702300D+01,0.2801900D+01,
    70      + 0.2887550D+01,0.2962050D+01,0.3025900D+01,0.3079700D+01,
    71      + 0.3125850D+01,0.3165850D+01,0.3199350D+01,0.3227400D+01,
    72      + 0.3251150D+01,0.3271500D+01,0.3288600D+01,0.3302300D+01,
    73      + 0.3312995D+01,0.3321445D+01,0.3328110D+01,0.3333270D+01,
    74      + 0.3337305D+01,0.3340500D+01,0.3343035D+01,0.1334510D+01,
    75      + 0.2435000D-02,0.4570000D-02,0.6575000D-02,0.8535000D-02,
    76      + 0.1050500D-01,0.1253000D-01,0.1481500D-01,0.1764500D-01,
    77      + 0.2128000D-01,0.2613500D-01,0.3272000D-01,0.4081000D-01,
    78      + 0.4957000D-01,0.5866000D-01,0.6827000D-01,0.7875500D-01,
    79      + 0.9065500D-01,0.1044050D+00,0.1200050D+00,0.1371050D+00,
    80      + 0.1548550D+00,0.1724050D+00,0.1891550D+00,0.2045550D+00,
    81      + 0.2182050D+00,0.2300550D+00,0.2403600D+00,0.2492200D+00,
    82      + 0.2568350D+00,0.2634550D+00,0.2691300D+00,0.2739150D+00,
    83      + 0.2780200D+00,0.2815750D+00,0.2845500D+00,0.2870400D+00,
    84      + 0.2891500D+00,0.2909600D+00,0.2924750D+00,0.2936900D+00,
    85      + 0.2946425D+00,0.2953940D+00,0.2959865D+00,0.2964455D+00,
    86      + 0.2968045D+00,0.2970885D+00,0.2973140D+00,0.2974975D+00,
    87      + 0.1740000D-03,0.3265000D-03,0.4695000D-03,0.6090000D-03,
    88      + 0.7495000D-03,0.8940000D-03,0.1057000D-02,0.1259000D-02,
    89      + 0.1518000D-02,0.1863500D-02,0.2332500D-02,0.2909000D-02,
    90      + 0.3533000D-02,0.4180500D-02,0.4865000D-02,0.5610500D-02,
    91      + 0.6455500D-02,0.7435000D-02,0.8550000D-02,0.9770000D-02,
    92      + 0.1103500D-01,0.1229000D-01,0.1348500D-01,0.1458000D-01,
    93      + 0.1555100D-01,0.1639550D-01,0.1713150D-01,0.1776300D-01,
    94      + 0.1830600D-01,0.1877800D-01,0.1918200D-01,0.1952250D-01,
    95      + 0.1981500D-01,0.2006850D-01,0.2028050D-01,0.2045801D-01,
    96      + 0.2060851D-01,0.2073750D-01,0.2084566D-01,0.2093241D-01,
    97      + 0.2100026D-01,0.2105381D-01,0.2109606D-01,0.2112876D-01,
    98      + 0.2115431D-01,0.2117456D-01,0.2119066D-01,0.2120376D-01,
    99      + 0.4885000D-05,0.9170000D-05,0.1319500D-04,0.1713000D-04,
    100      + 0.2108000D-04,0.2513500D-04,0.2971500D-04,0.3539500D-04,
    101      + 0.4268500D-04,0.5242500D-04,0.6562500D-04,0.8182500D-04,
    102      + 0.9937500D-04,0.1175750D-03,0.1368250D-03,0.1578250D-03,
    103      + 0.1816250D-03,0.2091750D-03,0.2404750D-03,0.2747750D-03,
    104      + 0.3103250D-03,0.3454750D-03,0.3790250D-03,0.4098750D-03,
    105      + 0.4372250D-03,0.4609750D-03,0.4816750D-03,0.4994750D-03,
    106      + 0.5147750D-03,0.5280750D-03,0.5394750D-03,0.5490700D-03,
    107      + 0.5572950D-03,0.5644250D-03,0.5703950D-03,0.5753899D-03,
    108      + 0.5796200D-03,0.5832500D-03,0.5862950D-03,0.5887350D-03,
    109      + 0.5906400D-03,0.5921450D-03,0.5933350D-03,0.5942565D-03,
    110      + 0.5949755D-03,0.5955440D-03,0.5959955D-03,0.5963635D-03,
    111      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    112      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    113      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    114      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    115      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    116      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    117      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    118      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    119      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    120      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    121      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    122      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    123      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    124      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    125      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    126      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    127      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    128      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    129      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    130      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    131      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    132      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    133      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    134      + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
    135      + 0.9525000D-05,0.1785500D-04,0.2567500D-04,0.3332000D-04,
    136      + 0.4100000D-04,0.4889000D-04,0.5779500D-04,0.6881000D-04,
    137      + 0.8296000D-04,0.1018600D-03,0.1275100D-03,0.1590600D-03,
    138      + 0.1932100D-03,0.2286100D-03,0.2660100D-03,0.3067600D-03,
    139      + 0.3529600D-03,0.4065100D-03,0.4674100D-03,0.5341100D-03,
    140      + 0.6032600D-03,0.6717100D-03,0.7370100D-03,0.7970100D-03,
    141      + 0.8501600D-03,0.8963600D-03,0.9366100D-03,0.9711101D-03,
    142      + 0.1000810D-02,0.1026660D-02,0.1048810D-02,0.1067460D-02,
    143      + 0.1083460D-02,0.1097310D-02,0.1108910D-02,0.1118640D-02,
    144      + 0.1126865D-02,0.1133915D-02,0.1139830D-02,0.1144570D-02,
    145      + 0.1148275D-02,0.1151200D-02,0.1153510D-02,0.1155300D-02,
    146      + 0.1156700D-02,0.1162205D-02,0.1170990D-02,0.1178145D-02,
    147      + 0.9360000D-04,0.1757000D-03,0.2528000D-03,0.3281500D-03,
    148      + 0.4038500D-03,0.4816500D-03,0.5694500D-03,0.6784000D-03,
    149      + 0.8184000D-03,0.1004900D-02,0.1257900D-02,0.1568900D-02,
    150      + 0.1905400D-02,0.2254400D-02,0.2623400D-02,0.3025400D-02,
    151      + 0.3480900D-02,0.4008900D-02,0.4609400D-02,0.5266900D-02,
    152      + 0.5948400D-02,0.6622900D-02,0.7266400D-02,0.7857900D-02,
    153      + 0.8381900D-02,0.8836901D-02,0.9233401D-02,0.9573901D-02,
    154      + 0.9866901D-02,0.1012140D-01,0.1033940D-01,0.1052340D-01,
    155      + 0.1068140D-01,0.1081840D-01,0.1093290D-01,0.1102855D-01,
    156      + 0.1110965D-01,0.1117920D-01,0.1123750D-01,0.1128425D-01,
    157      + 0.1132085D-01,0.1134975D-01,0.1137255D-01,0.1139020D-01,
    158      + 0.1140400D-01,0.1141492D-01,0.1142358D-01,0.1143063D-01,
    159      + 0.2500000D-03,0.4690000D-03,0.6745000D-03,0.8755000D-03,
    160      + 0.1077500D-02,0.1285000D-02,0.1519500D-02,0.1810000D-02,
    161      + 0.2182500D-02,0.2679500D-02,0.3353500D-02,0.4182000D-02,
    162      + 0.5079000D-02,0.6010000D-02,0.6994000D-02,0.8064000D-02,
    163      + 0.9279000D-02,0.1068900D-01,0.1228900D-01,0.1403900D-01,
    164      + 0.1585400D-01,0.1765400D-01,0.1937400D-01,0.2095400D-01,
    165      + 0.2235400D-01,0.2356900D-01,0.2462600D-01,0.2553350D-01,
    166      + 0.2631400D-01,0.2699250D-01,0.2757350D-01,0.2806300D-01,
    167      + 0.2848350D-01,0.2884800D-01,0.2915300D-01,0.2940850D-01,
    168      + 0.2962500D-01,0.2981050D-01,0.2996600D-01,0.3009050D-01,
    169      + 0.3018780D-01,0.3026480D-01,0.3032550D-01,0.3037250D-01,
    170      + 0.3040925D-01,0.3043835D-01,0.3046145D-01,0.3048025D-01,
    171      + 0.3585000D-03,0.6725000D-03,0.9675000D-03,0.1256000D-02,
    172      + 0.1545500D-02,0.1843000D-02,0.2179000D-02,0.2595500D-02,
    173      + 0.3130000D-02,0.3843500D-02,0.4813500D-02,0.6003500D-02,
    174      + 0.7288500D-02,0.8623499D-02,0.1003850D-01,0.1157850D-01,
    175      + 0.1331850D-01,0.1533350D-01,0.1762850D-01,0.2014350D-01,
    176      + 0.2274850D-01,0.2532850D-01,0.2779350D-01,0.3005850D-01,
    177      + 0.3206350D-01,0.3380350D-01,0.3531851D-01,0.3661850D-01,
    178      + 0.3773851D-01,0.3871351D-01,0.3954751D-01,0.4025051D-01,
    179      + 0.4085401D-01,0.4137701D-01,0.4181501D-01,0.4218151D-01,
    180      + 0.4249151D-01,0.4275751D-01,0.4298101D-01,0.4316001D-01,
    181      + 0.4330001D-01,0.4341061D-01,0.4349771D-01,0.4356516D-01,
    182      + 0.4361791D-01,0.4365961D-01,0.4369271D-01,0.4371971D-01,
    183      + 0.1685000D-03,0.3160000D-03,0.4545000D-03,0.5900000D-03,
    184      + 0.7260000D-03,0.8655000D-03,0.1023000D-02,0.1218500D-02,
    185      + 0.1469500D-02,0.1804500D-02,0.2259000D-02,0.2817500D-02,
    186      + 0.3422000D-02,0.4049500D-02,0.4712999D-02,0.5434999D-02,
    187      + 0.6252999D-02,0.7203000D-02,0.8283000D-02,0.9462999D-02,
    188      + 0.1068800D-01,0.1190300D-01,0.1306300D-01,0.1412800D-01,
    189      + 0.1507000D-01,0.1588850D-01,0.1660150D-01,0.1721300D-01,
    190      + 0.1773900D-01,0.1819650D-01,0.1858850D-01,0.1891850D-01,
    191      + 0.1920150D-01,0.1944700D-01,0.1965250D-01,0.1982450D-01,
    192      + 0.1997050D-01,0.2009550D-01,0.2020010D-01,0.2028410D-01,
    193      + 0.2034985D-01,0.2040175D-01,0.2044265D-01,0.2047435D-01,
    194      + 0.2049915D-01,0.2051875D-01,0.2053430D-01,0.2054695D-01 /
    195 
    196       DATA AE_ABI /
    197      + 0.3645000E-01,0.5211000E-01,0.5913000E-01,0.6203000E-01,
    198      + 0.6304000E-01,0.6340450E-01,0.6353275E-01,0.6358405E-01,
    199      + 0.6361780E-01,0.6364885E-01,0.6367990E-01,0.6371365E-01,
    200      + 0.6375955E-01,0.6383380E-01,0.6392965E-01,0.6403360E-01,
    201      + 0.6414810E-01,0.6426660E-01,0.6438010E-01,0.6448960E-01,
    202      + 0.6459510E-01,0.6468170E-01,0.6474110E-01,0.6478295E-01,
    203      + 0.6481670E-01,0.6484775E-01,0.6487610E-01,0.6490240E-01,
    204      + 0.3510000E-01,0.5018000E-01,0.5694000E-01,0.5973500E-01,
    205      + 0.6071000E-01,0.6106100E-01,0.6118450E-01,0.6123390E-01,
    206      + 0.6126640E-01,0.6129630E-01,0.6132620E-01,0.6135870E-01,
    207      + 0.6140290E-01,0.6147440E-01,0.6156670E-01,0.6166680E-01,
    208      + 0.6177730E-01,0.6189180E-01,0.6200130E-01,0.6210680E-01,
    209      + 0.6220820E-01,0.6229140E-01,0.6234860E-01,0.6238890E-01,
    210      + 0.6242140E-01,0.6245130E-01,0.6247860E-01,0.6250395E-01,
    211      + 0.3375000E-01,0.4825000E-01,0.5475000E-01,0.5744000E-01,
    212      + 0.5838000E-01,0.5871750E-01,0.5883625E-01,0.5888375E-01,
    213      + 0.5891500E-01,0.5894375E-01,0.5897250E-01,0.5900375E-01,
    214      + 0.5904625E-01,0.5911500E-01,0.5920375E-01,0.5930000E-01,
    215      + 0.5940650E-01,0.5951700E-01,0.5962200E-01,0.5972300E-01,
    216      + 0.5982050E-01,0.5990050E-01,0.5995550E-01,0.5999425E-01,
    217      + 0.6002550E-01,0.6005425E-01,0.6008050E-01,0.6010490E-01,
    218      + 0.3240000E-01,0.4632000E-01,0.5256000E-01,0.5514000E-01,
    219      + 0.5604000E-01,0.5636400E-01,0.5647800E-01,0.5652360E-01,
    220      + 0.5655360E-01,0.5658120E-01,0.5660880E-01,0.5663880E-01,
    221      + 0.5667960E-01,0.5674560E-01,0.5683080E-01,0.5692320E-01,
    222      + 0.5702520E-01,0.5713070E-01,0.5723140E-01,0.5732860E-01,
    223      + 0.5742220E-01,0.5749900E-01,0.5755180E-01,0.5758900E-01,
    224      + 0.5761900E-01,0.5764660E-01,0.5767180E-01,0.5769520E-01,
    225      + 0.3240000E-01,0.4632000E-01,0.5256000E-01,0.5514000E-01,
    226      + 0.5604000E-01,0.5636400E-01,0.5647800E-01,0.5652360E-01,
    227      + 0.5655360E-01,0.5658120E-01,0.5660880E-01,0.5663880E-01,
    228      + 0.5667960E-01,0.5674560E-01,0.5683080E-01,0.5692320E-01,
    229      + 0.5702520E-01,0.5713070E-01,0.5723140E-01,0.5732860E-01,
    230      + 0.5742220E-01,0.5749900E-01,0.5755180E-01,0.5758900E-01,
    231      + 0.5761900E-01,0.5764660E-01,0.5767180E-01,0.5769520E-01,
    232      + 0.3105000E-01,0.4439000E-01,0.5037000E-01,0.5284500E-01,
    233      + 0.5371000E-01,0.5402050E-01,0.5412975E-01,0.5417345E-01,
    234      + 0.5420220E-01,0.5422865E-01,0.5425510E-01,0.5428385E-01,
    235      + 0.5432295E-01,0.5438620E-01,0.5446785E-01,0.5455640E-01,
    236      + 0.5465390E-01,0.5475485E-01,0.5485145E-01,0.5494460E-01,
    237      + 0.5503430E-01,0.5510790E-01,0.5515850E-01,0.5519415E-01,
    238      + 0.5522290E-01,0.5524935E-01,0.5527350E-01,0.5529590E-01,
    239      + 0.2700000E-01,0.3860000E-01,0.4380000E-01,0.4595000E-01,
    240      + 0.4670000E-01,0.4697000E-01,0.4706500E-01,0.4710300E-01,
    241      + 0.4712800E-01,0.4715100E-01,0.4717400E-01,0.4719900E-01,
    242      + 0.4723300E-01,0.4728800E-01,0.4735900E-01,0.4743600E-01,
    243      + 0.4752100E-01,0.4760900E-01,0.4769300E-01,0.4777400E-01,
    244      + 0.4785200E-01,0.4791600E-01,0.4796000E-01,0.4799100E-01,
    245      + 0.4801600E-01,0.4803900E-01,0.4806000E-01,0.4807950E-01,
    246      + 0.2430000E-01,0.3474000E-01,0.3942000E-01,0.4135500E-01,
    247      + 0.4203000E-01,0.4227300E-01,0.4235850E-01,0.4239270E-01,
    248      + 0.4241520E-01,0.4243590E-01,0.4245660E-01,0.4247910E-01,
    249      + 0.4250970E-01,0.4255920E-01,0.4262310E-01,0.4269240E-01,
    250      + 0.4276890E-01,0.4284810E-01,0.4292370E-01,0.4299660E-01,
    251      + 0.4306680E-01,0.4312440E-01,0.4316400E-01,0.4319190E-01,
    252      + 0.4321440E-01,0.4323510E-01,0.4325400E-01,0.4327155E-01,
    253      + 0.2255000E-01,0.3225500E-01,0.3659500E-01,0.3838900E-01,
    254      + 0.3901500E-01,0.3924050E-01,0.3931985E-01,0.3935155E-01,
    255      + 0.3937240E-01,0.3939160E-01,0.3941080E-01,0.3943165E-01,
    256      + 0.3946005E-01,0.3950600E-01,0.3956530E-01,0.3962960E-01,
    257      + 0.3970055E-01,0.3977400E-01,0.3984415E-01,0.3991180E-01,
    258      + 0.3997695E-01,0.4003040E-01,0.4006715E-01,0.4009305E-01,
    259      + 0.4011390E-01,0.4013310E-01,0.4015065E-01,0.4016695E-01,
    260      + 0.2130000E-01,0.3044500E-01,0.3455500E-01,0.3625450E-01,
    261      + 0.3684700E-01,0.3706050E-01,0.3713575E-01,0.3716575E-01,
    262      + 0.3718550E-01,0.3720370E-01,0.3722190E-01,0.3724165E-01,
    263      + 0.3726850E-01,0.3731195E-01,0.3736805E-01,0.3742890E-01,
    264      + 0.3749605E-01,0.3756555E-01,0.3763190E-01,0.3769590E-01,
    265      + 0.3775750E-01,0.3780805E-01,0.3784280E-01,0.3786725E-01,
    266      + 0.3788700E-01,0.3790520E-01,0.3792180E-01,0.3793720E-01,
    267      + 0.2025000E-01,0.2895000E-01,0.3285000E-01,0.3446250E-01,
    268      + 0.3502500E-01,0.3522750E-01,0.3529875E-01,0.3532725E-01,
    269      + 0.3534600E-01,0.3536325E-01,0.3538050E-01,0.3539925E-01,
    270      + 0.3542475E-01,0.3546600E-01,0.3551925E-01,0.3557700E-01,
    271      + 0.3564075E-01,0.3570675E-01,0.3576975E-01,0.3583050E-01,
    272      + 0.3588900E-01,0.3593700E-01,0.3597000E-01,0.3599325E-01,
    273      + 0.3601200E-01,0.3602925E-01,0.3604500E-01,0.3605960E-01,
    274      + 0.1920000E-01,0.2745500E-01,0.3114500E-01,0.3267050E-01,
    275      + 0.3320300E-01,0.3339470E-01,0.3346215E-01,0.3348915E-01,
    276      + 0.3350690E-01,0.3352320E-01,0.3353950E-01,0.3355725E-01,
    277      + 0.3358140E-01,0.3362045E-01,0.3367085E-01,0.3372550E-01,
    278      + 0.3378585E-01,0.3384835E-01,0.3390800E-01,0.3396550E-01,
    279      + 0.3402090E-01,0.3406635E-01,0.3409760E-01,0.3411965E-01,
    280      + 0.3413740E-01,0.3415370E-01,0.3416860E-01,0.3418245E-01 /
    281 
    282       DATA BATM / 1222.6562D0,1144.9069D0,1305.5948D0,540.1778D0,0.D0  /
    283       DATA CATM / 994186.38D0,878153.55D0,636143.04D0,772170.16D0,1.D-9/
    284 
    285       DATA LAHG / 2000.0D2,4.0D5,1.0D6,4.0D6,1.0D7 /
     68      INTEGER I
     69
     70c--   BATM, CATM, LAHG:
     71c--   some parameters of the US standard atmosphere (see Corsika manual,
     72c--   appendix C). LAHG contains the limits of the 4 exponential layers,
     73c--   in which the mass overburden is given by T = AATM + BATM * exp(-h/CATM)
     74c--   The parameters AATM are not included in this code because they are not
     75c--   needed. The last layer of the US standard atmosphere (in which T varies
     76c--   linearly with h) is above 100 km and has not been included here for the
     77c--   same reason.
     78c--
     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
    28683      DATA Lm /3.70D5,3.85D5,4.0D5,4.17D5,4.17D5,4.35D5,5.0D5,5.56D5,
    28784     +         5.99D5,6.33D5, 6.67D5,7.04D5 /
    … …  
    307104      parameter (hscale = 7.4d5)
    308105
     106c--   Mean free path for scattering Rayleigh XR (g/cm^2)
     107      parameter (XR=2.970D3)
     108
     109c-- Set low limit of first atmospheric layer to the observation level
     110c-- so that the traversed atmospheric depth in the Rayleigh scattering
     111c-- will be calculated correctly.
     112
     113      LAHG(1) = obslev
     114
     115c-- 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
     124c--   AM, Jan 2002: now the argument height is directly the height above
     125c--   sea level, calculated in atm.c:
     126
     127      h = height
     128
    309129***********************************************************************
    310130*
    311 *     SCATTERING PARAMETERS FOR RAYLEIGH:
    312 *
    313 *     MEAN FREE PATH FOR SCATTERING RAYLEIGH  (g/cm^2)
    314       PARAMETER (XR=2.970D3)
    315 ***********************************************************************
    316        
    317 c-- Observation level at La Palma
    318       parameter (obslev = 2200.d2)
     131*     LARGE ZENITH ANGLE FACTOR (AIR MASS FACTOR):
    319132     
    320       T_Ray = 1.0
    321       T_Mie = 1.0
    322       T_Oz = 1.0
    323      
    324 c-- Height calculated using an obslev = H(LaPalma) <> 0
    325 c      H = -RT + SQRT(RT**2 + (height/COS(theta))**2 +
    326 c     +     (2.0D0*RT*height))
    327       h = -rt + sqrt((rt+obslev)**2 +
    328      +     ((height-obslev)/cos(theta))**2 +
    329      +     (2.0d0*(rt+obslev)*(height-obslev)))
    330 
    331       ROW = AINT(((H+1.)/1.0E5))     
    332              
    333 ***********************************************************************
    334 *
    335 *     LARGE ZENITH ANGLE FACTOR (AIR MASS):
    336      
    337 c fs : air mass factor = path_lenght(za) / path_lenght(vertical)
    338 c                        at point of emission of cherenkov photon
    339 c      => pure geometric correction on
    340 c         absorption lenght measured at vertical height [km^-1]   
    341 c--
    342 c      ma = (1.D0/H)*(SQRT((RT*COS(theta))**2+(2*RT*H)+H**2)-
    343 c     +     (RT*COS(theta)))
    344 c--
    345 c      mb = (1.d0/(h-obslev)) *
    346 c     &     ( -(rt+obslev)*cos(theta)
    347 c     &     +sqrt( (rt+h)**2 - ((rt+obslev)*sin(theta))**2) )
    348 c--
    349 
    350 c-- Air mass "m" calcualted using a one-exponential density
    351 c-- profile for the atmosphere, rho = rho_0 exp(-h/Hs)
    352 c-- with Hs = hscale = 7.4 km
     133c--   Air mass factor "m" calculated using a one-exponential density
     134c--   profile for the atmosphere, rho = rho_0 exp(-h/hscale) with
     135c--   hscale = 7.4 km. The air mass factor is defined as I(theta)/I(0),
     136c--   the ratio of the optical paths I (in g/cm2) traversed between two
     137c--   given heights, at theta and at 0 deg z.a.
    353138
    354139      Rcos2 = rt * cos(theta)**2
    355140      Rsin2 = rt * sin(theta)**2
    356      
    357       x1 = sqrt((2. * obslev + Rcos2) / (2. * hscale))
    358       x2 = sqrt((2. * h      + Rcos2) / (2. * hscale))
    359       x3 = sqrt((2. * obslev + rt)    / (2. * hscale))
    360       x4 = sqrt((2. * h      + rt)    / (2. * hscale))
    361 
    362 c--   AM Dec 2001, to avoid crash! A few photons seem to be "corrupted"
    363 c--    (have absurd value) in a cer file...
     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
     150c--   AM, Dec 2002: slightly changed the calculation of the air mass factor,
     151c--   for what I think is a better approximation (numerically the results are
     152c--   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
     159c--   AM Dec 2001, to avoid crash! Sometime a few photons seem to be
     160c--   "corrupted" (have absurd values) in cer files... Next lines avoid the
     161c--   crash. They will also return a -1 transmittance in the case the photon
     162c--   comes exactly from the observation level, because in that case the
     163c--   "air mass factor" would become infinity and the calculation is not valid!
    364164
    365165      if (abs(x3-x4) .lt. 1.e-10) then
    … …  
    368168      endif
    369169
    370 
    371170      e1 = derfc(x1)
    372171      e2 = derfc(x2)
    … …  
    376175      m = exp(-Rsin2 / (2. * hscale)) * ((e1 - e2) / (e3 - e4))
    377176
    378 
    379177**********************************************************************   
    380 *         
    381 *     RAYLEIGH SCATTERING
    382 
    383       RHOTOT = 0.0
    384       do 11 i=1,5
    385         RHOF(i) = 0
    386  11   continue
     178*     RAYLEIGH SCATTERING (by the molecules in the atmosphere)
     179*     SCATTERING M.F.P. FOR RAYLEIGH: XR = 2.970D3
     180**********************************************************************
     181
     182c--   Calculate the traversed "vertical thickness" of air using the
     183c--   US Standard atmosphere:
     184
     185      RHO_TOT = 0.0
    387186
    388187      DO 91 I=2,5
    389188        IF ( H .LT. LAHG(I) ) THEN
    390           RHOTOT = RHOTOT +
     189          RHO_TOT = RHO_TOT +
    391190     +        BATM(I-1)*(EXP(-LAHG(I-1)/CATM(I-1) ) -
    392191     +        EXP(-H/CATM(I-1)))
    393192          GOTO 92
    394193        ELSE
    395           RHOTOT = RHOTOT +
     194          RHO_TOT = RHO_TOT +
    396195     +        BATM(I-1)*(EXP(-LAHG(I-1)/CATM(I-1) ) -
    397196     +        EXP(-LAHG(I)/CATM(I-1)))
    … …  
    399198 91   CONTINUE
    400199 
    401 
    402  92   RHO_FI = m*RHOTOT
    403 
    404       T_Ray = EXP(-(RHO_FI/XR)*(400.D0/wavelength)**4)
     200c--   We now convert from "vertical thickness" to "slanted thickness"
     201c--   traversed by the photon on its way to the telescope, simply
     202c--   multiplying by the air mass factor m:
     203
     204 92   RHO_FI = m*RHO_TOT
     205
     206c--   Finally we calculate the transmission coefficient for the Rayleigh
     207c--   scattering:
     208c--   AM Dec 2002, introduced ABS below to account (in the future) for
     209c--   possible photons coming from below the observation level.
     210
     211      T_Ray = EXP(-ABS(RHO_FI/XR)*(400.D0/wavelength)**4)
    405212   
    406 
    407 ***************************************************************************
    408 *                                                                         *
    409 *     MIE ABSORPTION: WE USE FOR HEIGHTS LOWER THAN 10 Km AN EXACT FORMULA*
    410 *     FOR THE AEROSOL DENSITY, AND WE USE A TABLE FOR HIGHTS HIGHER THAN  *
    411 *     10 Km. WE TAKE THE TABLE FROM THE ALTERMAN ARTICLE.                 *
    412 ***************************************************************************
    413      
    414 
    415 ***************************************************************************
    416 *                                                                         *
    417 *     MIE SCATTERING PARAMETERS                                           *
    418       Hm = 1.2D5
    419 ***************************************************************************
    420 
    421       IF (ROW.GT.27) THEN
    422         ROW=28
    423       ENDIF
    424 
    425 
    426         CON_MI = 2
    427 
    428  1001   IF (ABS(LONG(CON_MI)-wavelength).LT.0.01) THEN
    429           TOT_AE =AE_ABI(ROW,CON_MI)
    430         ELSEIF (LONG(CON_MI).GT.wavelength) THEN
    431           A = (AE_ABI(ROW,CON_MI)-AE_ABI(ROW,CON_MI-1))/
    432      +        (LONG(CON_MI) - LONG(CON_MI-1))
    433           B = AE_ABI(ROW,CON_MI) - (A*LONG(CON_MI))
    434           TOT_AE = A*wavelength + B
    435         ELSE
    436           CON_MI = CON_MI +1
    437           GOTO 1001
    438         ENDIF
    439 
    440         T_Mie =  EXP(-(m*TOT_AE))           
    441 
    442 
    443 
    444 ***********************************************************************
    445 *                                                                     *   
    446 *     OZONE ABSORPTION: We used the Elterman table.                   *
    447 *                                                                     *
    448 ***********************************************************************
    449       IF (ROW.GT.47) THEN
    450         ROW = 47
    451       ENDIF
    452 
    453       CON_OZ = 2
    454 
    455  2001 IF (LONG(CON_OZ).GT.wavelength) THEN
    456              A = (OZ_ABI(ROW,CON_OZ)-OZ_ABI(ROW,CON_OZ-1))/
    457      +        (LONG(CON_OZ) - LONG(CON_OZ-1))
    458          B = OZ_ABI(ROW,CON_OZ) - (A*LONG(CON_OZ))
    459          TOT_OZ = (A* wavelength)+B
    460       ELSEIF (ABS(LONG(CON_OZ)-wavelength).LT.0.01) THEN
    461          TOT_OZ = OZ_ABI(ROW,CON_OZ)
    462       ELSE
    463         CON_OZ = CON_OZ + 1
    464         GOTO 2001
    465       ENDIF
    466      
    467       T_Oz = EXP(-(m*TOT_OZ))     
    468 
    469 ************************************************************************
    470 *                                                                      *
    471 *     TOTAL TRANSMISSION OF THE ATMOSPHERE: (RAYLEIGH + MIE + OZONE)   *
    472 ************************************************************************
    473 
    474       tr_atmos = T_Ray*T_Mie*T_Oz
     213      tr_atmos = T_Ray
    475214
    476215      RETURN
    477216
    478217      END
    479 
    480 c @endcode
    481 
    482 c EOF
Note: See TracChangeset for help on using the changeset viewer.