Changeset 1722 for trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c

Timestamp:

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

Author:

moralejo

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c (modified) (1 diff)

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  */
+
+
+