Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/Changelog =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/Changelog (revision 1720) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/Changelog (revision 1722) @@ -1,3 +1,33 @@ ** 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,5 +45,5 @@ 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. Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/Makefile =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/Makefile (revision 1720) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/Makefile (revision 1722) @@ -19,7 +19,7 @@ # # $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,5 +217,5 @@ 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,7 +229,11 @@ 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 Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c (revision 1720) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c (revision 1722) @@ -1,61 +1,255 @@ +/******************************************************************** +* * +* 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 #include #include + +#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 */ + + + Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.h (revision 1720) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.h (revision 1722) @@ -5,5 +5,4 @@ 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,14 +11,38 @@ #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 Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/attenu.f =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/attenu.f (revision 1720) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/attenu.f (revision 1722) @@ -7,281 +7,78 @@ * 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,59 +104,62 @@ 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,5 +168,4 @@ endif - e1 = derfc(x1) e2 = derfc(x2) @@ -376,22 +175,22 @@ 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,84 +198,20 @@ 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 Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/doc/Tdas0211.tex =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/doc/Tdas0211.tex (revision 1720) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/doc/Tdas0211.tex (revision 1722) @@ -33,4 +33,5 @@ \usepackage{magic-tdas} +\usepackage{amssymb} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -46,6 +47,6 @@ \author{A.Moralejo\\ \texttt{}} -\date{December 6, 2002\\} -\TDAScode{MAGIC-TDAS 02-11\\ 021206/AMoralejo} +\date{January 20, 2003\\} +\TDAScode{MAGIC-TDAS 02-11\\ 030120/AMoralejo} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -61,6 +62,17 @@ 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,5 +173,5 @@ 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,6 +194,6 @@ 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,6 +254,8 @@ 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,5 +302,5 @@ 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,5 +355,5 @@ 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,6 +397,6 @@ 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,7 +409,7 @@ 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,5 +427,43 @@ 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,41 +490,34 @@ 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,16 +560,4 @@ \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,4 +572,52 @@ 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,7 +638,11 @@ 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,7 +719,8 @@ \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,19 +762,22 @@ 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,5 +844,5 @@ 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,53 +885,50 @@ %------------------------------------------------------------ -\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,6 +937,26 @@ \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,51 +964,38 @@ \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,58 +1005,53 @@ % -\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,27 +1060,46 @@ \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,9 +1113,12 @@ & 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,17 +1127,343 @@ \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,20 +1496,28 @@ 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,4 +1540,22 @@ %%>>>> 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} Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/init.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/init.c (revision 1720) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/init.c (revision 1722) @@ -113,5 +113,8 @@ 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); Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/init.h =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/init.h (revision 1720) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/init.h (revision 1722) @@ -102,5 +102,5 @@ " - 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 */ \ Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.c (revision 1720) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/parms.c (revision 1722) @@ -12,4 +12,5 @@ char axisdev_filename[256], reflectivity_filename[256]; char geo_filename[256]; +char atmosphere[256]; /* atmospheric model */ /* Prototypes */ @@ -67,5 +68,4 @@ 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,5 +94,5 @@ break; case atm_model: - SetAtmModel(value_ptr); + strcpy(atmosphere, value_ptr); break; case verbose_level: Index: trunk/MagicSoft/Simulation/Detector/ReflectorII/reflector.c =================================================================== --- trunk/MagicSoft/Simulation/Detector/ReflectorII/reflector.c (revision 1720) +++ trunk/MagicSoft/Simulation/Detector/ReflectorII/reflector.c (revision 1722) @@ -47,4 +47,5 @@ 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,5 +471,7 @@ 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,9 +493,25 @@ 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]); + }