source: trunk/MagicSoft/Simulation/Detector/ReflectorII/atm.c@ 18477

/********************************************************************
*                                                                   *
*  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 "init.h"

/*  random numbers  */
#define RandomNumber  ranf()
#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(int model, float ol);
int absorption(float wlen, float height, float theta);
extern void attenu_(float *, float *, float *, float *, float *);  /* in Fortran        */
extern float ranf(void);

/* 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)
{
  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) */
  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;

}   /*  end of atm  */

int absorption(float wlen, float height, float theta)
{
  int ret = 0;  /*  0: passed, 1: absorbed  */
   
  if (RandomNumber > atm(wlen, height, theta)) ret=1;

  return ret; 
}   /*  end of absorption  */