Tools/ShowerUniversality/UnivParamNS/Atmosphere.cc
Go to the documentation of this file.
1 
19 #include "Atmosphere.h"
20 
21 #include <stdio.h>
22 #include <stdlib.h>
23 #include <math.h>
24 
25 using namespace AtmosphereNS;
26 
27 Atmosphere::Atmosphere(const bool UseGDAS, const std::string atmtype) : atmType(atmtype)
28 {
29  ok = false;
30 
31  if (atmType == "monthly") { //--- Montly (used mostly in MC simulations)
32 
33  //
34  // For explanation of the matrices a[i], b[i], c[i] and h[i]
35  // see CORSIKA_PHYSICS in CORSIKA documentation
36  // Eq. (2.1) and (2.2)
37  // T(h) = a + b exp(h / c) for i = 1...4
38  // T(h) = a - b * h/c for i = 5
39  // with [T]=g/cm**2
40  //
41 
42  float ai[12][5] = { { -129.86987635, -13.912578027, 1.137905219, -4.5501979957e-4, 0.01128292},
43  { -140.02703037, -32.154019677, 1.3324218584, -3.4759783013e-4, 0.01128292},
44  { -140.79368856, -36.408818366, 1.4042651261, -2.4583883619e-4, 0.01128292},
45  { -133.89496569, -47.089899137, 0.83278306486, -2.1054367331e-4, 0.01128292},
46  { -152.40981237, -12.095554799, 0.85136995339, -4.5928490669e-5, 0.01128292},
47  { -187.01662033, -25.280960312, 0.89432403304, 3.47276315e-5 , 0.01128292},
48  { -136.34855518, -20.914318979, 0.87280369409, 3.732759787e-5 , 0.01128292},
49  { -139.77721439, -25.948986634, 0.70618984827, -9.6699826826e-6, 0.01128292},
50  { -158.54644384, -95.110327312, 0.76408408672, -1.5130518952e-4, 0.01128292},
51  { -150.1396913, -58.140415154, 0.87315735073, -3.6221952628e-4, 0.01128292},
52  { -194.83705256, -68.148895317, 0.91697530761, -5.9280159541e-4, 0.01128292},
53  { -152.92783842, -37.51920314, 1.5800882678, -5.9051575509e-4, 0.01128292}
54  };
55 
56  float bi[12][5] = { {1170.0778448, 1310.6961298, 1490.6966654, 503.61356785, 1.},
57  {1177.9948918, 1226.5793800, 1544.3942757, 529.86754375, 1.},
58  {1177.4425334, 1189.0214700, 1556.0070888, 539.33390781, 1.},
59  {1174.1385323, 1215.7355043, 1425.6918596, 503.93169671, 1.},
60  {1190.3622315, 1265.7764387, 1437.9375049, 524.79662686, 1.},
61  {1227.6177212, 1178.8761669, 1445.4602411, 535.28753109, 1.},
62  {1182.9054740, 1213.2533581, 1413.9801198, 547.58686206, 1.},
63  {1184.5463072, 1205.4324091, 1366.0574680, 542.40761127, 1.},
64  {1200.3725926, 1165.8522743, 1373.2445748, 521.58648708, 1.},
65  {1189.9866918, 1193.1675932, 1422.7442647, 490.3837468, 1.},
66  {1231.0108397, 1164.4301252, 1412.7462127, 460.92737873, 1.},
67  {1190.2311741, 1189.4576063, 1575.7147090, 484.1030213, 1.}
68  };
69 
70  float ci[12][5] = { {971950.03798, 682326.92144, 615751.05849, 795110.76, 1.e9},
71  {990841.49847, 746882.89502, 608546.43122, 787929.8796, 1.e9},
72  {984028.35984, 766300.21844, 604631.06014, 782862.57672, 1.e9},
73  {967164.95864, 768155.10252, 621852.25389, 785601.62044, 1.e9},
74  {985619.00774, 681017.97599, 617989.22267, 776008.70875, 1.e9},
75  {1012299.9934, 732521.03651, 614783.70338, 771077.62332, 1.e9},
76  {947431.21788, 710893.29507, 619412.05748, 769605.79586, 1.e9},
77  {953640.44565, 723050.32076, 629158.20120, 772372.2386, 1.e9},
78  {986103.78389, 881190.80912, 628689.50007, 781005.59304, 1.e9},
79  {990579.96621, 800597.35355, 623333.46715, 793328.53093, 1.e9},
80  {1042817.5709, 832823.86466, 624149.14260, 805671.77828, 1.e9},
81  {997989.25368, 766606.85869, 599559.45014, 802421.45312, 1.e9}
82  };
83 
84  float hi[12][5] = { {0., 1070000., 1460000., 3660000., 10000000. },
85  {0., 980000., 1430000., 3550000., 10000000. },
86  {0., 900000., 1450000., 3500000., 10000000. },
87  {0., 990000., 1210000., 3800000., 10000000. },
88  {0., 950000., 1500000., 3750000., 10000000. },
89  {0., 750000., 1450000., 3700000., 10000000. },
90  {0., 890000., 1400000., 3700000., 10000000. },
91  {0., 880000., 1280000., 3830000., 10000000. },
92  {0., 690000., 1100000., 3820000., 10000000. },
93  {0., 900000., 1200000., 3800000., 10000000. },
94  {0., 700000., 1200000., 3810000., 10000000. },
95  {0., 870000., 1450000., 3480000., 10000000. }
96  };
97 
98  //--- GDAS
99  float ai_gdas[12][5] = {
100  { -136.72575606000001 , -31.63664304400000 , 1.88902340350000 , 0.00039201867984 , 0.01128000000000 },
101  { -137.25655861999999 , -31.79397889600000 , 2.06162275470000 , 0.00041243062289 , 0.01128000000000 },
102  { -132.36885161999999 , -29.07704662900000 , 2.09050150900000 , 0.00043534337925 , 0.01128000000000 },
103  { -129.99304119999999 , -21.84724843800000 , 1.52111364840000 , 0.00039559055121 , 0.01128000000000 },
104  { -125.11468467000000 , -14.59123562100000 , 0.93641128677000 , 0.00032475590985 , 0.01128000000000 },
105  { -126.17178851000000 , -7.72898528110000 , 0.81676828638000 , 0.00031947676891 , 0.01128000000000 },
106  { -126.17216789000000 , -8.61825375140000 , 0.74177836911000 , 0.00029350702097 , 0.01128000000000 },
107  { -123.27936204000000 , -10.05149304100000 , 0.84187346153000 , 0.00032422546759 , 0.01128000000000 },
108  { -126.94494665000001 , -9.55565369810000 , 0.74939405052000 , 0.00029823116961 , 0.01128000000000 },
109  { -133.13151124999999 , -13.97320926500000 , 0.83782634310000 , 0.00031117421760 , 0.01128000000000 },
110  { -134.72208165000001 , -18.17238290800000 , 1.11598068450000 , 0.00035217025515 , 0.01128000000000 },
111  { -135.40825208999999 , -22.83040902600000 , 1.42234534930000 , 0.00037512921774 , 0.01128000000000 }
112  };
113 
114  float bi_gdas[12][5] = {
115  { 1174.82983340000010 , 1204.82334530000003 , 1637.77035830000000 , 735.96095022999998 , 1.00000000000000 },
116  { 1176.09075650000000 , 1197.89511040000002 , 1646.46169549999991 , 755.18728656999997 , 1.00000000000000 },
117  { 1172.62277840000002 , 1215.39646769999990 , 1617.00992820000010 , 769.51991638000004 , 1.00000000000000 },
118  { 1172.32918780000000 , 1250.29227740000010 , 1542.62484130000007 , 713.10082850000003 , 1.00000000000000 },
119  { 1169.95113020000008 , 1277.67684880000002 , 1493.53037810000001 , 617.96607470000004 , 1.00000000000000 },
120  { 1171.09162760000004 , 1295.35164340000006 , 1455.30093439999996 , 595.11713507000002 , 1.00000000000000 },
121  { 1172.73406879999993 , 1258.91800790000002 , 1450.05371409999998 , 583.07727714999999 , 1.00000000000000 },
122  { 1169.76303600000006 , 1251.02198079999994 , 1436.64993720000007 , 627.42169844000000 , 1.00000000000000 },
123  { 1174.86764530000005 , 1251.55885289999992 , 1440.82575489999999 , 606.31473165000000 , 1.00000000000000 },
124  { 1176.98334730000010 , 1244.23453100000006 , 1464.01208550000001 , 622.11207419000004 , 1.00000000000000 },
125  { 1175.77379719999999 , 1238.95385039999996 , 1505.16143659999989 , 670.64752105000002 , 1.00000000000000 },
126  { 1174.64497099999994 , 1227.27536830000008 , 1585.71305619999998 , 691.23389637000002 , 1.00000000000000 }
127  };
128 
129  float ci_gdas[12][5] = {
130  { 982815.95247999997810 , 754029.87759000004735 , 594416.83822000003420 , 733974.36971999995876 , 1000000000.00000000000000 },
131  { 981369.61250000004657 , 756657.65382999996655 , 592969.89671000000089 , 731345.88332000002265 , 1000000000.00000000000000 },
132  { 972654.05630000005476 , 742769.21710000000894 , 595342.19851000001654 , 728921.61953999998514 , 1000000000.00000000000000 },
133  { 962396.55209999997169 , 711452.06672999996226 , 603480.61835000000428 , 735460.83741000003647 , 1000000000.00000000000000 },
134  { 947742.88769000000320 , 685089.57509000005666 , 609640.01931999996305 , 747555.95525999995880 , 1000000000.00000000000000 },
135  { 940102.98841999995057 , 661697.57542999996804 , 612702.06319999997504 , 749976.26832000003196 , 1000000000.00000000000000 },
136  { 934649.58886000001803 , 672975.82513000001200 , 614888.52457999996841 , 752631.28535999997985 , 1000000000.00000000000000 },
137  { 931569.97624999994878 , 678861.75135999999475 , 617363.34490999998525 , 746739.16140999994241 , 1000000000.00000000000000 },
138  { 936953.91919000004418 , 678906.60516000003554 , 618132.60560999996960 , 750154.67709000001196 , 1000000000.00000000000000 },
139  { 954151.40399999998044 , 692708.89815999998245 , 615439.43935999996029 , 747969.08132999995723 , 1000000000.00000000000000 },
140  { 964877.07765999995172 , 706199.57501999998931 , 610242.24563999997918 , 741412.74548000004143 , 1000000000.00000000000000 },
141  { 973884.44360999995843 , 723759.74682000000030 , 600308.13983000000007 , 738390.20524999999907 , 1000000000.00000000000000 }
142  };
143 
144  float hi_gdas[12][5] = {
145  { 0.00000000000000 , 940000.00000000000000 , 1530000.00000000000000 , 3160000.00000000000000 , 10000000.00000000000000 },
146  { 0.00000000000000 , 919999.99999999988358 , 1540000.00000000000000 , 3100000.00000000000000 , 10000000.00000000000000 },
147  { 0.00000000000000 , 960000.00000000000000 , 1520000.00000000000000 , 3070000.00000000000000 , 10000000.00000000000000 },
148  { 0.00000000000000 , 1000000.00000000000000 , 1490000.00000000000000 , 3260000.00000000000000 , 10000000.00000000000000 },
149  { 0.00000000000000 , 1019999.99999999988358 , 1510000.00000000000000 , 3590000.00000000000000 , 10000000.00000000000000 },
150  { 0.00000000000000 , 1010000.00000000000000 , 1600000.00000000000000 , 3670000.00000000046566 , 10000000.00000000000000 },
151  { 0.00000000000000 , 960000.00000000000000 , 1650000.00000000000000 , 3740000.00000000000000 , 10000000.00000000000000 },
152  { 0.00000000000000 , 960000.00000000000000 , 1590000.00000000000000 , 3629999.99999999953434 , 10000000.00000000000000 },
153  { 0.00000000000000 , 950000.00000000000000 , 1620000.00000000000000 , 3720000.00000000046566 , 10000000.00000000000000 },
154  { 0.00000000000000 , 950000.00000000000000 , 1550000.00000000000000 , 3650000.00000000000000 , 10000000.00000000000000 },
155  { 0.00000000000000 , 960000.00000000000000 , 1530000.00000000000000 , 3460000.00000000000000 , 10000000.00000000000000 },
156  { 0.00000000000000 , 960000.00000000000000 , 1560000.00000000000000 , 3329999.99999999953434 , 10000000.00000000000000 }
157  };
158 
159 
160  for (int imonth = 0; imonth < 12; imonth++) {
161  for (int i = 0; i < 5; i++) {
162  if (UseGDAS) {
163  amonth[imonth].ai[i] = ai_gdas[imonth][i];
164  amonth[imonth].bi[i] = bi_gdas[imonth][i];
165  amonth[imonth].ci[i] = ci_gdas[imonth][i];
166  amonth[imonth].hi[i] = hi_gdas[imonth][i];
167  } else {
168  amonth[imonth].ai[i] = ai[imonth][i];
169  amonth[imonth].bi[i] = bi[imonth][i];
170  amonth[imonth].ci[i] = ci[imonth][i];
171  amonth[imonth].hi[i] = hi[imonth][i];
172  }
173  }
174  }
175 
176  } else if (atmType == "seasonal") {
177 
178  float ai[5][5] = { {0, 0, 0, 0, 0},
179  {0, 0, 0, 0, 0},
180  {0, 0, 0, 0, 0},
181  {0, 0, 0, 0, 0},
182  {0, 0, 0, 0, 0}
183  };
184 
185  float bi[5][5] = { {0, 0, 0, 0, 0},
186  {0, 0, 0, 0, 0},
187  {0, 0, 0, 0, 0},
188  {0, 0, 0, 0, 0},
189  {0, 0, 0, 0, 0}
190  };
191 
192  float ci[5][5] = { {0, 0, 0, 0, 0},
193  {0, 0, 0, 0, 0},
194  {0, 0, 0, 0, 0},
195  {0, 0, 0, 0, 0},
196  {0, 0, 0, 0, 0}
197  };
198 
199  float hi[5][5] = { {0, 0, 0, 0, 0},
200  {0, 0, 0, 0, 0},
201  {0, 0, 0, 0, 0},
202  {0, 0, 0, 0, 0},
203  {0, 0, 0, 0, 0}
204  };
205 
206  float ai_gdas[5][5] = { {0, 0, 0, 0, 0},
207  {0, 0, 0, 0, 0},
208  {0, 0, 0, 0, 0},
209  {0, 0, 0, 0, 0},
210  {0, 0, 0, 0, 0}
211  };
212 
213  float bi_gdas[5][5] = { {0, 0, 0, 0, 0},
214  {0, 0, 0, 0, 0},
215  {0, 0, 0, 0, 0},
216  {0, 0, 0, 0, 0},
217  {0, 0, 0, 0, 0}
218  };
219 
220  float ci_gdas[5][5] = { {0, 0, 0, 0, 0},
221  {0, 0, 0, 0, 0},
222  {0, 0, 0, 0, 0},
223  {0, 0, 0, 0, 0},
224  {0, 0, 0, 0, 0}
225  };
226 
227  float hi_gdas[5][5] = { {0, 0, 0, 0, 0},
228  {0, 0, 0, 0, 0},
229  {0, 0, 0, 0, 0},
230  {0, 0, 0, 0, 0},
231  {0, 0, 0, 0, 0}
232  };
233 
234  for (int iseason = 0; iseason < 5; iseason++) {
235  for (int i = 0; i < 5; i++) {
236  if (UseGDAS) {
237  aseason[iseason].ai[i] = ai_gdas[iseason][i];
238  aseason[iseason].bi[i] = bi_gdas[iseason][i];
239  aseason[iseason].ci[i] = ci_gdas[iseason][i];
240  aseason[iseason].hi[i] = hi_gdas[iseason][i];
241  } else {
242  aseason[iseason].ai[i] = ai[iseason][i];
243  aseason[iseason].bi[i] = bi[iseason][i];
244  aseason[iseason].ci[i] = ci[iseason][i];
245  aseason[iseason].hi[i] = hi[iseason][i];
246  }
247  }
248  }
249  }
250 
251 }
252 
253 Atmosphere::~Atmosphere()
254 {
255 
256 }
257 
258 void Atmosphere::SetCurrentAtmosphere(AtmModel theAtm)
259 {
260 
261 
262  for (int i = 0; i < 5; i++) {
263  acurrent.ai[i] = theAtm.ai[i];
264  acurrent.bi[i] = theAtm.bi[i];
265  acurrent.ci[i] = theAtm.ci[i];
266  acurrent.hi[i] = theAtm.hi[i];
267  }
268 
269  acurrent.atmflag = 0;
270 
271  //Check if it is a monthly atmosphere
272  for (int imonth = 0; imonth < 12; imonth++) {
273  bool flag = true;
274  for (int j = 0; j < 5; j++) {
275  if (fabs(amonth[imonth].ai[j] / acurrent.ai[j] - 1) > 0.01) flag = false;
276  if (fabs(amonth[imonth].bi[j] / acurrent.bi[j] - 1) > 0.01) flag = false;
277  if (fabs(amonth[imonth].ci[j] / acurrent.ci[j] - 1) > 0.01) flag = false;
278  }
279 
280  if (flag) {
281  acurrent.atmflag = 1 + imonth;
282  break;
283  }
284  }
285  if (acurrent.atmflag == 0) {
286  ok = false;
287  printf("ERROR: I did not recognize the atmosphere %d \n", acurrent.atmflag);
288  exit(1);
289  } else {
290  ok = true;
291  //printf("Atmosphere Model: %d month \n",acurrent.atmflag);
292  }
293 
294 }
295 
296 
297 
298 void Atmosphere::SetCurrentAtmosphere(float* hi, float* ai, float* bi, float* ci)
299 {
300 
301  for (int i = 0; i < 5; i++) {
302  acurrent.ai[i] = ai[i];
303  acurrent.bi[i] = bi[i];
304  acurrent.ci[i] = ci[i];
305  acurrent.hi[i] = hi[i];
306  }
307 
308  acurrent.atmflag = 0;
309 
310  //Check if it is a monthly atmosphere
311  for (int imonth = 0; imonth < 12; imonth++) {
312  bool flag = true;
313  for (int j = 0; j < 5; j++) {
314  if (fabs(amonth[imonth].ai[j] / acurrent.ai[j] - 1) > 0.01) flag = false;
315  if (fabs(amonth[imonth].bi[j] / acurrent.bi[j] - 1) > 0.01) flag = false;
316  if (fabs(amonth[imonth].ci[j] / acurrent.ci[j] - 1) > 0.01) flag = false;
317  }
318 
319  if (flag) {
320  acurrent.atmflag = 1 + imonth;
321  break;
322  }
323  }
324  if (acurrent.atmflag == 0) {
325  ok = false;
326  printf("ERROR: I did not recognize the atmosphere %d \n", acurrent.atmflag);
327  exit(1);
328  } else {
329  ok = true;
330  printf("Atmosphere Model: %d month \n", acurrent.atmflag);
331  }
332 }
333 
334 void Atmosphere::SetCurrentAtmosphere(int imonth) // imonth=1,...,12
335 {
336  if (atmType == "monthly") {
337  if (imonth < 1 || imonth > 12) {
338  printf("Wrong atmosphere \n");
339  return;
340  }
341  for (int i = 0; i < 5; i++) {
342  acurrent.ai[i] = amonth[imonth - 1].ai[i];
343  acurrent.bi[i] = amonth[imonth - 1].bi[i];
344  acurrent.ci[i] = amonth[imonth - 1].ci[i];
345  acurrent.hi[i] = amonth[imonth - 1].hi[i];
346  }
347  } else if (atmType == "seasonal") {
348  if (imonth < 1 || imonth > 5) {
349  printf("Wrong atmosphere \n");
350  return;
351  }
352  for (int i = 0; i < 5; i++) {
353  acurrent.ai[i] = aseason[imonth - 1].ai[i];
354  acurrent.bi[i] = aseason[imonth - 1].bi[i];
355  acurrent.ci[i] = aseason[imonth - 1].ci[i];
356  acurrent.hi[i] = aseason[imonth - 1].hi[i];
357  }
358  }
359  acurrent.atmflag = imonth;
360  ok = true;
361 }
362 
363 void Atmosphere::ClearCurrentAtmosphere()
364 {
365  ok = false;
366  for (int i = 0; i < 5; i++) {
367  acurrent.ai[i] = 0.;
368  acurrent.bi[i] = 0.;
369  acurrent.ci[i] = 0.;
370  acurrent.hi[i] = 0.;
371  }
372 }
373 
374 
375 void Atmosphere::printCurrentAtmosphere(FILE* fp)
376 {
377  fprintf(fp, "---------------------\n");
378  fprintf(fp, "Current atmosphere %d (%s)\n", acurrent.atmflag, atmType.c_str());
379  fprintf(fp, "---------------------\n");
380 
381  fprintf(fp, " ai bi ci hi \n");
382  for (int i = 0; i < 5; i++)
383  fprintf(fp, "Layer=%d %f %f %f %f \n", i, acurrent.ai[i], acurrent.bi[i], acurrent.ci[i], acurrent.hi[i]);
384 
385 }
386 
387 float Atmosphere::VertDepthToHeight(float vdepth)
388 {
389  if (!ok) return 0.;
390  float h = 0.;
391 
392  float vdepth_min = acurrent.ai[4] - acurrent.bi[4] * 1.e7 / acurrent.ci[4];
393  float vdepth_max = acurrent.ai[0] + acurrent.bi[0];
394 
395  // float h0 = acurrent.ai[4] * acurrent.ci[4] / acurrent.bi[4];
396  if (vdepth < vdepth_min || vdepth > vdepth_max) {
397  //if (vdepth <= 0. ) h=h0;
398  if (vdepth < vdepth_min) h = (acurrent.ai[4] - vdepth) * acurrent.ci[4] / acurrent.bi[4];
399  else if (vdepth > vdepth_max) {
400  double density = acurrent.bi[0] / acurrent.ci[0];
401  h = (vdepth_max - vdepth) / density; // Constant density below h=0.
402  }
403  return h;
404  }
405 
406  //From Height in cm to Vertical Depth
407  for (int k = 0; k < 4; k++) {
408  h = acurrent.ci[k] * log(acurrent.bi[k] / (vdepth - acurrent.ai[k]));
409  if (h > acurrent.hi[k] && h < acurrent.hi[k + 1]) break;
410  }
411  return h;
412 }
413 
414 
415 float Atmosphere::SlantDepthToHeight(float sdepth, float theta)
416 {
417  if (!ok) return 0.;
418  float vdepth = sdepth * cos(theta);
419  return VertDepthToHeight(vdepth);
420 }
421 
422 
423 float Atmosphere::HeightToVertDepth(float h)
424 {
425  if (!ok) return 0.;
426  float vdepth = 0.;
427 
428  float h_min = 0., h_max = acurrent.hi[4];
429  // float h0 = acurrent.ai[4] * acurrent.ci[4] / acurrent.bi[4];
430  if (h > h_max || h < h_min) {
431  if (h < h_min) {
432  double density = acurrent.bi[0] / acurrent.ci[0];
433  vdepth = (acurrent.ai[0] + acurrent.bi[0]) + (density * (h_min - h)); //Constant density h<0.
434  } else if (h > h_max /*&& h < h0 */) vdepth = acurrent.ai[4] - acurrent.bi[4] * h / acurrent.ci[4];
435  //else if (h >h0) vdepth=0.;
436  return vdepth;
437  }
438 
439  for (int k = 0; k < 4; k++) {
440  if (h >= acurrent.hi[k] && h < acurrent.hi[k + 1]) {
441  vdepth = acurrent.ai[k] + acurrent.bi[k] * exp(-h / acurrent.ci[k]);
442  return vdepth;
443  }
444  }
445  return vdepth;
446 
447 }
448 
449 float Atmosphere::HeightToSlantDepth(float h, float theta)
450 {
451  if (!ok) return 0.;
452  return (HeightToVertDepth(h) / cos(theta));
453 }
454 
455 
456 
457 float Atmosphere::GetDensityAboveGround(float hground, float theta, float nRad)
458 {
459  //------
460  if (!ok) {
461  float mean = 0.;
462  for (int i = 1; i < 13; i++) {
463  SetCurrentAtmosphere(i);
464  mean += GetDensityAboveGround_i(hground, theta, nRad);
465  }
466  ok = false;
467  return mean / 12.;
468  }
469  //------
470  else return GetDensityAboveGround_i(hground, theta, nRad);
471 
472  return 0.;
473 }
474 
475 float Atmosphere::GetDensityAboveGround_i(float hground, float theta, float nRad)
476 {
477  if (!ok) return 0.;
478 
479  const float X0 = 37.;
480  float SlantDepth = HeightToSlantDepth(hground, theta) - nRad * X0;
481  float h = SlantDepthToHeight(SlantDepth, theta);
482  return GetDensity(h);
483 }
484 
485 float Atmosphere::GetDensity(float h) // g/cm^3
486 {
487  if (!ok) {
488  float mean = 0.;
489  for (int i = 1; i < 13; i++) {
490  SetCurrentAtmosphere(i);
491  mean += GetDensity_i(h);
492  }
493  ok = false;
494  return mean / 12.;
495  } else return GetDensity_i(h);
496 }
497 
498 float Atmosphere::GetDensity_i(float h) // g/cm^3
499 {
500  if (!ok) return 0.;
501  float density = 0.;
502 
503  float h_min = 0., h_max = acurrent.hi[4];
504  float h0 = acurrent.ai[4] * acurrent.ci[4] / acurrent.bi[4];
505  if (h > h_max || h < h_min) {
506  if (h < h_min) density = acurrent.bi[0] / acurrent.ci[0];
507  else if (h > h_max && h < h0) density = acurrent.ai[4] - acurrent.bi[4] / acurrent.ci[4];
508  else if (h > h0) density = 0.;
509  return density;
510  }
511 
512  for (int k = 0; k < 4; k++) {
513  if (h >= acurrent.hi[k] && h < acurrent.hi[k + 1]) {
514  density = acurrent.bi[k] * exp(-h / acurrent.ci[k]) / acurrent.ci[k];
515  return density;
516  }
517  }
518  return density;
519 }
520 
521 //-------------------------------------------------------------------
522 //-------------------------------------------------------------------
523 //-------------------------------------------------------------------
524 
525 
526 float Atmosphere::Get_DX_DL(float r, float psi, float SlantDepth, float theta, float hground, bool UseDL, bool IsDiffusive)
527 {
528  if (SlantDepth < 0.) return 0.;
529 
530  // ----------------
531  if (!ok) {
532  float mean = 0.;
533  for (int i = 1; i < 13; i++) {
534  SetCurrentAtmosphere(i);
535  mean += Get_DX_DL_i(r, psi, SlantDepth, theta, hground, UseDL, IsDiffusive);
536  }
537  ok = false;
538  return mean / 12.;
539  }
540  // ----------------
541  else return Get_DX_DL_i(r, psi, SlantDepth, theta, hground, UseDL, IsDiffusive);
542 
543  return 0.;
544 }
545 
546 
547 
548 float Atmosphere::Get_DX_DL_i(float r, float psi, float SlantDepth, float theta, float hground, bool UseDL, bool IsDiffusive)
549 {
550  if (!ok) return 0.;
551 
552  //---Diffusive (along the shower axis )
553  if (IsDiffusive) {
554  float dh = r * cos(psi) * sin(theta);
555  float h0 = hground + dh;
556 
557  if (UseDL) return (SlantDepthToHeight(SlantDepth, theta) - h0) / cos(theta) / 1000.; // [10 m]
558  else return HeightToSlantDepth(h0, theta) - SlantDepth; // [g/cm2]
559  }
560  //---Rectilinear ( from position to a point in shower axis )
561  else {
562  // X axis is shower axis, origin at h=hground
563 
564  // Maximum
565  double hmax = SlantDepthToHeight(SlantDepth, theta);
566  double rmax[3] = {(hmax - hground) / cos(theta), 0., 0.};
567 
568  // Ground position of the detector
569  double rdet[3] = {r * tan(theta)* cos(psi), r * sin(psi), -r * cos(psi)};
570  // Vector Detector->maximum
571  double rdetmax[3] = { rmax[0] - rdet[0], rmax[1] - rdet[1], rmax[2] - rdet[2] };
572 
573  // (Distance,alpha)
574  double l = sqrt(pow(rdetmax[0], 2) + pow(rdetmax[1], 2) + pow(rdetmax[2], 2));
575  double calpha = rdetmax[0] / l;
576 
577  if (calpha < 0.) return -100.;
578 
579  // Return DL
580  if (UseDL) return l / 1000.; // [10 m] , roughly [g/cm2]
581 
582  // Return DX
583  double ctheta = calpha * cos(theta) + sin(theta) * cos(psi) * sqrt(1. - calpha * calpha);
584  double DX = HeightToSlantDepth(hground, acos(ctheta)) - HeightToSlantDepth(hmax, acos(ctheta));
585 
586  return DX;
587  }
588 }
589 
590 
591 //-------------------------------------------------------------------
592 //-------------------------------------------------------------------
593 //-------------------------------------------------------------------
594 
595 
596 float Atmosphere::GetTimePF_delay(float x, float y, float z, float theta, float azi)
597 // Time origin: (0.,0.,0.)
598 // This routine returns the arrival time of plane front at (x,y,z)
599 {
600  float costheta = cos(theta);
601  float sintheta = sin(theta);
602  float cosazi = cos(azi);
603  float sinazi = sin(azi);
604 
605  //----Ground coordinate system with x-axis the projection of shower axis into Obs Level
606  double xr = x * cosazi + y * sinazi;
607  // double yr = -x * sinazi + y * cosazi;
608  double zr = z;
609 
610  float dT = -(xr * sintheta + zr * costheta) / c_cmperns;
611 
612  return dT;
613 }
614 
615 //-------------------------------------------------------------------
616 //-------------------------------------------------------------------
617 //-------------------------------------------------------------------
618 
619 float Atmosphere::GetTimePF(float r, float psi, float SlantDepth, float theta, float hground)
620 {
621  if (SlantDepth < 0.) return 0.;
622 
623  // ----------------
624  if (!ok) {
625  float mean = 0.;
626  for (int i = 1; i < 13; i++) {
627  SetCurrentAtmosphere(i);
628  mean += GetTimePF_i(r, psi, SlantDepth, theta, hground);
629  }
630  ok = false;
631  return mean / 12.;
632  }
633  // ----------------
634  else return GetTimePF_i(r, psi, SlantDepth, theta, hground);
635 
636  return 0.;
637 }
638 
639 float Atmosphere::GetTimePF_i(float r, float psi, float SlantDepth, float theta, float hground)
640 // Time origin: point in shower axis at Slantdepth
641 // This routine returns the time a plane front takes to propagate from the origin of times to a detector at (r,psi,hground)
642 {
643  if (!ok) return 0.;
644  float h2 = SlantDepthToHeight(SlantDepth, theta);
645  return GetTimePF_h(r, psi, h2, theta, hground);
646 }
647 
648 
649 float Atmosphere::GetTimePF_h(float r, float psi, float h2, float theta, float hground)
650 // Time origin: point in shower axis at h2
651 // This routine returns the time a plane front takes to propagate from the origin of times to a detector at (r,psi,hground)
652 {
653  if (!ok) return 0.;
654 
655  float dh = r * cos(psi) * sin(theta);
656  float h1 = hground + dh;
657  float t0 = (h2 - h1) / cos(theta) / c_cmperns;
658  return t0;
659 }
660 
661 //-------------------------------------------------------------------
662 //-------------------------------------------------------------------
663 //-------------------------------------------------------------------
664 
665 float Atmosphere::GetTimeCF(float r, float psi, float SlantDepth, float theta, float hground) //Time origin: plane front arrival at detector position
666 {
667  if (SlantDepth < 0.) return 0.;
668 
669  // ----------------
670  if (!ok) {
671  float mean = 0.;
672  for (int i = 1; i < 13; i++) {
673  SetCurrentAtmosphere(i);
674  mean += GetTimeCF_i(r, psi, SlantDepth, theta, hground);
675  }
676  ok = false;
677  return mean / 12.;
678  }
679  // ----------------
680  else return GetTimeCF_i(r, psi, SlantDepth, theta, hground);
681 
682  return 0.;
683 }
684 
685 float Atmosphere::GetTimeCF_i(float r, float psi, float SlantDepth, float theta, float hground) //Time origin: plane front arrival at detector position
686 {
687  if (!ok) return 0.;
688  double h1 = SlantDepthToHeight(SlantDepth, theta);
689  return GetTimeCF_h(r, psi, h1, theta, hground);
690 }
691 
692 
693 float Atmosphere::GetTimeCF_h(float r, float psi, float h1, float theta, float hground) //Time origin: plane front arrival at detector position
694 // Time origin: point in shower axis at h1
695 // This routine returns the time difference between the arrival time of a plane front and a particle propagating in straight line from the origin of times
696 {
697  double d = (h1 - hground) / cos(theta) - r * tan(theta) * cos(psi) ;
698  double t0 = (sqrt(d * d + r * r) - d) / c_cmperns;
699  return t0;
700 }
AtmosphereNS::Atmosphere::GetTimePF_i
float GetTimePF_i(float r, float psi, float SlantDepth, float theta, float hground)
Definition: Tools/ShowerUniversality/UnivParamNS/Atmosphere.cc:639
AtmosphereNS::Atmosphere::GetTimeCF
float GetTimeCF(float r, float psi, float SlantDepth, float theta, float hground)
Definition: Tools/ShowerUniversality/UnivParamNS/Atmosphere.cc:665
AtmosphereNS::Atmosphere::Get_DX_DL
float Get_DX_DL(float r, float psi, float SlantDepth, float theta, float hground, bool UseDL, bool IsDiffusive)
Definition: Tools/ShowerUniversality/UnivParamNS/Atmosphere.cc:526
AtmosphereNS::Atmosphere::GetTimeCF_i
float GetTimeCF_i(float r, float psi, float SlantDepth, float theta, float hground)
Definition: Tools/ShowerUniversality/UnivParamNS/Atmosphere.cc:685
std::pow
double pow(const double x, const unsigned int i)
Definition: MuonArrivalTime.h:17
exit
int exit
Definition: dump1090.h:237
AtmosphereNS::Atmosphere::GetTimeCF_h
float GetTimeCF_h(float r, float psi, float h1, float theta, float hground)
Definition: Tools/ShowerUniversality/UnivParamNS/Atmosphere.cc:693
AtmosphereNS::Atmosphere::Atmosphere
Atmosphere(const bool UseGDAS=false, const std::string atmtype="monthly")
Definition: Tools/ShowerUniversality/UnivParamNS/Atmosphere.cc:27
AtmosphereNS::Atmosphere::Get_DX_DL_i
float Get_DX_DL_i(float r, float psi, float SlantDepth, float theta, float hground, bool UseDL, bool IsDiffusive)
Definition: Tools/ShowerUniversality/UnivParamNS/Atmosphere.cc:548
AtmosphereNS::Atmosphere::GetTimePF
float GetTimePF(float r, float psi, float SlantDepth, float theta, float hground)
Definition: Tools/ShowerUniversality/UnivParamNS/Atmosphere.cc:619
det::if
if(dataRoot)
Definition: XXMLManager.h:1003
AtmosphereNS::Atmosphere::GetDensityAboveGround
float GetDensityAboveGround(float hground, float theta, float nRad)
Definition: Tools/ShowerUniversality/UnivParamNS/Atmosphere.cc:457
UnivParamNS::UseDL
const bool UseDL[4]
Definition: UnivParam.h:36
AtmosphereNS::Atmosphere::GetTimePF_delay
float GetTimePF_delay(float x, float y, float z, float theta, float azi)
Definition: Tools/ShowerUniversality/UnivParamNS/Atmosphere.cc:596
AtmosphereNS::Atmosphere::GetTimePF_h
float GetTimePF_h(float r, float psi, float hStart, float theta, float hground)
Definition: Tools/ShowerUniversality/UnivParamNS/Atmosphere.cc:649
AtmosphereNS::Atmosphere::GetDensityAboveGround_i
float GetDensityAboveGround_i(float hground, float theta, float nRad)
Definition: Tools/ShowerUniversality/UnivParamNS/Atmosphere.cc:475
, generated on Tue Sep 26 2023.