JBellhop.cc File Reference

Example program for sound ray tracing. More...

#include <iostream>
#include <iomanip>
#include <limits>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TH2D.h"
#include "JAcoustics/JSoundVelocity.hh"
#include "JLang/JException.hh"
#include "JTools/JQuantile.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Example program for sound ray tracing.

Author

mdejong

Definition in file JBellhop.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 317 of file JBellhop.cc.

{
  using namespace std;
  using namespace JPP;

  string                  outputFile;
  JSoundVelocity          V              = getSoundVelocity;    // default sound velocity
  double                  precision;
  int                     debug;

  try {

    JParser<> zap("Example program for sound ray tracing.");

    zap['o'] = make_field(outputFile)          = "";
    zap['V'] = make_field(V,            "sound velocity")                    = JPARSER::initialised();
    zap['e'] = make_field(precision)        = 1.0e-3;
    zap['d'] = make_field(debug)            = 2;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  V.set(-3500);

  const JBellhop bellhop(V, precision);

  if (outputFile == "") {

    double x0, z0;
    double x1, z1;
    /*
    x0 =    0.0;  z0 =  0.0;
    x1 = 2000.0;  z1 = 10.0;

    {
    */
    //
    for ( ; ; ) {

      cout << "> " << flush;
      cin  >> x0 >> z0 >> x1 >> z1;
      //

      try {
        
        const double     D      = hypot(x1 - x0, z1 - z0);
        const sin_t      s0     = bellhop(x0, z0, x1, z1); 
        const JBellhop_t result = bellhop(s0, x0, z0, x1, z1);
      
        cout << "s  " << SCIENTIFIC(12,5) << (x1 - x0) / D                    << endl;
        cout << "s0 " << SCIENTIFIC(12,5) << s0.value << ' ' << s0.option     << endl;
        cout << "s1 " << SCIENTIFIC(12,5) << result.s.value                   << endl;
        cout << "x' " << FIXED     (12,5) << result.x             << " [m]"   << endl;
        cout << "z' " << FIXED     (12,5) << result.z             << " [m]"   << endl;
        cout << "D  " << FIXED     (12,5) << D                    << " [m]"   << endl;
        cout << "D' " << FIXED     (12,5) << result.d             << " [m]"   << endl;
        cout << "t  " << FIXED     (12,5) << V.getTime(D, z0, z1) << " [s]"   << endl;
        cout << "t' " << FIXED     (12,5) << result.t             << " [s]"   << endl;


        const double v0 = V(z0);
        const double v1 = V(z1);
        const double b  = (v1 - v0) / (z1 - z0);
        const double c  =  v1 / v0;
        const double u  = (v0 + v1) / v1;
        const double v  = (v0 - v1) / v1;
        const double st = s0.value;

        cout << "v0 " << FIXED     (12,5) << v0                   << " [m/s]" << endl;
        cout << "v1 " << FIXED     (12,5) << v1                   << " [m/s]" << endl;

        double dx = 0.0;

        if (s0.option) {

          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;

        } else {

          dx = (sqrt(1.0 + st) * sqrt(1.0 - st) - sqrt(1.0   +   c * st) * sqrt(1.0   -   c * st)) * v0 / (   st   ) / b;
        }

        cout << "x' " << FIXED     (12,5) << dx                   << " [m]"   << endl;

        const double x =  s0.getc();
        const double y = (s0.option ? c * sqrt(u + st) * sqrt(v - st) : sqrt(1.0 + c * st) * sqrt(1.0 - c * st));
        
        cout << "t' " << FIXED     (12,5) << (atanh(x) - atanh(y)) / b        << endl;
      }
      catch(const exception& error) {
        cerr << error.what() << endl;
      }
    }

  } else {
    
    const double x0 =    0.0;
    const double z0 =    0.0;

    JQuantile Qx("x");
    JQuantile Qz("z");

    TFile out(outputFile.c_str(), "recreate");

    TH2D h2("h2", NULL, 50, 1.0, 3.5, 50, -1.0, 3.0);

    for (int ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
      for (int iy = 1; iy <= h2.GetYaxis()->GetNbins(); ++iy) {
      
        const double x1 = pow(10.0, h2.GetXaxis()->GetBinCenter(ix));
        const double z1 = pow(10.0, h2.GetYaxis()->GetBinCenter(iy));

        STATUS("x1 " << FIXED(12,6) << x1                   << " [m]" << ' ');  DEBUG(endl);
        STATUS("z1 " << FIXED(12,6) << z1                   << " [m]" << '\r'); DEBUG(endl);

        try {

          const double     D      = hypot(x1 - x0, z1 - z0);
          const sin_t      s0     = bellhop(x0, z0, x1, z1); 
          const JBellhop_t result = bellhop(s0, x0, z0, x1, z1);

          Qx.put(result.x - x1);
          Qz.put(result.z - z1);

          h2.SetBinContent(ix, iy, (result.t - V.getTime(D, z0, z1)) * 1.0e6);  // [us]
        }
        catch(const exception& error) {
          ERROR(error.what() << endl);
        }
      }
    }

    Qx.print(cout);
    Qz.print(cout);

    out.Write();
    out.Close();
  }
}