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 });
  }
}
 
 
/**
 * \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;
  int     debug;
 
  try {
 
    JParser<> zap("Example program to print geographical data for given detector.");
 
    zap['a'] = make_field(detectorFile);
    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         (isORCADetector(detector)) {
 
    name       = "ORCA";
    meridian   =  ORCA_MERIDIAN_CONVERGENCE_ANGLE_DEG;
    location   =  ORCA;
 
  } else if (isARCADetector(detector)) {
 
    name       = "ARCA";
    meridian   =  ARCA_MERIDIAN_CONVERGENCE_ANGLE_DEG;
    location   =  ARCA;
 
  } else {
 
    FATAL("No location" << endl);
  }
 
  cout << "Detector " << name << endl;
  cout << detector.getUTMGrid()     << endl;
  cout << detector.getUTMPosition() << endl;
 
  JGeographicalLocation geo = getGeographicalLocation(detector.getUTMGrid(), detector.getUTMPosition());
 
  cout << "Latitude  [deg]: " << FIXED(9,5) << getDegrees(location.getLatitude())  << " " << FIXED(9,5) << getDegrees(geo.getLatitude())  << endl;
  cout << "Longitude [deg]: " << FIXED(9,5) << getDegrees(location.getLongitude()) << " " << FIXED(9,5) << getDegrees(geo.getLongitude()) << endl;
  cout << "Latitude  [rad]: " << FIXED(9,5) << location.getLatitude()              << " " << FIXED(9,5) << geo.getLatitude()              << endl;
  cout << "Longitude [rad]: " << FIXED(9,5) << location.getLongitude()             << " " << FIXED(9,5) << geo.getLongitude()             << endl;
 
  cout << "Zone: " <<  detector.getUTMZone() << " " << getUTMZone(geo.getLongitude()) << endl;
  
  cout << "Meridian convergence angle [deg]: " << FIXED(9,5) << meridian << " " << FIXED(9,5) << getDegrees(getMeridianConvergenceAngle(detector.getUTMGrid(), detector.getUTMPosition())) << endl;
  cout << "Meridian convergence angle [deg]: " << FIXED(9,5) << getDegrees(getMeridianConvergenceAngle(location)) << " (approximate)" << endl;  
}