Use of tabulated atmospheric profiles and atmospheric refraction. More...
#include <stdio.h>#include <stdlib.h>#include <string.h>#include <math.h>#include "atmo.h"#include "fileopen.h"#include "rpolator.h"#include "rpolator.c"#include "mc_atmprof.c"
Macros | |
| #define | NLAYX 15 |
| #define | M_PI (3.14159265358979323846264338327950288) |
| / | |
| #define | WITH_RPOLATOR_CSPLINE 1 /* Take advantage of cubic splines if we have them available */ |
| #define | FAST_INTERPOLATION always |
| #define | FAST_INTERPOLATION2 always |
| #define | FAST_INTERPOLATION3 always |
| #define | WITH_THICKX_DIRECT |
| #define | UNUSED(x) UNUSED_ ## x |
| #define | MAX_PROFILE 120 |
| #define | MAX_FAST_PROFILE 40000 |
| #define | _XSTR_(s) _STR_(s) |
| Expand a macro first and then enclose in string. | |
| #define | _STR_(s) #s |
| Enclose in string without macro expansion. More... | |
| #define | SHOW_MACRO(s) |
| The following shows the same output after a "#define XYZ" and a "#define XYZ 1" (or compiled with "-DXYZ") More... | |
| int | atmosphere |
| The atmospheric profile number, 0 for built-in. More... | |
| static char * | atmprof_name = NULL |
| File name for atmospheric profile table. More... | |
| static int | num_prof |
| Number of levels in original table. More... | |
| static double | p_alt [MAX_PROFILE] |
| Altitude levels in given table (converted to [cm], upward order). More... | |
| static double | p_alt_rev [MAX_PROFILE] |
| Altitude levels in given table (converted to [cm], reversed order). More... | |
| static double | p_rho [MAX_PROFILE] |
| Density at given levels [g/cm^3] (upward height order) | |
| static double | p_rho_rev [MAX_PROFILE] |
| Density at given levels [g/cm^3], reversed order. | |
| static double | p_thick [MAX_PROFILE] |
| Atmospheric thickness at given levels [g/cm^2]. | |
| static double | p_log_rho [MAX_PROFILE] |
| Log of density at given levels (for better interpolation) | |
| static double | p_log_thick [MAX_PROFILE] |
| Log of atmospheric thickness at given levels (upward height order) | |
| static double | p_log_thick_rev [MAX_PROFILE] |
| Log of atmospheric thickness at given levels (reversed order) | |
| static double | p_log_n1 [MAX_PROFILE] |
| Log of (index of refraction minus 1.0) at given levels. More... | |
| static double | p_bend_ray_hori_a [MAX_PROFILE] |
| Coefficient for horizontal displacement refraction effect (reversed order) | |
| static double | p_bend_ray_time_a [MAX_PROFILE] |
| Coefficient for arrival time effect by refraction, density dependence (reversed order) | |
| static double | p_bend_ray_time0 [MAX_PROFILE] |
| Coefficient for arrival time effect by refraction, altitude dependence. | |
| static CsplinePar * | cs_alt_log_rho |
| Cubic spline parameters for log(density) versus altitude interpolation. | |
| static CsplinePar * | cs_alt_log_thick |
| Cubic spline parameters for log(thickness) versus altitude interpolation. | |
| static CsplinePar * | cs_alt_log_n1 |
| Cubic spline parameters for log(n-1) versus altitude interpolation. | |
| static CsplinePar * | cs_log_thick_alt |
| Cubic spline parameters for altitude versus log(thickness) reverse interpolation. | |
| static CsplinePar * | cs_bend_rho_hori_a |
| Cubic spline parameters for horizontal displacement refraction effect. | |
| static CsplinePar * | cs_bend_rho_time_a |
| Cubic spline parameters for arrival time effect by refraction, density dependence. | |
| static CsplinePar * | cs_bend_alt_time0 |
| Cubic spline parameters for arrival time effect by refraction, altitude dependence. | |
| static double | top_of_atmosphere = 112.83e5 |
| Height of top of atmosphere [cm], = p_alt[num_prof-1]. | |
| static double | top_of_atm_table |
| The heighest entry in the atmospheric profile table. More... | |
| static double | bottom_of_atmosphere = 0. |
| Bottom is normally at sea level but table could differ, = p_alt[0]. | |
| static double | bottom_log_thickness |
| Thickness at bottom, depending on profile, = p_log_thick[0]. | |
| static double | top_log_thickness |
| Thickness at top is zero, thus using an extrapolation from the second last value. More... | |
| static double | top_layer_2nd_altitude |
| Second highest altitude tabulated, using exponential density profile above that. More... | |
| static double | top_layer_rho0 |
| Sea-level altitude matching profile in top layer, without scaling for thickness. More... | |
| static double | top_layer_cfac |
| Scaling factor needed to match thickness at second last level. More... | |
| static double | top_layer_hscale |
| Height scale of exponention density profile in top layer. More... | |
| static double | top_layer_hscale_rho0_cfac |
| top_layer_rho0 * top_layer_cfac * top_layer_hscale | |
| static double | top_layer_hscale_rho0_cfac_inv |
| 1. More... | |
| static double | top_layer_exp_top |
| exp(-top_of_atmosphere/top_layer_hscale) | |
| static double * | fast_p_alt |
| Equidistant steps in altitude. | |
| static double * | fast_p_log_rho |
| Log(rho) at fast_p_alt steps. | |
| static double * | fast_p_log_thick |
| Log(thick) at fast_p_alt steps. | |
| static double * | fast_p_log_n1 |
| Log(n-1) at fast_p_alt steps. | |
| static double | fast_h_fac |
| Inverse of height difference per step. | |
| static double * | fast_p_thick |
| Direct interpolation of thickness at fast_p_alt steps. | |
| static double * | fast_p_rho |
| Direct interpolation of density at fast_p_alt steps. | |
| static double * | fast_p_n1 |
| Direct interpolation of n-1 at fast_p_alt steps. | |
| static double * | fast_p_eq_log_thick |
| Equidistant steps in log(thick) | |
| static double * | fast_p_heigh_rev |
| Altitude at given log(thick) steps. | |
| static double | fast_log_thick_fac |
| Inverse of log(thickness) difference per step. | |
| static double * | fast_p_bend_ray_hori_a |
| Coefficient for horizontal displacement refraction effect. | |
| static double * | fast_p_bend_ray_time_a |
| Coefficient for arrival time effect by refraction, density dependence. | |
| static double * | fast_p_bend_ray_time0 |
| Coefficient for arrival time effect by refraction, altitude dependence. | |
| static double * | fast_p_rho_rev |
| Density steps used in fast interpolation, in order of incr. More... | |
| static double | fast_rho_fac |
| static double | etadsn = 0.000283 * 994186.38 / 1222.656 |
| static double | observation_level |
| Altitude [cm] of observation level. | |
| static double | obs_level_refidx |
| static double | obs_level_thick |
| static void | init_refraction_tables () |
| Initialize tables needed for atmospheric refraction. More... | |
| static void | init_fast_interpolation () |
| An alternate interpolation method (which requires that the table is sufficiently fine-grained and equidistant) has to be initialized first. | |
| static void | init_corsika_atmosphere () |
| Take the atmospheric profile from CORSIKA built-in functions. More... | |
| static void | init_atmosphere_from_text_file () |
| Initialize atmospheric profiles. More... | |
| static void | init_atmosphere_from_atmprof (void) |
| static double | sum_log_dev_sq (double a, double b, double c, int np, double *h, double *t, double *rho) |
| static double | atm_exp_fit (double h1, double h2, double *ap, double *bp, double *cp, double *s0, int *npp) |
| Fit one atmosphere layer by an expontential density model. | |
| void | free_atmo_csplines (void) |
| static void | trace_ray_planar (double emlev, double olev, double za, double step, double *xo, double *to, double *so) |
| Trace a downward light ray in a planar atmosphere. More... | |
| static void | init_common_atmosphere (void) |
| Second-stage part of atmospheric profile initialization is common to CORSIKA default profile and tabulated profiles. | |
| void | atmset_ (int *iatmo, double *obslev) |
| Set number of atmospheric model profile to be used. More... | |
| void | atmnam_ (const char *aname, double *obslev) |
| Instead of setting the atmospheric profile by number, it gets set by name, also indicated by setting profile number to 99. More... | |
| void | atm_init (AtmProf *aprof) |
| This variant is not usable from the FORTRAN code side, thus no underscore at the end of the function name. More... | |
| double | rhofx_ (double *height) |
| Density of the atmosphere as a function of altitude. More... | |
| double | thickx_ (double *height) |
| Atmospheric thickness [g/cm**2] as a function of altitude. More... | |
| double | refim1x_ (double *height) |
| Index of refraction minus 1 as a function of altitude [cm]. More... | |
| double | refidx_ (double *height) |
| Index of refraction as a function of altitude [cm]. More... | |
| double | heighx_ (double *thick) |
| Altitude [cm] as a function of atmospheric thickness [g/cm**2]. More... | |
| void | raybnd_ (double *zem, cors_dbl_t *u, cors_dbl_t *v, double *w, cors_dbl_t *dx, cors_dbl_t *dy, cors_dbl_t *dt) |
| Calculate the bending of light due to atmospheric refraction. More... | |
| static double | sum_log_dev_sq (double a, double b, double c, int np, double *h, double *t, double *UNUSED(rho)) |
| Measure of deviation of model layers from tables. | |
| static double | fn_thick (double h, int nl, double *hl, double *a, double *b, double *c) |
| Corresponding to CORSIKA built-in function THICK; only used to show fit results. More... | |
| static double | fn_rhof (double h, int nl, double *hl, double *UNUSED(a), double *b, double *c) |
| Corresponding to CORSIKA built-in function RHOF; only used to show fit results. More... | |
| void | atmfit_ (int *nlp, double *hlay, double *aatm, double *batm, double *catm) |
| Fit the tabulated density profile for CORSIKA EGS part. More... | |
| void | show_atmo_macros (void) |
Detailed Description
- Date
CVS $Date: 2020/08/25 15:56:39 $
- Version
CVS $Revision: 1.51 $
This file provides code for use of external atmospheric models (in the form of text-format tables) with the CORSIKA program. Six atmospheric models as implemented in the MODTRAN program and as tabulated in MODTRAN documentation (F.X. Kneizys et al. 1996, 'The MODTRAN 2/3 Report and LOWTRAN 7 Model', Phillips Laboratory, Hanscom AFB, MA 01731-3010, U.S.A.) are provided as separate files (atmprof1.dat ... atmprof6.dat). User-provided atmospheric models should be given model numbers above 6. For model number 99 you may actually define a non-standard file name, rather than atmprof99.dat.
Note that for the Cherenkov part and the hadronic (and muon) part of CORSIKA the table values are directly interpolated but the electron/positron/gamma part (derived from EGS) uses special layers (at present 4 with exponential density decrease and the most upper layer with constant density). Parameters of these layers are fitted to tabulated values but not every possible atmospheric model fits very well with an exponential profile. You are adviced to check that the fit matches tabulated values to sufficient precision in the altitude ranges of interest to you. Try to adjust layer boundary altitudes in case of problems. The propagation of light without refraction (as implemented in CORSIKA, unless using the CURVED option) and with refraction (as implemented by this software) assumes a plane-parallel atmosphere. Instead of a single set of starting layer boundaries, newer version (Jan. 2019) try different sets of boundaries and use the most promising set for further fit improvement.
There are different ways to load and interpolate an atmospheric profile table. The usual procedure is that after loading the values from the file, the nearly exponential structure is taken advantage of and a relatively large number (order 10k) of equidistant support points are pre-interpolated from the original non-equidistant levels of the input point, saving binary-search of the matching interval for each interpolation. This FAST_INTERPOLATION option is activated by default but can be disabled by compiling with '-DNO_FAST_INTERPOLATION' (or '-DNO_FAST_INTERPOLATION2' if only fast interpolation for the reverse direction, from thickness to altitude, should be disabled.) The initial pre-interpolation used to be a linear interpolation in altitude versus log(density) or log(thickness) or log(n-1), respectively but with the 'rpolator' table interpolation code this can be turned into cubic spline interpolation. Unless problems with the CORSIKA build procedure prevent that the rpolator code and its cubic spline interpolation are intended to be the default scheme from now on (Jan. 2019). In either case for the default setting, the rpolator code and be enforced with '-DWITH_RPOLATOR' or disabled with '-DNO_RPOLATOR'. Cubic spline interpolation can be disabled with '-DNO_RPOLATOR_CSPLINE'. The fine-grained interpolation under the FAST_INTERPOLATION scheme continues to be linear in altitude versus log(density) etc.
Because the CORSIKA build procedure does not know yet about the additional source file, it gets included when compiling atmo.c, unless compiled with '-DWITH_RPOLATOR_SEPARATE'.
