Changeset 1722 for trunk

Timestamp:

01/21/03 15:02:35 (22 years ago)

Author:

moralejo

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Simulation/Detector/ReflectorII

Files:

• Changelog (modified) (2 diffs)
• Makefile (modified) (3 diffs)
• atm.c (modified) (1 diff)
• atm.h (modified) (2 diffs)
• attenu.f (modified) (5 diffs)
• doc/Tdas0211.ps.gz (modified) ( previous)
• doc/Tdas0211.tex (modified) (27 diffs)
• doc/eps/absorplot.eps.gz (added)
• doc/eps/aerosol.eps.gz (added)
• doc/eps/airmass.eps.gz (added)
• doc/eps/atmoscheme.eps.gz (added)
• doc/eps/mie.eps.gz (added)
• doc/eps/ozone.eps.gz (added)
• doc/eps/rayleigh.eps.gz (added)
• doc/eps/reflec.eps.gz (added)
• init.c (modified) (1 diff)
• init.h (modified) (1 diff)
• parms.c (modified) (3 diffs)
• reflector.c (modified) (3 diffs)

trunk/MagicSoft/Simulation/Detector/ReflectorII/Changelog

@@ -1 +1 @@
 ** Add changes at the beginning! **
+
+21/01/2003 A. Moralejo
+
+Removed from the output a null byte written right after the ascii label
+containing the program version number which is at the beginning of the rfl 
+file.
+
+19/12/2002 - 17/01/2003 A. Moralejo
+
+Lots of changes. Moved simulation of Mie scattering and Ozone absorption from
+attenu.f to atm.c, and changed the way it is done. Before it was not correct
+for large zenith angle (a variation like for Rayleigh scattering was assumed 
+for both Mie scattering and Ozone absorption).
+
+
+13/12/2002 A. Moralejo
+
+attenu.f: Found mistake in Mie absorption calculation (height from sea 
+level was taken as height above observation level). Changed optical depths
+table for Mie attenuation. Now they are no longer referred to 2 km height, 
+but to sea level. Detector level is taken into account later in calculation.
+
+10/12/2002 A. Moralejo
+
+attenu.f: Added comments, removed old/unnecessary code, corrected small
+mistake in Elterman's table for ozone optical depth.
+
+atm.c,h : removed atmospheric model "ATM_ISOTHERMAL" which actually was 
+not even implemented in the code, although it could be selected in the 
+input card.
 
 15/11/2002 A. Moralejo
@@ -15 +45 @@
 
 14/11/2002 A. Moralejo
-reflector.c, parms.c Added wobble mode option. Added wobble flag to the
+reflector.c, parms.c: Added wobble mode option. Added wobble flag to the
 output and also an atmospheric_model flag.
 

trunk/MagicSoft/Simulation/Detector/ReflectorII/Makefile

@@ -19 +19 @@
 #
 # $RCSfile: Makefile,v $
-# $Revision: 1.5 $
-# $Author: bigongia $ 
-# $Date: 2002-11-14 21:39:01 $
+# $Revision: 1.6 $
+# $Author: moralejo $ 
+# $Date: 2003-01-21 15:02:32 $
 #
 ##################################################################
@@ -217 +217 @@
 atm.o: /usr/include/sys/machine/machtime.h /usr/include/sys/rt_limits.h
 atm.o: /usr/include/string.h /usr/include/strings.h /usr/include/math.h
-atm.o: /usr/include/stdlib.h diag.h atm.h init.h
+atm.o: /usr/include/stdlib.h atm.h diag.h init.h
 ph2cph.o: /usr/include/stdio.h /usr/include/standards.h
 ph2cph.o: /usr/include/sys/seek.h /usr/include/va_list.h
@@ -229 +229 @@
 header.o: /usr/include/sys/types.h /usr/include/mach/machine/vm_types.h
 header.o: /usr/include/sys/select.h /usr/include/strings.h header.h
-attach.o: /usr/include/string.h /usr/include/standards.h
+attach.o: /usr/include/stdio.h /usr/include/standards.h
+attach.o: /usr/include/sys/seek.h /usr/include/va_list.h
 attach.o: /usr/include/sys/types.h /usr/include/mach/machine/vm_types.h
-attach.o: /usr/include/sys/select.h /usr/include/strings.h
+attach.o: /usr/include/sys/select.h /usr/include/getopt.h
+attach.o: /usr/include/sys/limits.h /usr/include/sys/machine/machlimits.h
+attach.o: /usr/include/sys/syslimits.h /usr/include/sys/machine/machtime.h
+attach.o: /usr/include/sys/rt_limits.h
 reflector.o: /usr/include/stdio.h /usr/include/standards.h
 reflector.o: /usr/include/sys/seek.h /usr/include/va_list.h

trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c

@@ +1 @@
+/********************************************************************
+*                                                                   *
+*  File: atm.c                                                      *
+*  Authors: J.C. Gonzalez, A. Moralejo                              *
+*                                                                   *
+* January 2002, A. Moralejo: lots of changes. Moved the code for    *
+*         the Mie scattering and ozone absorption from attenu.f to  *
+*         here, after some bugs were found. Now the implementation  *
+*         is different, we now precalculate the slant paths for the *
+*         aerosol and Ozone vertical profiles, and then do an       *
+*         interpolation in wavelength for every photon to get the   *
+*         optical depths. The parameters used, defined below,       *
+*         have been taken from "Atmospheric Optics", by L. Elterman *
+*         and R.B. Toolin, chapter 7 of the "Handbook of geophysics *
+*         and Space environments". (S.L. Valley, editor).           *
+*         McGraw-Hill, NY 1965.                                     *
+*                                                                   *
+*         WARNING: the Mie scattering and the Ozone absorption are  *
+*         implemented to work only with photons produced at a       *
+*         height a.s.l larger than the observation level. So this   *
+*         is not expected to work well for simulating the telescope *
+*         pointing at theta > 90 deg (for instance for neutrino     *
+*         studies. Rayleigh scattering works even for light coming  *
+*         from below.                                               *
+*                                                                   *
+*********************************************************************/
+
 #include <stdio.h>
 #include <string.h>
 #include <math.h>
+
+#include "atm.h"
 #include "diag.h"
-#include "atm.h"
 #include "init.h"
 
 /*  random numbers  */
 #define RandomNumber  ranf()
-
-/* AM Nov 2002: atmModel is now an external variable: */
-int atmModel=0;   /* current atm. model */
+#define STEPTHETA 1.74533e-2 /* aprox. 1 degree */
+
+#define MIN(x,y) ((x)<(y)? (x) : (y))
 
 /*  Function declarations  */
 static float atm(float wavelength, float height, float theta);
-void SetAtmModel(char *model);
+void SetAtmModel(int model, float ol);
 int absorption(float wlen, float height, float theta);
-extern void attenu_(float *, float *, float *, float *);  /* in Fortran */
+extern void attenu_(float *, float *, float *, float *, float *);  /* in Fortran        */
 extern float ranf(void);
 
-void SetAtmModel(char *model)
+/* aerosol_path contains the path integrals for the aerosol number
+ * density (relative to the number density at sea level) between the 
+ * observation level and a height h for different zenith angles. The 
+ * first index indicate height above sea level in units of 100m, the 
+ * second is the zenith angle in degrees.
+ */
+static float aerosol_path[301][90];
+
+/* ozone_path contains the path integrals for the ozone concentration
+ * between the observation level and a height h for different zenith 
+ * angles. The first index indicate height above sea level in units
+ * of 100m, the second is the zenith angle in degrees.
+ */
+static float ozone_path[501][90];
+
+static float obslev; /* observation level in cm */
+static double rt;    /* Earth radius in cm */
+static int atmModel;
+
+void SetAtmModel(int model, float ol)
 {
-    while (strcmp(model, AtmModelNames[atmModel]))
-        if (++atmModel == sizeof(AtmModelNames)/sizeof(AtmModelNames[0]))
-        {   atmModel = 0;
-            Error(ATM__NFND__ERR, model);
-            break;      }
-
-    Log(ATM__SET___LOG, AtmModelNames[atmModel]);
-
-    /*  'erfc' must be inited, so keep the following line.  */
-    Debug("Initing 'erfc': ret=%g\n", erfc(M_PI));
+  float Rcos2, sin2, rtsq, path_slant, h, dh, theta;
+  int j;
+
+  atmModel = model;
+  obslev = ol; 
+  rt= 6371315.E2; /* Earth radius (same as in Corsika) in cm */
+
+  if (atmModel == ATM_CORSIKA)
+    {
+      /* It follows a precalculation of the slant path integrals we need 
+       * for the estimate of the Mie scattering and Ozone absorption:
+       */
+
+      rtsq = sqrt(rt);
+      dh = 1.e3;
+
+      /* Mie (aerosol): */
+
+      for (j = 0; j < 90; j++)
+        {
+          theta = j * STEPTHETA;  /* aprox. steps of 1 deg */
+
+          path_slant = 0;
+          Rcos2 = rt * cos(theta)*cos(theta);
+          sin2 = sin(theta)*sin(theta);
+
+          for (h = obslev; h <= 30e5; h += dh)
+            {
+              if (fmod(h,1e4) == 0)
+                aerosol_path[(int)(h/1e4)][j] = path_slant;
+
+              path_slant += 
+                (aero_n[(int)floor(h/1.e5)] + (h/1.e5 - floor(h/1.e5))*
+                  (aero_n[(int)ceil(h/1.e5)]-aero_n[(int)floor(h/1.e5)]))
+                  /aero_n[0] * dh * (rt+h) /
+                    sqrt((rt+h)*(rt+h)-(rt+obslev)*(rt+obslev)*sin2);
+
+            }
+        }
+
+      /* Ozone absorption */
+
+      for (j = 0; j < 90; j++)
+        {
+          theta = j * STEPTHETA;  /* aprox. steps of 1 deg */
+          path_slant = 0;
+          Rcos2 = rt * cos(theta)*cos(theta);
+          sin2 = sin(theta)*sin(theta);
+
+          for (h = obslev; h <= 50e5; h += dh)
+            {
+              if (fmod(h,1e4) == 0)
+                ozone_path[(int)(h/1e4)][j] = path_slant;
+
+              path_slant += 
+                (oz_conc[(int)floor(h/1.e5)] + (h/1.e5 - floor(h/1.e5))*
+                  (oz_conc[(int)ceil(h/1.e5)]-oz_conc[(int)floor(h/1.e5)]))
+                  * dh * (rt+h) /
+                    sqrt((rt+h)*(rt+h)-(rt+obslev)*(rt+obslev)*sin2);
+            }
+        }
+
+    }
+
 }   /*  end of SetAtmModel  */
-            
+
 static float atm(float wavelength, float height, float theta)
-{   float transmittance = 1.0;  /* final atm transmittance (ret. value) */
-
-    switch(atmModel)
-    {   case ATM_NOATMOSPHERE:  /* no atm at all: transmittance = 100%  */
-             break;
-        case ATM_90PERCENT:     /* atm. with transmittance = 90%        */
-             transmittance = 0.9;
-             break;
-        case ATM_ISOTHERMAL:    /* isothermal approximation             */
-        /********************/
-             break;
-        case ATM_CORSIKA:       /* atmosphere as defined in CORSIKA     */
-             attenu_(&wavelength, &height, &theta, &transmittance);
-             break;
+{   
+  float transmittance = 1.0;  /* final atm transmittance (ret. value) */
+  float T_Ray, T_Mie, T_Oz;
+
+  float h; /* True height a.s.l. of the photon emission point in cm */
+  float tdist;
+  float beta0, path;
+
+  int index;
+
+  switch(atmModel)
+    {   
+    case ATM_NOATMOSPHERE:      /* no atm at all: transmittance = 100%  */
+      break;
+    case ATM_90PERCENT: /* atm. with transmittance = 90%        */
+      transmittance = 0.9;
+      break;
+    case ATM_CORSIKA:   /* atmosphere as defined in CORSIKA     */
+
+      /* Distance to telescope: */
+      tdist = (height-obslev)/cos(theta);
+
+      /* Avoid problems if photon is very close to telescope: */
+      if (fabs(tdist) < 1.)
+        {
+          transmittance = 1.;
+          break;
+        }
+
+      /*** We calculate h, the true emission height above sea level: ***/
+
+      h = -rt + sqrt((rt+obslev)*(rt+obslev) + tdist*tdist +
+                     (2*(rt+obslev)*(height-obslev)));
+
+      /******* Rayleigh scattering: *******/
+
+      attenu_(&wavelength, &h, &obslev, &theta, &T_Ray);
+
+
+      /******* Ozone absorption: *******/
+
+      if (h > 50.e5)
+        h = 50.e5;
+
+      /* First we get Vigroux Ozone absorption coefficient for the given 
+       * wavelength, through a linear interpolation:
+       */
+
+      for (index = 1; index < 11; index++)
+        if (wavelength < wl[index])
+          break;
+
+      beta0 = oz_vigroux[index-1]+(oz_vigroux[index]-oz_vigroux[index-1])*
+        (wavelength-wl[index-1])/(wl[index]-wl[index-1]);
+
+      /* from km^-1 to cm^-1 : */
+      beta0 *= 1e-5;
+      
+      /* Now use the pre-calculated values of the path integral 
+       * for h and theta: */
+
+      path = ozone_path[(int)floor(0.5+h/1e4)]
+        [(int)MIN(89,floor(0.5+theta/STEPTHETA))];
+
+      T_Oz = exp(-beta0*path);
+
+
+      /******* Mie (aerosol): *******/
+
+      if (h > 30.e5)
+        h = 30.e5;
+
+      /* First get Mie absorption coefficient at sea level for the given 
+       * wavelength, through a linear interpolation:
+       */
+
+      for (index = 1; index < 11; index++)
+        if (wavelength < wl[index])
+          break;
+
+      beta0 = aero_betap[index-1]+(aero_betap[index]-aero_betap[index-1])*
+        (wavelength-wl[index-1])/(wl[index]-wl[index-1]);
+
+      /* from km^-1 to cm^-1 : */
+      beta0 *= 1e-5;
+
+      /* Now use the pre-calculated values of the path integral 
+       * for h and theta: */
+
+      path = aerosol_path[(int)floor(0.5+h/1e4)]
+        [(int)MIN(89,floor(0.5+theta/STEPTHETA))];
+
+      T_Mie = exp(-beta0*path);
+
+
+      /* Calculate final transmission coefficient: */
+
+      transmittance = T_Ray * T_Oz * T_Mie;
+
+      break;
+
     }   /*  end of atm switch   */
 
-    return transmittance;
+
+  return transmittance;
+
 }   /*  end of atm  */
 
 int absorption(float wlen, float height, float theta)
-{   int ret = 0;        /*  0: passed, 1: absorbed  */
+{
+  int ret = 0;  /*  0: passed, 1: absorbed  */
    
-    if (RandomNumber > atm(wlen, height, theta)) ret=1;
-
-    return ret; 
+  if (RandomNumber > atm(wlen, height, theta)) ret=1;
+
+  return ret; 
 }   /*  end of absorption  */
+
+
+

trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.h

@@ -5 +5 @@
 T(ATM_NOATMOSPHERE),    /* no atmosphere at all: transmittance = 100% */ \
 T(ATM_90PERCENT),       /* atmosphere with transmittance = 90% */        \
-T(ATM_ISOTHERMAL),      /* isothermal approximation */                   \
 T(ATM_CORSIKA)          /* atmosphere as defined in CORSIKA */
   
@@ -12 +11 @@
 #undef T
 
-#define T(x)  #x        /* define T() as the string of x  */
-    const char *AtmModelNames[] = { ITEM_LIST };
-#undef T
 #undef ITEM_LIST
 
-/*  Strings  */
-#define ATM__NFND__ERR          /*  model name          */ \
-    " *** Atm model \"%s\" not found.\n"
-#define ATM__SET___LOG          /*  model name          */ \
-    "Setting atm model \"%s\".\n"
+/* A. Moralejo, January 2003: added some parameters for Mie scattering and 
+ * Ozone absorption derived from the clear standard atmosphere model in 
+ * "Atmospheric Optics", by L. Elterman and R.B. Toolin, chapter 7 of 
+ * the "Handbook of geophysics and Space environments". S.L. Valley, 
+ * editor. McGraw-Hill, NY 1965.
+ */
+
+/* Aerosol number density for 31 heights a.s.l., from 0 to 30 km, 
+ * in 1 km steps (units: cm^-3)
+ */
+
+float 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};
+
+/* Aerosol scattering coefficient at sea level for the following wavelengths:
+ * 280, 300, 320, 340, 360, 380, 400, 450, 500, 550, 600 and 650 nm
+ * (units: km^-1)
+ */
+
+float 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};
+
+float wl[12] = {280, 300, 320, 340, 360, 380, 400, 450, 500, 550, 600, 650};
+
+/* Ozone concentration for 51 heights a.s.l., from 0 to 50 km, 
+ * in 1 km steps (units: cm/km)
+ */
+
+float 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};
+
+/* Vigroux ozone absorption coefficients (cm-1) */
+
+float 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};
+
 
 #endif

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

@@ -7 +7 @@
 *             Jose Carlos Gonzalez                                  *
 *                                                                   *
+*     Modified December-2002, January-2003                          *
+*     Author: Abelardo Moralejo (moralejo@pd.infn.it)               *
+*                                                                   *
+*     Now this accounts only for Rayleigh scattering. Mie and       *
+*     Ozone absorption are now treated separatedly.                 *
+*                                                                   *
 *********************************************************************
-
-c @T \newpage
-
-
-c @section Source code of {\tt attenu.f}
-
-* @text 
-* Copyright $\copyright$ 1998, Aitor Ibarra Ibaibarriaga
+*
+* December 2002, Abelardo Moralejo: checked algorithms, removed 
+* old/unnecessary code, fixed some bugs, added comments.
+*
+* Fixed bugs (of small influence) in Mie absorption implementation: there were
+* errors in the optical depths table, as well as a confusion: height a.s.l.
+* was used as if it was height above the telescope level. The latter error was
+* also present in the Ozone aborption implementation.
+*
+* On the other hand, now we have the tables AE_ABI and OZ_ABI with optical 
+* depths between sea level and a height h (before it was between 2km a.s.l 
+* and a height h). So that we can simulate also in the future a different 
+* observation level. 
+*
+* AM: WARNING: IF VERY LARGE zenith angle simulations are to be done (say 
+* above 85 degrees, for neutrino primaries or any other purpose) this code 
+* will have to be adapted accordingly and checked, since in principle it has 
+* been written and tested to simulate the absorption of Cherenkov photons 
+* arriving at the telescope from above.
+*
+* AM: WARNING 2: not to be used for wavelengths outside the range 250-700 nm
+*
+* January 2003, Abelardo Moralejo: found error in Ozone absorption treatment.
+* At large zenith angles, the air mass used was the one calculated for 
+* Rayleigh scattering, but since the Ozone distribution is rather different
+* from the global distribution of air molecules, this is not a good 
+* approximation. Now I have left in this code only the Rayleigh scattering,
+* and moved to atm.c the Mie scattering and Ozone absorption calculated in
+* a better way.
+*
+* AM, Jan 2003: added obslev as an argument of the function. Changed the 
+* meaning of the argument height: now it is the true height above sea level 
+* at which a photon has been emitted, before it was the height given by
+* Corsika, measured in the vertical of the observer and not in the vertical 
+* of the emitting particle.
+*
 * @endtext
 c @code
 
-      SUBROUTINE attenu(wavelength, height, theta, tr_atmos)
-
-C----------------------------------------------------------------------C
-C  RHO (DENSITY) F(UNCTION)                                            C
-C                                                                      C
-C  CALCULATES DENSITY (G/CM**3) OF ATMOSPHERE DEPENDING ON HEIGHT (CM) C
-C  (US STANDARD ATMOSPHERE)                                            C
-C  THIS FUNCTION IS CALLED FROM ININKG, UPDATE, CERENE, CERENH         C
-C  ARGUMENT:                                                           C 
-C   ARG    = HEIGHT IN CM                                              C
-C----------------------------------------------------------------------C
-
-      COMMON /ATMOS/   AATM,BATM,CATM,LAHG,RHOF,LONG,OZ_ABI,AE_ABI
-      DOUBLE PRECISION AATM(5),BATM(5),CATM(5),LAHG(5),RHOF(5),Lm(12)
-      DOUBLE PRECISION H,RT, m,OZ_ABI(48,12),AE_ABI(28,12)
-      DOUBLE PRECISION XR, TOT_OZ, TOT_AE, T_Ray,T_Mie, T_Oz,
-     + RHO_TOT, RHO_FI, RHOFP
+      SUBROUTINE attenu(wavelength, height, obslev, theta, tr_atmos)
+
+      REAL wavelength, height, obslev, theta, tr_atmos
+      COMMON /ATMOS/   BATM,CATM,LAHG,LONG
+      DOUBLE PRECISION BATM(4),CATM(4),LAHG(5),Lm(12)
+      DOUBLE PRECISION H,RT, m
+      DOUBLE PRECISION XR, T_Ray, RHO_TOT, RHO_FI
       
-      DOUBLE PRECISION Rcos2, Rsin2
+      DOUBLE PRECISION Rcos2, Rsin2, pi, hscale
       DOUBLE PRECISION x1, x2, x3, x4
       DOUBLE PRECISION e1, e2, e3, e4
 
-      REAL wavelength, height, theta, tr_atmos
-      real trr, trm, tro
       REAL LONG(12)
-c fs: define obervation level
-      double precision obslev      
-      INTEGER I,CON_OZ,CON_MI J, ROW
-
-      DATA OZ_ABI / 
-     + 0.2880000D+00,0.5405000D+00,0.7775000D+00,0.1009000D+01,
-     + 0.1241500D+01,0.1480500D+01,0.1750500D+01,0.2085000D+01,
-     + 0.2514500D+01,0.3087500D+01,0.3864500D+01,0.4817500D+01,
-     + 0.5847500D+01,0.6917500D+01,0.8052500D+01,0.9287499D+01,
-     + 0.1068750D+02,0.1231250D+02,0.1415750D+02,0.1617750D+02,
-     + 0.1827250D+02,0.2034750D+02,0.2232750D+02,0.2414750D+02,
-     + 0.2575750D+02,0.2715250D+02,0.2836750D+02,0.2941100D+02,
-     + 0.3031000D+02,0.3109250D+02,0.3176300D+02,0.3232750D+02,
-     + 0.3281200D+02,0.3323200D+02,0.3358350D+02,0.3387750D+02,
-     + 0.3412650D+02,0.3434000D+02,0.3451900D+02,0.3466250D+02,
-     + 0.3477480D+02,0.3486355D+02,0.3493355D+02,0.3498775D+02,
-     + 0.3503010D+02,0.3506360D+02,0.3509020D+02,0.3511185D+02, 
-     + 0.2740000D-01,0.5140000D-01,0.7395000D-01,0.9600000D-01,
-     + 0.1181500D+00,0.1409000D+00,0.1666000D+00,0.1984500D+00,
-     + 0.2393500D+00,0.2939500D+00,0.3679500D+00,0.4589500D+00,
-     + 0.5573000D+00,0.6593000D+00,0.7673000D+00,0.8848000D+00,
-     + 0.1017800D+01,0.1172300D+01,0.1348300D+01,0.1540800D+01,
-     + 0.1740300D+01,0.1937800D+01,0.2126300D+01,0.2299800D+01,
-     + 0.2453300D+01,0.2586300D+01,0.2702300D+01,0.2801900D+01,
-     + 0.2887550D+01,0.2962050D+01,0.3025900D+01,0.3079700D+01,
-     + 0.3125850D+01,0.3165850D+01,0.3199350D+01,0.3227400D+01,
-     + 0.3251150D+01,0.3271500D+01,0.3288600D+01,0.3302300D+01,
-     + 0.3312995D+01,0.3321445D+01,0.3328110D+01,0.3333270D+01,
-     + 0.3337305D+01,0.3340500D+01,0.3343035D+01,0.1334510D+01,
-     + 0.2435000D-02,0.4570000D-02,0.6575000D-02,0.8535000D-02,
-     + 0.1050500D-01,0.1253000D-01,0.1481500D-01,0.1764500D-01,
-     + 0.2128000D-01,0.2613500D-01,0.3272000D-01,0.4081000D-01,
-     + 0.4957000D-01,0.5866000D-01,0.6827000D-01,0.7875500D-01,
-     + 0.9065500D-01,0.1044050D+00,0.1200050D+00,0.1371050D+00,
-     + 0.1548550D+00,0.1724050D+00,0.1891550D+00,0.2045550D+00,
-     + 0.2182050D+00,0.2300550D+00,0.2403600D+00,0.2492200D+00,
-     + 0.2568350D+00,0.2634550D+00,0.2691300D+00,0.2739150D+00,
-     + 0.2780200D+00,0.2815750D+00,0.2845500D+00,0.2870400D+00,
-     + 0.2891500D+00,0.2909600D+00,0.2924750D+00,0.2936900D+00,
-     + 0.2946425D+00,0.2953940D+00,0.2959865D+00,0.2964455D+00,
-     + 0.2968045D+00,0.2970885D+00,0.2973140D+00,0.2974975D+00,
-     + 0.1740000D-03,0.3265000D-03,0.4695000D-03,0.6090000D-03,
-     + 0.7495000D-03,0.8940000D-03,0.1057000D-02,0.1259000D-02,
-     + 0.1518000D-02,0.1863500D-02,0.2332500D-02,0.2909000D-02,
-     + 0.3533000D-02,0.4180500D-02,0.4865000D-02,0.5610500D-02,
-     + 0.6455500D-02,0.7435000D-02,0.8550000D-02,0.9770000D-02,
-     + 0.1103500D-01,0.1229000D-01,0.1348500D-01,0.1458000D-01,
-     + 0.1555100D-01,0.1639550D-01,0.1713150D-01,0.1776300D-01,
-     + 0.1830600D-01,0.1877800D-01,0.1918200D-01,0.1952250D-01,
-     + 0.1981500D-01,0.2006850D-01,0.2028050D-01,0.2045801D-01,
-     + 0.2060851D-01,0.2073750D-01,0.2084566D-01,0.2093241D-01,
-     + 0.2100026D-01,0.2105381D-01,0.2109606D-01,0.2112876D-01,
-     + 0.2115431D-01,0.2117456D-01,0.2119066D-01,0.2120376D-01,
-     + 0.4885000D-05,0.9170000D-05,0.1319500D-04,0.1713000D-04,
-     + 0.2108000D-04,0.2513500D-04,0.2971500D-04,0.3539500D-04,
-     + 0.4268500D-04,0.5242500D-04,0.6562500D-04,0.8182500D-04,
-     + 0.9937500D-04,0.1175750D-03,0.1368250D-03,0.1578250D-03,
-     + 0.1816250D-03,0.2091750D-03,0.2404750D-03,0.2747750D-03,
-     + 0.3103250D-03,0.3454750D-03,0.3790250D-03,0.4098750D-03,
-     + 0.4372250D-03,0.4609750D-03,0.4816750D-03,0.4994750D-03,
-     + 0.5147750D-03,0.5280750D-03,0.5394750D-03,0.5490700D-03,
-     + 0.5572950D-03,0.5644250D-03,0.5703950D-03,0.5753899D-03,
-     + 0.5796200D-03,0.5832500D-03,0.5862950D-03,0.5887350D-03,
-     + 0.5906400D-03,0.5921450D-03,0.5933350D-03,0.5942565D-03,
-     + 0.5949755D-03,0.5955440D-03,0.5959955D-03,0.5963635D-03,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.0000000D+00,0.0000000D+00,0.0000000D+00,0.0000000D+00,
-     + 0.9525000D-05,0.1785500D-04,0.2567500D-04,0.3332000D-04,
-     + 0.4100000D-04,0.4889000D-04,0.5779500D-04,0.6881000D-04,
-     + 0.8296000D-04,0.1018600D-03,0.1275100D-03,0.1590600D-03,
-     + 0.1932100D-03,0.2286100D-03,0.2660100D-03,0.3067600D-03,
-     + 0.3529600D-03,0.4065100D-03,0.4674100D-03,0.5341100D-03,
-     + 0.6032600D-03,0.6717100D-03,0.7370100D-03,0.7970100D-03,
-     + 0.8501600D-03,0.8963600D-03,0.9366100D-03,0.9711101D-03,
-     + 0.1000810D-02,0.1026660D-02,0.1048810D-02,0.1067460D-02,
-     + 0.1083460D-02,0.1097310D-02,0.1108910D-02,0.1118640D-02,
-     + 0.1126865D-02,0.1133915D-02,0.1139830D-02,0.1144570D-02,
-     + 0.1148275D-02,0.1151200D-02,0.1153510D-02,0.1155300D-02,
-     + 0.1156700D-02,0.1162205D-02,0.1170990D-02,0.1178145D-02,
-     + 0.9360000D-04,0.1757000D-03,0.2528000D-03,0.3281500D-03,
-     + 0.4038500D-03,0.4816500D-03,0.5694500D-03,0.6784000D-03,
-     + 0.8184000D-03,0.1004900D-02,0.1257900D-02,0.1568900D-02,
-     + 0.1905400D-02,0.2254400D-02,0.2623400D-02,0.3025400D-02,
-     + 0.3480900D-02,0.4008900D-02,0.4609400D-02,0.5266900D-02,
-     + 0.5948400D-02,0.6622900D-02,0.7266400D-02,0.7857900D-02,
-     + 0.8381900D-02,0.8836901D-02,0.9233401D-02,0.9573901D-02,
-     + 0.9866901D-02,0.1012140D-01,0.1033940D-01,0.1052340D-01,
-     + 0.1068140D-01,0.1081840D-01,0.1093290D-01,0.1102855D-01,
-     + 0.1110965D-01,0.1117920D-01,0.1123750D-01,0.1128425D-01,
-     + 0.1132085D-01,0.1134975D-01,0.1137255D-01,0.1139020D-01,
-     + 0.1140400D-01,0.1141492D-01,0.1142358D-01,0.1143063D-01,
-     + 0.2500000D-03,0.4690000D-03,0.6745000D-03,0.8755000D-03,
-     + 0.1077500D-02,0.1285000D-02,0.1519500D-02,0.1810000D-02,
-     + 0.2182500D-02,0.2679500D-02,0.3353500D-02,0.4182000D-02,
-     + 0.5079000D-02,0.6010000D-02,0.6994000D-02,0.8064000D-02,
-     + 0.9279000D-02,0.1068900D-01,0.1228900D-01,0.1403900D-01,
-     + 0.1585400D-01,0.1765400D-01,0.1937400D-01,0.2095400D-01,
-     + 0.2235400D-01,0.2356900D-01,0.2462600D-01,0.2553350D-01,
-     + 0.2631400D-01,0.2699250D-01,0.2757350D-01,0.2806300D-01,
-     + 0.2848350D-01,0.2884800D-01,0.2915300D-01,0.2940850D-01,
-     + 0.2962500D-01,0.2981050D-01,0.2996600D-01,0.3009050D-01,
-     + 0.3018780D-01,0.3026480D-01,0.3032550D-01,0.3037250D-01,
-     + 0.3040925D-01,0.3043835D-01,0.3046145D-01,0.3048025D-01,
-     + 0.3585000D-03,0.6725000D-03,0.9675000D-03,0.1256000D-02,
-     + 0.1545500D-02,0.1843000D-02,0.2179000D-02,0.2595500D-02,
-     + 0.3130000D-02,0.3843500D-02,0.4813500D-02,0.6003500D-02,
-     + 0.7288500D-02,0.8623499D-02,0.1003850D-01,0.1157850D-01,
-     + 0.1331850D-01,0.1533350D-01,0.1762850D-01,0.2014350D-01,
-     + 0.2274850D-01,0.2532850D-01,0.2779350D-01,0.3005850D-01,
-     + 0.3206350D-01,0.3380350D-01,0.3531851D-01,0.3661850D-01,
-     + 0.3773851D-01,0.3871351D-01,0.3954751D-01,0.4025051D-01,
-     + 0.4085401D-01,0.4137701D-01,0.4181501D-01,0.4218151D-01,
-     + 0.4249151D-01,0.4275751D-01,0.4298101D-01,0.4316001D-01,
-     + 0.4330001D-01,0.4341061D-01,0.4349771D-01,0.4356516D-01,
-     + 0.4361791D-01,0.4365961D-01,0.4369271D-01,0.4371971D-01,
-     + 0.1685000D-03,0.3160000D-03,0.4545000D-03,0.5900000D-03,
-     + 0.7260000D-03,0.8655000D-03,0.1023000D-02,0.1218500D-02,
-     + 0.1469500D-02,0.1804500D-02,0.2259000D-02,0.2817500D-02,
-     + 0.3422000D-02,0.4049500D-02,0.4712999D-02,0.5434999D-02,
-     + 0.6252999D-02,0.7203000D-02,0.8283000D-02,0.9462999D-02,
-     + 0.1068800D-01,0.1190300D-01,0.1306300D-01,0.1412800D-01,
-     + 0.1507000D-01,0.1588850D-01,0.1660150D-01,0.1721300D-01,
-     + 0.1773900D-01,0.1819650D-01,0.1858850D-01,0.1891850D-01,
-     + 0.1920150D-01,0.1944700D-01,0.1965250D-01,0.1982450D-01,
-     + 0.1997050D-01,0.2009550D-01,0.2020010D-01,0.2028410D-01,
-     + 0.2034985D-01,0.2040175D-01,0.2044265D-01,0.2047435D-01,
-     + 0.2049915D-01,0.2051875D-01,0.2053430D-01,0.2054695D-01 /
-
-      DATA AE_ABI /
-     + 0.3645000E-01,0.5211000E-01,0.5913000E-01,0.6203000E-01,
-     + 0.6304000E-01,0.6340450E-01,0.6353275E-01,0.6358405E-01,
-     + 0.6361780E-01,0.6364885E-01,0.6367990E-01,0.6371365E-01,
-     + 0.6375955E-01,0.6383380E-01,0.6392965E-01,0.6403360E-01,
-     + 0.6414810E-01,0.6426660E-01,0.6438010E-01,0.6448960E-01,
-     + 0.6459510E-01,0.6468170E-01,0.6474110E-01,0.6478295E-01,
-     + 0.6481670E-01,0.6484775E-01,0.6487610E-01,0.6490240E-01,
-     + 0.3510000E-01,0.5018000E-01,0.5694000E-01,0.5973500E-01,
-     + 0.6071000E-01,0.6106100E-01,0.6118450E-01,0.6123390E-01,
-     + 0.6126640E-01,0.6129630E-01,0.6132620E-01,0.6135870E-01,
-     + 0.6140290E-01,0.6147440E-01,0.6156670E-01,0.6166680E-01,
-     + 0.6177730E-01,0.6189180E-01,0.6200130E-01,0.6210680E-01,
-     + 0.6220820E-01,0.6229140E-01,0.6234860E-01,0.6238890E-01,
-     + 0.6242140E-01,0.6245130E-01,0.6247860E-01,0.6250395E-01,
-     + 0.3375000E-01,0.4825000E-01,0.5475000E-01,0.5744000E-01,
-     + 0.5838000E-01,0.5871750E-01,0.5883625E-01,0.5888375E-01,
-     + 0.5891500E-01,0.5894375E-01,0.5897250E-01,0.5900375E-01,
-     + 0.5904625E-01,0.5911500E-01,0.5920375E-01,0.5930000E-01,
-     + 0.5940650E-01,0.5951700E-01,0.5962200E-01,0.5972300E-01,
-     + 0.5982050E-01,0.5990050E-01,0.5995550E-01,0.5999425E-01,
-     + 0.6002550E-01,0.6005425E-01,0.6008050E-01,0.6010490E-01,
-     + 0.3240000E-01,0.4632000E-01,0.5256000E-01,0.5514000E-01,
-     + 0.5604000E-01,0.5636400E-01,0.5647800E-01,0.5652360E-01,
-     + 0.5655360E-01,0.5658120E-01,0.5660880E-01,0.5663880E-01,
-     + 0.5667960E-01,0.5674560E-01,0.5683080E-01,0.5692320E-01,
-     + 0.5702520E-01,0.5713070E-01,0.5723140E-01,0.5732860E-01,
-     + 0.5742220E-01,0.5749900E-01,0.5755180E-01,0.5758900E-01,
-     + 0.5761900E-01,0.5764660E-01,0.5767180E-01,0.5769520E-01,
-     + 0.3240000E-01,0.4632000E-01,0.5256000E-01,0.5514000E-01,
-     + 0.5604000E-01,0.5636400E-01,0.5647800E-01,0.5652360E-01,
-     + 0.5655360E-01,0.5658120E-01,0.5660880E-01,0.5663880E-01,
-     + 0.5667960E-01,0.5674560E-01,0.5683080E-01,0.5692320E-01,
-     + 0.5702520E-01,0.5713070E-01,0.5723140E-01,0.5732860E-01,
-     + 0.5742220E-01,0.5749900E-01,0.5755180E-01,0.5758900E-01,
-     + 0.5761900E-01,0.5764660E-01,0.5767180E-01,0.5769520E-01,
-     + 0.3105000E-01,0.4439000E-01,0.5037000E-01,0.5284500E-01,
-     + 0.5371000E-01,0.5402050E-01,0.5412975E-01,0.5417345E-01,
-     + 0.5420220E-01,0.5422865E-01,0.5425510E-01,0.5428385E-01,
-     + 0.5432295E-01,0.5438620E-01,0.5446785E-01,0.5455640E-01,
-     + 0.5465390E-01,0.5475485E-01,0.5485145E-01,0.5494460E-01,
-     + 0.5503430E-01,0.5510790E-01,0.5515850E-01,0.5519415E-01,
-     + 0.5522290E-01,0.5524935E-01,0.5527350E-01,0.5529590E-01,
-     + 0.2700000E-01,0.3860000E-01,0.4380000E-01,0.4595000E-01,
-     + 0.4670000E-01,0.4697000E-01,0.4706500E-01,0.4710300E-01,
-     + 0.4712800E-01,0.4715100E-01,0.4717400E-01,0.4719900E-01,
-     + 0.4723300E-01,0.4728800E-01,0.4735900E-01,0.4743600E-01,
-     + 0.4752100E-01,0.4760900E-01,0.4769300E-01,0.4777400E-01,
-     + 0.4785200E-01,0.4791600E-01,0.4796000E-01,0.4799100E-01,
-     + 0.4801600E-01,0.4803900E-01,0.4806000E-01,0.4807950E-01,
-     + 0.2430000E-01,0.3474000E-01,0.3942000E-01,0.4135500E-01,
-     + 0.4203000E-01,0.4227300E-01,0.4235850E-01,0.4239270E-01,
-     + 0.4241520E-01,0.4243590E-01,0.4245660E-01,0.4247910E-01,
-     + 0.4250970E-01,0.4255920E-01,0.4262310E-01,0.4269240E-01,
-     + 0.4276890E-01,0.4284810E-01,0.4292370E-01,0.4299660E-01,
-     + 0.4306680E-01,0.4312440E-01,0.4316400E-01,0.4319190E-01,
-     + 0.4321440E-01,0.4323510E-01,0.4325400E-01,0.4327155E-01,
-     + 0.2255000E-01,0.3225500E-01,0.3659500E-01,0.3838900E-01,
-     + 0.3901500E-01,0.3924050E-01,0.3931985E-01,0.3935155E-01,
-     + 0.3937240E-01,0.3939160E-01,0.3941080E-01,0.3943165E-01,
-     + 0.3946005E-01,0.3950600E-01,0.3956530E-01,0.3962960E-01,
-     + 0.3970055E-01,0.3977400E-01,0.3984415E-01,0.3991180E-01,
-     + 0.3997695E-01,0.4003040E-01,0.4006715E-01,0.4009305E-01,
-     + 0.4011390E-01,0.4013310E-01,0.4015065E-01,0.4016695E-01,
-     + 0.2130000E-01,0.3044500E-01,0.3455500E-01,0.3625450E-01,
-     + 0.3684700E-01,0.3706050E-01,0.3713575E-01,0.3716575E-01,
-     + 0.3718550E-01,0.3720370E-01,0.3722190E-01,0.3724165E-01,
-     + 0.3726850E-01,0.3731195E-01,0.3736805E-01,0.3742890E-01,
-     + 0.3749605E-01,0.3756555E-01,0.3763190E-01,0.3769590E-01,
-     + 0.3775750E-01,0.3780805E-01,0.3784280E-01,0.3786725E-01,
-     + 0.3788700E-01,0.3790520E-01,0.3792180E-01,0.3793720E-01,
-     + 0.2025000E-01,0.2895000E-01,0.3285000E-01,0.3446250E-01,
-     + 0.3502500E-01,0.3522750E-01,0.3529875E-01,0.3532725E-01,
-     + 0.3534600E-01,0.3536325E-01,0.3538050E-01,0.3539925E-01,
-     + 0.3542475E-01,0.3546600E-01,0.3551925E-01,0.3557700E-01,
-     + 0.3564075E-01,0.3570675E-01,0.3576975E-01,0.3583050E-01,
-     + 0.3588900E-01,0.3593700E-01,0.3597000E-01,0.3599325E-01,
-     + 0.3601200E-01,0.3602925E-01,0.3604500E-01,0.3605960E-01,
-     + 0.1920000E-01,0.2745500E-01,0.3114500E-01,0.3267050E-01,
-     + 0.3320300E-01,0.3339470E-01,0.3346215E-01,0.3348915E-01,
-     + 0.3350690E-01,0.3352320E-01,0.3353950E-01,0.3355725E-01,
-     + 0.3358140E-01,0.3362045E-01,0.3367085E-01,0.3372550E-01,
-     + 0.3378585E-01,0.3384835E-01,0.3390800E-01,0.3396550E-01,
-     + 0.3402090E-01,0.3406635E-01,0.3409760E-01,0.3411965E-01,
-     + 0.3413740E-01,0.3415370E-01,0.3416860E-01,0.3418245E-01 /
-
-      DATA BATM / 1222.6562D0,1144.9069D0,1305.5948D0,540.1778D0,0.D0  /
-      DATA CATM / 994186.38D0,878153.55D0,636143.04D0,772170.16D0,1.D-9/
-
-      DATA LAHG / 2000.0D2,4.0D5,1.0D6,4.0D6,1.0D7 /
+      INTEGER I
+
+c--   BATM, CATM, LAHG:
+c--   some parameters of the US standard atmosphere (see Corsika manual, 
+c--   appendix C). LAHG contains the limits of the 4 exponential layers, 
+c--   in which the mass overburden is given by T = AATM + BATM * exp(-h/CATM)
+c--   The parameters AATM are not included in this code because they are not 
+c--   needed. The last layer of the US standard atmosphere (in which T varies
+c--   linearly with h) is above 100 km and has not been included here for the
+c--   same reason.
+c--
+      DATA BATM / 1222.6562D0,1144.9069D0,1305.5948D0,540.1778D0 /
+      DATA CATM / 994186.38D0,878153.55D0,636143.04D0,772170.16D0 /
+      DATA LAHG / 0.,4.0D5,1.0D6,4.0D6,1.0D7 /
+
       DATA Lm /3.70D5,3.85D5,4.0D5,4.17D5,4.17D5,4.35D5,5.0D5,5.56D5,
      +         5.99D5,6.33D5, 6.67D5,7.04D5 /
@@ -307 +104 @@
       parameter (hscale = 7.4d5)
 
+c--   Mean free path for scattering Rayleigh XR (g/cm^2)
+      parameter (XR=2.970D3)
+
+c-- Set low limit of first atmospheric layer to the observation level
+c-- so that the traversed atmospheric depth in the Rayleigh scattering 
+c-- will be calculated correctly.
+
+      LAHG(1) = obslev
+
+c-- For the case of simulating a telescope higher than 4 km...
+
+      if (obslev .gt. LAHG(2)) then
+         LAHG(2) = obslev
+      end if
+
+      T_Ray = 1.0
+
+
+c--   AM, Jan 2002: now the argument height is directly the height above 
+c--   sea level, calculated in atm.c:
+
+      h = height
+
 ***********************************************************************
 *
-*     SCATTERING PARAMETERS FOR RAYLEIGH:
-*
-*     MEAN FREE PATH FOR SCATTERING RAYLEIGH  (g/cm^2)
-      PARAMETER (XR=2.970D3)
-***********************************************************************
-        
-c-- Observation level at La Palma
-      parameter (obslev = 2200.d2)
+*     LARGE ZENITH ANGLE FACTOR (AIR MASS FACTOR):
       
-      T_Ray = 1.0
-      T_Mie = 1.0
-      T_Oz = 1.0
-      
-c-- Height calculated using an obslev = H(LaPalma) <> 0
-c      H = -RT + SQRT(RT**2 + (height/COS(theta))**2 +
-c     +     (2.0D0*RT*height))
-      h = -rt + sqrt((rt+obslev)**2 +
-     +     ((height-obslev)/cos(theta))**2 +
-     +     (2.0d0*(rt+obslev)*(height-obslev)))
-
-      ROW = AINT(((H+1.)/1.0E5))      
-              
-***********************************************************************
-*
-*     LARGE ZENITH ANGLE FACTOR (AIR MASS):
-      
-c fs : air mass factor = path_lenght(za) / path_lenght(vertical) 
-c                        at point of emission of cherenkov photon
-c      => pure geometric correction on 
-c         absorption lenght measured at vertical height [km^-1]   
-c--
-c      ma = (1.D0/H)*(SQRT((RT*COS(theta))**2+(2*RT*H)+H**2)-
-c     +     (RT*COS(theta)))
-c--
-c      mb = (1.d0/(h-obslev)) * 
-c     &     ( -(rt+obslev)*cos(theta) 
-c     &     +sqrt( (rt+h)**2 - ((rt+obslev)*sin(theta))**2) )
-c--
-
-c-- Air mass "m" calcualted using a one-exponential density
-c-- profile for the atmosphere, rho = rho_0 exp(-h/Hs)
-c-- with Hs = hscale = 7.4 km
+c--   Air mass factor "m" calculated using a one-exponential density
+c--   profile for the atmosphere, rho = rho_0 exp(-h/hscale) with 
+c--   hscale = 7.4 km. The air mass factor is defined as I(theta)/I(0), 
+c--   the ratio of the optical paths I (in g/cm2) traversed between two 
+c--   given heights, at theta and at 0 deg z.a. 
 
       Rcos2 = rt * cos(theta)**2
       Rsin2 = rt * sin(theta)**2
-      
-      x1 = sqrt((2. * obslev + Rcos2) / (2. * hscale))
-      x2 = sqrt((2. * h      + Rcos2) / (2. * hscale))
-      x3 = sqrt((2. * obslev + rt)    / (2. * hscale))
-      x4 = sqrt((2. * h      + rt)    / (2. * hscale))
-
-c--   AM Dec 2001, to avoid crash! A few photons seem to be "corrupted" 
-c--    (have absurd value) in a cer file...
+
+*     
+*     Obsolete lines:
+*     x1 = sqrt((2. * obslev + Rcos2) / (2. * hscale))
+*     x2 = sqrt((2. * h      + Rcos2) / (2. * hscale))
+*     x3 = sqrt((2. * obslev + rt)    / (2. * hscale))
+*     x4 = sqrt((2. * h      + rt)    / (2. * hscale))
+*     
+
+c--   AM, Dec 2002: slightly changed the calculation of the air mass factor,
+c--   for what I think is a better approximation (numerically the results are 
+c--   almost exactly the same, this change is irrelevant!):
+
+      x1 = sqrt(Rcos2 / (2*hscale))
+      x2 = sqrt((2*(h-obslev) + Rcos2) / (2*hscale))
+      x3 = sqrt(rt / (2*hscale))
+      x4 = sqrt((2*(h-obslev) + rt) / (2*hscale))
+
+c--   AM Dec 2001, to avoid crash! Sometime a few photons seem to be 
+c--   "corrupted" (have absurd values) in cer files... Next lines avoid the
+c--   crash. They will also return a -1 transmittance in the case the photon 
+c--   comes exactly from the observation level, because in that case the 
+c--   "air mass factor" would become infinity and the calculation is not valid!
 
       if (abs(x3-x4) .lt. 1.e-10) then
@@ -368 +168 @@
       endif
 
-
       e1 = derfc(x1)
       e2 = derfc(x2)
@@ -376 +175 @@
       m = exp(-Rsin2 / (2. * hscale)) * ((e1 - e2) / (e3 - e4))
 
-
 **********************************************************************   
-*          
-*     RAYLEIGH SCATTERING
-
-      RHOTOT = 0.0
-      do 11 i=1,5
-        RHOF(i) = 0
- 11   continue
+*     RAYLEIGH SCATTERING (by the molecules in the atmosphere)
+*     SCATTERING M.F.P. FOR RAYLEIGH: XR = 2.970D3
+**********************************************************************
+
+c--   Calculate the traversed "vertical thickness" of air using the
+c--   US Standard atmosphere:
+
+      RHO_TOT = 0.0
 
       DO 91 I=2,5
         IF ( H .LT. LAHG(I) ) THEN
-          RHOTOT = RHOTOT +
+          RHO_TOT = RHO_TOT +
      +        BATM(I-1)*(EXP(-LAHG(I-1)/CATM(I-1) ) -
      +        EXP(-H/CATM(I-1)))
           GOTO 92
         ELSE
-          RHOTOT = RHOTOT +
+          RHO_TOT = RHO_TOT +
      +        BATM(I-1)*(EXP(-LAHG(I-1)/CATM(I-1) ) -
      +        EXP(-LAHG(I)/CATM(I-1)))
@@ -399 +198 @@
  91   CONTINUE
  
-
- 92   RHO_FI = m*RHOTOT
-
-      T_Ray = EXP(-(RHO_FI/XR)*(400.D0/wavelength)**4)
+c--   We now convert from "vertical thickness" to "slanted thickness"
+c--   traversed by the photon on its way to the telescope, simply
+c--   multiplying by the air mass factor m:
+
+ 92   RHO_FI = m*RHO_TOT
+
+c--   Finally we calculate the transmission coefficient for the Rayleigh
+c--   scattering:
+c--   AM Dec 2002, introduced ABS below to account (in the future) for 
+c--   possible photons coming from below the observation level.
+
+      T_Ray = EXP(-ABS(RHO_FI/XR)*(400.D0/wavelength)**4)
    
-
-***************************************************************************
-*                                                                         *
-*     MIE ABSORPTION: WE USE FOR HEIGHTS LOWER THAN 10 Km AN EXACT FORMULA*
-*     FOR THE AEROSOL DENSITY, AND WE USE A TABLE FOR HIGHTS HIGHER THAN  *
-*     10 Km. WE TAKE THE TABLE FROM THE ALTERMAN ARTICLE.                 *
-***************************************************************************
-      
-
-***************************************************************************
-*                                                                         *
-*     MIE SCATTERING PARAMETERS                                           *
-      Hm = 1.2D5
-***************************************************************************
-
-      IF (ROW.GT.27) THEN
-        ROW=28
-      ENDIF
-
-
-        CON_MI = 2
-
- 1001   IF (ABS(LONG(CON_MI)-wavelength).LT.0.01) THEN
-          TOT_AE =AE_ABI(ROW,CON_MI)
-        ELSEIF (LONG(CON_MI).GT.wavelength) THEN
-          A = (AE_ABI(ROW,CON_MI)-AE_ABI(ROW,CON_MI-1))/
-     +        (LONG(CON_MI) - LONG(CON_MI-1))
-          B = AE_ABI(ROW,CON_MI) - (A*LONG(CON_MI))
-          TOT_AE = A*wavelength + B
-        ELSE
-          CON_MI = CON_MI +1
-          GOTO 1001
-        ENDIF
-
-        T_Mie =  EXP(-(m*TOT_AE))           
-
-
-
-***********************************************************************
-*                                                                     *    
-*     OZONE ABSORPTION: We used the Elterman table.                   *
-*                                                                     *
-***********************************************************************
-      IF (ROW.GT.47) THEN
-        ROW = 47
-      ENDIF
-
-      CON_OZ = 2
-
- 2001 IF (LONG(CON_OZ).GT.wavelength) THEN
-             A = (OZ_ABI(ROW,CON_OZ)-OZ_ABI(ROW,CON_OZ-1))/
-     +        (LONG(CON_OZ) - LONG(CON_OZ-1))
-         B = OZ_ABI(ROW,CON_OZ) - (A*LONG(CON_OZ))
-         TOT_OZ = (A* wavelength)+B
-      ELSEIF (ABS(LONG(CON_OZ)-wavelength).LT.0.01) THEN
-         TOT_OZ = OZ_ABI(ROW,CON_OZ)
-      ELSE
-        CON_OZ = CON_OZ + 1
-        GOTO 2001
-      ENDIF
-      
-      T_Oz = EXP(-(m*TOT_OZ))      
-
-************************************************************************
-*                                                                      *
-*     TOTAL TRANSMISSION OF THE ATMOSPHERE: (RAYLEIGH + MIE + OZONE)   *
-************************************************************************
-
-      tr_atmos = T_Ray*T_Mie*T_Oz
+      tr_atmos = T_Ray
 
       RETURN
 
       END
-
-c @endcode
-
-c EOF

trunk/MagicSoft/Simulation/Detector/ReflectorII/doc/Tdas0211.tex

@@ -33 +33 @@
 
 \usepackage{magic-tdas}
+\usepackage{amssymb}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -46 +47 @@
 \author{A.Moralejo\\ 
         \texttt{<moralejo@pd.infn.it>}} 
-\date{December 6, 2002\\}
-\TDAScode{MAGIC-TDAS 02-11\\ 021206/AMoralejo}
+\date{January 20, 2003\\}
+\TDAScode{MAGIC-TDAS 02-11\\ 030120/AMoralejo}
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -61 +62 @@
 been included here for clarity. Two important bugs regarding the
 ray-tracing routine have been corrected in this new version.
+\par
+NOTE: In December 2002, a first release of Reflector 0.6 was made,
+together with a first version of the present TDAS note. Immediately
+after that (too short a time to justify a new version number), some
+other changes were implemented in the program to improve the
+atmospheric absorption routines, and the dependence of mirror
+reflectivity with wavelength was introduced. This document is the
+manual of the final 0.6 version of the reflector, including
+explanations of these last modifications.
+
 \end{abstract}
 
+\newpage
 %% contents %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 \thetableofcontents
@@ -161 +173 @@
 see the resulting spot in fig. \ref{spot_inf_f1697}. A completely
 independent ray-tracing program was used to verify that this is the
-spot that our 17 m tesellated paraboloid should produce.
+spot that a f/1 16.97 m $\varnothing$ tesellated paraboloid should produce.
 %
 \begin{figure}[h]
@@ -182 +194 @@
 possible mirror misalignments and surface irregularities (by the way,
 a feature which somehow optimistically was not included in the
-simulations shown in the proposal). In this way we can check just the
-ray tracing.
+simulations shown in the MAGIC proposal). In this way we can check
+just the ray tracing.
 %
 \begin{figure}[h!]
@@ -242 +254 @@
 intended, and looks more like an error, but anyway we checked that this
 was not the reason for the problems in the ray tracing: a test
-magic.def file was created with all parameters calculated as for a f =
-1697 cm paraboloid, and no significant difference could be seen. That
+magic.def file was created with {\it all} the mirror parameters
+(positions, orientations and radii of curvature)
+calculated as for a f = 1697 cm paraboloid, and no significant
+difference could be seen. That 
 is: i) the individual mirror orientations are the dominant factor, and
 the overall dish in the old reflector behaved indeed like a f = 1697
@@ -288 +302 @@
 determine the particle's trajectory as long as one knows that they
 refer to a downgoing versor, in which case one gets the third
-direction cosine as $w = -\sqrt{u^2+v^2}$, with a minus sign. However,
+direction cosine as $w = -\sqrt{1-u^2-v^2}$, with a minus sign. However,
 in {\it ph2cph.c} we found exactly the opposite (lines 116 to 118 in
 v.05):
@@ -341 +355 @@
 The general transformation between both is a simple rotation,
 since for the sake of simplicity we assume in the simulation that the
-origins always coincide. In Reflector v.05 or older the rotation
+origins always coincide. As we have said, in Reflector v.05 or older the rotation
 matrix was wrong: it had been written assuming that ($\Phi$, $\Theta$)
 indicated the direction towards which the telescope pointed. Actually,
@@ -383 +397 @@
 The bug was certainly present in versions 0.4 and 0.5, but may be even
 older. Nevertheless, there is no doubt that the reflector program used
-for the simulation shown in the MAGIC design report (which must have been
-an early version of the present one) was working fine. Extensive proof
+for the simulation shown in the MAGIC design report was working
+fine. Extensive proof 
 of this is provided in an appendix of the design report. A plausible
 explanation could be that, up to some date, the data being read in by
@@ -395 +409 @@
 file shows no record of any change in this respect, but given that we
 have always used a slightly modified Corsika, it would not be
-surprising if the Cherenkov output was modified at some point. There
-is no documentation on this, so if anyone has any relevant information,
-please make it public.
+surprising if the Cherenkov output was modified at some point in the
+MAGIC version of Corsika (MMCS). There is no documentation on this, so
+if anyone has any relevant information, please make it public.
 %
 \paragraph {Influence of the bug on image analysis\\}
@@ -413 +427 @@
 dramatic effect than the defocusing which the bug was producing.
 %
+\subsection {Performance of ray-tracing in the new version}
+In figure \ref{coma} we show the images of a point-like source at 10
+km from the telescope, produced with the Reflector version 0.6,
+and using the new version of the magic.def file (see see next
+section). No noise has been introduced in the reflection, the observed
+spots are just the result of optical aberrations. The light source has
+been put at slightly different viewing angles
+from the telescope. The results are comparable to those in the design
+report, actually these are a bit better, the difference probably being
+that the focal lengths of the mirror tiles in the older magic.def file
+were discretized in only eight values, while now they change rather
+continuously. Some images of a point source at infinity (a star) can be
+seen in fig. \ref{coma_star}. We can see that for any incidence angle,
+the area within which 50$\%$ of the light is concentrated is smaller
+than that of a small pixel.
+\par
+In figure \ref{refl06images} the images of three gamma events ($\theta
+= 10^\circ$, E = 16, 46, 232 GeV), the same of fig. \ref{evcompare}
+are shown. They have been produced with Reflector 0.6 assuming perfect
+spherical mirrors (left) and realistic ones (right). The images look
+reasonable, much sharper than with the older versions, even when the
+mirror imperfections are taken into account.
+%
+\begin{figure}[h!]
+  \begin{center}
+    \epsfig{file=eps/timing.eps,width=\textwidth}
+    \caption{Test of reflector isochrony. The arrival time
+distributions of photons in the camera are shown for (buggy) Reflector
+0.5 and for Reflector 0.6. The sketch in the center shows the test for
+the case in which the light beam is paralel to the telescope axis
+(left plot). On the right, the same test has been made with light
+arriving 1 degree off axis.
+ \label{timing}}
+  \end{center} 
+\vspace*{-1cm}
+\end{figure}
+%
 \subsection {Second bug: photon timing}
+%
 After a first release of Reflector 0.6 had been made public, another
 important bug was found. We tried to check whether the simulated
@@ -438 +490 @@
 from scratch!).
 %
-\begin{figure}[h]
+\subsection{The new magic.def file}
+A new magic.def file (see sect. \ref{neededfiles}) has been created
+and included in the Reflector 0.6 package. Now the number of
+individual mirror tiles is 956, matching 
+the number and distribution of the final MAGIC design. The mirror
+centers and orientations are those corresponding to a paraboloid of
+1697 cm focal (hence the camera plane is placed at 1700 cm from the
+dish). The focal lengths have been calculated by R. Mirzoyan taking
+into account the so called ``shortening effect'' (see design report).
+A new axisdev.dat file (se again \ref{neededfiles}) with data for the
+956 mirrors is also included.
+%
+\subsection{The new reflectivity.dat file}
+Up to version 0.5 of the program, the reflectivity of the mirrors
+(which the program reads in from the file reflectivity.dat) was
+considered to be equal to $90\%$ for all wavelengths. Following
+measurements performed in Padua of a mirror sample, the reflectivity
+has been found to be between 80 and 90$\%$ in the range from 280 to 650 nm,
+with a dependence on wavelength which is shown in figure \ref{reflec}.
+%
+\begin{figure}[h!]
   \begin{center}
-    \epsfig{file=eps/timing.eps,width=\textwidth}
-    \caption{Test of reflector isochrony. The arrival time
-distributions of photons in the camera are shown for (buggy) Reflector
-0.5 and for Reflector 0.6. The sketch in the center shows the test for
-the case in which the light beam is paralel to the telescope axis
-(left plot). On the right, the same test has been made with light
-arriving 1 degree off axis.
- \label{timing}}
+    \epsfig{file=eps/reflec.eps,width=0.7\textwidth,height=0.4\textwidth}
+    \caption{Reflectivity of the MAGIC telescope mirrors as a
+function of the wavelength of the incident light. The dip around 350
+nm is due to destructive interference of the light reflected on the
+aluminum plate with that reflected on the protective coating of the
+mirror. \label{reflec}}
   \end{center} 
-\vspace*{-1cm}
-\end{figure}
-
-\subsection {Performance of the new version}
-In figure \ref{coma} we show the images of a point-like source at 10
-km from the telescope, produced with the Reflector version 0.6,
-and using the new version of the magic.def file (see see next
-section). No noise has been introduced in the reflection, the observed
-spots are just the result of optical aberrations. The light source has
-been put at slightly different viewing angles
-from the telescope. The results are comparable to those in the design
-report, actually these are a bit better, the difference probably being
-that the focal lengths of the mirror tiles in the older magic.def file
-were discretized in only eight values, while now they change rather
-continuously. Some images of a point source at infinity (a star) can be
-seen in fig. \ref{coma_star}. We can see that for any incidence angle,
-the area within which 50$\%$ of the light is concentrated is smaller
-than that of a small pixel.
-\par
-In figure \ref{refl06images} the images of three gamma events ($\theta
-= 10^\circ$, E = 16, 46, 232 GeV), the same of fig. \ref{evcompare}
-are shown. They have been produced with Reflector 0.6 assuming perfect
-spherical mirrors (left) and realistic ones (right). The images look
-reasonable, much sharper than with the older versions, even when the
-mirror imperfections are taken into account.
-\par
+\end{figure}
 %
 \begin{figure}[p]
@@ -515 +560 @@
 \end{figure}
 %
-\subsection{The new magic.def file}
-A new magic.def file (see sect. \ref{neededfiles}) has been created
-and included in the Reflector 0.6 package. Now the number of
-individual mirror tiles is 956, matching 
-the number and distribution of the final MAGIC design. The mirror
-centers and orientations are those corresponding to a paraboloid of
-1697 cm focal (hence the camera plane is placed at 1700 cm from the
-dish). The focal lengths have been calculated by R. Mirzoyan taking
-into account the so called ``shortening effect'' (see design report).
-A new axisdev.dat file (se again \ref{neededfiles}) with data for the
-956 mirrors is also included.
-%
 \subsection{The {\itshape cermaker} program}
 A test program to produce cer files (input for the reflector)
@@ -539 +572 @@
 file is called {\it cer000001}, and can be read by the reflector program.
 %
+\subsection{Changes in the atmospheric absorption routines}
+%
+In the present version of the program, the simulation of atmospheric
+absorption has undergone major changes with respect to the original
+implementation by J.C. Gonz\'alez and Aitor Ibarra. Three options are
+available in the simulation: ATM\_NOATMOSPHERE, ATM\_90PERCENT and
+ATM\_CORSIKA (see section \ref{commands}). In the 
+first two, 100$\%$ and 90$\%$ of the emitted light respectively
+reaches the telescope mirror, regardless of the emission height. The
+changes in the present version affect only the third option, which is
+the recommended one for running reflector in the standard MAGIC MC
+production. In this case, a detailed estimate of the atmospheric 
+transmission is done. 
+\par
+Three contributions to the atmospheric opacity
+are considered: Rayleigh scattering, Mie scattering by aerosols, and
+absorption by ozone. Details on how the effect of these three
+components is calculated are given in appendix A. In reflector
+versions older than 0.6,
+due to a bug, the contribution of Mie scattering and Ozone absorption
+was slightly overestimated because the vertical height of the emission
+point above sea level was interpreted as height above the
+telescope. The difference is small, since the density of the aerosols
+responsible for Mie scattering decreases very fast with height and
+therefore the absorption is hardly increased by going up 2.2 km in the
+region were most of the Cherenkov light is produced. In the case of
+ozone, its contribution is only important in the very low end of the
+Cherenkov light spectrum, and so the effect of the bug in the total
+amount of light reaching the telescope is negligible. Another change
+is that now the observation level is read in from the cer file, instead
+of being a fixed parameter in the routines. In this way the reflector
+program has become more flexible, and can now be used to process cer
+files produced for a detector at a different observation level.
+\par
+Another problem in the old implementation of absorption was that the
+variation of the optical depth of the atmosphere with the zenith angle
+$\theta$ was assumed to be the same for Mie scattering, ozone
+absorption and Rayleigh scattering: the so called {\it air
+mass} (see appendix A for details) was therefore calculated only once,
+using the overall density profile of air molecules (which does not
+match that of aerosols or ozone), and used to account for the
+variation of all three effects with $\theta$, while rigorously
+speaking it is only valid for Rayleigh scattering. Now, the simulation
+of Mie scattering and ozone absorption has been moved from {\it
+attenu.f} to {\it atm.c} and is done in a more accurate way which
+takes into account the vertical profile of aerosols and ozone (see
+appendix A for details).
+%
 \subsection{Other changes in Reflector 0.6}
 
@@ -557 +638 @@
 section \ref{opt}).
 
-\item New output format (see sect. \ref{out}): added a Run header,
-which is like that of Corsika, plus a couple of variables concerning
-the reflector parameters: the wobble mode and the atmospheric model
+\item New output format (see sect. \ref{out}): we have removed a null
+byte that was written immediately after the ascii label containing the
+reflector program version number at the begining of the file. This
+byte was there for historical reasons and had no function
+whatsoever. Then we have added a Run header, which is like that of
+Corsika, plus a couple of variables concerning the reflector
+parameters: the wobble mode and the atmospheric model 
 used for the simulation. The event header has also been changed to
 include all the information present in the Corsika event
@@ -634 +719 @@
 \end{itemize}
 
-All these files are usually in the {\bf
-MagicProgs/Simulation/Detector/Data/} directory and {\it in principle} you
-should {\bf not} make any change in them to run the program.
+All these files (included in the reflector package) are usually kept
+in the {\bf MagicProgs/Simulation/Detector/Data/} directory and {\it
+in principle} you should {\bf not} make any change in them to run the
+program.
 
 %------------------------------------------------------------
@@ -676 +762 @@
         Wobble-observation mode (see TDAS 01-05 by W. Wittek).
 
-\item[ct\_file ../Data/magic.def]
+\item[ct\_file /path/magic.def]
 
         The ct\_file statement defines where the program can find the 
-        telescope characteristics. The path in the example above is
-        correct to run reflector in
-        MagicProgs/Simu\-la\-tion/De\-tector/Reflector\_0.6/ directory. 
-        If you want to run it in a different directory you have to modify the 
-        path accordingly. 
+        telescope characteristics. Usually, the magic.def file is kept in
+        MagicProgs/Simulation/Detector/Data
 
 \item[atm\_model ATM\_CORSIKA]
-        The atm\_model statement says to the program what kind of
+
+        The atm\_model statement tells the program what kind of
         atmospheric absorption model to use. Possible choices are: 
-        ATM\_CORSIKA, ATM\_ISO\-THERMAL, ATM\_90\-PER\-CENT and
-        ATM\_NO\-ATMO\-SPHE\-RE.
-        
+        ATM\_NOATMOSPHERE,\\ ATM\_90PERCENT and ATM\_CORSIKA,
+        corresponding respectively to no absorption, a 10$\%$
+        absorption and a model using the US Standard atmosphere (see
+        Corsika manual, appendix C) for the Rayleigh scattering and a
+        model by L. Elterman \cite{elterman64,elterman65} for the Mie
+        scattering and ozone absorption (see appendix A). The third
+        model should be chosen for the standard MC MAGIC production.
+
 \item[cer\_files] 
 
@@ -755 +844 @@
         File containing mirror reflectivity as a function of
         wavelength (see section \ref{neededfiles}). If this option is
-        not supplied, the program will look for
+        not supplied, the program will look for \\
         ``../Data/reflectivity.dat'' as previous versions of
         Reflector did.
@@ -796 +885 @@
 %------------------------------------------------------------
 
-\newpage
 \section{Output file \label{out}} 
-The output file begins with two ascii lines:\\
+The reflector output file begins with two ascii lines, the first of
+which informs us of the program version with which it has been
+produced (NOTE: in the following, the dollar symbol \verb"$" stands
+for a carriage return):\\
 \\
-\verb"reflector 0.6" \\
-\verb"START---RUN" \\ 
-After the \verb"START---RUN" flag there is a carriage return, and then
-the run header which is basically the one from Corsika with two added
-variables, {\it wobble\_mode} and {\it atmospheric\_model}. Check the
-Corsika manual for the meaning and units of the rest of them. All of the
-variables 4-byte real numbers except the first, which is a 4 character
-string containing the run header ascii label from Corsika:
+\verb"reflector 0.6$START---RUN$" \\ 
+\\
+Then there is run header which is basically the one from Corsika with
+two added variables, {\it wobble\_mode} and {\it
+atmospheric\_model}. Check the Corsika manual for the meaning and
+units of the rest of them. All of the variables are 4-byte real numbers
+except the first, which is a 4 character string containing the run
+header ascii label from Corsika:
 \vspace*{0.5cm}
 \\
 %
-\begin{tabular}{ll}
-\parbox{5cm}{Variable} & Description \\
+\begin{tabular}{lll}
+\multicolumn{2}{c}{Variable} & Description \\
 \hline 
-& \\
-ASCII Label & 'RUNH' \\
-RunNumber & \\
-date & \\
-Corsika\_version & \\
-NumObsLev & \\
-HeightLev[10] & \\
-SlopeSpec & \\
-ELowLim & \\
-EUppLim & \\
-EGS4\_flag & \\
-NKG\_flag & \\
-Ecutoffh & \\
-Ecutoffe & \\
-Ecutoffg & \\
-C[50] & \\
-wobble\_mode & Wobble mode with which the reflector was run (TDAS
+&& \\
+1 & ASCII Label & 'RUNH' \\
+2 & RunNumber & \\
+3 & date & \\
+4 & Corsika\_version & \\
+5 & NumObsLev & \\
+6 to 15 & HeightLev[10] & \\
+16 & SlopeSpec & \\
+17 & ELowLim & \\
+18 & EUppLim & \\
+19 & EGS4\_flag & \\
+20 & NKG\_flag & \\
+21 & Ecutoffh & \\
+22 & Ecutoffm & \\
+23 & Ecutoffe & \\
+24 & Ecutoffg & \\
+25 to 74 & C[50] & \\
+75 & wobble\_mode & Wobble mode with which the reflector was run (TDAS
 01-05) \\
-atmospheric\_model & Atmospheric model used for the absorption
+76 & atmospheric\_model & Atmospheric model used for the absorption
 simulation \\
-& 0 = no atmosphere;  1 = atm\_90percent;  \\
-& 2 = atm\_isothermal;  3 = atm\_corsika. \\
-dummy1[18] & not used \\
-CKA[40] & \\
-CETA[5] & \\
-CSTRBA[11] & \\
-dummy2[104] & not used \\
-AATM[5] & \\
-BATM[5] & \\
-CATM[5] & \\
-NFL[4] & \\
+&& 0 = no atmosphere;  1 = atm\_90percent;  \\
+&& 2 = atm\_corsika. \\
+77 to 94 & dummy1[18] & not used \\
+95 to 134 & CKA[40] & \\
+135 to 139 & CETA[5] & \\
 &\\
 \hline 
@@ -851 +937 @@
 \newpage
 
-Then there comes a ``\verb"START-EVENT"'' flag, followed by a carriage
-return and then the binary event header. Each variable is a 4-byte
+\begin{tabular}{lll}
+\multicolumn{2}{c}{Variable} & \parbox{11cm}{Description} \\
+\hline 
+&& \\
+140 to 140 & CSTRBA[11] & \\
+151 to 254 & dummy2[104] & not used \\
+255 to 259 & AATM[5] & \\
+260 to 264 & BATM[5] & \\
+265 to 269 & CATM[5] & \\
+270 to 273 & NFL[4] & \\
+&\\
+\hline 
+\end{tabular}
+\vspace*{0.5cm}
+\\
+%
+Then there comes a carriage return followed by the ascii flag which
+indicates the start of an event, and again a carriage return:\\
+\\
+\verb"$START-EVENT$"\\
+\\
+and then the binary event header. Each variable is a 4-byte
 float number except for the first one which is the event header label
 from Corsika (a string of 4 characters). Some of of the variables from
@@ -858 +964 @@
 \vspace*{0.5cm}
 \\
-\begin{tabular}{ll}
-\parbox{5cm}{Variable} & Description \\
+\begin{tabular}{lll}
+\multicolumn{2}{c}{Variable} & Description \\
 \hline
-& \\ 
-ASCII label     & 'EVTH' \\
-EvtNumber       & Event Number \\
-PrimaryID       & Primary particle identification code \\
-Etotal          & Primary particle total energy (GeV) \\
-Thick0          & CORSIKA's starting altitude in g/cm2 \\
-FirstTarget     & CORSIKA's number of first target if fixed \\
-zFirstInt       & Height of first interaction in cm \\
-p[3]            & Primary particle momentum in x,y,-z directions (GeV) \\
-Theta           & Primary particle zenith angle (rad) \\
-Phi             & Primary particle azimuth angle (rad) \\       
-
-NumRndSeq       & Number of different CORSIKA random sequences (max. 10) \\
-RndData[10][3]  & RndData[i][0]: integer seed of sequence i \\
-                & RndData[i][1]: number of offset random calls (mod $10^6$) of sequence i. \\ 
-                & RndData[i][2]: number of offset random calls ($/10^6$) of sequence i. \\ 
-
-RunNumber       & Run number \\
-DateRun         & Date of run yymmdd \\
-Corsika\_version & Version of {\it CORSIKA} \\
-
-NumObsLev       & Number of observation levels (should be always 1 for
+&& \\ 
+1 & ASCII label     & 'EVTH' \\
+2 & EvtNumber       & Event Number \\
+3 & PrimaryID   & Primary particle identification code \\
+4 & Etotal              & Primary particle total energy (GeV) \\
+5 & Thick0              & CORSIKA's starting altitude in g/cm2 \\
+6 & FirstTarget         & CORSIKA's number of first target if fixed \\
+7 & zFirstInt   & Height of first interaction in cm \\
+8 to 10 & p[3]          & Primary particle momentum in x,y,-z directions (GeV) \\
+11 & Theta              & Primary particle zenith angle (rad) \\
+12 & Phi                & Primary particle azimuth angle (rad) \\       
+
+13 & NumRndSeq  & Number of different CORSIKA random sequences (max. 10) \\
+14 to 43 & RndData[10][3]       & RndData[i][0]: integer seed of sequence i \\
+&               & RndData[i][1]: number of offset random calls (mod $10^6$) of sequence i. \\ 
+&               & RndData[i][2]: number of offset random calls ($/10^6$) of sequence i. \\ 
+
+44 & RunNumber  & Run number \\
+45 & DateRun    & Date of run yymmdd \\
+46 & Corsika\_version & Version of {\it CORSIKA} \\
+
+47 & NumObsLev  & Number of observation levels (should be always 1 for
 us) \\
-HeightLev[10]   & Height of observation levels in cm \\
-
-SlopeSpec       & Energy spectrum slope \\
-ELowLim         & Energy lower limit (GeV) \\
-EUppLim         & Energy upper limit (GeV) \\
-Ecutoffh        & \\
-Ecutoffm        & \\
-Ecutoffe        & \\
-Ecutoffg        & \\
-NFLAIN          & \\
-NFLDIF          & \\
-NFLPI0          & \\
-NFLPIF          & \\
-NFLCHE          & \\
-NFRAGM          & \\
-Bx              & \\
-By              & \\
-EGS4yn          & \\
-NKGyn           & \\
-GHEISHAyn       & \\
-VENUSyn         & \\
+48 to 57 & HeightLev[10]        & Height of observation levels in cm \\
+
+58 & SlopeSpec  & Energy spectrum slope \\
+59 & ELowLim    & Energy lower limit (GeV) \\
+60 & EUppLim    & Energy upper limit (GeV) \\
+61 & Ecutoffh        & \\
+62 & Ecutoffm        & \\
+63 & Ecutoffe        & \\
 & \\
 \hline 
@@ -912 +1005 @@
 %
 
-\begin{tabular}{ll}
-
-\parbox{5cm}{Variable} & Description \\
+\begin{tabular}{lll}
+
+\multicolumn{2}{c}{Variable} & \parbox{11cm}{Description} \\
 \hline 
-& \\
-CERENKOVyn      & \\
-NEUTRINOyn      & \\
-HORIZONTyn      & \\
-COMPUTER        & \\
-
-ThetaMin        & Minimum Theta of primaries (deg) \\
-ThetaMax        & Maximum Theta of primaries (deg) \\
-PhiMin          & Minimum Phi of primaries (deg) \\
-PhiMax          & Maximum Phi of primaries (deg) \\
-CBunchSize      & \\
-CDetInX         & \\
-CDetInY         & \\
-CSpacInX        & \\
-CSpacInY        & \\
-CLenInX         & \\
-CLenInY         & \\
-COutput         & \\
-AngleNorthX     & \\
-MuonInfo        & \\
-StepLength      & \\
-CWaveLower      & Wavelength lower limit (nm) \\
-CWaveUpper      & Wavelength upper limit (nm) \\
-Multipl         & \\
-CorePos[2][20]  & Core positions of randomized shower \\
-SIBYLL[2]       & \\
-QGSJET[2]       & \\
-DPMJET[2]       & \\
-VENUS\_cross    & \\
-mu\_mult\_scat  & \\
-NKG\_range      & \\
-EFRCTHN[2]      & \\
-WMAX[2]         & \\
-rthin\_rmax     & \\
-viewcone\_angles[2] & Inner and outer angles of Corsika's VIEWCONE
-option. \\
-telescopePhi    & Telescope azimuth (rad). Measured from South, counter-clockwise \\
-telescopeTheta  & Telescope zenith angle (rad) \\
-TimeFirst       & Arrival time on camera of first photon (ns) \\
-TimeLast        & Arrival time on camera of last photon (ns) \\
-
-& 6 next variables: CORSIKA longitudinal particle fit parameters \\
-& \hspace{0.5cm} (see CORSIKA manual for precise meaning and units)\\
-longi\_Nmax        & Numer of charged particles at maximum \\
-longi\_t0          & Atmospheric depth of shower starting point (N=0) \\
-longi\_tmax        & Atmospheric depth of shower maximum (g/cm$^2$) \\
-longi\_a           & \\
-longi\_b           & For {\bf longi\_a}, {\bf longi\_b}, {\bf longi\_c}, see CORSIKA manual  \\
-longi\_c           &  \\
-longi\_chi2        & $\chi^2/dof$ of the fit\\
+&& \\
+64 & Ecutoffg        & \\
+65 & NFLAIN             & \\
+66 & NFLDIF             & \\
+67 & NFLPI0             & \\
+68 & NFLPIF             & \\
+69 & NFLCHE             & \\
+70 & NFRAGM             & \\
+71 & Bx         & \\
+72 & By         & \\
+73 & EGS4yn     & \\
+74 & NKGyn              & \\
+75 & GHEISHAyn  & \\
+76 & VENUSyn            & \\
+77 & CERENKOVyn & \\
+78 & NEUTRINOyn & \\
+79 & HORIZONTyn & \\
+80 & COMPUTER   & \\
+81 & ThetaMin   & Minimum Theta of primaries (deg) \\
+82 & ThetaMax   & Maximum Theta of primaries (deg) \\
+83 & PhiMin     & Minimum Phi of primaries (deg) \\
+84 & PhiMax     & Maximum Phi of primaries (deg) \\
+85 & CBunchSize & \\
+86 & CDetInX    & \\
+87 & CDetInY    & \\
+88 & CSpacInX   & \\
+89 & CSpacInY   & \\
+90 & CLenInX    & \\
+91 & CLenInY    & \\
+92 & COutput    & \\
+93 & AngleNorthX& \\
+94 & MuonInfo   & \\
+95 & StepLength & \\
+96 & CWaveLower & Wavelength lower limit (nm) \\
+97 & CWaveUpper & Wavelength upper limit (nm) \\
+98 & Multipl    & \\
+99 to 138 & CorePos[2][20] & Core positions of randomized shower \\
+139 to 140 & SIBYLL[2]  & \\
+141 to 142 & QGSJET[2]  & \\
+143 to 144 & DPMJET[2]  & \\
+145 & VENUS\_cross      & \\
+146 & mu\_mult\_scat    & \\
+147 & NKG\_range        & \\
+148 to 149 & EFRCTHN[2] & \\
+150 to 151 & WMAX[2]    & \\
 & \\
 \hline 
@@ -972 +1060 @@
 \newpage
 
-\begin{tabular}{ll}
-\parbox{5cm}{Variable} & Description \\
+\begin{tabular}{lll}
+\multicolumn{2}{c}{Variable} & \parbox{11cm}{Description} \\
 \hline 
-& \\
-CORSIKAPhs      & Original photons written by {\it CORSIKA} \\ 
-AtmAbsPhs       & Photons absorbed by the atmosphere  \\
-MirrAbsPhs      & Photons absorbed by the mirror      \\ 
-OutOfMirrPhs    & Photons outside the mirror          \\ 
-BlackSpotPhs    & Photons lost in the "black spot"    \\ 
-OutOfChamPhs    & Photons outside the camera         \\ 
-CPhotons        & Photons reaching the camera        \\ 
-
-elec\_cph\_fraction & Fraction of C-photons produced by electrons \\
-muon\_cph\_fraction & Fraction of C-photons produced by muons \\
-other\_cph\_fraction & Fraction of C-photons produced by electrons \\
+&& \\
+152 & rthin\_rmax       & \\
+153 to 154 & viewcone\_angles[2] & Inner and outer angles of Corsika's VIEWCONE
+option. \\
+155 & telescopePhi      & Telescope azimuth (rad). Measured from South, counter-clockwise \\
+156 & telescopeTheta    & Telescope zenith angle (rad) \\
+157 & TimeFirst         & Arrival time on camera of first photon (ns) \\
+158 & TimeLast          & Arrival time on camera of last photon (ns) \\
+
+&& 6 next variables: CORSIKA longitudinal particle fit parameters \\
+&& \hspace{0.5cm} (see CORSIKA manual for precise meaning and units)\\
+159 & longi\_Nmax        & Numer of charged particles at maximum \\
+160 & longi\_t0          & Atmospheric depth of shower starting point (N=0) \\
+161 & longi\_tmax        & Atmospheric depth of shower maximum (g/cm$^2$) \\
+162 & longi\_a           & \\
+163 & longi\_b           & For {\bf longi\_a}, {\bf longi\_b}, {\bf longi\_c}, see CORSIKA manual  \\
+164 & longi\_c           &  \\
+165 & longi\_chi2        & $\chi^2/dof$ of the fit\\
+166 & CORSIKAPhs      & Original photons written by {\it CORSIKA} \\ 
+167 & AtmAbsPhs       & Photons absorbed by the atmosphere  \\
+168 & MirrAbsPhs      & Photons absorbed by the mirror      \\ 
+169 & OutOfMirrPhs    & Photons outside the mirror          \\ 
+170 & BlackSpotPhs    & Photons lost in the "black spot"    \\ 
+171 & OutOfChamPhs    & Photons outside the camera         \\ 
+172 & CPhotons        & Photons reaching the camera        \\ 
+
+173 & elec\_cph\_fraction & Fraction of C-photons produced by electrons \\
+174 & muon\_cph\_fraction & Fraction of C-photons produced by muons \\
+175 & other\_cph\_fraction & Fraction of C-photons produced by electrons \\
+176 to 182 & dummy[7] & not used \\
 & \\
 \hline 
 \end{tabular}
-
-\vspace*{1cm}
-The event header is followed by 8-word blocks, one for each photon
-that reaches the camera. A photon block contains the following
-variables:
+%
+\vspace*{0.5cm}
+\\
+The event header is followed by 8-word blocks, one for each Cherenkov
+photon that reaches the camera. A photon block contains the following
+variables (as 4-byte float numbers):
 \vspace*{0.5cm}
 \\
@@ -1006 +1113 @@
         & n is the index from 1 to ICERML (see Corsika manual) for the case in which each Corsika \\
         & shower is used more than once (normally, in MMCS will be just 1). \\
-        & \\
 x, y    & Impact point in camera coordinates (cm) \\
 u, v    & Director cosines of down-going versor indicating the photon direction \\
-t       & Arrival time on camera (ns) \\
-h       & Production height (cm) \\
+t       & Arrival time on camera (ns), measured from the time of first
+interaction of the primary \\
+h       & Production height (cm), measured above sea level on the
+vertical of the telescope location \\
+        & (it is not the {\it true} height which would be measured on
+the vertical of the emitting particle!)  \\ 
 phi     & Incidence angle with respect to camera plane (rad) \\
 & \\
@@ -1017 +1127 @@
 \vspace*{0.5cm}
 \\
-After the last event photon block there is a blank line, an \verb$END---EVENT$
-flag, another blank line and then the following event. After the last
-event in a run it appears the flag ``\verb$END-----RUN$'', while after all
-the processed runs, a ``\verb$END----FILE$'' flag is written. Finally,
-after this flag an ``ascii tail'' has been added to the file: it
-consists of the ascii files {\it magic.def}, {\it axisdev.dat} and {\it
-reflectivity.dat} one after the other separated by blank lines. In
-this way all the relevant parameters used to produce the output are
+After the last photon block of an event there is a carriage return
+followed by the ascii flag indicating the event end, and then two more
+carriage returns before the ascii flag of the beginning of the next
+event, and so on:\\
+\\
+\verb"$END---EVENT$$START-EVENT$Event_header....$END---EVENT$$END-----RUN$$START---RUN$..."\\
+\\
+The flag ``\verb$END-----RUN$'' appears after the last event in a run
+(that is, the last event processed of each of the input cer
+files). After the last processed run, an end of file flag is
+written:\\
+\\
+\verb"...$END---EVENT$$END-----RUN$$END----FILE$"\\
+\\
+Finally, after this flag an ``ascii tail'' has been attached to the file:
+it consists of the ascii files {\it magic.def}, {\it axisdev.dat} and
+{\it reflectivity.dat} one after the other, separated by carriage returns:
+\\
+\\ 
+\verb"$magic.def$axisdev.dat$reflectivity.dat"\\
+\\
+In this way all the relevant parameters used to produce the output are
 kept together with the reflected events.
-
 %------------------------------------------------------------
-
-\section{Appendix A}
+\newpage
+\renewcommand{\thesubsection}{A.\arabic{subsection}}
+\section*{Appendix A : atmospheric absorption}
+\addcontentsline{toc}{section}{Appendix A : atmospheric absorption}
+%
+The simulation of the absorption of Cherenkov light in the atmosphere
+has been included in the {\it Reflector} program because this feature
+was not yet available in the first versions of CORSIKA used within the
+MAGIC collaboration. In the latest CORSIKA versions, the atmospheric
+absorption has been included as an option, but it is not
+compatible with the simulation of a curved atmosphere \cite{cor02},
+and hence we have kept this step as a part of our reflector
+simulation. This appendix describes how the atmospheric absorption is
+implemented in the program when the option ATM\_CORSIKA (see section
+\ref{commands}) is chosen.
+\par
+The geometry of the problem is sketched in figure
+\ref{fig:atmoscheme}. A Cherenkov photon is emitted in point A and
+travels towards the telescope placed at B. At any moment, the height
+$h$ of the photon above sea level is related to the distance $L$
+between the photon and the telescope through
+%
+\begin{equation}
+(R+h)^2 = (R+h_1)^2 + L^2 + 2 L \; (R+h_1) \; \cos \theta
+\label{eq:height}
+\end{equation}
+%
+, where $R$ is the Earth radius, $h_1$ the height (a.s.l.) of the
+observation level and $\theta$ is the zenith angle of the photon
+trajectory measured at the telescope site. The Cherenkov output of
+CORSIKA contains for each photon the height $h_C$ of the emission
+point A, measured along the vertical of the observer. The {\it true
+vertical height} $h_2$ of the emission point can be obtained by
+replacing $L$ by $(h_C-h_1)/\cos \theta$ in equation
+(\ref{eq:height}).
+%
+\begin{figure}[ht]
+\begin{center}
+\mbox{ \epsfig{file=eps/atmoscheme.eps,width=0.8\textwidth} }
+\end{center}
+\caption[]
+{Calculation of the true vertical height $h_2$ of the emission point
+of a Cherenkov photon (point A), and the optical path traversed down to the
+telescope (point B).}
+\label{fig:atmoscheme}
+\end{figure}
+%
+\par
+The optical path $I(\theta, h_1, h_2)$ traversed by the photon can be
+calculated integrating the air density along the trajectory
+$\overline{\text{AB}}$. For the case $h/R \ll 1$, we can drop in
+(\ref{eq:height}) the terms in $(h/R)^2$ and smaller, and then solve
+for L. Deriving now with respect to $h$, we have:
+%
+\begin{equation}
+\frac{dL}{dh} \simeq \sqrt{\frac{R}{2\;(h-h_1)+R\;\cos^2 \theta}}
+\qquad \text{for} \qquad \frac{h}{R} \ll  1
+\label{eq:dldh}
+\end{equation}
+%
+\subsection{Rayleigh scattering}
+Rayleigh scattering is the scattering of light by particles smaller
+than its wavelength. These are in our case the air molecules. The
+transmission coefficient due to Rayleigh scattering is a strong
+function of the wavelength $\lambda$:
+%
+\begin{equation}
+T_{\text{Rayl}} (\lambda) = \exp \Biggl[ - \frac{I(\theta, h_1, h_2)}{x_R} \; \Biggl(\frac{400 \;
+\text{nm}}{\lambda}\Biggl)^4 \; \Biggl] 
+\label{eq:rayleigh}
+\end{equation}
+%
+Here $I(\theta, h_1, h_2)$ is the optical path (in g/cm$^2$) traversed
+between points A and B, and $x_R = 2970$ g/cm$^2$ is the mean free path of
+the Rayleigh scattering at $\lambda = 400$ nm. 
+\par
+A convenient way of expressing the optical path is the following:
+%
+\begin{equation}
+I\;(\theta, h_1, h_2) = (x_1 - x_2) \cdot \mathcal{AM}\;(\theta, h_1, h_2) 
+\end{equation}
+Here $x_{i=1,2}$ is the mass overburden of the atmosphere above a
+height $h_i$ (in the 
+vertical direction) and $\mathcal{AM}$ is the so called {\it air
+mass}\footnote{If we set $h_2 = \infty$, we have the usual definition
+of {\it air mass} in optical astronomy.}, defined as
+%
+\begin{equation}
+\mathcal{AM} \equiv \frac{I\;(\theta,h_1,h_2)}{I\;(0^\circ,h_1,h_2)}
+\label{eq:airmass}
+\end{equation}
+%
+, which is a function mainly of the zenith angle $\theta$ (see
+fig. \ref{fig:airmass}). In our simulation, for the calculation of the
+mass overburden $x_i$ we have used the U.S. standard atmosphere
+parametrized by J. Linsley \cite{cor02}, the same we used in Corsika
+for the shower development simulation. It consists of five layers: in
+the lower four the density decreases exponentially with height, and in
+the upper one the mass overburden cecreases linearly until it vanishes
+at $h = 112.8$ km.
+%
+\begin{figure}[ht]
+\begin{center}
+\mbox{ \epsfig{file=eps/airmass.eps,width=0.8\textwidth} }
+\end{center}
+\caption[]
+{Dependence with zenith angle of the air mass as defined in the
+text. The air mass has been calculated for an exponential atmosphere
+with scale height $H = 7.4$ km, for the observation level of MAGIC
+(2.2 km a.s.l.), and for light coming from $h_2 = 10$ and at 100 km a.s.l. As
+we can see the dependence with the emission height $h_2$ is very small.}
+\label{fig:airmass}
+\end{figure}
+%
+\par
+For the estimate of $\mathcal{AM}$, a simpler atmospheric model has
+been used, in which the vertical density profile is described by a
+single exponential: $\rho = \rho_0 \; e^{-h/H}$  with scale height $H
+= 7.4$ km. This simplifies the calculations, and is accurate enough
+for our purposes. Using (\ref{eq:dldh}) the optical path $I(\theta,
+h_1, h_2)$ can then be obtained approximately as:
+%
+\begin{equation}
+I(\theta, h_1, h_2) = \int_A^B \rho\;(h)\; \frac{dL}{dh} \; dh \simeq
+\sqrt{\frac{R}{2}} \;
+\int_{h_1}^{h_2} \frac{\rho_0
+\;e^{-h/H}}{\sqrt{h-h_1+\frac{1}{2}\;R\;\cos^2 \theta}} \; dh 
+\label{eq:optpath}
+\end{equation}
+%
+and finally:
+% 
+\begin{equation}
+\mathcal{AM} \simeq e^{-\frac{R \sin^2 \theta}{2H}}
+\cdot 
+\frac{
+\text{erfc}\;(\sqrt{\frac{R \cos^2 \theta}{2H}})
+\;-\;
+\text{erfc}\;(\sqrt{\frac{2 (h_2 - h_1) + R \cos^2 \theta}{2H}})
+}
+{
+\text{erfc}\;(\sqrt{\frac{R}{2H}})
+\;-\;
+\text{erfc}\;(\sqrt{\frac{2(h_2 - h_1) + R}{2 H}})
+}
+\label{eq:airmass2}
+\end{equation}
+%
+where we have used the complementary error function $\text{erfc}\;(x)
+= \frac{2}{\sqrt{\pi}}\int_x^\infty e^{-t^2} dt$. From 
+(\ref{eq:airmass2}), $\mathcal{AM}$ can be readily evaluated for any
+value of $\theta$, $h_1$ and $h_2$. In fig. \ref{fig:rayleigh} the
+resulting Rayleigh transmission coefficient $T_{\text{Rayl}}$ finally
+obtained from (\ref{eq:rayleigh}) is plotted versus the zenith angle
+for three wavelengths.
+%
+\begin{figure}[ht]
+\begin{center}
+\mbox{ \epsfig{file=eps/rayleigh.eps,width=\textwidth} }
+\end{center}
+\caption[]
+{Rayleigh transmission coefficient as a function of zenith angle for three
+different wavelengths. The solid, dashed and dotted lines correspond
+respectively to light coming from 5, 10 and 20 km distance from the
+telescope.}
+\label{fig:rayleigh}
+\end{figure}
+%
+\subsection{Mie scattering}
+Cherenkov light also suffers Mie scattering through interaction with
+small dust particles suspended in the air (aerosols), whose size is
+comparable to the wavelength of the light. In our simulation of the
+attenuation due to aerosols we have used the model proposed by
+Elterman \cite{elterman64,elterman65}, which considers an aerosol
+number density $N_p$ which (roughly) decreases exponentially up to 10
+km a.s.l. with scale height $H \simeq 1.2$ km, followed by a more
+tenuous layer between 10 and 30 km (see fig. \ref{fig:aerosol}a). In
+this model, the aerosol size distribution is considered to be
+unchanged with altitude.
+%
+\begin{figure}[ht]
+\begin{center}
+\mbox{ \epsfig{file=eps/aerosol.eps,width=0.9\textwidth} }
+\end{center}
+\caption[]
+{Aerosol model by Elterman: in (a), number density of
+aerosols as a function of height above sea level; in (b), aerosol
+attenuation coefficient at sea level as a function of wavelength.}
+\label{fig:aerosol}
+\end{figure}
+%
+\par
+Measured values of the aerosol attenuation coefficients at sea level
+$\beta_p(0)$ for different wavelengths \cite{elterman65} are shown in
+figure \ref{fig:aerosol}b. To obtain the attenuation coefficient at a
+given height $h$, we simply do $\beta_p(h, \lambda) = \beta_p(0,
+\lambda) \cdot N_p(h)/N_p(0)$. In the Elterman model the aerosol
+transmission coefficient for the trajectory from A to B depicted in
+figure \ref{fig:atmoscheme} would then be: 
+%
+\begin{equation}
+T_{\text{Mie}}(\lambda) = e^{-\tau_{mie}} \quad \text{, with}\quad
+\tau_{mie}(h_1, h_2, \theta, \lambda) = \frac{\beta_p(0, \lambda)}{N_p(0)} \;
+\int_{h_1}^{h_2} 
+\; N_p(h) \;\; \frac{dL}{dh}\; dh
+\label{eq:aerosoltau}
+\end{equation}
+%
+Here $\tau_{mie}$ is the aerosol optical depth of the path from A to
+B. Given that the aerosol density distribution is not a simple
+exponential, we have to do the integral in (\ref{eq:aerosoltau})
+numerically. The integral depends on $h_2$ and on $\theta$, through
+$dL/dh$ (it also 
+depends on $h_1$, but this is the observation level and therefore
+fixed). In this case we use the exact expression for $\;dL/dh\;$ which can be
+obtained from (\ref{eq:height}). At the beginning of each simulation
+run we calculate and store the results of the integral for values of
+$\theta$ between 0 and 90$^\circ$ (in steps of $1^\circ$), and of
+$h_2$ from $h_1$ up to 30 km (in steps of 100 m). To do the integral a
+linear interpolation has been used to obtain the value of $N_p$ for
+any height. During the simulation of each Cherenkov photon, we get the
+corresponding precalculated value of the integral and deduce
+$T_{\text{Mie}}$ from expression (\ref{eq:aerosoltau}).
+%
+\begin{figure}[ht]
+\begin{center}
+\mbox{ \epsfig{file=eps/mie.eps,width=\textwidth} }
+\end{center}
+\caption[]
+{Aerosol transmission coefficient for three different wavelengths as a
+function of the distance between the photon emission point and the
+telescope. Plots for four different zenith angles between 0 and 80
+degrees are shown.}
+\label{fig:mie}
+\end{figure}
+%
+\par
+Since in this model the aerosols are concentrated mainly at very low
+altitude, the transmission coefficient is more or less constant above
+a certain height (which depends on $\theta$), as can be seen in
+fig. \ref{fig:mie}. For instance, for vertically incident 300 nm light
+emitted higher than 4 km above the telescope, the Mie transmission is
+about 0.95.
+%
+\subsection{Absorption by Ozone}
+%
+The absorption of Cherenkov light by Ozone has been implemented
+following also the Elterman standard atmosphere \cite{elterman65}. The
+coefficient for ozone absorption is given by
+%
+\begin{equation}
+\beta_3(h,\lambda) = A_v(\lambda) \cdot D_3(h)
+\end{equation}
+%
+, where $A_v(\lambda)$ is the Vigroux \cite{vigroux53} ozone
+absorption coefficient (cm$^{-1}$) and $D_3(h)$ is the ozone
+concentration (cm km$^{-1}$) according to Elterman. The transmission
+coefficient of Ozone in the path $\overline{\text{AB}}$ is then:
+%
+\begin{equation}
+T_{\text{Ozone}}(\lambda) = e^{-\tau_{oz}} \quad \text{, with}\quad
+\tau_{oz}(h_1, h_2, \theta, \lambda) = A_v(\lambda) \; \int_{h_1}^{h_2} 
+\; D_3(h) \;\; \frac{dL}{dh}\; dh
+\label{eq:ozonetau}
+\end{equation}
+%
+\begin{figure}[ht]
+\begin{center}
+\mbox{ \epsfig{file=eps/ozone.eps,width=\textwidth} }
+\end{center}
+\caption[]
+{Ozone concentration vertical profile (a) and Vigroux coefficients for
+Ozone absorption (b). The Vigroux coefficients for $\lambda =$ 380 and
+400 nm are zero.}
+\label{fig:ozone}
+\end{figure}
+%
+\par
+Once again, like in the case of Mie scattering, the optical depth
+$\tau_{oz}$ is the product of a factor which depends on $\lambda$ and
+an integral which depends on $h_1$, $h_2$ and $\theta$. We proceed in
+the same way as before, precalculating the values of the integral in
+steps of $\Delta\theta = 1^\circ$ and $\Delta h = 100$ m, up to a
+height of 50 km a.s.l. (where the ozone concentration becomes
+negligible), and then reading the appropriate values for every
+simulated photon.
+\par
+Finally, the overall atmospheric transmission coefficient is
+calculated as
+%
+\begin{equation}
+T_{total} = T_{Ray} \cdot T_{Mie} \cdot T_{Ozone}
+\end{equation}
+%
+In figure \ref{fig:absorplot} the atmospheric transmission as a
+function of the distance to the telescope for $\theta = 0^\circ$ is
+shown. Ozone absorption turns out to be dominant below 320 nm, while
+Rayleigh scattering is the main cause of the loss of photons at longer
+wavelengths.
+%
+\begin{figure}[ht]
+\begin{center}
+\mbox{ \epsfig{file=eps/absorplot.eps,width=\textwidth} }
+\end{center}
+\caption[]
+{Total transmission coefficient for vertically incident light as a
+function of the distance between the emission point and the
+telescope. The contributions of absorption by ozone and of Rayleigh
+and Mie scattering are also shown for comparison.}
+\label{fig:absorplot}
+\end{figure}
+%
+\newpage
+\section*{Appendix B : files in the reflector package}
+\addcontentsline{toc}{section}{Appendix B : files in the reflector package}
 
 The list of all Reflector files follows.
@@ -1060 +1496 @@
 MagicProgs/Simulation/Detector/Reflector_0.6/version.h
 
-MagicProgs/Simulation/Detector/Reflector_0.6/doc/Tdas0211.ps.gz
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/Tdas0211.ps
 MagicProgs/Simulation/Detector/Reflector_0.6/doc/Tdas0211.tex
 MagicProgs/Simulation/Detector/Reflector_0.6/doc/magic-tdas.sty
 MagicProgs/Simulation/Detector/Reflector_0.6/doc/magiclogo.eps
-MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/colimation.eps.gz
-MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/coma.eps.gz
-MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/coma_star.eps.gz
-MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/coorsystems.eps.gz
-MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/evcompare.eps.gz
-MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/parabola.eps.gz
-MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/refl06images.eps.gz
-MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/spot10kmf1694.eps.gz
-MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/spot10kmf1700.eps.gz
-MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/spot_inf_f1697.eps.gz
-MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/telecoor.eps.gz
-MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/timing.eps.gz
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/absorplot.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/aerosol.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/airmass.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/atmoscheme.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/colimation.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/coma.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/coma_star.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/coorsystems.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/evcompare.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/mie.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/ozone.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/parabola.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/rayleigh.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/refl06images.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/reflec.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/spot10kmf1694.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/spot10kmf1700.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/spot_inf_f1697.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/telecoor.eps
+MagicProgs/Simulation/Detector/Reflector_0.6/doc/eps/timing.eps
 
 MagicProgs/Simulation/Detector/Reflector_0.6/tester/Makefile
@@ -1096 +1540 @@
 %%>>>> Or the following if you include here by hand your 
 %%>>>> bibliographic entries
+
+\begin{thebibliography}{00}
+
+\bibitem{elterman64}
+L. Elterman, Applied Optics Vol. 3, No. 6 (1964) 745.
+
+\bibitem{elterman65}
+L. Elterman, R.B. Toolin, S.L. Valley (editor), Handbook of
+geophysics and space environments, McGraw-Hill, N.Y. (1965).
+
+\bibitem{cor02}D. Heck and J. Knapp, EAS Simulation with CORSIKA: A User's
+Manual, 2002.
+
+\bibitem{vigroux53}
+E. Vigroux, Contributions \`a l'étude expérimentale de l'absorption
+de l'ozone, Annales de Physique, v. 8 (1953) 709.
+
+\end{thebibliography}
 
 \end{document}

trunk/MagicSoft/Simulation/Detector/ReflectorII/init.c

@@ -113 +113 @@
     fseek(rflfile, 0L, SEEK_SET);
     fprintf(rflfile, "%s %s", QUOTE(PROGRAM), QUOTE(VERSION));
-    fputc(0, rflfile);
+
+    /* Removed by A.M. January 2003: */
+    /*    fputc(0, rflfile); */
+
     fwrite(FLAG_START_OF_RUN, SIZE_OF_FLAGS, 1, rflfile);
 

trunk/MagicSoft/Simulation/Detector/ReflectorII/init.h

@@ -102 +102 @@
     "        - D. Bastieri & C. Bigongiari - Jan 2000\n" \
     "    Version 0.6\n" \
-    "        - A. Moralejo - October 2002\n" \
+    "        - A. Moralejo - January 2003\n" \
     "  ################################################\n\n"
 #define SIGN_ERROR_FTL          /*  program, version    */ \

trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.c

@@ -12 +12 @@
 char axisdev_filename[256], reflectivity_filename[256];
 char geo_filename[256];
+char atmosphere[256]; /* atmospheric model */
 
 /*  Prototypes  */
@@ -67 +68 @@
     extern FILE *geofile;       /*  geo file (init)     */
     extern void SetVerbose(int vlevel); /*  from diag.c */
-    extern void SetAtmModel(char *model); /*  from atm.c        */
     extern int ParseLine(FILE *parfile, /*  FILE with parms     */
                          const char *token_list[], /*  array w/tokens   */
@@ -94 +94 @@
                 break;
               case atm_model:
-                SetAtmModel(value_ptr);
+                strcpy(atmosphere, value_ptr);
                 break;
               case verbose_level:

trunk/MagicSoft/Simulation/Detector/ReflectorII/reflector.c

@@ -47 +47 @@
 extern CerEventHeader *cheadp;  /*  var inited in header.c      */
 extern RflEventHeader *rheadp;  /*  var inited in header.c      */
+extern void SetAtmModel(int model, float ol); /*  from atm.c    */
 
 /*  Prototypes  */
@@ -470 +471 @@
   RflRunHeader RflRunHead;
 
-  extern int atmModel;      /* current atm. model */
+  const char *AtmModelNames[]={"ATM_NOATMOSPHERE","ATM_90PERCENT","ATM_CORSIKA"};
+  extern char atmosphere[256];  /* current atm. model */
+  int model=0;
 
   do
@@ -490 +493 @@
           if (newFile)
             {
+              while (strcmp(atmosphere, AtmModelNames[model]))
+                if (++model == sizeof(AtmModelNames)/sizeof(AtmModelNames[0]))
+                  {   
+                    model = 0;
+                    Error(" *** Atm model \"%s\" not found.\n", atmosphere);
+                    break;      
+                  }
+
               /* Write Reflector "run header" (one per cer file!): */
               memcpy(&RflRunHead, cheadp, sizeof(RflRunHead));
               RflRunHead.wobble_mode = wobble_position;
-              RflRunHead.atmospheric_model = atmModel;
+              RflRunHead.atmospheric_model = model;
               fwrite(&RflRunHead, sizeof(RflRunHead), 1, rflfile);
+
+              /* Set up atmosphere (do some necessary calculations) once we 
+               * read in the observation level from the run header: 
+               */
+
+              Log("Setting atm model \"%s\" and observation level %.0f cm.\n", AtmModelNames[model], RflRunHead.HeightLev[0]);
+              SetAtmModel(model,RflRunHead.HeightLev[0]);
+
             }