getGeographicalLocation.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <sstream>
#include <iomanip>
 
#include "JLang/JManip.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
 
#include "JCompass/JCompassSupportkit.hh"
 
#include "JLang/JException.hh"
 
#include "JAstronomy/JAstronomy.hh"
 
#include "LatLong-UTMconversion.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
namespace {
 
  using JUTM::JUTMGrid;
  using JUTM::JUTMPosition;
  using JLANG::JParseError;
  using JASTRONOMY::JGeographicalLocation;
 
 
  /**
   * Get Earth's ellipsoid.
   *
   * \param  grid               UTM grid
   * \return                    ellipsoid
   */
  inline const Ellipsoid& getEllipsoid(const JUTMGrid& grid)
  {
    using namespace std;
 
    string buffer = grid.getWGS();
    size_t pos    = buffer.find_first_of("0123456789");
 
    const string wgs = buffer.substr(0, pos);
 
    int zone;
 
    istringstream(buffer.substr(pos)) >> zone;
  
    for (int i = 0; i != sizeof(ellipsoid)/sizeof(ellipsoid[0]); ++i) {
 
      buffer = ellipsoid[i].ellipsoidName;
      pos    = buffer.find(wgs);
 
      if (pos != string::npos) {
 
        pos += wgs.size() + 1;
 
        int val;
      
        istringstream(buffer.substr(pos)) >> val;
 
        if (zone == val) {
          return ellipsoid[i];
        }
      }
    }
 
    THROW(JParseError, "Invalid UTM grid " << grid);
  }
                                                  
 
  /**
   * Get geographical location.
   *
   * \param  grid               UTM grid
   * \param  position           UTM position
   * \return                    geographical location
   */
  inline JGeographicalLocation getGeographicalLocation(const JUTMGrid&     grid,
                                                       const JUTMPosition& position)
  {
    using namespace JPP;
 
    const Ellipsoid& el = getEllipsoid(grid);
 
    double lat, lon;
        
    UTMtoLL(el.id, position.getUTMNorth(),  position.getUTMEast(), grid.getUTMZone().c_str(), lat, lon);
 
    return JGeographicalLocation(getRadians(lat), getRadians(lon));
  }
 
 
  /**
   * Get Meridian convergence angle
   *
   * \param  grid               UTM grid
   * \param  position           UTM position
   * \return                    meridian convergence angle [rad]
   */
  inline double getMeridianConvergenceAngle(const JUTMGrid&     grid,
                                            const JUTMPosition& position)
  {
    const JGeographicalLocation location = getGeographicalLocation(grid, position);
 
    const Ellipsoid& el = getEllipsoid(grid);
 
    return getMeridianConvergenceAngle(location, { el.EquatorialRadius, el.eccentricitySquared });
  }
 
 
  const char* const calculated_t   =  "calculated";       //!< calulated values
  const char* const fixed_t        =  "fixed";            //!< fixed     values
}
 
 
/**
 * \file
 *
 * Example program to print geographical data for given detector.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  string  detectorFile;
  string  option;
  int     debug;
 
  try {
 
    JParser<> zap("Example program to print geographical data for given detector.");
 
    zap['a'] = make_field(detectorFile);
    zap['O'] = make_field(option)        = calculated_t, fixed_t;
    zap['d'] = make_field(debug)         = 2;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  JDetector detector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
 
  string                name;
  double                meridian;
  JGeographicalLocation location;
 
  if        (option == calculated_t) {
 
    meridian = getDegrees(getMeridianConvergenceAngle(detector.getUTMGrid(), detector.getUTMPosition()));
    location = getGeographicalLocation(detector.getUTMGrid(), detector.getUTMPosition());
 
  } else if (option == fixed_t) {
 
    if         (isORCADetector(detector)) {
 
      meridian   =  ORCA_MERIDIAN_CONVERGENCE_ANGLE_DEG;
      location   =  ORCA;
 
    } else if (isARCADetector(detector)) {
 
      meridian   =  ARCA_MERIDIAN_CONVERGENCE_ANGLE_DEG;
      location   =  ARCA;
 
    } else {
 
      FATAL("No location" << endl);
    }
 
  } else {
 
    FATAL("Invalid option " << option << endl);
  }
  
  if (debug >= debug_t) {
 
    cout << "grid: " << detector.getUTMGrid()     << endl;
    cout << "UTM:  " << detector.getUTMPosition() << endl;
    cout << "zone: " << detector.getUTMZone()     << ' ' << getUTMZone(location.getLongitude()) << endl;
 
    cout << "Meridian convergence angle [deg]: " << FIXED(9,5) << meridian                                          << " (" << option << ")" << endl;
    cout << "Meridian convergence angle [deg]: " << FIXED(9,5) << getDegrees(getMeridianConvergenceAngle(location)) << " (calculated)"       << endl;
  }
 
  cout << "Latitude  [deg]: " << FIXED(9,5) << getDegrees(location.getLatitude())    << endl;
  cout << "Longitude [deg]: " << FIXED(9,5) << getDegrees(location.getLongitude())   << endl;
  cout << "Latitude  [rad]: " << FIXED(9,5) << location.getLatitude()                << endl;
  cout << "Longitude [rad]: " << FIXED(9,5) << location.getLongitude()               << endl;
  cout << "Meridian convergence angle [deg]: " << FIXED(9,5) << meridian             << endl;
  cout << "Meridian convergence angle [rad]: " << FIXED(9,5) << getRadians(meridian) << endl;
}