Changeset 1722


Ignore:
Timestamp:
01/21/03 15:02:35 (22 years ago)
Author:
moralejo
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Simulation/Detector/ReflectorII
Files:
8 added
11 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Simulation/Detector/ReflectorII/Changelog

    r1624 r1722  
    11** Add changes at the beginning! **
     2
     321/01/2003 A. Moralejo
     4
     5Removed from the output a null byte written right after the ascii label
     6containing the program version number which is at the beginning of the rfl
     7file.
     8
     919/12/2002 - 17/01/2003 A. Moralejo
     10
     11Lots of changes. Moved simulation of Mie scattering and Ozone absorption from
     12attenu.f to atm.c, and changed the way it is done. Before it was not correct
     13for large zenith angle (a variation like for Rayleigh scattering was assumed
     14for both Mie scattering and Ozone absorption).
     15
     16
     1713/12/2002 A. Moralejo
     18
     19attenu.f: Found mistake in Mie absorption calculation (height from sea
     20level was taken as height above observation level). Changed optical depths
     21table for Mie attenuation. Now they are no longer referred to 2 km height,
     22but to sea level. Detector level is taken into account later in calculation.
     23
     2410/12/2002 A. Moralejo
     25
     26attenu.f: Added comments, removed old/unnecessary code, corrected small
     27mistake in Elterman's table for ozone optical depth.
     28
     29atm.c,h : removed atmospheric model "ATM_ISOTHERMAL" which actually was
     30not even implemented in the code, although it could be selected in the
     31input card.
    232
    33315/11/2002 A. Moralejo
     
    1545
    164614/11/2002 A. Moralejo
    17 reflector.c, parms.c Added wobble mode option. Added wobble flag to the
     47reflector.c, parms.c: Added wobble mode option. Added wobble flag to the
    1848output and also an atmospheric_model flag.
    1949
  • trunk/MagicSoft/Simulation/Detector/ReflectorII/Makefile

    r1614 r1722  
    1919#
    2020# $RCSfile: Makefile,v $
    21 # $Revision: 1.5 $
    22 # $Author: bigongia $
    23 # $Date: 2002-11-14 21:39:01 $
     21# $Revision: 1.6 $
     22# $Author: moralejo $
     23# $Date: 2003-01-21 15:02:32 $
    2424#
    2525##################################################################
     
    217217atm.o: /usr/include/sys/machine/machtime.h /usr/include/sys/rt_limits.h
    218218atm.o: /usr/include/string.h /usr/include/strings.h /usr/include/math.h
    219 atm.o: /usr/include/stdlib.h diag.h atm.h init.h
     219atm.o: /usr/include/stdlib.h atm.h diag.h init.h
    220220ph2cph.o: /usr/include/stdio.h /usr/include/standards.h
    221221ph2cph.o: /usr/include/sys/seek.h /usr/include/va_list.h
     
    229229header.o: /usr/include/sys/types.h /usr/include/mach/machine/vm_types.h
    230230header.o: /usr/include/sys/select.h /usr/include/strings.h header.h
    231 attach.o: /usr/include/string.h /usr/include/standards.h
     231attach.o: /usr/include/stdio.h /usr/include/standards.h
     232attach.o: /usr/include/sys/seek.h /usr/include/va_list.h
    232233attach.o: /usr/include/sys/types.h /usr/include/mach/machine/vm_types.h
    233 attach.o: /usr/include/sys/select.h /usr/include/strings.h
     234attach.o: /usr/include/sys/select.h /usr/include/getopt.h
     235attach.o: /usr/include/sys/limits.h /usr/include/sys/machine/machlimits.h
     236attach.o: /usr/include/sys/syslimits.h /usr/include/sys/machine/machtime.h
     237attach.o: /usr/include/sys/rt_limits.h
    234238reflector.o: /usr/include/stdio.h /usr/include/standards.h
    235239reflector.o: /usr/include/sys/seek.h /usr/include/va_list.h
  • trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c

    r1614 r1722  
     1/********************************************************************
     2*                                                                   *
     3*  File: atm.c                                                      *
     4*  Authors: J.C. Gonzalez, A. Moralejo                              *
     5*                                                                   *
     6* January 2002, A. Moralejo: lots of changes. Moved the code for    *
     7*         the Mie scattering and ozone absorption from attenu.f to  *
     8*         here, after some bugs were found. Now the implementation  *
     9*         is different, we now precalculate the slant paths for the *
     10*         aerosol and Ozone vertical profiles, and then do an       *
     11*         interpolation in wavelength for every photon to get the   *
     12*         optical depths. The parameters used, defined below,       *
     13*         have been taken from "Atmospheric Optics", by L. Elterman *
     14*         and R.B. Toolin, chapter 7 of the "Handbook of geophysics *
     15*         and Space environments". (S.L. Valley, editor).           *
     16*         McGraw-Hill, NY 1965.                                     *
     17*                                                                   *
     18*         WARNING: the Mie scattering and the Ozone absorption are  *
     19*         implemented to work only with photons produced at a       *
     20*         height a.s.l larger than the observation level. So this   *
     21*         is not expected to work well for simulating the telescope *
     22*         pointing at theta > 90 deg (for instance for neutrino     *
     23*         studies. Rayleigh scattering works even for light coming  *
     24*         from below.                                               *
     25*                                                                   *
     26*********************************************************************/
     27
    128#include <stdio.h>
    229#include <string.h>
    330#include <math.h>
     31
     32#include "atm.h"
    433#include "diag.h"
    5 #include "atm.h"
    634#include "init.h"
    735
    836/*  random numbers  */
    937#define RandomNumber  ranf()
    10 
    11 /* AM Nov 2002: atmModel is now an external variable: */
    12 int atmModel=0;   /* current atm. model */
     38#define STEPTHETA 1.74533e-2 /* aprox. 1 degree */
     39
     40#define MIN(x,y) ((x)<(y)? (x) : (y))
    1341
    1442/*  Function declarations  */
    1543static float atm(float wavelength, float height, float theta);
    16 void SetAtmModel(char *model);
     44void SetAtmModel(int model, float ol);
    1745int absorption(float wlen, float height, float theta);
    18 extern void attenu_(float *, float *, float *, float *);  /* in Fortran */
     46extern void attenu_(float *, float *, float *, float *, float *);  /* in Fortran        */
    1947extern float ranf(void);
    2048
    21 void SetAtmModel(char *model)
     49/* aerosol_path contains the path integrals for the aerosol number
     50 * density (relative to the number density at sea level) between the
     51 * observation level and a height h for different zenith angles. The
     52 * first index indicate height above sea level in units of 100m, the
     53 * second is the zenith angle in degrees.
     54 */
     55static float aerosol_path[301][90];
     56
     57/* ozone_path contains the path integrals for the ozone concentration
     58 * between the observation level and a height h for different zenith
     59 * angles. The first index indicate height above sea level in units
     60 * of 100m, the second is the zenith angle in degrees.
     61 */
     62static float ozone_path[501][90];
     63
     64static float obslev; /* observation level in cm */
     65static double rt;    /* Earth radius in cm */
     66static int atmModel;
     67
     68void SetAtmModel(int model, float ol)
    2269{
    23     while (strcmp(model, AtmModelNames[atmModel]))
    24         if (++atmModel == sizeof(AtmModelNames)/sizeof(AtmModelNames[0]))
    25         {   atmModel = 0;
    26             Error(ATM__NFND__ERR, model);
    27             break;      }
    28 
    29     Log(ATM__SET___LOG, AtmModelNames[atmModel]);
    30 
    31     /*  'erfc' must be inited, so keep the following line.  */
    32     Debug("Initing 'erfc': ret=%g\n", erfc(M_PI));
     70  float Rcos2, sin2, rtsq, path_slant, h, dh, theta;
     71  int j;
     72
     73  atmModel = model;
     74  obslev = ol;
     75  rt= 6371315.E2; /* Earth radius (same as in Corsika) in cm */
     76
     77  if (atmModel == ATM_CORSIKA)
     78    {
     79      /* It follows a precalculation of the slant path integrals we need
     80       * for the estimate of the Mie scattering and Ozone absorption:
     81       */
     82
     83      rtsq = sqrt(rt);
     84      dh = 1.e3;
     85
     86      /* Mie (aerosol): */
     87
     88      for (j = 0; j < 90; j++)
     89        {
     90          theta = j * STEPTHETA;  /* aprox. steps of 1 deg */
     91
     92          path_slant = 0;
     93          Rcos2 = rt * cos(theta)*cos(theta);
     94          sin2 = sin(theta)*sin(theta);
     95
     96          for (h = obslev; h <= 30e5; h += dh)
     97            {
     98              if (fmod(h,1e4) == 0)
     99                aerosol_path[(int)(h/1e4)][j] = path_slant;
     100
     101              path_slant +=
     102                (aero_n[(int)floor(h/1.e5)] + (h/1.e5 - floor(h/1.e5))*
     103                  (aero_n[(int)ceil(h/1.e5)]-aero_n[(int)floor(h/1.e5)]))
     104                  /aero_n[0] * dh * (rt+h) /
     105                    sqrt((rt+h)*(rt+h)-(rt+obslev)*(rt+obslev)*sin2);
     106
     107            }
     108        }
     109
     110      /* Ozone absorption */
     111
     112      for (j = 0; j < 90; j++)
     113        {
     114          theta = j * STEPTHETA;  /* aprox. steps of 1 deg */
     115          path_slant = 0;
     116          Rcos2 = rt * cos(theta)*cos(theta);
     117          sin2 = sin(theta)*sin(theta);
     118
     119          for (h = obslev; h <= 50e5; h += dh)
     120            {
     121              if (fmod(h,1e4) == 0)
     122                ozone_path[(int)(h/1e4)][j] = path_slant;
     123
     124              path_slant +=
     125                (oz_conc[(int)floor(h/1.e5)] + (h/1.e5 - floor(h/1.e5))*
     126                  (oz_conc[(int)ceil(h/1.e5)]-oz_conc[(int)floor(h/1.e5)]))
     127                  * dh * (rt+h) /
     128                    sqrt((rt+h)*(rt+h)-(rt+obslev)*(rt+obslev)*sin2);
     129            }
     130        }
     131
     132    }
     133
    33134}   /*  end of SetAtmModel  */
    34            
     135
    35136static float atm(float wavelength, float height, float theta)
    36 {   float transmittance = 1.0;  /* final atm transmittance (ret. value) */
    37 
    38     switch(atmModel)
    39     {   case ATM_NOATMOSPHERE:  /* no atm at all: transmittance = 100%  */
    40              break;
    41         case ATM_90PERCENT:     /* atm. with transmittance = 90%        */
    42              transmittance = 0.9;
    43              break;
    44         case ATM_ISOTHERMAL:    /* isothermal approximation             */
    45         /********************/
    46              break;
    47         case ATM_CORSIKA:       /* atmosphere as defined in CORSIKA     */
    48              attenu_(&wavelength, &height, &theta, &transmittance);
    49              break;
     137{   
     138  float transmittance = 1.0;  /* final atm transmittance (ret. value) */
     139  float T_Ray, T_Mie, T_Oz;
     140
     141  float h; /* True height a.s.l. of the photon emission point in cm */
     142  float tdist;
     143  float beta0, path;
     144
     145  int index;
     146
     147  switch(atmModel)
     148    {   
     149    case ATM_NOATMOSPHERE:      /* no atm at all: transmittance = 100%  */
     150      break;
     151    case ATM_90PERCENT: /* atm. with transmittance = 90%        */
     152      transmittance = 0.9;
     153      break;
     154    case ATM_CORSIKA:   /* atmosphere as defined in CORSIKA     */
     155
     156      /* Distance to telescope: */
     157      tdist = (height-obslev)/cos(theta);
     158
     159      /* Avoid problems if photon is very close to telescope: */
     160      if (fabs(tdist) < 1.)
     161        {
     162          transmittance = 1.;
     163          break;
     164        }
     165
     166      /*** We calculate h, the true emission height above sea level: ***/
     167
     168      h = -rt + sqrt((rt+obslev)*(rt+obslev) + tdist*tdist +
     169                     (2*(rt+obslev)*(height-obslev)));
     170
     171      /******* Rayleigh scattering: *******/
     172
     173      attenu_(&wavelength, &h, &obslev, &theta, &T_Ray);
     174
     175
     176      /******* Ozone absorption: *******/
     177
     178      if (h > 50.e5)
     179        h = 50.e5;
     180
     181      /* First we get Vigroux Ozone absorption coefficient for the given
     182       * wavelength, through a linear interpolation:
     183       */
     184
     185      for (index = 1; index < 11; index++)
     186        if (wavelength < wl[index])
     187          break;
     188
     189      beta0 = oz_vigroux[index-1]+(oz_vigroux[index]-oz_vigroux[index-1])*
     190        (wavelength-wl[index-1])/(wl[index]-wl[index-1]);
     191
     192      /* from km^-1 to cm^-1 : */
     193      beta0 *= 1e-5;
     194     
     195      /* Now use the pre-calculated values of the path integral
     196       * for h and theta: */
     197
     198      path = ozone_path[(int)floor(0.5+h/1e4)]
     199        [(int)MIN(89,floor(0.5+theta/STEPTHETA))];
     200
     201      T_Oz = exp(-beta0*path);
     202
     203
     204      /******* Mie (aerosol): *******/
     205
     206      if (h > 30.e5)
     207        h = 30.e5;
     208
     209      /* First get Mie absorption coefficient at sea level for the given
     210       * wavelength, through a linear interpolation:
     211       */
     212
     213      for (index = 1; index < 11; index++)
     214        if (wavelength < wl[index])
     215          break;
     216
     217      beta0 = aero_betap[index-1]+(aero_betap[index]-aero_betap[index-1])*
     218        (wavelength-wl[index-1])/(wl[index]-wl[index-1]);
     219
     220      /* from km^-1 to cm^-1 : */
     221      beta0 *= 1e-5;
     222
     223      /* Now use the pre-calculated values of the path integral
     224       * for h and theta: */
     225
     226      path = aerosol_path[(int)floor(0.5+h/1e4)]
     227        [(int)MIN(89,floor(0.5+theta/STEPTHETA))];
     228
     229      T_Mie = exp(-beta0*path);
     230
     231
     232      /* Calculate final transmission coefficient: */
     233
     234      transmittance = T_Ray * T_Oz * T_Mie;
     235
     236      break;
     237
    50238    }   /*  end of atm switch   */
    51239
    52     return transmittance;
     240
     241  return transmittance;
     242
    53243}   /*  end of atm  */
    54244
    55245int absorption(float wlen, float height, float theta)
    56 {   int ret = 0;        /*  0: passed, 1: absorbed  */
     246{
     247  int ret = 0;  /*  0: passed, 1: absorbed  */
    57248   
    58     if (RandomNumber > atm(wlen, height, theta)) ret=1;
    59 
    60     return ret;
     249  if (RandomNumber > atm(wlen, height, theta)) ret=1;
     250
     251  return ret;
    61252}   /*  end of absorption  */
     253
     254
     255
  • trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.h

    r725 r1722  
    55T(ATM_NOATMOSPHERE),    /* no atmosphere at all: transmittance = 100% */ \
    66T(ATM_90PERCENT),       /* atmosphere with transmittance = 90% */        \
    7 T(ATM_ISOTHERMAL),      /* isothermal approximation */                   \
    87T(ATM_CORSIKA)          /* atmosphere as defined in CORSIKA */
    98 
     
    1211#undef T
    1312
    14 #define T(x)  #x        /* define T() as the string of x  */
    15     const char *AtmModelNames[] = { ITEM_LIST };
    16 #undef T
    1713#undef ITEM_LIST
    1814
    19 /*  Strings  */
    20 #define ATM__NFND__ERR          /*  model name          */ \
    21     " *** Atm model \"%s\" not found.\n"
    22 #define ATM__SET___LOG          /*  model name          */ \
    23     "Setting atm model \"%s\".\n"
     15/* A. Moralejo, January 2003: added some parameters for Mie scattering and
     16 * Ozone absorption derived from the clear standard atmosphere model in
     17 * "Atmospheric Optics", by L. Elterman and R.B. Toolin, chapter 7 of
     18 * the "Handbook of geophysics and Space environments". S.L. Valley,
     19 * editor. McGraw-Hill, NY 1965.
     20 */
     21
     22/* Aerosol number density for 31 heights a.s.l., from 0 to 30 km,
     23 * in 1 km steps (units: cm^-3)
     24 */
     25
     26float aero_n[31] = {2.0e2, 8.7e1, 3.8e1, 1.6e1, 7.2, 3.1, 1.1, 4.0e-1, 1.4e-1, 5.0e-2, 2.6e-2, 2.3e-2, 2.1e-2, 2.3e-2, 2.5e-2, 4.1e-2, 6.7e-2, 7.3e-2, 8.0e-2, 9.0e-2, 8.6e-2, 8.2e-2, 8.0e-2, 7.6e-2, 5.2e-2, 3.6e-2, 2.5e-2, 2.4e-2, 2.2e-2, 2.0e-2, 1.9e-2};
     27
     28/* Aerosol scattering coefficient at sea level for the following wavelengths:
     29 * 280, 300, 320, 340, 360, 380, 400, 450, 500, 550, 600 and 650 nm
     30 * (units: km^-1)
     31 */
     32
     33float aero_betap[12] = {0.27, 0.26, 0.25, 0.24, 0.24, 0.23, 0.20, 0.180, 0.167, 0.158, 0.150, 0.142};
     34
     35float wl[12] = {280, 300, 320, 340, 360, 380, 400, 450, 500, 550, 600, 650};
     36
     37/* Ozone concentration for 51 heights a.s.l., from 0 to 50 km,
     38 * in 1 km steps (units: cm/km)
     39 */
     40
     41float oz_conc[51]={0.3556603E-02, 0.3264150E-02, 0.2933961E-02, 0.2499999E-02, 0.2264150E-02, 0.2207546E-02, 0.2160377E-02, 0.2226414E-02, 0.2283018E-02, 0.2811320E-02, 0.3499999E-02, 0.4603772E-02, 0.6207545E-02, 0.8452828E-02, 0.9528299E-02, 0.9905657E-02, 0.1028302E-01, 0.1113207E-01, 0.1216981E-01, 0.1424528E-01, 0.1641509E-01, 0.1839622E-01, 0.1971697E-01, 0.1981131E-01, 0.1933962E-01, 0.1801886E-01, 0.1632075E-01, 0.1405660E-01, 0.1226415E-01, 0.1066037E-01, 0.9028300E-02, 0.7933960E-02, 0.6830187E-02, 0.5820753E-02, 0.4830188E-02, 0.4311319E-02, 0.3613206E-02, 0.3018867E-02, 0.2528301E-02, 0.2169811E-02, 0.1858490E-02, 0.1518867E-02, 0.1188679E-02, 0.9301884E-03, 0.7443394E-03, 0.5764149E-03, 0.4462263E-03, 0.3528301E-03, 0.2792452E-03, 0.2226415E-03, 0.1858490E-03};
     42
     43/* Vigroux ozone absorption coefficients (cm-1) */
     44
     45float oz_vigroux[12]= {1.06e2, 1.01e1, 8.98e-1, 6.40e-2, 1.80e-3, 0, 0, 3.50e-3, 3.45e-2, 9.20e-2, 1.32e-1, 6.20e-2};
     46
    2447
    2548#endif
  • 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
  • trunk/MagicSoft/Simulation/Detector/ReflectorII/doc/Tdas0211.tex

    r1672 r1722  
    3333
    3434\usepackage{magic-tdas}
     35\usepackage{amssymb}
    3536
    3637%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
    4647\author{A.Moralejo\\
    4748        \texttt{<moralejo@pd.infn.it>}}
    48 \date{December 6, 2002\\}
    49 \TDAScode{MAGIC-TDAS 02-11\\ 021206/AMoralejo}
     49\date{January 20, 2003\\}
     50\TDAScode{MAGIC-TDAS 02-11\\ 030120/AMoralejo}
    5051%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    5152
     
    6162been included here for clarity. Two important bugs regarding the
    6263ray-tracing routine have been corrected in this new version.
     64\par
     65NOTE: In December 2002, a first release of Reflector 0.6 was made,
     66together with a first version of the present TDAS note. Immediately
     67after that (too short a time to justify a new version number), some
     68other changes were implemented in the program to improve the
     69atmospheric absorption routines, and the dependence of mirror
     70reflectivity with wavelength was introduced. This document is the
     71manual of the final 0.6 version of the reflector, including
     72explanations of these last modifications.
     73
    6374\end{abstract}
    6475
     76\newpage
    6577%% contents %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    6678\thetableofcontents
     
    161173see the resulting spot in fig. \ref{spot_inf_f1697}. A completely
    162174independent ray-tracing program was used to verify that this is the
    163 spot that our 17 m tesellated paraboloid should produce.
     175spot that a f/1 16.97 m $\varnothing$ tesellated paraboloid should produce.
    164176%
    165177\begin{figure}[h]
     
    182194possible mirror misalignments and surface irregularities (by the way,
    183195a feature which somehow optimistically was not included in the
    184 simulations shown in the proposal). In this way we can check just the
    185 ray tracing.
     196simulations shown in the MAGIC proposal). In this way we can check
     197just the ray tracing.
    186198%
    187199\begin{figure}[h!]
     
    242254intended, and looks more like an error, but anyway we checked that this
    243255was not the reason for the problems in the ray tracing: a test
    244 magic.def file was created with all parameters calculated as for a f =
    245 1697 cm paraboloid, and no significant difference could be seen. That
     256magic.def file was created with {\it all} the mirror parameters
     257(positions, orientations and radii of curvature)
     258calculated as for a f = 1697 cm paraboloid, and no significant
     259difference could be seen. That
    246260is: i) the individual mirror orientations are the dominant factor, and
    247261the overall dish in the old reflector behaved indeed like a f = 1697
     
    288302determine the particle's trajectory as long as one knows that they
    289303refer to a downgoing versor, in which case one gets the third
    290 direction cosine as $w = -\sqrt{u^2+v^2}$, with a minus sign. However,
     304direction cosine as $w = -\sqrt{1-u^2-v^2}$, with a minus sign. However,
    291305in {\it ph2cph.c} we found exactly the opposite (lines 116 to 118 in
    292306v.05):
     
    341355The general transformation between both is a simple rotation,
    342356since for the sake of simplicity we assume in the simulation that the
    343 origins always coincide. In Reflector v.05 or older the rotation
     357origins always coincide. As we have said, in Reflector v.05 or older the rotation
    344358matrix was wrong: it had been written assuming that ($\Phi$, $\Theta$)
    345359indicated the direction towards which the telescope pointed. Actually,
     
    383397The bug was certainly present in versions 0.4 and 0.5, but may be even
    384398older. Nevertheless, there is no doubt that the reflector program used
    385 for the simulation shown in the MAGIC design report (which must have been
    386 an early version of the present one) was working fine. Extensive proof
     399for the simulation shown in the MAGIC design report was working
     400fine. Extensive proof
    387401of this is provided in an appendix of the design report. A plausible
    388402explanation could be that, up to some date, the data being read in by
     
    395409file shows no record of any change in this respect, but given that we
    396410have always used a slightly modified Corsika, it would not be
    397 surprising if the Cherenkov output was modified at some point. There
    398 is no documentation on this, so if anyone has any relevant information,
    399 please make it public.
     411surprising if the Cherenkov output was modified at some point in the
     412MAGIC version of Corsika (MMCS). There is no documentation on this, so
     413if anyone has any relevant information, please make it public.
    400414%
    401415\paragraph {Influence of the bug on image analysis\\}
     
    413427dramatic effect than the defocusing which the bug was producing.
    414428%
     429\subsection {Performance of ray-tracing in the new version}
     430In figure \ref{coma} we show the images of a point-like source at 10
     431km from the telescope, produced with the Reflector version 0.6,
     432and using the new version of the magic.def file (see see next
     433section). No noise has been introduced in the reflection, the observed
     434spots are just the result of optical aberrations. The light source has
     435been put at slightly different viewing angles
     436from the telescope. The results are comparable to those in the design
     437report, actually these are a bit better, the difference probably being
     438that the focal lengths of the mirror tiles in the older magic.def file
     439were discretized in only eight values, while now they change rather
     440continuously. Some images of a point source at infinity (a star) can be
     441seen in fig. \ref{coma_star}. We can see that for any incidence angle,
     442the area within which 50$\%$ of the light is concentrated is smaller
     443than that of a small pixel.
     444\par
     445In figure \ref{refl06images} the images of three gamma events ($\theta
     446= 10^\circ$, E = 16, 46, 232 GeV), the same of fig. \ref{evcompare}
     447are shown. They have been produced with Reflector 0.6 assuming perfect
     448spherical mirrors (left) and realistic ones (right). The images look
     449reasonable, much sharper than with the older versions, even when the
     450mirror imperfections are taken into account.
     451%
     452\begin{figure}[h!]
     453  \begin{center}
     454    \epsfig{file=eps/timing.eps,width=\textwidth}
     455    \caption{Test of reflector isochrony. The arrival time
     456distributions of photons in the camera are shown for (buggy) Reflector
     4570.5 and for Reflector 0.6. The sketch in the center shows the test for
     458the case in which the light beam is paralel to the telescope axis
     459(left plot). On the right, the same test has been made with light
     460arriving 1 degree off axis.
     461 \label{timing}}
     462  \end{center}
     463\vspace*{-1cm}
     464\end{figure}
     465%
    415466\subsection {Second bug: photon timing}
     467%
    416468After a first release of Reflector 0.6 had been made public, another
    417469important bug was found. We tried to check whether the simulated
     
    438490from scratch!).
    439491%
    440 \begin{figure}[h]
     492\subsection{The new magic.def file}
     493A new magic.def file (see sect. \ref{neededfiles}) has been created
     494and included in the Reflector 0.6 package. Now the number of
     495individual mirror tiles is 956, matching
     496the number and distribution of the final MAGIC design. The mirror
     497centers and orientations are those corresponding to a paraboloid of
     4981697 cm focal (hence the camera plane is placed at 1700 cm from the
     499dish). The focal lengths have been calculated by R. Mirzoyan taking
     500into account the so called ``shortening effect'' (see design report).
     501A new axisdev.dat file (se again \ref{neededfiles}) with data for the
     502956 mirrors is also included.
     503%
     504\subsection{The new reflectivity.dat file}
     505Up to version 0.5 of the program, the reflectivity of the mirrors
     506(which the program reads in from the file reflectivity.dat) was
     507considered to be equal to $90\%$ for all wavelengths. Following
     508measurements performed in Padua of a mirror sample, the reflectivity
     509has been found to be between 80 and 90$\%$ in the range from 280 to 650 nm,
     510with a dependence on wavelength which is shown in figure \ref{reflec}.
     511%
     512\begin{figure}[h!]
    441513  \begin{center}
    442     \epsfig{file=eps/timing.eps,width=\textwidth}
    443     \caption{Test of reflector isochrony. The arrival time
    444 distributions of photons in the camera are shown for (buggy) Reflector
    445 0.5 and for Reflector 0.6. The sketch in the center shows the test for
    446 the case in which the light beam is paralel to the telescope axis
    447 (left plot). On the right, the same test has been made with light
    448 arriving 1 degree off axis.
    449  \label{timing}}
     514    \epsfig{file=eps/reflec.eps,width=0.7\textwidth,height=0.4\textwidth}
     515    \caption{Reflectivity of the MAGIC telescope mirrors as a
     516function of the wavelength of the incident light. The dip around 350
     517nm is due to destructive interference of the light reflected on the
     518aluminum plate with that reflected on the protective coating of the
     519mirror. \label{reflec}}
    450520  \end{center}
    451 \vspace*{-1cm}
    452 \end{figure}
    453 
    454 \subsection {Performance of the new version}
    455 In figure \ref{coma} we show the images of a point-like source at 10
    456 km from the telescope, produced with the Reflector version 0.6,
    457 and using the new version of the magic.def file (see see next
    458 section). No noise has been introduced in the reflection, the observed
    459 spots are just the result of optical aberrations. The light source has
    460 been put at slightly different viewing angles
    461 from the telescope. The results are comparable to those in the design
    462 report, actually these are a bit better, the difference probably being
    463 that the focal lengths of the mirror tiles in the older magic.def file
    464 were discretized in only eight values, while now they change rather
    465 continuously. Some images of a point source at infinity (a star) can be
    466 seen in fig. \ref{coma_star}. We can see that for any incidence angle,
    467 the area within which 50$\%$ of the light is concentrated is smaller
    468 than that of a small pixel.
    469 \par
    470 In figure \ref{refl06images} the images of three gamma events ($\theta
    471 = 10^\circ$, E = 16, 46, 232 GeV), the same of fig. \ref{evcompare}
    472 are shown. They have been produced with Reflector 0.6 assuming perfect
    473 spherical mirrors (left) and realistic ones (right). The images look
    474 reasonable, much sharper than with the older versions, even when the
    475 mirror imperfections are taken into account.
    476 \par
     521\end{figure}
    477522%
    478523\begin{figure}[p]
     
    515560\end{figure}
    516561%
    517 \subsection{The new magic.def file}
    518 A new magic.def file (see sect. \ref{neededfiles}) has been created
    519 and included in the Reflector 0.6 package. Now the number of
    520 individual mirror tiles is 956, matching
    521 the number and distribution of the final MAGIC design. The mirror
    522 centers and orientations are those corresponding to a paraboloid of
    523 1697 cm focal (hence the camera plane is placed at 1700 cm from the
    524 dish). The focal lengths have been calculated by R. Mirzoyan taking
    525 into account the so called ``shortening effect'' (see design report).
    526 A new axisdev.dat file (se again \ref{neededfiles}) with data for the
    527 956 mirrors is also included.
    528 %
    529562\subsection{The {\itshape cermaker} program}
    530563A test program to produce cer files (input for the reflector)
     
    539572file is called {\it cer000001}, and can be read by the reflector program.
    540573%
     574\subsection{Changes in the atmospheric absorption routines}
     575%
     576In the present version of the program, the simulation of atmospheric
     577absorption has undergone major changes with respect to the original
     578implementation by J.C. Gonz\'alez and Aitor Ibarra. Three options are
     579available in the simulation: ATM\_NOATMOSPHERE, ATM\_90PERCENT and
     580ATM\_CORSIKA (see section \ref{commands}). In the
     581first two, 100$\%$ and 90$\%$ of the emitted light respectively
     582reaches the telescope mirror, regardless of the emission height. The
     583changes in the present version affect only the third option, which is
     584the recommended one for running reflector in the standard MAGIC MC
     585production. In this case, a detailed estimate of the atmospheric
     586transmission is done.
     587\par
     588Three contributions to the atmospheric opacity
     589are considered: Rayleigh scattering, Mie scattering by aerosols, and
     590absorption by ozone. Details on how the effect of these three
     591components is calculated are given in appendix A. In reflector
     592versions older than 0.6,
     593due to a bug, the contribution of Mie scattering and Ozone absorption
     594was slightly overestimated because the vertical height of the emission
     595point above sea level was interpreted as height above the
     596telescope. The difference is small, since the density of the aerosols
     597responsible for Mie scattering decreases very fast with height and
     598therefore the absorption is hardly increased by going up 2.2 km in the
     599region were most of the Cherenkov light is produced. In the case of
     600ozone, its contribution is only important in the very low end of the
     601Cherenkov light spectrum, and so the effect of the bug in the total
     602amount of light reaching the telescope is negligible. Another change
     603is that now the observation level is read in from the cer file, instead
     604of being a fixed parameter in the routines. In this way the reflector
     605program has become more flexible, and can now be used to process cer
     606files produced for a detector at a different observation level.
     607\par
     608Another problem in the old implementation of absorption was that the
     609variation of the optical depth of the atmosphere with the zenith angle
     610$\theta$ was assumed to be the same for Mie scattering, ozone
     611absorption and Rayleigh scattering: the so called {\it air
     612mass} (see appendix A for details) was therefore calculated only once,
     613using the overall density profile of air molecules (which does not
     614match that of aerosols or ozone), and used to account for the
     615variation of all three effects with $\theta$, while rigorously
     616speaking it is only valid for Rayleigh scattering. Now, the simulation
     617of Mie scattering and ozone absorption has been moved from {\it
     618attenu.f} to {\it atm.c} and is done in a more accurate way which
     619takes into account the vertical profile of aerosols and ozone (see
     620appendix A for details).
     621%
    541622\subsection{Other changes in Reflector 0.6}
    542623
     
    557638section \ref{opt}).
    558639
    559 \item New output format (see sect. \ref{out}): added a Run header,
    560 which is like that of Corsika, plus a couple of variables concerning
    561 the reflector parameters: the wobble mode and the atmospheric model
     640\item New output format (see sect. \ref{out}): we have removed a null
     641byte that was written immediately after the ascii label containing the
     642reflector program version number at the begining of the file. This
     643byte was there for historical reasons and had no function
     644whatsoever. Then we have added a Run header, which is like that of
     645Corsika, plus a couple of variables concerning the reflector
     646parameters: the wobble mode and the atmospheric model
    562647used for the simulation. The event header has also been changed to
    563648include all the information present in the Corsika event
     
    634719\end{itemize}
    635720
    636 All these files are usually in the {\bf
    637 MagicProgs/Simulation/Detector/Data/} directory and {\it in principle} you
    638 should {\bf not} make any change in them to run the program.
     721All these files (included in the reflector package) are usually kept
     722in the {\bf MagicProgs/Simulation/Detector/Data/} directory and {\it
     723in principle} you should {\bf not} make any change in them to run the
     724program.
    639725
    640726%------------------------------------------------------------
     
    676762        Wobble-observation mode (see TDAS 01-05 by W. Wittek).
    677763
    678 \item[ct\_file ../Data/magic.def]
     764\item[ct\_file /path/magic.def]
    679765
    680766        The ct\_file statement defines where the program can find the
    681         telescope characteristics. The path in the example above is
    682         correct to run reflector in
    683         MagicProgs/Simu\-la\-tion/De\-tector/Reflector\_0.6/ directory.
    684         If you want to run it in a different directory you have to modify the
    685         path accordingly.
     767        telescope characteristics. Usually, the magic.def file is kept in
     768        MagicProgs/Simulation/Detector/Data
    686769
    687770\item[atm\_model ATM\_CORSIKA]
    688         The atm\_model statement says to the program what kind of
     771
     772        The atm\_model statement tells the program what kind of
    689773        atmospheric absorption model to use. Possible choices are:
    690         ATM\_CORSIKA, ATM\_ISO\-THERMAL, ATM\_90\-PER\-CENT and
    691         ATM\_NO\-ATMO\-SPHE\-RE.
    692        
     774        ATM\_NOATMOSPHERE,\\ ATM\_90PERCENT and ATM\_CORSIKA,
     775        corresponding respectively to no absorption, a 10$\%$
     776        absorption and a model using the US Standard atmosphere (see
     777        Corsika manual, appendix C) for the Rayleigh scattering and a
     778        model by L. Elterman \cite{elterman64,elterman65} for the Mie
     779        scattering and ozone absorption (see appendix A). The third
     780        model should be chosen for the standard MC MAGIC production.
     781
    693782\item[cer\_files]
    694783
     
    755844        File containing mirror reflectivity as a function of
    756845        wavelength (see section \ref{neededfiles}). If this option is
    757         not supplied, the program will look for
     846        not supplied, the program will look for \\
    758847        ``../Data/reflectivity.dat'' as previous versions of
    759848        Reflector did.
     
    796885%------------------------------------------------------------
    797886
    798 \newpage
    799887\section{Output file \label{out}}
    800 The output file begins with two ascii lines:\\
     888The reflector output file begins with two ascii lines, the first of
     889which informs us of the program version with which it has been
     890produced (NOTE: in the following, the dollar symbol \verb"$" stands
     891for a carriage return):\\
    801892\\
    802 \verb"reflector 0.6" \\
    803 \verb"START---RUN" \\
    804 After the \verb"START---RUN" flag there is a carriage return, and then
    805 the run header which is basically the one from Corsika with two added
    806 variables, {\it wobble\_mode} and {\it atmospheric\_model}. Check the
    807 Corsika manual for the meaning and units of the rest of them. All of the
    808 variables 4-byte real numbers except the first, which is a 4 character
    809 string containing the run header ascii label from Corsika:
     893\verb"reflector 0.6$START---RUN$" \\
     894\\
     895Then there is run header which is basically the one from Corsika with
     896two added variables, {\it wobble\_mode} and {\it
     897atmospheric\_model}. Check the Corsika manual for the meaning and
     898units of the rest of them. All of the variables are 4-byte real numbers
     899except the first, which is a 4 character string containing the run
     900header ascii label from Corsika:
    810901\vspace*{0.5cm}
    811902\\
    812903%
    813 \begin{tabular}{ll}
    814 \parbox{5cm}{Variable} & Description \\
     904\begin{tabular}{lll}
     905\multicolumn{2}{c}{Variable} & Description \\
    815906\hline
    816 & \\
    817 ASCII Label & 'RUNH' \\
    818 RunNumber & \\
    819 date & \\
    820 Corsika\_version & \\
    821 NumObsLev & \\
    822 HeightLev[10] & \\
    823 SlopeSpec & \\
    824 ELowLim & \\
    825 EUppLim & \\
    826 EGS4\_flag & \\
    827 NKG\_flag & \\
    828 Ecutoffh & \\
    829 Ecutoffe & \\
    830 Ecutoffg & \\
    831 C[50] & \\
    832 wobble\_mode & Wobble mode with which the reflector was run (TDAS
     907&& \\
     9081 & ASCII Label & 'RUNH' \\
     9092 & RunNumber & \\
     9103 & date & \\
     9114 & Corsika\_version & \\
     9125 & NumObsLev & \\
     9136 to 15 & HeightLev[10] & \\
     91416 & SlopeSpec & \\
     91517 & ELowLim & \\
     91618 & EUppLim & \\
     91719 & EGS4\_flag & \\
     91820 & NKG\_flag & \\
     91921 & Ecutoffh & \\
     92022 & Ecutoffm & \\
     92123 & Ecutoffe & \\
     92224 & Ecutoffg & \\
     92325 to 74 & C[50] & \\
     92475 & wobble\_mode & Wobble mode with which the reflector was run (TDAS
    83392501-05) \\
    834 atmospheric\_model & Atmospheric model used for the absorption
     92676 & atmospheric\_model & Atmospheric model used for the absorption
    835927simulation \\
    836 & 0 = no atmosphere;  1 = atm\_90percent;  \\
    837 & 2 = atm\_isothermal;  3 = atm\_corsika. \\
    838 dummy1[18] & not used \\
    839 CKA[40] & \\
    840 CETA[5] & \\
    841 CSTRBA[11] & \\
    842 dummy2[104] & not used \\
    843 AATM[5] & \\
    844 BATM[5] & \\
    845 CATM[5] & \\
    846 NFL[4] & \\
     928&& 0 = no atmosphere;  1 = atm\_90percent;  \\
     929&& 2 = atm\_corsika. \\
     93077 to 94 & dummy1[18] & not used \\
     93195 to 134 & CKA[40] & \\
     932135 to 139 & CETA[5] & \\
    847933&\\
    848934\hline
     
    851937\newpage
    852938
    853 Then there comes a ``\verb"START-EVENT"'' flag, followed by a carriage
    854 return and then the binary event header. Each variable is a 4-byte
     939\begin{tabular}{lll}
     940\multicolumn{2}{c}{Variable} & \parbox{11cm}{Description} \\
     941\hline
     942&& \\
     943140 to 140 & CSTRBA[11] & \\
     944151 to 254 & dummy2[104] & not used \\
     945255 to 259 & AATM[5] & \\
     946260 to 264 & BATM[5] & \\
     947265 to 269 & CATM[5] & \\
     948270 to 273 & NFL[4] & \\
     949&\\
     950\hline
     951\end{tabular}
     952\vspace*{0.5cm}
     953\\
     954%
     955Then there comes a carriage return followed by the ascii flag which
     956indicates the start of an event, and again a carriage return:\\
     957\\
     958\verb"$START-EVENT$"\\
     959\\
     960and then the binary event header. Each variable is a 4-byte
    855961float number except for the first one which is the event header label
    856962from Corsika (a string of 4 characters). Some of of the variables from
     
    858964\vspace*{0.5cm}
    859965\\
    860 \begin{tabular}{ll}
    861 \parbox{5cm}{Variable} & Description \\
     966\begin{tabular}{lll}
     967\multicolumn{2}{c}{Variable} & Description \\
    862968\hline
    863 & \\
    864 ASCII label     & 'EVTH' \\
    865 EvtNumber       & Event Number \\
    866 PrimaryID       & Primary particle identification code \\
    867 Etotal          & Primary particle total energy (GeV) \\
    868 Thick0          & CORSIKA's starting altitude in g/cm2 \\
    869 FirstTarget     & CORSIKA's number of first target if fixed \\
    870 zFirstInt       & Height of first interaction in cm \\
    871 p[3]            & Primary particle momentum in x,y,-z directions (GeV) \\
    872 Theta           & Primary particle zenith angle (rad) \\
    873 Phi             & Primary particle azimuth angle (rad) \\       
    874 
    875 NumRndSeq       & Number of different CORSIKA random sequences (max. 10) \\
    876 RndData[10][3]  & RndData[i][0]: integer seed of sequence i \\
    877                 & RndData[i][1]: number of offset random calls (mod $10^6$) of sequence i. \\
    878                 & RndData[i][2]: number of offset random calls ($/10^6$) of sequence i. \\
    879 
    880 RunNumber       & Run number \\
    881 DateRun         & Date of run yymmdd \\
    882 Corsika\_version & Version of {\it CORSIKA} \\
    883 
    884 NumObsLev       & Number of observation levels (should be always 1 for
     969&& \\
     9701 & ASCII label     & 'EVTH' \\
     9712 & EvtNumber       & Event Number \\
     9723 & PrimaryID   & Primary particle identification code \\
     9734 & Etotal              & Primary particle total energy (GeV) \\
     9745 & Thick0              & CORSIKA's starting altitude in g/cm2 \\
     9756 & FirstTarget         & CORSIKA's number of first target if fixed \\
     9767 & zFirstInt   & Height of first interaction in cm \\
     9778 to 10 & p[3]          & Primary particle momentum in x,y,-z directions (GeV) \\
     97811 & Theta              & Primary particle zenith angle (rad) \\
     97912 & Phi                & Primary particle azimuth angle (rad) \\       
     980
     98113 & NumRndSeq  & Number of different CORSIKA random sequences (max. 10) \\
     98214 to 43 & RndData[10][3]       & RndData[i][0]: integer seed of sequence i \\
     983&               & RndData[i][1]: number of offset random calls (mod $10^6$) of sequence i. \\
     984&               & RndData[i][2]: number of offset random calls ($/10^6$) of sequence i. \\
     985
     98644 & RunNumber  & Run number \\
     98745 & DateRun    & Date of run yymmdd \\
     98846 & Corsika\_version & Version of {\it CORSIKA} \\
     989
     99047 & NumObsLev  & Number of observation levels (should be always 1 for
    885991us) \\
    886 HeightLev[10]   & Height of observation levels in cm \\
    887 
    888 SlopeSpec       & Energy spectrum slope \\
    889 ELowLim         & Energy lower limit (GeV) \\
    890 EUppLim         & Energy upper limit (GeV) \\
    891 Ecutoffh        & \\
    892 Ecutoffm        & \\
    893 Ecutoffe        & \\
    894 Ecutoffg        & \\
    895 NFLAIN          & \\
    896 NFLDIF          & \\
    897 NFLPI0          & \\
    898 NFLPIF          & \\
    899 NFLCHE          & \\
    900 NFRAGM          & \\
    901 Bx              & \\
    902 By              & \\
    903 EGS4yn          & \\
    904 NKGyn           & \\
    905 GHEISHAyn       & \\
    906 VENUSyn         & \\
     99248 to 57 & HeightLev[10]        & Height of observation levels in cm \\
     993
     99458 & SlopeSpec  & Energy spectrum slope \\
     99559 & ELowLim    & Energy lower limit (GeV) \\
     99660 & EUppLim    & Energy upper limit (GeV) \\
     99761 & Ecutoffh        & \\
     99862 & Ecutoffm        & \\
     99963 & Ecutoffe        & \\
    9071000& \\
    9081001\hline
     
    9121005%
    9131006
    914 \begin{tabular}{ll}
    915 
    916 \parbox{5cm}{Variable} & Description \\
     1007\begin{tabular}{lll}
     1008
     1009\multicolumn{2}{c}{Variable} & \parbox{11cm}{Description} \\
    9171010\hline
    918 & \\
    919 CERENKOVyn      & \\
    920 NEUTRINOyn      & \\
    921 HORIZONTyn      & \\
    922 COMPUTER        & \\
    923 
    924 ThetaMin        & Minimum Theta of primaries (deg) \\
    925 ThetaMax        & Maximum Theta of primaries (deg) \\
    926 PhiMin          & Minimum Phi of primaries (deg) \\
    927 PhiMax          & Maximum Phi of primaries (deg) \\
    928 CBunchSize      & \\
    929 CDetInX         & \\
    930 CDetInY         & \\
    931 CSpacInX        & \\
    932 CSpacInY        & \\
    933 CLenInX         & \\
    934 CLenInY         & \\
    935 COutput         & \\
    936 AngleNorthX     & \\
    937 MuonInfo        & \\
    938 StepLength      & \\
    939 CWaveLower      & Wavelength lower limit (nm) \\
    940 CWaveUpper      & Wavelength upper limit (nm) \\
    941 Multipl         & \\
    942 CorePos[2][20]  & Core positions of randomized shower \\
    943 SIBYLL[2]       & \\
    944 QGSJET[2]       & \\
    945 DPMJET[2]       & \\
    946 VENUS\_cross    & \\
    947 mu\_mult\_scat  & \\
    948 NKG\_range      & \\
    949 EFRCTHN[2]      & \\
    950 WMAX[2]         & \\
    951 rthin\_rmax     & \\
    952 viewcone\_angles[2] & Inner and outer angles of Corsika's VIEWCONE
    953 option. \\
    954 telescopePhi    & Telescope azimuth (rad). Measured from South, counter-clockwise \\
    955 telescopeTheta  & Telescope zenith angle (rad) \\
    956 TimeFirst       & Arrival time on camera of first photon (ns) \\
    957 TimeLast        & Arrival time on camera of last photon (ns) \\
    958 
    959 & 6 next variables: CORSIKA longitudinal particle fit parameters \\
    960 & \hspace{0.5cm} (see CORSIKA manual for precise meaning and units)\\
    961 longi\_Nmax        & Numer of charged particles at maximum \\
    962 longi\_t0          & Atmospheric depth of shower starting point (N=0) \\
    963 longi\_tmax        & Atmospheric depth of shower maximum (g/cm$^2$) \\
    964 longi\_a           & \\
    965 longi\_b           & For {\bf longi\_a}, {\bf longi\_b}, {\bf longi\_c}, see CORSIKA manual  \\
    966 longi\_c           &  \\
    967 longi\_chi2        & $\chi^2/dof$ of the fit\\
     1011&& \\
     101264 & Ecutoffg        & \\
     101365 & NFLAIN             & \\
     101466 & NFLDIF             & \\
     101567 & NFLPI0             & \\
     101668 & NFLPIF             & \\
     101769 & NFLCHE             & \\
     101870 & NFRAGM             & \\
     101971 & Bx         & \\
     102072 & By         & \\
     102173 & EGS4yn     & \\
     102274 & NKGyn              & \\
     102375 & GHEISHAyn  & \\
     102476 & VENUSyn            & \\
     102577 & CERENKOVyn & \\
     102678 & NEUTRINOyn & \\
     102779 & HORIZONTyn & \\
     102880 & COMPUTER   & \\
     102981 & ThetaMin   & Minimum Theta of primaries (deg) \\
     103082 & ThetaMax   & Maximum Theta of primaries (deg) \\
     103183 & PhiMin     & Minimum Phi of primaries (deg) \\
     103284 & PhiMax     & Maximum Phi of primaries (deg) \\
     103385 & CBunchSize & \\
     103486 & CDetInX    & \\
     103587 & CDetInY    & \\
     103688 & CSpacInX   & \\
     103789 & CSpacInY   & \\
     103890 & CLenInX    & \\
     103991 & CLenInY    & \\
     104092 & COutput    & \\
     104193 & AngleNorthX& \\
     104294 & MuonInfo   & \\
     104395 & StepLength & \\
     104496 & CWaveLower & Wavelength lower limit (nm) \\
     104597 & CWaveUpper & Wavelength upper limit (nm) \\
     104698 & Multipl    & \\
     104799 to 138 & CorePos[2][20] & Core positions of randomized shower \\
     1048139 to 140 & SIBYLL[2]  & \\
     1049141 to 142 & QGSJET[2]  & \\
     1050143 to 144 & DPMJET[2]  & \\
     1051145 & VENUS\_cross      & \\
     1052146 & mu\_mult\_scat    & \\
     1053147 & NKG\_range        & \\
     1054148 to 149 & EFRCTHN[2] & \\
     1055150 to 151 & WMAX[2]    & \\
    9681056& \\
    9691057\hline
     
    9721060\newpage
    9731061
    974 \begin{tabular}{ll}
    975 \parbox{5cm}{Variable} & Description \\
     1062\begin{tabular}{lll}
     1063\multicolumn{2}{c}{Variable} & \parbox{11cm}{Description} \\
    9761064\hline
    977 & \\
    978 CORSIKAPhs      & Original photons written by {\it CORSIKA} \\
    979 AtmAbsPhs       & Photons absorbed by the atmosphere  \\
    980 MirrAbsPhs      & Photons absorbed by the mirror      \\
    981 OutOfMirrPhs    & Photons outside the mirror          \\
    982 BlackSpotPhs    & Photons lost in the "black spot"    \\
    983 OutOfChamPhs    & Photons outside the camera         \\
    984 CPhotons        & Photons reaching the camera        \\
    985 
    986 elec\_cph\_fraction & Fraction of C-photons produced by electrons \\
    987 muon\_cph\_fraction & Fraction of C-photons produced by muons \\
    988 other\_cph\_fraction & Fraction of C-photons produced by electrons \\
     1065&& \\
     1066152 & rthin\_rmax       & \\
     1067153 to 154 & viewcone\_angles[2] & Inner and outer angles of Corsika's VIEWCONE
     1068option. \\
     1069155 & telescopePhi      & Telescope azimuth (rad). Measured from South, counter-clockwise \\
     1070156 & telescopeTheta    & Telescope zenith angle (rad) \\
     1071157 & TimeFirst         & Arrival time on camera of first photon (ns) \\
     1072158 & TimeLast          & Arrival time on camera of last photon (ns) \\
     1073
     1074&& 6 next variables: CORSIKA longitudinal particle fit parameters \\
     1075&& \hspace{0.5cm} (see CORSIKA manual for precise meaning and units)\\
     1076159 & longi\_Nmax        & Numer of charged particles at maximum \\
     1077160 & longi\_t0          & Atmospheric depth of shower starting point (N=0) \\
     1078161 & longi\_tmax        & Atmospheric depth of shower maximum (g/cm$^2$) \\
     1079162 & longi\_a           & \\
     1080163 & longi\_b           & For {\bf longi\_a}, {\bf longi\_b}, {\bf longi\_c}, see CORSIKA manual  \\
     1081164 & longi\_c           &  \\
     1082165 & longi\_chi2        & $\chi^2/dof$ of the fit\\
     1083166 & CORSIKAPhs      & Original photons written by {\it CORSIKA} \\
     1084167 & AtmAbsPhs       & Photons absorbed by the atmosphere  \\
     1085168 & MirrAbsPhs      & Photons absorbed by the mirror      \\
     1086169 & OutOfMirrPhs    & Photons outside the mirror          \\
     1087170 & BlackSpotPhs    & Photons lost in the "black spot"    \\
     1088171 & OutOfChamPhs    & Photons outside the camera         \\
     1089172 & CPhotons        & Photons reaching the camera        \\
     1090
     1091173 & elec\_cph\_fraction & Fraction of C-photons produced by electrons \\
     1092174 & muon\_cph\_fraction & Fraction of C-photons produced by muons \\
     1093175 & other\_cph\_fraction & Fraction of C-photons produced by electrons \\
     1094176 to 182 & dummy[7] & not used \\
    9891095& \\
    9901096\hline
    9911097\end{tabular}
    992 
    993 \vspace*{1cm}
    994 The event header is followed by 8-word blocks, one for each photon
    995 that reaches the camera. A photon block contains the following
    996 variables:
     1098%
     1099\vspace*{0.5cm}
     1100\\
     1101The event header is followed by 8-word blocks, one for each Cherenkov
     1102photon that reaches the camera. A photon block contains the following
     1103variables (as 4-byte float numbers):
    9971104\vspace*{0.5cm}
    9981105\\
     
    10061113        & n is the index from 1 to ICERML (see Corsika manual) for the case in which each Corsika \\
    10071114        & shower is used more than once (normally, in MMCS will be just 1). \\
    1008         & \\
    10091115x, y    & Impact point in camera coordinates (cm) \\
    10101116u, v    & Director cosines of down-going versor indicating the photon direction \\
    1011 t       & Arrival time on camera (ns) \\
    1012 h       & Production height (cm) \\
     1117t       & Arrival time on camera (ns), measured from the time of first
     1118interaction of the primary \\
     1119h       & Production height (cm), measured above sea level on the
     1120vertical of the telescope location \\
     1121        & (it is not the {\it true} height which would be measured on
     1122the vertical of the emitting particle!)  \\
    10131123phi     & Incidence angle with respect to camera plane (rad) \\
    10141124& \\
     
    10171127\vspace*{0.5cm}
    10181128\\
    1019 After the last event photon block there is a blank line, an \verb$END---EVENT$
    1020 flag, another blank line and then the following event. After the last
    1021 event in a run it appears the flag ``\verb$END-----RUN$'', while after all
    1022 the processed runs, a ``\verb$END----FILE$'' flag is written. Finally,
    1023 after this flag an ``ascii tail'' has been added to the file: it
    1024 consists of the ascii files {\it magic.def}, {\it axisdev.dat} and {\it
    1025 reflectivity.dat} one after the other separated by blank lines. In
    1026 this way all the relevant parameters used to produce the output are
     1129After the last photon block of an event there is a carriage return
     1130followed by the ascii flag indicating the event end, and then two more
     1131carriage returns before the ascii flag of the beginning of the next
     1132event, and so on:\\
     1133\\
     1134\verb"$END---EVENT$$START-EVENT$Event_header....$END---EVENT$$END-----RUN$$START---RUN$..."\\
     1135\\
     1136The flag ``\verb$END-----RUN$'' appears after the last event in a run
     1137(that is, the last event processed of each of the input cer
     1138files). After the last processed run, an end of file flag is
     1139written:\\
     1140\\
     1141\verb"...$END---EVENT$$END-----RUN$$END----FILE$"\\
     1142\\
     1143Finally, after this flag an ``ascii tail'' has been attached to the file:
     1144it consists of the ascii files {\it magic.def}, {\it axisdev.dat} and
     1145{\it reflectivity.dat} one after the other, separated by carriage returns:
     1146\\
     1147\\
     1148\verb"$magic.def$axisdev.dat$reflectivity.dat"\\
     1149\\
     1150In this way all the relevant parameters used to produce the output are
    10271151kept together with the reflected events.
    1028 
    10291152%------------------------------------------------------------
    1030 
    1031 \section{Appendix A}
     1153\newpage
     1154\renewcommand{\thesubsection}{A.\arabic{subsection}}
     1155\section*{Appendix A : atmospheric absorption}
     1156\addcontentsline{toc}{section}{Appendix A : atmospheric absorption}
     1157%
     1158The simulation of the absorption of Cherenkov light in the atmosphere
     1159has been included in the {\it Reflector} program because this feature
     1160was not yet available in the first versions of CORSIKA used within the
     1161MAGIC collaboration. In the latest CORSIKA versions, the atmospheric
     1162absorption has been included as an option, but it is not
     1163compatible with the simulation of a curved atmosphere \cite{cor02},
     1164and hence we have kept this step as a part of our reflector
     1165simulation. This appendix describes how the atmospheric absorption is
     1166implemented in the program when the option ATM\_CORSIKA (see section
     1167\ref{commands}) is chosen.
     1168\par
     1169The geometry of the problem is sketched in figure
     1170\ref{fig:atmoscheme}. A Cherenkov photon is emitted in point A and
     1171travels towards the telescope placed at B. At any moment, the height
     1172$h$ of the photon above sea level is related to the distance $L$
     1173between the photon and the telescope through
     1174%
     1175\begin{equation}
     1176(R+h)^2 = (R+h_1)^2 + L^2 + 2 L \; (R+h_1) \; \cos \theta
     1177\label{eq:height}
     1178\end{equation}
     1179%
     1180, where $R$ is the Earth radius, $h_1$ the height (a.s.l.) of the
     1181observation level and $\theta$ is the zenith angle of the photon
     1182trajectory measured at the telescope site. The Cherenkov output of
     1183CORSIKA contains for each photon the height $h_C$ of the emission
     1184point A, measured along the vertical of the observer. The {\it true
     1185vertical height} $h_2$ of the emission point can be obtained by
     1186replacing $L$ by $(h_C-h_1)/\cos \theta$ in equation
     1187(\ref{eq:height}).
     1188%
     1189\begin{figure}[ht]
     1190\begin{center}
     1191\mbox{ \epsfig{file=eps/atmoscheme.eps,width=0.8\textwidth} }
     1192\end{center}
     1193\caption[]
     1194{Calculation of the true vertical height $h_2$ of the emission point
     1195of a Cherenkov photon (point A), and the optical path traversed down to the
     1196telescope (point B).}
     1197\label{fig:atmoscheme}
     1198\end{figure}
     1199%
     1200\par
     1201The optical path $I(\theta, h_1, h_2)$ traversed by the photon can be
     1202calculated integrating the air density along the trajectory
     1203$\overline{\text{AB}}$. For the case $h/R \ll 1$, we can drop in
     1204(\ref{eq:height}) the terms in $(h/R)^2$ and smaller, and then solve
     1205for L. Deriving now with respect to $h$, we have:
     1206%
     1207\begin{equation}
     1208\frac{dL}{dh} \simeq \sqrt{\frac{R}{2\;(h-h_1)+R\;\cos^2 \theta}}
     1209\qquad \text{for} \qquad \frac{h}{R} \ll  1
     1210\label{eq:dldh}
     1211\end{equation}
     1212%
     1213\subsection{Rayleigh scattering}
     1214Rayleigh scattering is the scattering of light by particles smaller
     1215than its wavelength. These are in our case the air molecules. The
     1216transmission coefficient due to Rayleigh scattering is a strong
     1217function of the wavelength $\lambda$:
     1218%
     1219\begin{equation}
     1220T_{\text{Rayl}} (\lambda) = \exp \Biggl[ - \frac{I(\theta, h_1, h_2)}{x_R} \; \Biggl(\frac{400 \;
     1221\text{nm}}{\lambda}\Biggl)^4 \; \Biggl]
     1222\label{eq:rayleigh}
     1223\end{equation}
     1224%
     1225Here $I(\theta, h_1, h_2)$ is the optical path (in g/cm$^2$) traversed
     1226between points A and B, and $x_R = 2970$ g/cm$^2$ is the mean free path of
     1227the Rayleigh scattering at $\lambda = 400$ nm.
     1228\par
     1229A convenient way of expressing the optical path is the following:
     1230%
     1231\begin{equation}
     1232I\;(\theta, h_1, h_2) = (x_1 - x_2) \cdot \mathcal{AM}\;(\theta, h_1, h_2)
     1233\end{equation}
     1234Here $x_{i=1,2}$ is the mass overburden of the atmosphere above a
     1235height $h_i$ (in the
     1236vertical direction) and $\mathcal{AM}$ is the so called {\it air
     1237mass}\footnote{If we set $h_2 = \infty$, we have the usual definition
     1238of {\it air mass} in optical astronomy.}, defined as
     1239%
     1240\begin{equation}
     1241\mathcal{AM} \equiv \frac{I\;(\theta,h_1,h_2)}{I\;(0^\circ,h_1,h_2)}
     1242\label{eq:airmass}
     1243\end{equation}
     1244%
     1245, which is a function mainly of the zenith angle $\theta$ (see
     1246fig. \ref{fig:airmass}). In our simulation, for the calculation of the
     1247mass overburden $x_i$ we have used the U.S. standard atmosphere
     1248parametrized by J. Linsley \cite{cor02}, the same we used in Corsika
     1249for the shower development simulation. It consists of five layers: in
     1250the lower four the density decreases exponentially with height, and in
     1251the upper one the mass overburden cecreases linearly until it vanishes
     1252at $h = 112.8$ km.
     1253%
     1254\begin{figure}[ht]
     1255\begin{center}
     1256\mbox{ \epsfig{file=eps/airmass.eps,width=0.8\textwidth} }
     1257\end{center}
     1258\caption[]
     1259{Dependence with zenith angle of the air mass as defined in the
     1260text. The air mass has been calculated for an exponential atmosphere
     1261with scale height $H = 7.4$ km, for the observation level of MAGIC
     1262(2.2 km a.s.l.), and for light coming from $h_2 = 10$ and at 100 km a.s.l. As
     1263we can see the dependence with the emission height $h_2$ is very small.}
     1264\label{fig:airmass}
     1265\end{figure}
     1266%
     1267\par
     1268For the estimate of $\mathcal{AM}$, a simpler atmospheric model has
     1269been used, in which the vertical density profile is described by a
     1270single exponential: $\rho = \rho_0 \; e^{-h/H}$  with scale height $H
     1271= 7.4$ km. This simplifies the calculations, and is accurate enough
     1272for our purposes. Using (\ref{eq:dldh}) the optical path $I(\theta,
     1273h_1, h_2)$ can then be obtained approximately as:
     1274%
     1275\begin{equation}
     1276I(\theta, h_1, h_2) = \int_A^B \rho\;(h)\; \frac{dL}{dh} \; dh \simeq
     1277\sqrt{\frac{R}{2}} \;
     1278\int_{h_1}^{h_2} \frac{\rho_0
     1279\;e^{-h/H}}{\sqrt{h-h_1+\frac{1}{2}\;R\;\cos^2 \theta}} \; dh
     1280\label{eq:optpath}
     1281\end{equation}
     1282%
     1283and finally:
     1284%
     1285\begin{equation}
     1286\mathcal{AM} \simeq e^{-\frac{R \sin^2 \theta}{2H}}
     1287\cdot
     1288\frac{
     1289\text{erfc}\;(\sqrt{\frac{R \cos^2 \theta}{2H}})
     1290\;-\;
     1291\text{erfc}\;(\sqrt{\frac{2 (h_2 - h_1) + R \cos^2 \theta}{2H}})
     1292}
     1293{
     1294\text{erfc}\;(\sqrt{\frac{R}{2H}})
     1295\;-\;
     1296\text{erfc}\;(\sqrt{\frac{2(h_2 - h_1) + R}{2 H}})
     1297}
     1298\label{eq:airmass2}
     1299\end{equation}
     1300%
     1301where we have used the complementary error function $\text{erfc}\;(x)
     1302= \frac{2}{\sqrt{\pi}}\int_x^\infty e^{-t^2} dt$. From
     1303(\ref{eq:airmass2}), $\mathcal{AM}$ can be readily evaluated for any
     1304value of $\theta$, $h_1$ and $h_2$. In fig. \ref{fig:rayleigh} the
     1305resulting Rayleigh transmission coefficient $T_{\text{Rayl}}$ finally
     1306obtained from (\ref{eq:rayleigh}) is plotted versus the zenith angle
     1307for three wavelengths.
     1308%
     1309\begin{figure}[ht]
     1310\begin{center}
     1311\mbox{ \epsfig{file=eps/rayleigh.eps,width=\textwidth} }
     1312\end{center}
     1313\caption[]
     1314{Rayleigh transmission coefficient as a function of zenith angle for three
     1315different wavelengths. The solid, dashed and dotted lines correspond
     1316respectively to light coming from 5, 10 and 20 km distance from the
     1317telescope.}
     1318\label{fig:rayleigh}
     1319\end{figure}
     1320%
     1321\subsection{Mie scattering}
     1322Cherenkov light also suffers Mie scattering through interaction with
     1323small dust particles suspended in the air (aerosols), whose size is
     1324comparable to the wavelength of the light. In our simulation of the
     1325attenuation due to aerosols we have used the model proposed by
     1326Elterman \cite{elterman64,elterman65}, which considers an aerosol
     1327number density $N_p$ which (roughly) decreases exponentially up to 10
     1328km a.s.l. with scale height $H \simeq 1.2$ km, followed by a more
     1329tenuous layer between 10 and 30 km (see fig. \ref{fig:aerosol}a). In
     1330this model, the aerosol size distribution is considered to be
     1331unchanged with altitude.
     1332%
     1333\begin{figure}[ht]
     1334\begin{center}
     1335\mbox{ \epsfig{file=eps/aerosol.eps,width=0.9\textwidth} }
     1336\end{center}
     1337\caption[]
     1338{Aerosol model by Elterman: in (a), number density of
     1339aerosols as a function of height above sea level; in (b), aerosol
     1340attenuation coefficient at sea level as a function of wavelength.}
     1341\label{fig:aerosol}
     1342\end{figure}
     1343%
     1344\par
     1345Measured values of the aerosol attenuation coefficients at sea level
     1346$\beta_p(0)$ for different wavelengths \cite{elterman65} are shown in
     1347figure \ref{fig:aerosol}b. To obtain the attenuation coefficient at a
     1348given height $h$, we simply do $\beta_p(h, \lambda) = \beta_p(0,
     1349\lambda) \cdot N_p(h)/N_p(0)$. In the Elterman model the aerosol
     1350transmission coefficient for the trajectory from A to B depicted in
     1351figure \ref{fig:atmoscheme} would then be:
     1352%
     1353\begin{equation}
     1354T_{\text{Mie}}(\lambda) = e^{-\tau_{mie}} \quad \text{, with}\quad
     1355\tau_{mie}(h_1, h_2, \theta, \lambda) = \frac{\beta_p(0, \lambda)}{N_p(0)} \;
     1356\int_{h_1}^{h_2}
     1357\; N_p(h) \;\; \frac{dL}{dh}\; dh
     1358\label{eq:aerosoltau}
     1359\end{equation}
     1360%
     1361Here $\tau_{mie}$ is the aerosol optical depth of the path from A to
     1362B. Given that the aerosol density distribution is not a simple
     1363exponential, we have to do the integral in (\ref{eq:aerosoltau})
     1364numerically. The integral depends on $h_2$ and on $\theta$, through
     1365$dL/dh$ (it also
     1366depends on $h_1$, but this is the observation level and therefore
     1367fixed). In this case we use the exact expression for $\;dL/dh\;$ which can be
     1368obtained from (\ref{eq:height}). At the beginning of each simulation
     1369run we calculate and store the results of the integral for values of
     1370$\theta$ between 0 and 90$^\circ$ (in steps of $1^\circ$), and of
     1371$h_2$ from $h_1$ up to 30 km (in steps of 100 m). To do the integral a
     1372linear interpolation has been used to obtain the value of $N_p$ for
     1373any height. During the simulation of each Cherenkov photon, we get the
     1374corresponding precalculated value of the integral and deduce
     1375$T_{\text{Mie}}$ from expression (\ref{eq:aerosoltau}).
     1376%
     1377\begin{figure}[ht]
     1378\begin{center}
     1379\mbox{ \epsfig{file=eps/mie.eps,width=\textwidth} }
     1380\end{center}
     1381\caption[]
     1382{Aerosol transmission coefficient for three different wavelengths as a
     1383function of the distance between the photon emission point and the
     1384telescope. Plots for four different zenith angles between 0 and 80
     1385degrees are shown.}
     1386\label{fig:mie}
     1387\end{figure}
     1388%
     1389\par
     1390Since in this model the aerosols are concentrated mainly at very low
     1391altitude, the transmission coefficient is more or less constant above
     1392a certain height (which depends on $\theta$), as can be seen in
     1393fig. \ref{fig:mie}. For instance, for vertically incident 300 nm light
     1394emitted higher than 4 km above the telescope, the Mie transmission is
     1395about 0.95.
     1396%
     1397\subsection{Absorption by Ozone}
     1398%
     1399The absorption of Cherenkov light by Ozone has been implemented
     1400following also the Elterman standard atmosphere \cite{elterman65}. The
     1401coefficient for ozone absorption is given by
     1402%
     1403\begin{equation}
     1404\beta_3(h,\lambda) = A_v(\lambda) \cdot D_3(h)
     1405\end{equation}
     1406%
     1407, where $A_v(\lambda)$ is the Vigroux \cite{vigroux53} ozone
     1408absorption coefficient (cm$^{-1}$) and $D_3(h)$ is the ozone
     1409concentration (cm km$^{-1}$) according to Elterman. The transmission
     1410coefficient of Ozone in the path $\overline{\text{AB}}$ is then:
     1411%
     1412\begin{equation}
     1413T_{\text{Ozone}}(\lambda) = e^{-\tau_{oz}} \quad \text{, with}\quad
     1414\tau_{oz}(h_1, h_2, \theta, \lambda) = A_v(\lambda) \; \int_{h_1}^{h_2}
     1415\; D_3(h) \;\; \frac{dL}{dh}\; dh
     1416\label{eq:ozonetau}
     1417\end{equation}
     1418%
     1419\begin{figure}[ht]
     1420\begin{center}
     1421\mbox{ \epsfig{file=eps/ozone.eps,width=\textwidth} }
     1422\end{center}
     1423\caption[]
     1424{Ozone concentration vertical profile (a) and Vigroux coefficients for
     1425Ozone absorption (b). The Vigroux coefficients for $\lambda =$ 380 and
     1426400 nm are zero.}
     1427\label{fig:ozone}
     1428\end{figure}
     1429%
     1430\par
     1431Once again, like in the case of Mie scattering, the optical depth
     1432$\tau_{oz}$ is the product of a factor which depends on $\lambda$ and
     1433an integral which depends on $h_1$, $h_2$ and $\theta$. We proceed in
     1434the same way as before, precalculating the values of the integral in
     1435steps of $\Delta\theta = 1^\circ$ and $\Delta h = 100$ m, up to a
     1436height of 50 km a.s.l. (where the ozone concentration becomes
     1437negligible), and then reading the appropriate values for every
     1438simulated photon.
     1439\par
     1440Finally, the overall atmospheric transmission coefficient is
     1441calculated as
     1442%
     1443\begin{equation}
     1444T_{total} = T_{Ray} \cdot T_{Mie} \cdot T_{Ozone}
     1445\end{equation}
     1446%
     1447In figure \ref{fig:absorplot} the atmospheric transmission as a
     1448function of the distance to the telescope for $\theta = 0^\circ$ is
     1449shown. Ozone absorption turns out to be dominant below 320 nm, while
     1450Rayleigh scattering is the main cause of the loss of photons at longer
     1451wavelengths.
     1452%
     1453\begin{figure}[ht]
     1454\begin{center}
     1455\mbox{ \epsfig{file=eps/absorplot.eps,width=\textwidth} }
     1456\end{center}
     1457\caption[]
     1458{Total transmission coefficient for vertically incident light as a
     1459function of the distance between the emission point and the
     1460telescope. The contributions of absorption by ozone and of Rayleigh
     1461and Mie scattering are also shown for comparison.}
     1462\label{fig:absorplot}
     1463\end{figure}
     1464%
     1465\newpage
     1466\section*{Appendix B : files in the reflector package}
     1467\addcontentsline{toc}{section}{Appendix B : files in the reflector package}
    10321468
    10331469The list of all Reflector files follows.
     
    10601496MagicProgs/Simulation/Detector/Reflector_0.6/version.h
    10611497
    1062 MagicProgs/Simulation/Detector/Reflector_0.6/doc/Tdas0211.ps.gz
     1498MagicProgs/Simulation/Detector/Reflector_0.6/doc/Tdas0211.ps
    10631499MagicProgs/Simulation/Detector/Reflector_0.6/doc/Tdas0211.tex
    10641500MagicProgs/Simulation/Detector/Reflector_0.6/doc/magic-tdas.sty
    10651501MagicProgs/Simulation/Detector/Reflector_0.6/doc/magiclogo.eps
    1066 MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/colimation.eps.gz
    1067 MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/coma.eps.gz
    1068 MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/coma_star.eps.gz
    1069 MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/coorsystems.eps.gz
    1070 MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/evcompare.eps.gz
    1071 MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/parabola.eps.gz
    1072 MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/refl06images.eps.gz
    1073 MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/spot10kmf1694.eps.gz
    1074 MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/spot10kmf1700.eps.gz
    1075 MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/spot_inf_f1697.eps.gz
    1076 MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/telecoor.eps.gz
    1077 MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/timing.eps.gz
     1502MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/absorplot.eps
     1503MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/aerosol.eps
     1504MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/airmass.eps
     1505MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/atmoscheme.eps
     1506MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/colimation.eps
     1507MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/coma.eps
     1508MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/coma_star.eps
     1509MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/coorsystems.eps
     1510MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/evcompare.eps
     1511MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/mie.eps
     1512MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/ozone.eps
     1513MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/parabola.eps
     1514MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/rayleigh.eps
     1515MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/refl06images.eps
     1516MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/reflec.eps
     1517MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/spot10kmf1694.eps
     1518MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/spot10kmf1700.eps
     1519MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/spot_inf_f1697.eps
     1520MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/telecoor.eps
     1521MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/timing.eps
    10781522
    10791523MagicProgs/Simulation/Detector/Reflector_0.6/tester/Makefile
     
    10961540%%>>>> Or the following if you include here by hand your
    10971541%%>>>> bibliographic entries
     1542
     1543\begin{thebibliography}{00}
     1544
     1545\bibitem{elterman64}
     1546L. Elterman, Applied Optics Vol. 3, No. 6 (1964) 745.
     1547
     1548\bibitem{elterman65}
     1549L. Elterman, R.B. Toolin, S.L. Valley (editor), Handbook of
     1550geophysics and space environments, McGraw-Hill, N.Y. (1965).
     1551
     1552\bibitem{cor02}D. Heck and J. Knapp, EAS Simulation with CORSIKA: A User's
     1553Manual, 2002.
     1554
     1555\bibitem{vigroux53}
     1556E. Vigroux, Contributions \`a l'étude expérimentale de l'absorption
     1557de l'ozone, Annales de Physique, v. 8 (1953) 709.
     1558
     1559\end{thebibliography}
    10981560
    10991561\end{document}
  • trunk/MagicSoft/Simulation/Detector/ReflectorII/init.c

    r1614 r1722  
    113113    fseek(rflfile, 0L, SEEK_SET);
    114114    fprintf(rflfile, "%s %s", QUOTE(PROGRAM), QUOTE(VERSION));
    115     fputc(0, rflfile);
     115
     116    /* Removed by A.M. January 2003: */
     117    /*    fputc(0, rflfile); */
     118
    116119    fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile);
    117120
  • trunk/MagicSoft/Simulation/Detector/ReflectorII/init.h

    r1614 r1722  
    102102    "        - D. Bastieri & C. Bigongiari - Jan 2000\n" \
    103103    "    Version 0.6\n" \
    104     "        - A. Moralejo - October 2002\n" \
     104    "        - A. Moralejo - January 2003\n" \
    105105    "  ################################################\n\n"
    106106#define SIGN_ERROR_FTL          /*  program, version    */ \
  • trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.c

    r1614 r1722  
    1212char axisdev_filename[256], reflectivity_filename[256];
    1313char geo_filename[256];
     14char atmosphere[256]; /* atmospheric model */
    1415
    1516/*  Prototypes  */
     
    6768    extern FILE *geofile;       /*  geo file (init)     */
    6869    extern void SetVerbose(int vlevel); /*  from diag.c */
    69     extern void SetAtmModel(char *model); /*  from atm.c        */
    7070    extern int ParseLine(FILE *parfile, /*  FILE with parms     */
    7171                         const char *token_list[], /*  array w/tokens   */
     
    9494                break;
    9595              case atm_model:
    96                 SetAtmModel(value_ptr);
     96                strcpy(atmosphere, value_ptr);
    9797                break;
    9898              case verbose_level:
  • trunk/MagicSoft/Simulation/Detector/ReflectorII/reflector.c

    r1622 r1722  
    4747extern CerEventHeader *cheadp;  /*  var inited in header.c      */
    4848extern RflEventHeader *rheadp;  /*  var inited in header.c      */
     49extern void SetAtmModel(int model, float ol); /*  from atm.c    */
    4950
    5051/*  Prototypes  */
     
    470471  RflRunHeader RflRunHead;
    471472
    472   extern int atmModel;      /* current atm. model */
     473  const char *AtmModelNames[]={"ATM_NOATMOSPHERE","ATM_90PERCENT","ATM_CORSIKA"};
     474  extern char atmosphere[256];  /* current atm. model */
     475  int model=0;
    473476
    474477  do
     
    490493          if (newFile)
    491494            {
     495              while (strcmp(atmosphere, AtmModelNames[model]))
     496                if (++model == sizeof(AtmModelNames)/sizeof(AtmModelNames[0]))
     497                  {   
     498                    model = 0;
     499                    Error(" *** Atm model \"%s\" not found.\n", atmosphere);
     500                    break;     
     501                  }
     502
    492503              /* Write Reflector "run header" (one per cer file!): */
    493504              memcpy(&RflRunHead, cheadp, sizeof(RflRunHead));
    494505              RflRunHead.wobble_mode = wobble_position;
    495               RflRunHead.atmospheric_model = atmModel;
     506              RflRunHead.atmospheric_model = model;
    496507              fwrite(&RflRunHead, sizeof(RflRunHead), 1, rflfile);
     508
     509              /* Set up atmosphere (do some necessary calculations) once we
     510               * read in the observation level from the run header:
     511               */
     512
     513              Log("Setting atm model \"%s\" and observation level %.0f cm.\n", AtmModelNames[model], RflRunHead.HeightLev[0]);
     514              SetAtmModel(model,RflRunHead.HeightLev[0]);
     515
    497516            }
    498517
Note: See TracChangeset for help on using the changeset viewer.