Jpp 20.0.0-195-g190c9e876
the software that should make you happy
Loading...
Searching...
No Matches
JPhysicsSupportkit.hh
Go to the documentation of this file.
1#ifndef __JPHYSICS_JPHYSICSSUPPORTKIT__
2#define __JPHYSICS_JPHYSICSSUPPORTKIT__
3
4#include <cmath>
5
6#include "JTools/JFunction1D_t.hh"
7
8#include "JPhysics/JConstants.hh"
9
10
11/**
12 * \author mdejong
13 */
14
15namespace JPHYSICS {}
16namespace JPP { using namespace JPHYSICS; }
17
18/**
19 * Auxiliary methods for light properties of deep-sea water.
20 */
21namespace JPHYSICS {
22
23 using JTOOLS::JPolint1Function1D_t;
24 using JTOOLS::JGridSplineFunction1D_t;
25
26
27 /**
28 * Get multiple Coulomb scattering angle.
29 *
30 * The formula is taken from reference:
31 * Particle Data Book, formula 27.14.
32 *
33 * \param E Energy [GeV]
34 * \param x distance [m]
35 * \param X0 radiation length [m]
36 * \param M mass [GeV]
37 * \param Q charge [unit]
38 * \return angle [rad]
39 */
40 inline double getThetaMCS(const double E,
41 const double x,
42 const double X0,
43 const double M,
44 const double Q)
45 {
46 if (E > M) {
47
48 const double p = sqrt((E + M) * (E - M));
49 const double beta = p / E;
50
51 return THETA_MCS * Q * sqrt(x/X0) * (1.0 + 0.038*log(x/X0)) / (beta*p);
52 }
53
54 return 0.0;
55 }
56
57
58 /**
59 * Get multiple Coulomb scattering angle for muon.
60 *
61 * \param E Energy [GeV]
62 * \param x distance [m]
63 * \return angle [rad]
64 */
65 inline double getThetaMCS(const double E, const double x)
66 {
67 return getThetaMCS(E, x, X0_WATER_M, MASS_MUON, 1.0);
68 }
69
70
71 /**
72 * Auxiliary method to describe light scattering in water (Henyey-Greenstein).
73 *
74 * \param g angular dependence parameter
75 * \param x cosine scattering angle
76 * \return probability
77 */
78 inline double henyey_greenstein(const double g,
79 const double x)
80
81 {
82 const double a0 = (1.0 - g*g) / (4*PI);
83 const double y = 1.0 + g*g - 2.0*g*x;
84
85 return a0 / (y*sqrt(y));
86 }
87
88
89 /**
90 * Auxiliary method to describe light scattering in water (Heneyey-Greenstein).
91 *
92 * \param x cosine scattering angle
93 * \return probability
94 */
95 inline double henyey_greenstein(const double x)
96 {
97 static const double g = 0.924;
98
99 return henyey_greenstein(g, x);
100 }
101
102
103 /**
104 * Auxiliary method to describe light scattering in water (Rayleigh).
105 *
106 * \param a angular dependence parameter
107 * \param x cosine scattering angle
108 * \return probability
109 */
110 inline double rayleigh(const double a,
111 const double x)
112
113 {
114 const double a0 = 1.0 / (1.0 + a/3.0) / (4*PI);
115
116 return a0 * (1.0 + a*x*x);
117 }
118
119
120 /**
121 * Auxiliary method to describe light scattering in water (Rayleigh).
122 *
123 * \param x cosine scattering angle
124 * \return probability
125 */
126 inline double rayleigh(const double x)
127 {
128 return rayleigh(0.835, x);
129 }
130
131
132 /**
133 * Model specific function to describe light scattering in water (f4).
134 *
135 * \param x cosine scattering angle
136 * \return probability
137 */
138 inline double f4(const double x)
139 {
140 static const double g1 = 0.77;
141 static const double g2 = 0.00;
142 static const double f = 1.00;
143
144 return f * henyey_greenstein(g1,x) + (1.0 - f) * henyey_greenstein(g2,x);
145 }
146
147
148 /**
149 * Model specific function to describe light scattering in water (p00075).
150 *
151 * \param x cosine scattering angle
152 * \return probability
153 */
154 inline double p00075(const double x)
155 {
156 static const double g = 0.924;
157 static const double f = 0.25;
158
159 return f * rayleigh(x) + (1.0 - f) * henyey_greenstein(g,x);
160 }
161
162
163 /**
164 * Model specific function to describe light scattering in water (Petzhold).
165 *
166 * Values are taken from reference C.D. Mobley "Light and Water", ISBN 0-12-502750-8, pag. 111, Table 3.10 (right column).
167 *
168 * \param x cosine scattering angle
169 * \return probability
170 */
171 inline double petzhold(const double x)
172 {
173 using namespace JPP;
174
175 JPolint1Function1D_t f1;
176
177 if (f1.empty()) {
178 // angle probability
179 // [deg]
180 f1[ 0.100] = 1.767e+03;
181 f1[ 0.126] = 1.296e+03;
182 f1[ 0.158] = 9.502e+02;
183 f1[ 0.200] = 6.991e+02;
184 f1[ 0.251] = 5.140e+02;
185 f1[ 0.316] = 3.764e+02;
186 f1[ 0.398] = 2.763e+02;
187 f1[ 0.501] = 2.188e+02;
188 f1[ 0.631] = 1.444e+02;
189 f1[ 0.794] = 1.022e+02;
190 f1[ 1.000] = 7.161e+01;
191 f1[ 1.259] = 4.958e+01;
192 f1[ 1.585] = 3.395e+01;
193 f1[ 1.995] = 2.281e+01;
194 f1[ 2.512] = 1.516e+01;
195 f1[ 3.162] = 1.002e+01;
196 f1[ 3.981] = 6.580e+00;
197 f1[ 5.012] = 4.295e+00;
198 f1[ 6.310] = 2.807e+00;
199 f1[ 7.943] = 1.819e+00;
200 f1[ 10.000] = 1.153e+00;
201 f1[ 15.000] = 4.893e-01;
202 f1[ 20.000] = 2.444e-01;
203 f1[ 25.000] = 1.472e-01;
204 f1[ 30.000] = 8.609e-02;
205 f1[ 35.000] = 5.931e-02;
206 f1[ 40.000] = 4.210e-02;
207 f1[ 45.000] = 3.067e-02;
208 f1[ 50.000] = 2.275e-02;
209 f1[ 55.000] = 1.699e-02;
210 f1[ 60.000] = 1.313e-02;
211 f1[ 65.000] = 1.046e-02;
212 f1[ 70.000] = 8.488e-03;
213 f1[ 75.000] = 6.976e-03;
214 f1[ 80.000] = 5.842e-03;
215 f1[ 85.000] = 4.953e-03;
216 f1[ 90.000] = 4.292e-03;
217 f1[ 95.000] = 3.782e-03;
218 f1[100.000] = 3.404e-03;
219 f1[105.000] = 3.116e-03;
220 f1[110.000] = 2.912e-03;
221 f1[115.000] = 2.797e-03;
222 f1[120.000] = 2.686e-03;
223 f1[125.000] = 2.571e-03;
224 f1[130.000] = 2.476e-03;
225 f1[135.000] = 2.377e-03;
226 f1[140.000] = 2.329e-03;
227 f1[145.000] = 2.313e-03;
228 f1[150.000] = 2.365e-03;
229 f1[155.000] = 2.506e-03;
230 f1[160.000] = 2.662e-03;
231 f1[165.000] = 2.835e-03;
232 f1[170.000] = 3.031e-03;
233 f1[175.000] = 3.092e-03;
234 f1[180.000] = 3.154e-03;
235
236 f1.mul(1.166); // normalisation
237
238 f1.compile();
239 }
240
241 if (x >= +1.0)
242 return f1. begin()->getY();
243 else if (x <= -1.0)
244 return f1.rbegin()->getY();
245
246 const double a = acos(x) * 180.0 / PI;
247
248 if (a < f1. begin()->getX())
249 return f1. begin()->getY();
250 else if (a > f1.rbegin()->getX())
251 return f1.rbegin()->getY();
252 else
253 return f1(a);
254 }
255
256
257 /**
258 * Auxiliary data structure for scattering lengths of deep-sea water.
259 *
260 * Use the Kopelevich model for spectral volume scattering functions.\n
261 * The model separates the contributions by:
262 * - pure: pure sea water;
263 * - small: 'small' particles (size < 1 um);
264 * - large: 'large' particles (size > 1 um).
265 *
266 * Values are taken from reference C.D. Mobley "Light and Water", ISBN 0-12-502750-8, pag. 119.
267 */
268 struct JMobley {
269 /**
270 * Constructor.
271 *
272 * \param lambda wavelength of light [nm]
273 */
274 JMobley(const double lambda)
275 {
276 static const double Vs = 0.0075;
277 static const double Vl = 0.0075;
278 static const double bw = 0.0017;
279 static const double bs = 1.340;
280 static const double bl = 0.312;
281
282 const double x = 550.0/lambda;
283
284 Lis_pure = bw * pow(x, 4.3);
285 Lis_small = bs * Vs * pow(x, 1.7);
286 Lis_large = bl * Vl * pow(x, 0.3);
287 }
288
289 double Lis_pure; //!< inverse scattering length due to pure water [m^-1]
290 double Lis_small; //!< inverse scattering length due to small particles [m^-1]
291 double Lis_large; //!< inverse scattering length due to large particles [m^-1]
292 };
293
294
295 /**
296 * Absorption length of deep-sea water according Smith & Baker.
297 *
298 * Raymond C. Smith and Karen S. Baker, "Optical properties of the clearest natural waters (200800 nm)," Appl. Opt. 20, 177-184 (1981)
299 *
300 * \param lambda wavelength of light [nm]
301 * \return absorption length [m]
302 */
303 inline double getAbsorptionLengthSmithAndBaker(const double lambda)
304 {
305 static JGridSplineFunction1D_t f1;
306
307 if (f1.empty()) {
308 // wave- absorption
309 // length coefficient
310 // [nm] [1/m]
311 f1[200] = 3.0700;
312 f1[210] = 1.9900;
313 f1[220] = 1.3100;
314 f1[230] = 0.9270;
315 f1[240] = 0.7200;
316 f1[250] = 0.5590;
317 f1[260] = 0.4570;
318 f1[270] = 0.3730;
319 f1[280] = 0.2880;
320 f1[290] = 0.2150;
321 f1[300] = 0.1410;
322 f1[310] = 0.1050;
323 f1[320] = 0.0844;
324 f1[330] = 0.0678;
325 f1[340] = 0.0561;
326 f1[350] = 0.0463;
327 f1[360] = 0.0379;
328 f1[370] = 0.0300;
329 f1[380] = 0.0220;
330 f1[390] = 0.0191;
331 f1[400] = 0.0171;
332 f1[410] = 0.0162;
333 f1[420] = 0.0153;
334 f1[430] = 0.0144;
335 f1[440] = 0.0145;
336 f1[450] = 0.0145;
337 f1[460] = 0.0156;
338 f1[470] = 0.0156;
339 f1[480] = 0.0176;
340 f1[490] = 0.0196;
341 f1[500] = 0.0257;
342 f1[510] = 0.0357;
343 f1[520] = 0.0477;
344 f1[530] = 0.0507;
345 f1[540] = 0.0558;
346 f1[550] = 0.0638;
347 f1[560] = 0.0708;
348 f1[570] = 0.0799;
349 f1[580] = 0.1080;
350 f1[590] = 0.1570;
351 f1[600] = 0.2440;
352 f1[610] = 0.2890;
353 f1[620] = 0.3090;
354 f1[630] = 0.3190;
355 f1[640] = 0.3290;
356 f1[650] = 0.3490;
357 f1[660] = 0.4000;
358 f1[670] = 0.4300;
359 f1[680] = 0.4500;
360 f1[690] = 0.5000;
361 f1[700] = 0.6500;
362 f1[710] = 0.8390;
363 f1[720] = 1.1690;
364 f1[730] = 1.7990;
365 f1[740] = 2.3800;
366 f1[750] = 2.4700;
367 f1[760] = 2.5500;
368 f1[770] = 2.5100;
369 f1[780] = 2.3600;
370 f1[790] = 2.1600;
371 f1[800] = 2.0700;
372
373 f1.compile();
374 }
375
376 if (lambda > f1.getXmin() && lambda < f1.getXmax())
377 return 1.0 / f1(lambda);
378 else
379 return 0.0;
380 }
381
382
383 /**
384 * Absorption length of pure water.
385 *
386 * CITATION:
387 * Jonasz M. 2007. Absorption coefficient of water: Data sources (www.tpdsci.com/Tpc/AbsCfOfWaterDat.php).
388 * In: Top. Part. Disp. Sci. (www.tpdsci.com).
389 *
390 * DATA FROM:
391 * Wozniak B., Wozniak S. B., Tyszka K., Ostrowska M., Majchrowski R., Ficek D., Dera J. 2005.
392 * Modelling the light absorption properties of particulate matter forming organic particles suspended in seawater. Part 2.
393 * Modelling results. Oceanologia 47, 621-662.
394 * see also
395 * Wozniak B., Dera J. 2007.
396 * Light absorption in sea water. Springer, Berlin, 456 pp. (see p. 62)
397 *
398 * NOTES:
399 * As stated by the data authors, the data are based on measurement results obtained by various authors
400 * (interpolated by a linear approximation where applicable):
401 * Wavelength Reference
402 * - 200-335 nm Smith R.C., Baker K.S. 1981. Optical properties of the clearest natural waters (200-800 nm). Appl. Opt. 20, 177-184.
403 * - 340-370 nm Sogandares F.M., Fry, E.S. 1997. Absorption spectrum (340 -640 nm) of pure water. I. Photothermal measurements Appl. Opt. 36, 8699-8709.
404 * - 380-700 nm Pope R.M., Fry E.S. 1997. Absorption spectrum (380 -700 nm) of pure water. II. Integrating cavity measurements. Appl. Opt. 36, 8710-8723
405 *
406 * \param lambda wavelength of light [nm]
407 * \return absorption length [m]
408 */
409 inline double getAbsorptionLengthOfPureWater(const double lambda)
410 {
411 static JGridSplineFunction1D_t f1;
412
413 if (f1.empty()) {
414 // wave- absorption
415 // length coefficient
416 // [nm] [1/m]
417 f1[0.200e3] = 3.07;
418 f1[0.205e3] = 2.53;
419 f1[0.210e3] = 1.99;
420 f1[0.215e3] = 1.65;
421 f1[0.220e3] = 1.31;
422 f1[0.225e3] = 1.1185;
423 f1[0.230e3] = 0.927;
424 f1[0.235e3] = 0.8235;
425 f1[0.240e3] = 0.72;
426 f1[0.245e3] = 0.6395;
427 f1[0.250e3] = 0.559;
428 f1[0.255e3] = 0.508;
429 f1[0.260e3] = 0.457;
430 f1[0.265e3] = 0.415;
431 f1[0.270e3] = 0.373;
432 f1[0.275e3] = 0.3305;
433 f1[0.280e3] = 0.288;
434 f1[0.285e3] = 0.2515;
435 f1[0.290e3] = 0.215;
436 f1[0.295e3] = 0.178;
437 f1[0.300e3] = 0.141;
438 f1[0.305e3] = 0.123;
439 f1[0.310e3] = 0.105;
440 f1[0.315e3] = 0.0947;
441 f1[0.320e3] = 0.0844;
442 f1[0.325e3] = 0.0761;
443 f1[0.330e3] = 0.0678;
444 f1[0.335e3] = 0.06195;
445 f1[0.340e3] = 0.0325;
446 f1[0.345e3] = 0.02645;
447 f1[0.350e3] = 0.0204;
448 f1[0.355e3] = 0.018;
449 f1[0.360e3] = 0.0156;
450 f1[0.365e3] = 0.0135;
451 f1[0.370e3] = 0.0114;
452 f1[0.375e3] = 0.011385;
453 f1[0.380e3] = 0.01137;
454 f1[0.385e3] = 0.00941;
455 f1[0.390e3] = 0.00851;
456 f1[0.395e3] = 0.00813;
457 f1[0.400e3] = 0.00663;
458 f1[0.405e3] = 0.0053;
459 f1[0.410e3] = 0.00473;
460 f1[0.415e3] = 0.00444;
461 f1[0.420e3] = 0.00454;
462 f1[0.425e3] = 0.00478;
463 f1[0.430e3] = 0.00495;
464 f1[0.435e3] = 0.0053;
465 f1[0.440e3] = 0.00635;
466 f1[0.445e3] = 0.00751;
467 f1[0.450e3] = 0.00922;
468 f1[0.455e3] = 0.00962;
469 f1[0.460e3] = 0.00979;
470 f1[0.465e3] = 0.01011;
471 f1[0.470e3] = 0.0106;
472 f1[0.475e3] = 0.0114;
473 f1[0.480e3] = 0.0127;
474 f1[0.485e3] = 0.0136;
475 f1[0.490e3] = 0.015;
476 f1[0.495e3] = 0.0173;
477 f1[0.500e3] = 0.0204;
478 f1[0.505e3] = 0.0256;
479 f1[0.510e3] = 0.0325;
480 f1[0.515e3] = 0.0396;
481 f1[0.520e3] = 0.0409;
482 f1[0.525e3] = 0.0417;
483 f1[0.530e3] = 0.0434;
484 f1[0.535e3] = 0.0452;
485 f1[0.540e3] = 0.0474;
486 f1[0.545e3] = 0.0511;
487 f1[0.550e3] = 0.0565;
488 f1[0.555e3] = 0.0596;
489 f1[0.560e3] = 0.0619;
490 f1[0.565e3] = 0.0642;
491 f1[0.570e3] = 0.0695;
492 f1[0.575e3] = 0.0772;
493 f1[0.580e3] = 0.0896;
494 f1[0.585e3] = 0.11;
495 f1[0.590e3] = 0.1351;
496 f1[0.595e3] = 0.1672;
497 f1[0.600e3] = 0.2224;
498 f1[0.605e3] = 0.2577;
499 f1[0.610e3] = 0.2644;
500 f1[0.615e3] = 0.2678;
501 f1[0.620e3] = 0.2755;
502 f1[0.625e3] = 0.2834;
503 f1[0.630e3] = 0.2916;
504 f1[0.635e3] = 0.3012;
505 f1[0.640e3] = 0.3108;
506 f1[0.645e3] = 0.325;
507 f1[0.650e3] = 0.34;
508 f1[0.655e3] = 0.371;
509 f1[0.660e3] = 0.41;
510 f1[0.665e3] = 0.429;
511 f1[0.670e3] = 0.439;
512 f1[0.675e3] = 0.448;
513 f1[0.680e3] = 0.465;
514 f1[0.685e3] = 0.486;
515 f1[0.690e3] = 0.516;
516 f1[0.695e3] = 0.559;
517 f1[0.700e3] = 0.624;
518
519 f1.compile();
520 }
521
522 if (lambda > f1.getXmin() && lambda < f1.getXmax())
523 return 1.0 / f1(lambda);
524 else
525 return 0.0;
526 }
527}
528
529#endif
JFunction1D_t.hh
JConstants.hh
Physics constants.
g1
Double_t g1(const Double_t x)
Function.
Definition JQuantiles.cc:25
JPHYSICS
Auxiliary methods for light properties of deep-sea water.
Definition JAbstractMedium.hh:9
JPHYSICS::henyey_greenstein
double henyey_greenstein(const double g, const double x)
Auxiliary method to describe light scattering in water (Henyey-Greenstein).
Definition JPhysicsSupportkit.hh:78
JPHYSICS::p00075
double p00075(const double x)
Model specific function to describe light scattering in water (p00075).
Definition JPhysicsSupportkit.hh:154
JPHYSICS::getThetaMCS
double getThetaMCS(const double E, const double x, const double X0, const double M, const double Q)
Get multiple Coulomb scattering angle.
Definition JPhysicsSupportkit.hh:40
JPHYSICS::petzhold
double petzhold(const double x)
Model specific function to describe light scattering in water (Petzhold).
Definition JPhysicsSupportkit.hh:171
JPHYSICS::rayleigh
double rayleigh(const double a, const double x)
Auxiliary method to describe light scattering in water (Rayleigh).
Definition JPhysicsSupportkit.hh:110
JPHYSICS::X0_WATER_M
static const double X0_WATER_M
Radiation length pure water [m].
Definition JPhysics/JConstants.hh:48
JPHYSICS::getAbsorptionLengthSmithAndBaker
double getAbsorptionLengthSmithAndBaker(const double lambda)
Absorption length of deep-sea water according Smith & Baker.
Definition JPhysicsSupportkit.hh:303
JPHYSICS::f4
double f4(const double x)
Model specific function to describe light scattering in water (f4).
Definition JPhysicsSupportkit.hh:138
JPHYSICS::getAbsorptionLengthOfPureWater
double getAbsorptionLengthOfPureWater(const double lambda)
Absorption length of pure water.
Definition JPhysicsSupportkit.hh:409
JPHYSICS::THETA_MCS
static const double THETA_MCS
Multiple Coulomb scattering constant [GeV].
Definition JPhysics/JConstants.hh:33
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition JAAnetToolkit.hh:48
JPHYSICS::JMobley
Auxiliary data structure for scattering lengths of deep-sea water.
Definition JPhysicsSupportkit.hh:268
JPHYSICS::JMobley::Lis_large
double Lis_large
inverse scattering length due to large particles [m^-1]
Definition JPhysicsSupportkit.hh:291
JPHYSICS::JMobley::JMobley
JMobley(const double lambda)
Constructor.
Definition JPhysicsSupportkit.hh:274
JPHYSICS::JMobley::Lis_small
double Lis_small
inverse scattering length due to small particles [m^-1]
Definition JPhysicsSupportkit.hh:290
JPHYSICS::JMobley::Lis_pure
double Lis_pure
inverse scattering length due to pure water [m^-1]
Definition JPhysicsSupportkit.hh:289
JTOOLS::JGridSplineFunction1D_t
Type definition of a spline interpolation based on a JGridCollection with result type double.
Definition JFunction1D_t.hh:61
JTOOLS::JPolint1Function1D_t
Type definition of a 1st degree polynomial interpolation with result type double.
Definition JFunction1D_t.hh:187
Generated by doxygen 1.12.0