Jpp  18.4.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JBellhop.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <iomanip>
3 #include <limits>
4 #include <cmath>
5 
6 #include "TROOT.h"
7 #include "TFile.h"
8 #include "TH2D.h"
9 
11 #include "JLang/JException.hh"
12 #include "JTools/JQuantile.hh"
13 
14 #include "Jeep/JParser.hh"
15 #include "Jeep/JMessage.hh"
16 
17 
18 namespace {
19 
21  using JLANG::JException;
22  using JMATH::JMath;
23 
24 
25  /**
26  * Auxiliary fata structure for high-precision sinus of angle.
27  *
28  * If sin_t::option is <tt>true</tt>, value is offset by one; else value is offset by zero.
29  */
30  struct sin_t {
31  /**
32  * Get sin.
33  *
34  * \return sin
35  */
36  double gets() const
37  {
38  return (option ? 1.0 + value: value);
39  }
40 
41  /**
42  * Get cos.
43  *
44  * \return cos
45  */
46  double getc() const
47  {
48  return (option ? sqrt(2.0 + value) * sqrt(0.0 - value) : sqrt(1.0 + value) * sqrt(1.0 - value));
49  }
50 
51  bool option; //!< option
52  double value; //!< value
53  };
54 
55 
56  /**
57  * Auxiliary data structure for tracing of sound ray.
58  */
59  struct JBellhop_t {
60  /**
61  * Default constructor.
62  */
63  JBellhop_t() :
64  s({false,0.0}),
65  x(0.0),
66  z(0.0),
67  t(0.0),
68  d(0.0)
69  {}
70 
71  /**
72  * Constructor.
73  *
74  * \param s sin angle with respect to vertical
75  * \param x x-position [m]
76  * \param z z-position [m]
77  * \param t time [s]
78  * \param d distance [m]
79  */
80  JBellhop_t(const sin_t& s,
81  const double x,
82  const double z,
83  const double t,
84  const double d) :
85  s(s),
86  x(x),
87  z(z),
88  t(t),
89  d(d)
90  {}
91 
92  sin_t s; //!< sinus angle with respect to vertical
93  double x; //!< x-position [m]
94  double z; //!< z-position [m]
95  double t; //!< time [s]
96  double d; //!< distance [m]
97  };
98 
99 
100  /**
101  * Data structure for tracing of sound ray.
102  */
103  struct JBellhop {
104 
105  static constexpr double STEP_SIZE_M = 0.001; //!< step size [m]
106  static constexpr size_t NUMBER_OF_ITERATIONS = 1000; //!< maximum number of iterations
107 
108 
109  /**
110  * Constructor.
111  *
112  * \param V sound velocity
113  * \param precision precision
114  */
115  JBellhop(const JAbstractSoundVelocity& V, const double precision) :
116  V(V),
117  precision(precision)
118  {}
119 
120 
121  /**
122  * Estimate sinus of angle of sound ray with respect to vertical at start position to reach end position.
123  *
124  * \param x0 x-position at start
125  * \param z0 z-position at start
126  * \param x1 x-position at end
127  * \param z1 z-position at end
128  * \return sinus of angle of sound ray
129  */
130  sin_t operator()(const double x0,
131  const double z0,
132  const double x1,
133  const double z1) const
134  {
135  using namespace std;
136 
137  const double st = fabs(x1 - x0) / hypot(x1 - x0, z1 - z0);
138 
139  double st_min = (z1 > z0 ? st : 0.0);
140  double st_max = (z1 > z0 ? 1.0 : st);
141 
142  if (fabs(x1 - x0) > fabs(z1 - z0)) {
143 
144  st_min -= 1.0; // precision
145  st_max -= 1.0; //
146 
147  for (size_t i = 0; i != NUMBER_OF_ITERATIONS; ++i) {
148 
149  double s0 = 0.5*(st_min + st_max);
150 
151  double s = s0;
152  double x = x0;
153  double z = z0;
154 
155  while ((x - x1) * copysign(1.0, x1 - x0) < -0.5 * STEP_SIZE_M) {
156 
157  const double c = sqrt(2.0 + s) * sqrt(-s);
158  const double dx = copysign(STEP_SIZE_M, x1 - x0);
159  const double dz = copysign(STEP_SIZE_M * c / (1.0 + s), z1 - z0);
160 
161  s += (1.0 + s) * (V(z + dz) - V(z)) / V(z + 0.5*dz);
162  x += dx;
163  z += dz;
164 
165  if (s > 0.0) {
166 
167  s = 0.0;
168 
169  } else if (s <= -1.0) {
170 
171  z = numeric_limits<double>::max();
172 
173  break;
174  }
175  }
176 
177  if (fabs(z - z1) <= precision) {
178  return { true, s0 };
179  }
180 
181  if ((z - z1) * copysign(1.0, z1 - z0) > 0.0)
182  st_min = s0;
183  else
184  st_max = s0;
185  }
186 
187  } else {
188 
189  for (size_t i = 0; i != NUMBER_OF_ITERATIONS; ++i) {
190 
191  double s0 = 0.5*(st_min + st_max);
192 
193  double s = s0;
194  double x = x0;
195  double z = z0;
196 
197  while ((z - z1) * copysign(1.0, z1 - z0) < -0.5 * STEP_SIZE_M) {
198 
199  const double c = sqrt((1.0 + s) * (1.0 - s));
200  const double dx = copysign(STEP_SIZE_M * s / c, x1 - x0);
201  const double dz = copysign(STEP_SIZE_M, z1 - z0);
202 
203  s += s * (V(z + dz) - V(z)) / V(z + 0.5*dz);
204  x += dx;
205  z += dz;
206 
207  if (s >= 1.0) {
208 
209  x = numeric_limits<double>::max();
210 
211  break;
212 
213  } else if (s < 0.0) {
214 
215  s = 0.0;
216  }
217  }
218 
219  if (fabs(x - x1) <= precision) {
220  return { false, s0 };
221  }
222 
223  if ((x - x1) * copysign(1.0, x1 - x0) < 0.0)
224  st_min = s0;
225  else
226  st_max = s0;
227  }
228  }
229 
230  THROW(JException, "Reached maximum number of iterations " << NUMBER_OF_ITERATIONS);
231  }
232 
233 
234  /**
235  * Sound ray tracing.
236  *
237  * \param s0 sinus angle of sound ray
238  * \param x0 x-position at start
239  * \param z0 z-position at start
240  * \param x1 x-position at end
241  * \param z1 z-position at end
242  * \return sound ray trace
243  */
244  JBellhop_t operator()(const sin_t& s0,
245  const double x0,
246  const double z0,
247  const double x1,
248  const double z1) const
249  {
250  using namespace std;
251 
252  double x = x0;
253  double z = z0;
254  double s = s0.value;
255  double t = 0.0;
256  double d = 0.0;
257 
258  if (s0.option) {
259 
260  while ((x - x1) * copysign(1.0, x1 - x0) < -0.5 * STEP_SIZE_M) {
261 
262  const double c = sqrt(2.0 + s) * sqrt(-s);
263  const double dx = copysign(STEP_SIZE_M, x1 - x0);
264  const double dz = copysign(STEP_SIZE_M * c / (1.0 + s), z1 - z0);
265  const double v = V(z + 0.5*dz);
266 
267  s += (1.0 + s) * (V(z + dz) - V(z)) / v;
268  x += dx;
269  z += dz;
270  d += sqrt(dx*dx + dz*dz);
271  t += sqrt(dx*dx + dz*dz) / v;
272 
273  if (s > 0.0)
274  s = 0.0;
275  else if (s < -1.0)
276  s = -1.0;
277  }
278 
279  } else {
280 
281  while ((z - z1) * copysign(1.0, z1 - z0) < -0.5 * STEP_SIZE_M) {
282 
283  const double c = sqrt((1.0 + s) * (1.0 - s));
284  const double dx = copysign(STEP_SIZE_M * s / c, x1 - x0);
285  const double dz = copysign(STEP_SIZE_M, z1 - z0);
286  const double v = V(z + 0.5*dz);
287 
288  s += s * (V(z + dz) - V(z)) / v;
289  x += dx;
290  z += dz;
291  d += sqrt(dx*dx + dz*dz);
292  t += sqrt(dx*dx + dz*dz) / v;
293 
294  if (s > 1.0)
295  s = 1.0;
296  else if (s < 0.0)
297  s = 0.0;
298  }
299  }
300 
301  return JBellhop_t({s0.option,s}, x, z, t, d);
302  }
303 
304 
305  const JAbstractSoundVelocity& V;
306  double precision;
307  };
308 }
309 
310 
311 /**
312  * \file
313  *
314  * Example program for sound ray tracing.
315  * \author mdejong
316  */
317 int main(int argc, char **argv)
318 {
319  using namespace std;
320  using namespace JPP;
321 
322  string outputFile;
323  JSoundVelocity V = getSoundVelocity; // default sound velocity
324  double precision;
325  int debug;
326 
327  try {
328 
329  JParser<> zap("Example program for sound ray tracing.");
330 
331  zap['o'] = make_field(outputFile) = "";
332  zap['V'] = make_field(V, "sound velocity") = JPARSER::initialised();
333  zap['e'] = make_field(precision) = 1.0e-3;
334  zap['d'] = make_field(debug) = 2;
335 
336  zap(argc, argv);
337  }
338  catch(const exception &error) {
339  FATAL(error.what() << endl);
340  }
341 
342 
343  V.set(-3500);
344 
345  const JBellhop bellhop(V, precision);
346 
347  if (outputFile == "") {
348 
349  double x0, z0;
350  double x1, z1;
351  /*
352  x0 = 0.0; z0 = 0.0;
353  x1 = 2000.0; z1 = 10.0;
354 
355  {
356  */
357  //
358  for ( ; ; ) {
359 
360  cout << "> " << flush;
361  cin >> x0 >> z0 >> x1 >> z1;
362  //
363 
364  try {
365 
366  const double D = hypot(x1 - x0, z1 - z0);
367  const sin_t s0 = bellhop(x0, z0, x1, z1);
368  const JBellhop_t result = bellhop(s0, x0, z0, x1, z1);
369 
370  cout << "s " << SCIENTIFIC(12,5) << (x1 - x0) / D << endl;
371  cout << "s0 " << SCIENTIFIC(12,5) << s0.value << ' ' << s0.option << endl;
372  cout << "s1 " << SCIENTIFIC(12,5) << result.s.value << endl;
373  cout << "x' " << FIXED (12,5) << result.x << " [m]" << endl;
374  cout << "z' " << FIXED (12,5) << result.z << " [m]" << endl;
375  cout << "D " << FIXED (12,5) << D << " [m]" << endl;
376  cout << "D' " << FIXED (12,5) << result.d << " [m]" << endl;
377  cout << "t " << FIXED (12,5) << V.getTime(D, z0, z1) << " [s]" << endl;
378  cout << "t' " << FIXED (12,5) << result.t << " [s]" << endl;
379 
380 
381  const double v0 = V(z0);
382  const double v1 = V(z1);
383  const double b = (v1 - v0) / (z1 - z0);
384  const double c = v1 / v0;
385  const double u = (v0 + v1) / v1;
386  const double v = (v0 - v1) / v1;
387  const double st = s0.value;
388 
389  cout << "v0 " << FIXED (12,5) << v0 << " [m/s]" << endl;
390  cout << "v1 " << FIXED (12,5) << v1 << " [m/s]" << endl;
391 
392  double dx = 0.0;
393 
394  if (s0.option) {
395 
396  dx = (sqrt(2.0 + st) * sqrt(0.0 - st) - sqrt(1.0 + c + c * st) * sqrt(1.0 - c - c * st)) * v0 / (1.0 + st) / b;
397 
398  } else {
399 
400  dx = (sqrt(1.0 + st) * sqrt(1.0 - st) - sqrt(1.0 + c * st) * sqrt(1.0 - c * st)) * v0 / ( st ) / b;
401  }
402 
403  cout << "x' " << FIXED (12,5) << dx << " [m]" << endl;
404 
405  const double x = s0.getc();
406  const double y = (s0.option ? c * sqrt(u + st) * sqrt(v - st) : sqrt(1.0 + c * st) * sqrt(1.0 - c * st));
407 
408  cout << "t' " << FIXED (12,5) << (atanh(x) - atanh(y)) / b << endl;
409  }
410  catch(const exception& error) {
411  cerr << error.what() << endl;
412  }
413  }
414 
415  } else {
416 
417  const double x0 = 0.0;
418  const double z0 = 0.0;
419 
420  JQuantile Qx("x");
421  JQuantile Qz("z");
422 
423  TFile out(outputFile.c_str(), "recreate");
424 
425  TH2D h2("h2", NULL, 50, 1.0, 3.5, 50, -1.0, 3.0);
426 
427  for (int ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
428  for (int iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {
429 
430  const double x1 = pow(10.0, h2.GetXaxis()->GetBinCenter(ix));
431  const double z1 = pow(10.0, h2.GetYaxis()->GetBinCenter(iy));
432 
433  STATUS("x1 " << FIXED(12,6) << x1 << " [m]" << ' '); DEBUG(endl);
434  STATUS("z1 " << FIXED(12,6) << z1 << " [m]" << '\r'); DEBUG(endl);
435 
436  try {
437 
438  const double D = hypot(x1 - x0, z1 - z0);
439  const sin_t s0 = bellhop(x0, z0, x1, z1);
440  const JBellhop_t result = bellhop(s0, x0, z0, x1, z1);
441 
442  Qx.put(result.x - x1);
443  Qz.put(result.z - z1);
444 
445  h2.SetBinContent(ix, iy, (result.t - V.getTime(D, z0, z1)) * 1.0e6); // [us]
446  }
447  catch(const exception& error) {
448  ERROR(error.what() << endl);
449  }
450  }
451  }
452 
453  Qx.print(cout);
454  Qz.print(cout);
455 
456  out.Write();
457  out.Close();
458  }
459 }
Utility class to parse command line options.
Definition: JParser.hh:1514
General exception.
Definition: JException.hh:24
Exceptions.
int main(int argc, char *argv[])
Definition: Main.cc:15
Auxiliary base class for aritmetic operations of derived class types.
Definition: JMath.hh:109
Sound velocity.
#define STATUS(A)
Definition: JMessage.hh:63
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Definition: JParser.hh:83
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:446
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
string outputFile
static const JSoundVelocity getSoundVelocity(1541.0,-17.0e-3,-2000.0)
Function object for velocity of sound.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:1989
#define ERROR(A)
Definition: JMessage.hh:66
T pow(const T &x, const double y)
Power .
Definition: JMath.hh:97
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
$WORKDIR ev_configure_dqsimulator txt echo process $DQ_SIMULATOR $i $SOURCE_HOST[$index] csh c(setenv ROOTSYS $ROOTSYS &&source $JPP_DIR/setenv.csh $JPP_DIR &&($DQ_SIMULATOR\-u\$NAME\$\-H\$SERVER\$\-M\$LOGGER\$\-d $DEBUG</dev/null > &/dev/null &))'
then JMuonMCEvt f $INPUT_FILE o $INTERMEDIATE_FILE d
Definition: JMuonPath.sh:47
Utility class to parse command line options.
Interface for depth dependend velocity of sound.
data_type v[N+1][M+1]
Definition: JPolint.hh:866
double u[N+1]
Definition: JPolint.hh:865
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:486
do echo Generating $dir eval D
Definition: JDrawLED.sh:53
int debug
debug level
#define DEBUG(A)
Message macros.
Definition: JMessage.hh:62