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