Latitude.cpp

Go to the documentation of this file.

#include "Latitude.h"

#include <cmath>

#include "Constants.h"
#include "Distance.h"
#include "IException.h"
#include "IString.h"
#include "PvlGroup.h"
#include "SpecialPixel.h"
#include "TProjection.h"

namespace Isis {
  Latitude::Latitude() : Angle() {
    m_equatorialRadius = NULL;
    m_polarRadius = NULL;

    m_errors = AllowPastPole;
  }


  Latitude::Latitude(double latitude, Angle::Units latitudeUnits,
                     ErrorChecking errors) : Angle() {
    m_equatorialRadius = NULL;
    m_polarRadius = NULL;

    m_errors = errors;

    setPlanetocentric(latitude, latitudeUnits);
  }


  Latitude::Latitude(Angle latitude, ErrorChecking errors) : Angle() {
    m_equatorialRadius = NULL;
    m_polarRadius = NULL;

    m_errors = errors;

    setPlanetocentric(latitude.radians(), Radians);
  }


  Latitude::Latitude(Angle latitude, PvlGroup mapping,
            ErrorChecking errors) : Angle(latitude) {
    m_equatorialRadius = NULL;
    m_polarRadius = NULL;

    if (mapping.hasKeyword("EquatorialRadius") && mapping.hasKeyword("PolarRadius")) {
      m_equatorialRadius = new Distance(toDouble(mapping["EquatorialRadius"][0]),
          Distance::Meters);
      m_polarRadius = new Distance(toDouble(mapping["PolarRadius"][0]),
          Distance::Meters);
    }
    else {
      PvlGroup radiiGrp = TProjection::TargetRadii(mapping["TargetName"]);

      m_equatorialRadius = new Distance(toDouble(radiiGrp["EquatorialRadius"][0]),
          Distance::Meters);
      m_polarRadius = new Distance(toDouble(radiiGrp["PolarRadius"][0]),
          Distance::Meters);
    }

    m_errors = errors;

    if(mapping["LatitudeType"][0] == "Planetographic") {
      setPlanetographic(latitude.radians(), Radians);
    }
    else if(mapping["LatitudeType"][0] == "Planetocentric") {
      setPlanetocentric(latitude.radians(), Radians);
    }
    else {
      QString msg = "Latitude type [" + mapping["LatitudeType"][0] +
        "] is not recognized";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
  }


  Latitude::Latitude(double latitude,
            PvlGroup mapping,
            Angle::Units latitudeUnits,
            ErrorChecking errors) : Angle(latitude, latitudeUnits) {
    m_equatorialRadius = NULL;
    m_polarRadius = NULL;

    if (mapping.hasKeyword("EquatorialRadius") && mapping.hasKeyword("PolarRadius")) {
      m_equatorialRadius = new Distance(toDouble(mapping["EquatorialRadius"][0]),
          Distance::Meters);
      m_polarRadius = new Distance(toDouble(mapping["PolarRadius"][0]),
          Distance::Meters);
    }
    else {
      PvlGroup radiiGrp = TProjection::TargetRadii(mapping["TargetName"]);

      m_equatorialRadius = new Distance(toDouble(radiiGrp["EquatorialRadius"][0]),
          Distance::Meters);
      m_polarRadius = new Distance(toDouble(radiiGrp["PolarRadius"][0]),
          Distance::Meters);
    }

    m_errors = errors;

    if(mapping["LatitudeType"][0] == "Planetographic") {
      setPlanetographic(latitude, latitudeUnits);
    }
    else if(mapping["LatitudeType"][0] == "Planetocentric") {
      setPlanetocentric(latitude, latitudeUnits);
    }
    else {
      IString msg = "Latitude type [" + IString(mapping["LatitudeType"][0]) +
        "] is not recognized";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
  }


  Latitude::Latitude(double latitude,
                     Distance equatorialRadius, Distance polarRadius,
                     CoordinateType latType,
                     Angle::Units latitudeUnits,
                     ErrorChecking errors) : Angle(latitude, latitudeUnits) {
    m_equatorialRadius = NULL;
    m_polarRadius = NULL;

    m_equatorialRadius = new Distance(equatorialRadius);
    m_polarRadius = new Distance(polarRadius);

    m_errors = errors;

    if (latType == Planetocentric) {
      setPlanetocentric(latitude, latitudeUnits);
    }
    else if (latType == Planetographic) {
      setPlanetographic(latitude, latitudeUnits);
    }
    else {
      IString msg = "Enumeration value [" + IString(latType) + "] is not a "
        "valid CoordinateType";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
  }


  Latitude::Latitude(const Latitude &latitudeToCopy) : Angle(latitudeToCopy) {
    m_equatorialRadius = NULL;
    m_polarRadius = NULL;

    m_errors = latitudeToCopy.m_errors;

    if(latitudeToCopy.m_equatorialRadius) {
      m_equatorialRadius = new Distance(*latitudeToCopy.m_equatorialRadius);
    }

    if(latitudeToCopy.m_polarRadius) {
      m_polarRadius = new Distance(*latitudeToCopy.m_polarRadius);
    }
  }


  Latitude::~Latitude() {
    if (m_equatorialRadius) {
      delete m_equatorialRadius;
      m_equatorialRadius = NULL;
    }

    if (m_polarRadius) {
      delete m_polarRadius;
      m_polarRadius = NULL;
    }
  }


  double Latitude::planetocentric(Angle::Units units) const {
    return angle(units);
  }


  void Latitude::setPlanetocentric(double latitude, Angle::Units units) {
    setAngle(latitude, units);
  }


  double Latitude::planetographic(Angle::Units units) const {
    if (m_equatorialRadius == NULL || m_polarRadius == NULL) {
      IString msg = "Latitude [" + IString(degrees()) + " degrees] cannot "
          "be converted to Planetographic without the planetary radii, please "
          "use the other Latitude constructor";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    if (*this > Angle(90.0, Angle::Degrees) ||
        *this < Angle(-90.0, Angle::Degrees)) {
      IString msg = "Latitudes outside of the -90/90 range cannot be converted "
          "between Planetographic and Planetocentric";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    if(!isValid()) {
      IString msg = "Invalid planetographic latitudes are not currently "
          "supported";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    double ographicLatitude = atan(tan(radians()) *
        (*m_equatorialRadius / *m_polarRadius) *
        (*m_equatorialRadius / *m_polarRadius));

    // This theoretically should just be an angle, but make it a Latitude so
    //   we can access angle
    return Latitude(ographicLatitude, Angle::Radians).angle(units);
  }


  void Latitude::setPlanetographic(double latitude, Angle::Units units) {
    if (m_equatorialRadius == NULL || m_polarRadius == NULL) {
      IString msg = "Latitude [" + IString(latitude) + "] cannot be "
          "converted to Planetocentic without the planetary radii, please use "
          "the other Latitude constructor";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    Angle inputAngle(latitude, units);

    if (inputAngle > Angle(90.0, Angle::Degrees) ||
        inputAngle < Angle(-90.0, Angle::Degrees)) {
      IString msg = "Latitudes outside of the -90/90 range cannot be converted "
          "between Planetographic and Planetocentric";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    if(IsSpecial(latitude)) {
      IString msg = "Invalid planetographic latitudes are not currently "
          "supported";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    double ocentricLatitude = atan(tan(inputAngle.radians()) *
        (*m_polarRadius / *m_equatorialRadius) *
        (*m_polarRadius / *m_equatorialRadius));

    // Sometimes the trig functions return the negative of the expected value at the pole.
    if ((ocentricLatitude > 0) != (inputAngle.radians() > 0)) {
      ocentricLatitude *= -1;
    }

    setAngle(ocentricLatitude, Angle::Radians);
  }


  bool Latitude::inRange(Latitude min, Latitude max) const {
    // Validity check on the range
    if (min > max) {
      IString msg = "Minimum latitude [" + IString(min.degrees()) +
        "] degrees is greater than maximum latitude [" +
        IString(max.degrees()) + "] degrees";
      throw IException(IException::User, msg, _FILEINFO_);
    }

    // Provide a little wriggle room for precision problems
    Angle epsilon(DBL_EPSILON, Angle::Degrees);
    Latitude adjustedMin = min - epsilon;
    Latitude adjustedMax = max + epsilon;

    // Is this latitude between the min and the max
    return *this >= adjustedMin && *this <= adjustedMax;
  }


  Latitude& Latitude::operator=(const Latitude & latitudeToCopy) {
    if(this == &latitudeToCopy) return *this;

    m_equatorialRadius = NULL;
    m_polarRadius = NULL;

    m_errors = latitudeToCopy.m_errors;

    if(latitudeToCopy.m_equatorialRadius) {
      m_equatorialRadius = new Distance(*latitudeToCopy.m_equatorialRadius);
    }

    if(latitudeToCopy.m_polarRadius) {
      m_polarRadius = new Distance(*latitudeToCopy.m_polarRadius);
    }

    setPlanetocentric(latitudeToCopy.planetocentric());

    return *this;
  }
  

  Latitude Latitude::add(Angle angleToAdd, PvlGroup mapping) {

    CoordinateType latType;

    Distance equatorialRadius;
    Distance polarRadius;
    if (mapping.hasKeyword("EquatorialRadius") && mapping.hasKeyword("PolarRadius")) {
      equatorialRadius = Distance(toDouble(mapping["EquatorialRadius"][0]),
          Distance::Meters);
      polarRadius = Distance(toDouble(mapping["PolarRadius"][0]),
          Distance::Meters);
    }
    else {
      PvlGroup radiiGrp = TProjection::TargetRadii(mapping["TargetName"]);

      equatorialRadius = Distance(toDouble(radiiGrp["EquatorialRadius"][0]),
          Distance::Meters);
      polarRadius = Distance(toDouble(radiiGrp["PolarRadius"][0]),
          Distance::Meters);
    }

    if(mapping["LatitudeType"][0] == "Planetocentric")
      latType = Planetocentric;
    else if (mapping["LatitudeType"][0] == "Planetographic")
      latType = Planetographic;
    else {
      IString msg = "Latitude type [" + IString(mapping["LatitudeType"][0]) +
        "] is not recognized";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    return add(angleToAdd, equatorialRadius, polarRadius, latType);
  }

  
  Latitude Latitude::add(Angle angleToAdd, Distance equatorialRadius, Distance polarRadius,
                         CoordinateType latType) {
    Latitude result;
    
    switch (latType) {
      case Planetocentric:
        result = Latitude(planetocentric() + angleToAdd.radians(), equatorialRadius, polarRadius,
                          latType, Radians, m_errors);
        break;

      case Planetographic:
        result = Latitude(planetographic() + angleToAdd.radians(), equatorialRadius, polarRadius,
                          latType, Radians, m_errors);
        break;
    }
    
    return result;
  }

 
  void Latitude::setAngle(double angle, const Angle::Units &units) {
    // Check for passing 90 degrees if that error checking is on
    if (!IsSpecial(angle) && (m_errors & AllowPastPole) != AllowPastPole) {
      Angle tmpAngle(angle, units);
      if(tmpAngle > Angle(90, Angle::Degrees) ||
         tmpAngle < Angle(-90, Angle::Degrees)) {
        IString msg = "Latitudes past 90 degrees are not valid. The latitude "
            "[" + IString(tmpAngle.degrees()) + " degrees] is not allowed";
        throw IException(IException::Programmer, msg, _FILEINFO_);
      }
    }

    Angle::setAngle(angle, units);
  }
}