38 return (option ? 1.0 + value: value);
48 return (option ? sqrt(2.0 + value) * sqrt(0.0 - value) : sqrt(1.0 + value) * sqrt(1.0 - value));
80 JBellhop_t(
const sin_t& s,
105 static constexpr
double STEP_SIZE_M = 0.001;
115 JBellhop(
const JAbstractSoundVelocity&
V,
const double precision) :
130 sin_t operator()(
const double x0,
133 const double z1)
const
137 const double st = fabs(x1 - x0) / hypot(x1 - x0, z1 - z0);
139 double st_min = (z1 > z0 ? st : 0.0);
140 double st_max = (z1 > z0 ? 1.0 : st);
142 if (fabs(x1 - x0) > fabs(z1 - z0)) {
149 double s0 = 0.5*(st_min + st_max);
155 while ((x - x1) * copysign(1.0, x1 - x0) < -0.5 * STEP_SIZE_M) {
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);
161 s += (1.0 + s) * (
V(z + dz) -
V(z)) /
V(z + 0.5*dz);
169 }
else if (s <= -1.0) {
171 z = numeric_limits<double>::max();
177 if (fabs(z - z1) <= precision) {
181 if ((z - z1) * copysign(1.0, z1 - z0) > 0.0)
191 double s0 = 0.5*(st_min + st_max);
197 while ((z - z1) * copysign(1.0, z1 - z0) < -0.5 * STEP_SIZE_M) {
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);
203 s += s * (
V(z + dz) -
V(z)) /
V(z + 0.5*dz);
209 x = numeric_limits<double>::max();
213 }
else if (s < 0.0) {
219 if (fabs(x - x1) <= precision) {
220 return {
false, s0 };
223 if ((x - x1) * copysign(1.0, x1 - x0) < 0.0)
244 JBellhop_t operator()(
const sin_t& s0,
248 const double z1)
const
260 while ((x - x1) * copysign(1.0, x1 - x0) < -0.5 * STEP_SIZE_M) {
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);
267 s += (1.0 + s) * (
V(z + dz) -
V(z)) / v;
270 d += sqrt(dx*dx + dz*dz);
271 t += sqrt(dx*dx + dz*dz) /
v;
281 while ((z - z1) * copysign(1.0, z1 - z0) < -0.5 * STEP_SIZE_M) {
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);
288 s += s * (
V(z + dz) -
V(z)) / v;
291 d += sqrt(dx*dx + dz*dz);
292 t += sqrt(dx*dx + dz*dz) /
v;
301 return JBellhop_t({s0.option,s},
x, z, t,
d);
305 const JAbstractSoundVelocity&
V;
317 int main(
int argc,
char **argv)
329 JParser<> zap(
"Example program for sound ray tracing.");
338 catch(
const exception &error) {
339 FATAL(error.what() << endl);
345 const JBellhop bellhop(
V, precision);
360 cout <<
"> " << flush;
361 cin >> x0 >> z0 >> x1 >> z1;
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);
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;
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;
389 cout <<
"v0 " <<
FIXED (12,5) << v0 <<
" [m/s]" << endl;
390 cout <<
"v1 " <<
FIXED (12,5) << v1 <<
" [m/s]" << endl;
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;
400 dx = (sqrt(1.0 + st) * sqrt(1.0 - st) - sqrt(1.0 + c * st) * sqrt(1.0 - c * st)) * v0 / ( st ) / b;
403 cout <<
"x' " <<
FIXED (12,5) << dx <<
" [m]" << endl;
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));
408 cout <<
"t' " <<
FIXED (12,5) << (atanh(x) - atanh(y)) / b << endl;
410 catch(
const exception& error) {
411 cerr << error.what() << endl;
417 const double x0 = 0.0;
418 const double z0 = 0.0;
425 TH2D h2(
"h2", NULL, 50, 1.0, 3.5, 50, -1.0, 3.0);
427 for (
int ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
428 for (
int iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {
430 const double x1 =
pow(10.0, h2.GetXaxis()->GetBinCenter(ix));
431 const double z1 =
pow(10.0, h2.GetYaxis()->GetBinCenter(iy));
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);
442 Qx.put(result.x - x1);
443 Qz.put(result.z - z1);
445 h2.SetBinContent(ix, iy, (result.t -
V.getTime(D, z0, z1)) * 1.0e6);
447 catch(
const exception& error) {
448 ERROR(error.what() << endl);
Utility class to parse command line options.
int main(int argc, char *argv[])
Auxiliary base class for aritmetic operations of derived class types.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Empty structure for specification of parser element that is initialised (i.e. does not require input)...
Auxiliary data structure for floating point format specification.
V(JDAQEvent-JTriggerReprocessor)*1.0/(JDAQEvent+1.0e-10)
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
T pow(const T &x, const double y)
Power .
General purpose messaging.
$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
Utility class to parse command line options.
Interface for depth dependend velocity of sound.
Auxiliary data structure for floating point format specification.
do echo Generating $dir eval D
#define DEBUG(A)
Message macros.