LunarAzimuthalEqualArea.cpp

Go to the documentation of this file.

#include "LunarAzimuthalEqualArea.h"

#include <cmath>
#include <cfloat>

#include "Constants.h"
#include "IException.h"
#include "TProjection.h"
#include "Pvl.h"
#include "PvlGroup.h"

using namespace std;
namespace Isis {
  LunarAzimuthalEqualArea::LunarAzimuthalEqualArea(
    Pvl &label) : TProjection::TProjection(label) {
    try {
      // Try to read the mapping group
      PvlGroup &mapGroup = label.findGroup("Mapping", Pvl::Traverse);

      // Get the max libration
      m_maxLibration = mapGroup["MaximumLibration"];
      m_maxLibration *= PI / 180.0;
    }
    catch(IException &e) {
      QString message = "Invalid label group [Mapping]";
      throw IException(IException::Unknown, message, _FILEINFO_);
    }
  }

  LunarAzimuthalEqualArea::~LunarAzimuthalEqualArea() {
  }

   bool LunarAzimuthalEqualArea::operator== (const TProjection &proj) {
     if (!Projection::operator==(proj))
       return false;
     // dont use != (it is a recusive plunge)
     //  if (Projection::operator!=(proj)) return false;
   
     LunarAzimuthalEqualArea *LKAEA = (LunarAzimuthalEqualArea *) &proj;
     if (LKAEA->m_maxLibration != m_maxLibration)
       return false;
     return true;
   }

   QString LunarAzimuthalEqualArea::Name() const {
     return "LunarAzimuthalEqualArea";
   }

   QString LunarAzimuthalEqualArea::Version() const {
     return "0.1";
   }

  bool LunarAzimuthalEqualArea::SetGround(const double lat,
                                          const double lon) {
    // Convert longitude to radians
    m_longitude = lon;
    double lonRadians = lon * Isis::PI / 180.0;
    if (m_longitudeDirection == PositiveWest)
      lonRadians *= -1.0;

    // Now convert latitude to radians... it must be planetographic
    m_latitude = lat;
    double latRadians = lat;
    if (IsPlanetocentric())
      latRadians = ToPlanetographic(latRadians);
    latRadians *= Isis::PI / 180.0;

    double x, y;
    if (lonRadians == 0.0 && latRadians == 0.0) {
      x = 0.0;
      y = 0.0;
      SetComputedXY(x, y);
      m_good = true;
      return true;
    }

    double E = acos(cos(latRadians) * cos(lonRadians));
    double test = (sin(lonRadians) * cos(latRadians)) / sin(E);

    if (test > 1.0) test = 1.0;
    else if (test < -1.0) test = -1.0;

    double D = HALFPI - asin(test);
    if (latRadians < 0.0)
      D = -D;

    double radius = m_equatorialRadius;
    double PFAC = (HALFPI + m_maxLibration) / HALFPI;
    double RP = radius * sin(E / PFAC);

    x = RP * cos(D);
    y = RP * sin(D);

    SetComputedXY(x, y);
    m_good = true;
    return true;

  } // of SetGround


  bool LunarAzimuthalEqualArea::SetCoordinate(const double x,
      const double y) {
    // Save the coordinate
    SetXY(x, y);

    double RP = sqrt((x * x) + (y * y));

    double lat, lon;
    if (y == 0.0 && x == 0.0) {
      lat = 0.0;
      lon = 0.0;
      return true;
    }

    double radius = m_equatorialRadius;

    double D = atan2(y, x);
    double test = RP / radius;
    if (abs(test) > 1.0) {
      return false;
    }

    double EPSILON = 0.0000000001;
    double PFAC = (HALFPI + m_maxLibration) / HALFPI;
    double E = PFAC * asin(RP / radius);

    lat = HALFPI - (acos(sin(D) * sin(E)));

    if (abs(HALFPI - abs(lat)) <= EPSILON) {
      lon = 0.0;
    }
    else {
      test = sin(E) * cos(D) / sin(HALFPI - lat);
      if (test > 1.0) test = 1.0;
      else if (test < -1.0) test = -1.0;

      lon = asin(test);
    }

    if (E >= HALFPI) {
      if (lon <= 0.0) lon = -PI - lon;
      else lon = PI - lon;
    }

    // Convert to degrees
    m_latitude = lat * 180.0 / Isis::PI;
    m_longitude = lon * 180.0 / Isis::PI;

    // Cleanup the latitude
    if (IsPlanetocentric())
      m_latitude = ToPlanetocentric(m_latitude);

    m_good = true;
    return m_good;
  }


  bool LunarAzimuthalEqualArea::XYRange(double &minX, double &maxX,
                                        double &minY, double &maxY) {
    // Check the corners of the lat/lon range
    XYRangeCheck(m_minimumLatitude, m_minimumLongitude);
    XYRangeCheck(m_maximumLatitude, m_minimumLongitude);
    XYRangeCheck(m_minimumLatitude, m_maximumLongitude);
    XYRangeCheck(m_maximumLatitude, m_maximumLongitude);

    // If the latitude range contains 0
    if ((m_minimumLatitude < 0.0) && (m_maximumLatitude > 0.0)) {
      XYRangeCheck(0.0, m_minimumLongitude);
      XYRangeCheck(0.0, m_maximumLongitude);
    }

    // If the longitude range contains 0
    if ((m_minimumLongitude < 0.0) && (m_maximumLongitude > 0.0)) {
      XYRangeCheck(m_minimumLatitude, 0.0);
      XYRangeCheck(m_maximumLatitude, 0.0);
    }

    // Make sure everything is ordered
    if (m_minimumX >= m_maximumX) return false;
    if (m_minimumY >= m_maximumY) return false;

    // Return X/Y min/maxs
    minX = m_minimumX;
    maxX = m_maximumX;
    minY = m_minimumY;
    maxY = m_maximumY;
    return true;
  }

  PvlGroup LunarAzimuthalEqualArea::Mapping() {
    PvlGroup mapping = TProjection::Mapping();
    mapping += m_mappingGrp["MaximumLibration"];
    return mapping;
  }

} // end namespace isis

extern "C" Isis::TProjection *LunarAzimuthalEqualAreaPlugin(Isis::Pvl &lab,
    bool allowDefaults) {
  return new Isis::LunarAzimuthalEqualArea(lab);
}