UNESCO.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <string>

#include "TROOT.h"
#include "TFile.h"
#include "TH2D.h"

#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

/**
 * Get sound velocity.
 *
 * \param  D            Depth [m]
 * \param  S            Salinity [ppk]
 * \param  T            Temperature [C]
 * \return              sound velocity [m/s]
 */
inline double getVelocity(const double D, const double S, const double T)
{
  return (1448.96
          + 4.591 * T
          - 5.304e-2  * T*T
          + 2.374e-4  * T*T*T
          + 1.340 * (S - 35)
          + 1.630e-2  * D
          + 1.675e-7  * D*D
          - 1.025e-2  * T * (S - 35)
          - 7.139e-13 * T*D*D*D);
}

/**
 * \file
 *
 * Example program to plot UNESCO sound velocity.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string   outputFile;
  double   D;
  int      debug;

  try {

    JParser<> zap("Example program to plot UNESCO sound velocity.");

    zap['o'] = make_field(outputFile);
    zap['D'] = make_field(D,  "Depth [m]");
    zap['d'] = make_field(debug)            = 2;

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

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

  TH2D h2("h2", NULL, 1000, 10.0, 20.0, 1000, 25.0, 40.0);

  for (Int_t ix = 1; ix <= h2.GetXaxis()->GetNbins(); ++ix) {
    for (Int_t iy = 1; iy <= h2.GetXaxis()->GetNbins(); ++iy) {

      const Double_t T = h2.GetXaxis()->GetBinCenter(ix);
      const Double_t S = h2.GetYaxis()->GetBinCenter(iy);

      const double V   = getVelocity(D, S, T);
      const double vs  = getVelocity(D, S, T)  -  getVelocity(D + 1, S, T);
      
      h2.SetBinContent(ix, iy, V);
    }
  }

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