TProjection.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#include "TProjection.h"
 
#include <QObject>
 
#include <cfloat>
#include <cmath>
#include <iomanip>
#include <sstream>
#include <vector>
 
#include <SpiceUsr.h>
 
#include "Constants.h"
#include "Displacement.h"
#include "FileName.h"
#include "IException.h"
#include "IString.h"
#include "Longitude.h"
#include "NaifStatus.h"
#include "Pvl.h"
#include "PvlGroup.h"
#include "PvlKeyword.h"
#include "SpecialPixel.h"
#include "Target.h"
#include "WorldMapper.h"
 
using namespace std;
namespace Isis {
  TProjection::TProjection(Pvl &label) : Projection::Projection(label) {
    try {
      // Mapping group is read by the parent
      // Get the radii from the EquatorialRadius and PolarRadius keywords
      if ((m_mappingGrp.hasKeyword("EquatorialRadius")) &&
          (m_mappingGrp.hasKeyword("PolarRadius"))) {
        m_equatorialRadius = m_mappingGrp["EquatorialRadius"];
        m_polarRadius = m_mappingGrp["PolarRadius"];
      }
      // Get the radii
      try {
         PvlGroup radii = Target::radiiGroup(label, m_mappingGrp);
         m_equatorialRadius = radii["EquatorialRadius"];
         m_polarRadius = radii["PolarRadius"];
      }
      catch (IException &e) {
        QString msg = "Projection failed. No target radii are available "
                      "through keywords [EquatorialRadius and PolarRadius] "
                      "or [TargetName].";
        throw IException(e, IException::Unknown, msg, _FILEINFO_);
      }
 
      // Check the radii for validity
      if (m_equatorialRadius <= 0.0) {
        QString msg = "Projection failed. Invalid value for keyword "
                      "[EquatorialRadius]. It must be greater than zero";
        throw IException(IException::Unknown, msg, _FILEINFO_);
      }
      if (m_polarRadius <= 0.0) {
        QString msg = "Projection failed. Invalid value for keyword "
                      "[PolarRadius]. It must be greater than zero";
        throw IException(IException::Unknown, msg, _FILEINFO_);
      }
 
      // Get the LatitudeType
      if ((QString) m_mappingGrp["LatitudeType"] == "Planetographic") {
        m_latitudeType = Planetographic;
      }
      else if ((QString) m_mappingGrp["LatitudeType"] == "Planetocentric") {
        m_latitudeType = Planetocentric;
      }
      else {
        QString msg = "Projection failed. Invalid value for keyword "
                      "[LatitudeType] must be "
                      "[Planetographic or Planetocentric]";
        throw IException(IException::Unknown, msg, _FILEINFO_);
      }
 
      // Get the LongitudeDirection
      if ((QString) m_mappingGrp["LongitudeDirection"] == "PositiveWest") {
        m_longitudeDirection = PositiveWest;
      }
      else if ((QString) m_mappingGrp["LongitudeDirection"] == "PositiveEast") {
        m_longitudeDirection = PositiveEast;
      }
      else {
        QString msg = "Projection failed. Invalid value for keyword "
                      "[LongitudeDirection] must be "
                      "[PositiveWest or PositiveEast]";
        throw IException(IException::Unknown, msg, _FILEINFO_);
      }
 
      // Get the LongitudeDomain
      if ((QString) m_mappingGrp["LongitudeDomain"] == "360") {
        m_longitudeDomain = 360;
      }
      else if ((QString) m_mappingGrp["LongitudeDomain"] == "180") {
        m_longitudeDomain = 180;
      }
      else {
        QString msg = "Projection failed. Invalid value for keyword "
                      "[LongitudeDomain] must be [180 or 360]";
        throw IException(IException::Unknown, msg, _FILEINFO_);
      }
 
      // Get the ground range if it exists
      m_groundRangeGood = false;
      if ((m_mappingGrp.hasKeyword("MinimumLatitude")) &&
          (m_mappingGrp.hasKeyword("MaximumLatitude")) &&
          (m_mappingGrp.hasKeyword("MinimumLongitude")) &&
          (m_mappingGrp.hasKeyword("MaximumLongitude"))) {
        m_minimumLatitude  = m_mappingGrp["MinimumLatitude"];
        m_maximumLatitude  = m_mappingGrp["MaximumLatitude"];
        m_minimumLongitude = m_mappingGrp["MinimumLongitude"];
        m_maximumLongitude = m_mappingGrp["MaximumLongitude"];
 
        if ((m_minimumLatitude < -90.0) || (m_minimumLatitude > 90.0)) {
          QString msg = "Projection failed. "
                        "[MinimumLatitude] of [" + toString(m_minimumLatitude)
                        + "] is outside the range of [-90:90]";
          throw IException(IException::Unknown, msg, _FILEINFO_);
        }
 
        if ((m_maximumLatitude < -90.0) || (m_maximumLatitude > 90.0)) {
          QString msg = "Projection failed. "
                        "[MaximumLatitude] of [" + toString(m_maximumLatitude)
                        + "] is outside the range of [-90:90]";
          throw IException(IException::Unknown, msg, _FILEINFO_);
        }
 
        if (m_minimumLatitude >= m_maximumLatitude) {
          QString msg = "Projection failed. "
                        "[MinimumLatitude,MaximumLatitude] of ["
                        + toString(m_minimumLatitude) + ","
                        + toString(m_maximumLatitude) + "] are not "
                        + "properly ordered";
          throw IException(IException::Unknown, msg, _FILEINFO_);
        }
 
        if (m_minimumLongitude >= m_maximumLongitude) {
          QString msg = "Projection failed. "
                        "[MinimumLongitude,MaximumLongitude] of ["
                        + toString(m_minimumLongitude) + "," 
                        + toString(m_maximumLongitude) + "] are not "
                        + "properly ordered";
          throw IException(IException::Unknown, msg, _FILEINFO_);
        }
 
        m_groundRangeGood = true;
      }
      else {
        // if no ground range is given, initialize the min/max lat/lon to 0
        m_minimumLatitude  = 0.0;
        m_maximumLatitude  = 0.0;
        m_minimumLongitude = 0.0;
        m_maximumLongitude = 0.0;
      }
 
      // Initialize miscellaneous protected data elements
      if (m_equatorialRadius < m_polarRadius) {
        QString msg = "Projection failed. Invalid keyword value(s). "
                      "[EquatorialRadius] = " + toString(m_equatorialRadius)
                      + " must be greater than or equal to [PolarRadius] = "
                      + toString(m_polarRadius);
        throw IException(IException::Unknown, msg, _FILEINFO_);
      }
      else {
        m_eccentricity = 1.0 -
                         (m_polarRadius * m_polarRadius) /
                         (m_equatorialRadius * m_equatorialRadius);
        m_eccentricity = sqrt(m_eccentricity);
      }
 
      // initialize the rest of the x,y,lat,lon member variables
      m_latitude = Null;
      m_longitude = Null;
 
      // If we made it to this point, we have what we need for a triaxial projection
      setProjectionType(Triaxial);
    }
    catch (IException &e) {
      QString msg = "Projection failed.  Invalid label group [Mapping]";
      throw IException(e, IException::Unknown, msg, _FILEINFO_);
    }
  }
 
  TProjection::~TProjection() {
  }
 
  bool TProjection::operator== (const Projection &proj) {
    if (!Projection::operator==(proj)) return false;
    TProjection *tproj = (TProjection *) &proj;
    if (EquatorialRadius() != tproj->EquatorialRadius()) return false;
    if (PolarRadius() != tproj->PolarRadius()) return false;
    if (IsPlanetocentric() != tproj->IsPlanetocentric()) return false;
    if (IsPositiveWest() != tproj->IsPositiveWest()) return false;
    return true;
  }
 
 
  double TProjection::EquatorialRadius() const {
    return m_equatorialRadius;
  }
 
  double TProjection::PolarRadius() const {
    return m_polarRadius;
  }
 
  double TProjection::Eccentricity() const {
    return m_eccentricity;
  }
 
  double TProjection::LocalRadius(double latitude) const {
    if (latitude == Null) {
      throw IException(IException::Unknown, 
                       "Unable to calculate local radius. The given latitude value [" 
                       + toString(latitude) + "] is invalid.", 
                       _FILEINFO_);
    }
    double a = m_equatorialRadius;
    double c = m_polarRadius;
    // to save calculations, if the target is spherical, return the eq. rad
    if (a - c < DBL_EPSILON) {
      return a;
    }
    else {
      double lat = latitude * PI / 180.0;
      return  a * c / sqrt(pow(c * cos(lat), 2) + pow(a * sin(lat), 2));
    }
  }
 
  double TProjection::LocalRadius() const {
    return LocalRadius(m_latitude);
  }
 
 
  double TProjection::TrueScaleLatitude() const {
    return 0.0;
  }
 
  bool TProjection::IsEquatorialCylindrical() {
    return false;
  }
 
  bool TProjection::IsPlanetocentric() const {
    return m_latitudeType == Planetocentric;
  }
 
  bool TProjection::IsPlanetographic() const {
    return m_latitudeType == Planetographic;
  }
 
  double TProjection::ToPlanetocentric(const double lat) const {
    return ToPlanetocentric(lat, m_equatorialRadius, m_polarRadius);
  }
 
  double TProjection::ToPlanetocentric(const double lat,
                                      double eRadius, double pRadius) {
    if (lat == Null || abs(lat) > 90.0) {
      throw IException(IException::Unknown, 
                       "Unable to convert to Planetocentric. The given latitude value [" 
                       + toString(lat) + "] is invalid.", 
                       _FILEINFO_);
    }
    double mylat = lat;
    if (abs(mylat) < 90.0) {  // So tan doesn't fail
      mylat *= PI / 180.0;
      mylat = atan(tan(mylat) * (pRadius / eRadius) *
                   (pRadius / eRadius));
      mylat *= 180.0 / PI;
    }
    return mylat;
  }
 
  double TProjection::ToPlanetographic(const double lat) const {
    return ToPlanetographic(lat,  m_equatorialRadius, m_polarRadius);
  }
 
  double TProjection::ToPlanetographic(double lat,
                                      double eRadius, double pRadius) {
    //Account for double rounding error.
    if (qFuzzyCompare(fabs(lat), 90.0)) {
      lat = qRound(lat);
    }
    if (lat == Null || fabs(lat) > 90.0) {
      throw IException(IException::Unknown, 
                       "Unable to convert to Planetographic. The given latitude value [" 
                       + toString(lat) + "] is invalid.", 
                       _FILEINFO_);
    }
    double mylat = lat;
    if (fabs(mylat) < 90.0) {  // So tan doesn't fail
      mylat *= PI / 180.0;
      mylat = atan(tan(mylat) * (eRadius / pRadius) *
                   (eRadius / pRadius));
      mylat *= 180.0 / PI;
    }
    return mylat;
  }
 
  QString TProjection::LatitudeTypeString() const {
    if (m_latitudeType == Planetographic) return "Planetographic";
    return "Planetocentric";
  }
 
  bool TProjection::IsPositiveEast() const {
    return m_longitudeDirection == PositiveEast;
  }
 
  bool TProjection::IsPositiveWest() const {
    return m_longitudeDirection == PositiveWest;
  }
 
  double TProjection::ToPositiveEast(const double lon, const int domain) {
    if (lon == Null) {
      throw IException(IException::Unknown, 
                       "Unable to convert to PositiveEast. The given longitude value [" 
                       + toString(lon) + "] is invalid.", 
                       _FILEINFO_);
    }
    double mylon = lon;
 
    mylon *= -1;
 
    if (domain == 360) {
      mylon = To360Domain(mylon);
    }
    else if (domain == 180) {
      mylon = To180Domain(mylon);
    }
    else {
      QString msg = "Unable to convert longitude.  Domain [" + toString(domain) 
                    + "] is not 180 or 360.";
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
 
    return mylon;
  }
 
  double TProjection::ToPositiveWest(const double lon, const int domain) {
    if (lon == Null) {
      throw IException(IException::Unknown, 
                       "Unable to convert to PositiveWest. The given longitude value [" 
                       + toString(lon) + "] is invalid.", 
                       _FILEINFO_);
    }
    double mylon = lon;
 
    mylon *= -1;
 
    if (domain == 360) {
      mylon = To360Domain(mylon);
    }
    else if (domain == 180) {
      mylon = To180Domain(mylon);
    }
    else {
      QString msg = "Unable to convert longitude.  Domain [" + toString(domain)
                    + "] is not 180 or 360.";
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
 
    return mylon;
  }
 
  QString TProjection::LongitudeDirectionString() const {
    if (m_longitudeDirection == PositiveEast) return "PositiveEast";
    return "PositiveWest";
  }
 
  bool TProjection::Has180Domain() const {
    return m_longitudeDomain == 180;
  }
 
  bool TProjection::Has360Domain() const {
    return m_longitudeDomain == 360;
  }
 
  double TProjection::To180Domain(const double lon) {
    if (lon == Null) {
      throw IException(IException::Unknown, 
                       "Unable to convert to 180 degree domain. The given longitude value [" 
                       + toString(lon) + "] is invalid.", 
                       _FILEINFO_);
    }
    return Isis::Longitude(lon, Angle::Degrees).force180Domain().degrees();
  }
 
  double TProjection::To360Domain(const double lon) {
    if (lon == Null) {
      throw IException(IException::Unknown, 
                       "Unable to convert to 360 degree domain. The given longitude value [" 
                       + toString(lon) + "] is invalid.", 
                       _FILEINFO_);
    }
    double result = lon;
 
    if ( (lon < 0.0 || lon > 360.0) &&
        !qFuzzyCompare(lon, 0.0) && !qFuzzyCompare(lon, 360.0)) {
     result = Isis::Longitude(lon, Angle::Degrees).force360Domain().degrees();
    }
 
    return result;
  }
 
  QString TProjection::LongitudeDomainString() const {
    if (m_longitudeDomain == 360) return "360";
    return "180";
  }
 
  double TProjection::MinimumLatitude() const {
    return m_minimumLatitude;
  }
 
  double TProjection::MaximumLatitude() const {
    return m_maximumLatitude;
  }
 
  double TProjection::MinimumLongitude() const {
    return m_minimumLongitude;
  }
 
  double TProjection::MaximumLongitude() const {
    return m_maximumLongitude;
  }
 
  bool TProjection::SetGround(const double lat, const double lon) {
    if (lat == Null || lon == Null) {
      m_good = false;
      return m_good;
    }
    else {
      m_latitude = lat;
      m_longitude = lon;
      m_good = true;
      SetComputedXY(lon, lat);
    }
    return m_good;
  }
 
  bool TProjection::SetCoordinate(const double x, const double y) {
    if (x == Null || y == Null) {
      m_good = false;
    }
    else {
      m_good = true;
      SetXY(x, y);
      m_latitude = XCoord();
      m_longitude = YCoord();
    }
    return m_good;
  }
 
 
  double TProjection::Latitude() const {
    return m_latitude;
  }
 
  double TProjection::Longitude() const {
    return m_longitude;
  }
 
 
  bool TProjection::SetUniversalGround(const double lat, const double lon) {
    if (lat == Null || lon == Null) {
      m_good = false;
      return m_good;
    }
    // Deal with the longitude first
    m_longitude = lon;
    if (m_longitudeDirection == PositiveWest) m_longitude = -lon;
    if (m_longitudeDomain == 180) {
      m_longitude = To180Domain(m_longitude);
    }
    else {
      // Do this because longitudeDirection could cause (-360,0)
      m_longitude = To360Domain(m_longitude);
    }
 
    // Deal with the latitude
    if (m_latitudeType == Planetographic) {
      m_latitude = ToPlanetographic(lat);
    }
    else {
      m_latitude = lat;
    }
 
    // Now the lat/lon are in user defined coordinates so set them
    return SetGround(m_latitude, m_longitude);
  }
 
 
  bool TProjection::SetUnboundUniversalGround(const double lat, const double lon) {
    if (lat == Null || lon == Null) {
      m_good = false;
      return m_good;
    }
    // Deal with the longitude first
    m_longitude = lon;
    if (m_longitudeDirection == PositiveWest) m_longitude = -lon;
 
    // Deal with the latitude
    if (m_latitudeType == Planetographic) {
      m_latitude = ToPlanetographic(lat);
    }
    else {
      m_latitude = lat;
    }
 
    // Now the lat/lon are in user defined coordinates so set them
    return SetGround(m_latitude, m_longitude);
  }
 
 
  double TProjection::UniversalLatitude() {
    double lat = m_latitude;
    if (m_latitudeType == Planetographic) lat = ToPlanetocentric(lat);
    return lat;
  }
 
  double TProjection::UniversalLongitude() {
    double lon = m_longitude;
    if (m_longitudeDirection == PositiveWest) lon = -lon;
    lon = To360Domain(lon);
    return lon;
  }
 
 
  double TProjection::Scale() const {
    if (m_mapper != NULL) {
      double lat = TrueScaleLatitude() * PI / 180.0;
      double a = m_polarRadius * cos(lat);
      double b = m_equatorialRadius * sin(lat);
      double localRadius = m_equatorialRadius * m_polarRadius /
                           sqrt(a * a + b * b);
 
      return localRadius / m_mapper->Resolution();
    }
    else {
      return 1.0;
    }
  }
 
 
  bool TProjection::XYRange(double &minX, double &maxX, 
                           double &minY, double &maxY) {
    if (minX == Null || maxX == Null || minY == Null || maxY == Null) {
      return false;
    }
    if (m_groundRangeGood) {
      minX = m_minimumLongitude;
      maxX = m_maximumLongitude;
      minY = m_minimumLatitude;
      maxY = m_maximumLatitude;
      return true;
    }
    return false;
  }
 
  void TProjection::XYRangeCheck(const double latitude, const double longitude) {
 
    if (latitude == Null || longitude == Null) {
      m_good = false;
      return;
    }
 
    SetGround(latitude, longitude);
    if (!IsGood()) return;
 
    if (XCoord() < m_minimumX) m_minimumX = XCoord();
    if (XCoord() > m_maximumX) m_maximumX = XCoord();
    if (YCoord() < m_minimumY) m_minimumY = YCoord();
    if (YCoord() > m_maximumY) m_maximumY = YCoord();
    return;
  }
 
 
  bool TProjection::inLongitudeRange(double minLon, 
                                     double maxLon, 
                                     double longitude) {
    // get the min/max range closest to 0.0 lon
    double adjustedLon =    To360Domain(longitude);
    double adjustedMinLon = To360Domain(minLon);
    double adjustedMaxLon = To360Domain(maxLon); 
 
    if (adjustedMinLon > adjustedMaxLon) {
      if (adjustedLon > adjustedMinLon) {
        adjustedLon -= 360;
      }
      adjustedMinLon -= 360;
    }
 
    // if this range covers all longitudes, then the given longitude is clearly in range
    if (qFuzzyCompare(maxLon - minLon, 360.0)) {
      return true;
    }
    else if (adjustedMinLon <= adjustedLon && adjustedLon <= adjustedMaxLon) { 
      return true;
    }
    else {
      return false;
    }
  }
 
 
  bool TProjection::inLongitudeRange(double longitude) {
    return inLongitudeRange(MinimumLongitude(), MaximumLongitude(), longitude);
  }
 
 
  bool TProjection::inLatitudeRange(double latitude) {
    if (MaximumLatitude() - MinimumLatitude() == 180) {
      return true;
    }
    else if (MinimumLatitude() <= latitude && latitude <= MaximumLatitude()) { 
      return true;
    }
    else {
      return false;
    }
  }
 
 
  bool TProjection::xyRangeOblique(double &minX, double &maxX, 
                                  double &minY, double &maxY) {
    if (minX == Null || maxX == Null || minY == Null || maxY == Null) {
      return false;
    }
    //For oblique, we'll have to walk all 4 sides to find out min/max x/y values
    if (!HasGroundRange()) return false; // Don't have min/max lat/lon, 
                                        //can't continue
 
    m_specialLatCases.clear();
    m_specialLonCases.clear();
 
    // First, search longitude for min X/Y
    double minFoundX1, minFoundX2;
    double minFoundY1, minFoundY2;
 
    // Search for minX between minlat and maxlat along minlon
    doSearch(MinimumLatitude(), MaximumLatitude(), 
             minFoundX1, MinimumLongitude(), true, true, true);
    // Search for minX between minlat and maxlat along maxlon
    doSearch(MinimumLatitude(), MaximumLatitude(), 
             minFoundX2, MaximumLongitude(), true, true, true);
    // Search for minY between minlat and maxlat along minlon
    doSearch(MinimumLatitude(), MaximumLatitude(), 
             minFoundY1, MinimumLongitude(), false, true, true);
    // Search for minY between minlat and maxlat along maxlon
    doSearch(MinimumLatitude(), MaximumLatitude(), 
             minFoundY2, MaximumLongitude(), false, true, true);
 
    // Second, search latitude for min X/Y
    double minFoundX3, minFoundX4;
    double minFoundY3, minFoundY4;
 
    // Search for minX between minlon and maxlon along minlat
    doSearch(MinimumLongitude(), MaximumLongitude(), 
             minFoundX3, MinimumLatitude(), true, false, true);
    // Search for minX between minlon and maxlon along maxlat
    doSearch(MinimumLongitude(), MaximumLongitude(), 
             minFoundX4, MaximumLatitude(), true, false, true);
    // Search for minY between minlon and maxlon along minlat
    doSearch(MinimumLongitude(), MaximumLongitude(), 
             minFoundY3, MinimumLatitude(), false, false, true);
    // Search for minY between minlon and maxlon along maxlat
    doSearch(MinimumLongitude(), MaximumLongitude(), 
             minFoundY4, MaximumLatitude(), false, false, true);
 
    // We've searched all possible minimums, go ahead and store the lowest
    double minFoundX5 = min(minFoundX1, minFoundX2);
    double minFoundX6 = min(minFoundX3, minFoundX4);
    m_minimumX = min(minFoundX5, minFoundX6);
 
    double minFoundY5 = min(minFoundY1, minFoundY2);
    double minFoundY6 = min(minFoundY3, minFoundY4);
    m_minimumY = min(minFoundY5, minFoundY6);
 
    // Search longitude for max X/Y
    double maxFoundX1, maxFoundX2;
    double maxFoundY1, maxFoundY2;
 
    // Search for maxX between minlat and maxlat along minlon
    doSearch(MinimumLatitude(), MaximumLatitude(), 
             maxFoundX1, MinimumLongitude(), true, true, false);
    // Search for maxX between minlat and maxlat along maxlon
    doSearch(MinimumLatitude(), MaximumLatitude(), 
             maxFoundX2, MaximumLongitude(), true, true, false);
    // Search for maxY between minlat and maxlat along minlon
    doSearch(MinimumLatitude(), MaximumLatitude(), 
             maxFoundY1, MinimumLongitude(), false, true, false);
    // Search for maxY between minlat and maxlat along maxlon
    doSearch(MinimumLatitude(), MaximumLatitude(), 
             maxFoundY2, MaximumLongitude(), false, true, false);
 
    // Search latitude for max X/Y
    double maxFoundX3, maxFoundX4;
    double maxFoundY3, maxFoundY4;
 
    // Search for maxX between minlon and maxlon along minlat
    doSearch(MinimumLongitude(), MaximumLongitude(), 
             maxFoundX3, MinimumLatitude(), true, false, false);
    // Search for maxX between minlon and maxlon along maxlat
    doSearch(MinimumLongitude(), MaximumLongitude(), 
             maxFoundX4, MaximumLatitude(), true, false, false);
    // Search for maxY between minlon and maxlon along minlat
    doSearch(MinimumLongitude(), MaximumLongitude(), 
             maxFoundY3, MinimumLatitude(), false, false, false);
    // Search for maxY between minlon and maxlon along maxlat
    doSearch(MinimumLongitude(), MaximumLongitude(), 
             maxFoundY4, MaximumLatitude(), false, false, false);
 
    // We've searched all possible maximums, go ahead and store the highest
    double maxFoundX5 = max(maxFoundX1, maxFoundX2);
    double maxFoundX6 = max(maxFoundX3, maxFoundX4);
    m_maximumX = max(maxFoundX5, maxFoundX6);
 
    double maxFoundY5 = max(maxFoundY1, maxFoundY2);
    double maxFoundY6 = max(maxFoundY3, maxFoundY4);
    m_maximumY = max(maxFoundY5, maxFoundY6);
 
    // Look along discontinuities for more extremes
    vector<double> specialLatCases = m_specialLatCases;
    for (unsigned int specialLatCase = 0; 
        specialLatCase < specialLatCases.size(); 
        specialLatCase ++) {
      double minX, maxX, minY, maxY;
 
      // Search for minX between minlon and maxlon along latitude discontinuities
      doSearch(MinimumLongitude(), MaximumLongitude(), 
               minX, specialLatCases[specialLatCase], true,  false, true);
      // Search for minY between minlon and maxlon along latitude discontinuities
      doSearch(MinimumLongitude(), MaximumLongitude(), 
               minY, specialLatCases[specialLatCase], false, false, true);
      // Search for maxX between minlon and maxlon along latitude discontinuities
      doSearch(MinimumLongitude(), MaximumLongitude(), 
               maxX, specialLatCases[specialLatCase], true,  false, false);
      // Search for maxX between minlon and maxlon along latitude discontinuities
      doSearch(MinimumLongitude(), MaximumLongitude(), 
               maxY, specialLatCases[specialLatCase], false, false, false);
 
      m_minimumX = min(minX, m_minimumX);
      m_maximumX = max(maxX, m_maximumX);
      m_minimumY = min(minY, m_minimumY);
      m_maximumY = max(maxY, m_maximumY);
    }
 
    vector<double> specialLonCases = m_specialLonCases;
    for (unsigned int specialLonCase = 0; 
        specialLonCase < specialLonCases.size(); 
        specialLonCase ++) {
      double minX, maxX, minY, maxY;
 
      // Search for minX between minlat and maxlat along longitude discontinuities
      doSearch(MinimumLatitude(), MaximumLatitude(), 
               minX, specialLonCases[specialLonCase], true,  true, true);
      // Search for minY between minlat and maxlat along longitude discontinuities
      doSearch(MinimumLatitude(), MaximumLatitude(), 
               minY, specialLonCases[specialLonCase], false, true, true);
      // Search for maxX between minlat and maxlat along longitude discontinuities
      doSearch(MinimumLatitude(), MaximumLatitude(), 
               maxX, specialLonCases[specialLonCase], true,  true, false);
      // Search for maxY between minlat and maxlat along longitude discontinuities
      doSearch(MinimumLatitude(), MaximumLatitude(), 
               maxY, specialLonCases[specialLonCase], false, true, false);
 
      m_minimumX = min(minX, m_minimumX);
      m_maximumX = max(maxX, m_maximumX);
      m_minimumY = min(minY, m_minimumY);
      m_maximumY = max(maxY, m_maximumY);
    }
 
    m_specialLatCases.clear();
    m_specialLonCases.clear();
 
    // 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;
  }
 
  void TProjection::doSearch(double minBorder, double maxBorder, 
                            double &extremeVal, const double constBorder, 
                            bool searchX, bool searchLongitude, bool findMin) {
    if (minBorder == Null || maxBorder == Null || constBorder == Null) {
      return;
    }
    const double TOLERANCE = PixelResolution()/2;
    const int NUM_ATTEMPTS = (unsigned int)DBL_DIG; // It's unsafe to go past 
                                                    // this precision
 
    double minBorderX, minBorderY, maxBorderX, maxBorderY;
    int attempts = 0;
 
    do {
      findExtreme(minBorder, maxBorder, minBorderX, minBorderY, maxBorderX, 
                  maxBorderY, constBorder, searchX, searchLongitude, findMin);
      if (minBorderX == Null && maxBorderX == Null 
          && minBorderY == Null && maxBorderY == Null ) {
        attempts = NUM_ATTEMPTS;
        continue;
      }
      attempts ++;
    }
    while ((fabs(minBorderX - maxBorderX) > TOLERANCE 
           || fabs(minBorderY - maxBorderY) > TOLERANCE)
           && (attempts < NUM_ATTEMPTS)); 
    // check both x and y distance in case symmetry of map
    // For example, if minBorderX = maxBorderX but minBorderY = -maxBorderY,
    // these points may not be close enough.
 
    if (attempts >= NUM_ATTEMPTS) {
      // We zoomed in on a discontinuity because our range never shrank, this
      // will need to be rechecked later. 
      // *min and max border should be nearly identical, so it doesn't matter
      //  which is used here
      if (searchLongitude) {
        m_specialLatCases.push_back(minBorder);
      }
      else {
        m_specialLonCases.push_back(minBorder);
      }
    }
 
    // These values will always be accurate, even over a discontinuity
    if (findMin) {
      if (searchX) extremeVal = min(minBorderX, maxBorderX);
      else         extremeVal = min(minBorderY, maxBorderY);
    }
    else {
      if (searchX) extremeVal = max(minBorderX, maxBorderX);
      else         extremeVal = max(minBorderY, maxBorderY);
    }
    return;
  }
 
  void TProjection::findExtreme(double &minBorder,  double &maxBorder,
                               double &minBorderX, double &minBorderY,
                               double &maxBorderX, double &maxBorderY, 
                               const double constBorder, bool searchX, 
                               bool searchLongitude, bool findMin) {
    if (minBorder == Null || maxBorder == Null || constBorder == Null) {
      minBorderX = Null;
      minBorderY = minBorderX;
      minBorderY = minBorderX;
      return;
    }
    if (!searchLongitude && (fabs(fabs(constBorder) - 90.0) < DBL_EPSILON)) {
      // it is impossible to search "along" a pole
      setSearchGround(minBorder, constBorder, searchLongitude);
      minBorderX = XCoord();
      minBorderY = YCoord();
      maxBorderY = minBorderY;
      return;
    }
    // Always do 10 steps
    const double STEP_SIZE = (maxBorder - minBorder) / 10.0;
    const double LOOP_END = maxBorder + (STEP_SIZE / 2.0); // This ensures we do 
                                                           // all of the steps
                                                           // properly
    double currBorderVal = minBorder;
    setSearchGround(minBorder, constBorder, searchLongitude);
 
    // this makes sure that the initial currBorderVal is valid before entering
    // the loop below
    if (!m_good){
      // minBorder = currBorderVal+STEP_SIZE < LOOP_END until setGround is good?
      // then, if still not good return?
      while (!m_good && currBorderVal <= LOOP_END) {
        currBorderVal+=STEP_SIZE;
        if (searchLongitude && (currBorderVal - 90.0 > DBL_EPSILON)) {
          currBorderVal = 90.0;
        }
        setSearchGround(currBorderVal, constBorder, searchLongitude);
      }
      if (!m_good) {
        minBorderX = Null;
        minBorderY = Null;
        return;
      }
    }
 
    // save the values of three consecutive steps from the minBorder towards
    // the maxBorder along the constBorder. initialize these three border
    // values (the non-constant lat or lon)
    double border1 = currBorderVal;
    double border2 = currBorderVal;
    double border3 = currBorderVal;
 
    // save the coordinate (x or y) values that correspond to the first
    // two borders that are being saved.
    // initialize these two coordinate values (x or y)
    double value1 = (searchX) ? XCoord() : YCoord();
    double value2 = value1;
 
    // initialize the extreme coordinate value 
    // -- this is the largest coordinate value found so far
    double extremeVal2 = value2;
 
    // initialize the extreme border values
    // -- these are the borders on either side of the extreme coordinate value
    double extremeBorder1 = minBorder;
    double extremeBorder3 = minBorder;
 
    while (currBorderVal <= LOOP_END) {
 
      // this conditional was added to prevent trying to SetGround with an
      // invalid latitude greater than 90 degrees. There is no need check for
      // latitude less than -90 since we start at the minBorder (already
      // assumed to be valid) and step forward toward (and possibly past)
      // maxBorder
      if (searchLongitude && (currBorderVal - 90.0 > DBL_EPSILON)) {
        currBorderVal = 90.0;
      }
 
      // update the current border value along constBorder
      currBorderVal += STEP_SIZE;
      setSearchGround(currBorderVal, constBorder, searchLongitude);
      if (!m_good){ 
        continue;
      } 
                     
      // update the border and coordinate values 
      border3 = border2;
      border2 = border1;
      border1 = currBorderVal;                                             
      value2 = value1;
      value1 = (searchX) ? XCoord() : YCoord();
 
      if ((findMin && value2 < extremeVal2) 
          || (!findMin && value2 > extremeVal2)) {
        // Compare the coordinate value associated with the center border with 
        // the current extreme. If the updated coordinate value is more extreme
        // (smaller or larger, depending on findMin), then we update the
        // extremeVal and it's borders.
        extremeVal2 = value2;
 
        extremeBorder3 = border3;
        extremeBorder1 = border1;
      }
    }
 
    // update min/max border values to the values on either side of the most 
    // extreme coordinate found in this call to this method
    
    minBorder = extremeBorder3; // Border 3 is lagging and thus smaller
 
    // since the loop steps past the original maxBorder, we want to retain 
    // the original maxBorder value so we don't go outside of the original
    // min/max range given
    if (extremeBorder1 <= maxBorder ) {
      maxBorder = extremeBorder1; // Border 1 is leading and thus larger
    }
 
    // update minBorder coordinate values
    setSearchGround(minBorder, constBorder, searchLongitude);
    // if (!m_good){
    //   this should not happen since minBorder has already been verified in 
    //   the while loop above
    // }
 
    minBorderX = XCoord();
    minBorderY = YCoord();
 
    // update maxBorder coordinate values
    setSearchGround(maxBorder, constBorder, searchLongitude);
    // if (!m_good){
    //   this should not happen since maxBorder has already been verified in
    //   the while loop above
    // }
 
    maxBorderX = XCoord();
    maxBorderY = YCoord();
    return;
  }
 
  void TProjection::setSearchGround(const double variableBorder, 
                                   const double constBorder, 
                                   bool variableIsLat) {
    if (variableBorder == Null || constBorder == Null) {
      return;
    }
    double lat, lon;
    if (variableIsLat) {
      lat = variableBorder;
      lon = constBorder;
    }
    else {
      lat = constBorder;
      lon = variableBorder;
    }
    SetGround(lat, lon);
    return;
  }
 
 
  PvlGroup TProjection::Mapping() {
    PvlGroup mapping("Mapping");
 
    QStringList keyNames;
    keyNames << "TargetName" << "ProjectionName" << "EquatorialRadius" << "PolarRadius"
             << "LatitudeType" << "LongitudeDirection" << "LongitudeDomain"
             << "PixelResolution" << "Scale" << "UpperLeftCornerX" << "UpperLeftCornerY"
             << "MinimumLatitude" << "MaximumLatitude" << "MinimumLongitude" << "MaximumLongitude"
             << "Rotation";
 
    foreach (QString keyName, keyNames) {
      if (m_mappingGrp.hasKeyword(keyName)) {
        mapping += m_mappingGrp[keyName];
      }
    }
 
    return mapping;
  }
 
 
  PvlGroup TProjection::MappingLatitudes() {
    PvlGroup mapping("Mapping");
 
    if (HasGroundRange()) {
      mapping += m_mappingGrp["MinimumLatitude"];
      mapping += m_mappingGrp["MaximumLatitude"];
    }
 
    return mapping;
  }
 
  PvlGroup TProjection::MappingLongitudes() {
    PvlGroup mapping("Mapping");
 
    if (HasGroundRange()) {
      mapping += m_mappingGrp["MinimumLongitude"];
      mapping += m_mappingGrp["MaximumLongitude"];
    }
 
    return mapping;
  }
 
  double TProjection::qCompute(const double sinPhi) const {
    if (m_eccentricity < DBL_EPSILON) {
      QString msg = "Snyder's q variable should only be computed for "
                    "ellipsoidal projections.";
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
    double eSinPhi = m_eccentricity * sinPhi;
    return (1 - m_eccentricity * m_eccentricity)
           * (sinPhi / (1 - eSinPhi * eSinPhi) 
               - 1 / (2 * m_eccentricity) * log( (1 - eSinPhi) / (1 + eSinPhi) ));
           // Note: We know that q is well defined since 
           //   0 < e < 1 and -1 <= sin(phi) <= 1 
           //   implies that -1 < e*sin(phi) < 1
           //   Thus, there are no 0 denominators and the log domain is 
           //   satisfied, (1-e*sin(phi))/(1+e*sin(phi)) > 0
  }
 
  double TProjection::phi2Compute(const double t) const {
    double localPhi = HALFPI - 2.0 * atan(t);
    double halfEcc = 0.5 * Eccentricity();
    double difference = DBL_MAX;
    int iteration = 0;
    // a failure in this loop means an exception will be thrown, which is
    //   an expensive operation. So letting let the loop iterate quite a bit
    //   is not a big deal. Also, the user will be unable to project at all
    //   if this fails so more reason to let this loop iterate: better to
    //   function slow than not at all.
    const int MAX_ITERATIONS = 45;
 
    while ((iteration < MAX_ITERATIONS) && (difference > 0.0000000001))  {
      double eccTimesSinphi = Eccentricity() * sin(localPhi);
      double newPhi = HALFPI -
                      2.0 * atan(t * pow((1.0 - eccTimesSinphi) /
                                         (1.0 + eccTimesSinphi), halfEcc));
      difference = fabs(newPhi - localPhi);
      localPhi = newPhi;
      iteration++;
    }
 
    if (iteration >= MAX_ITERATIONS) {
      QString msg = "Failed to converge in TProjection::phi2Compute()";
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
 
    return localPhi;
  }
 
  double TProjection::mCompute(const double sinphi, const double cosphi) const {
    double eccTimesSinphi = Eccentricity() * sinphi;
    double denominator = sqrt(1.0 - eccTimesSinphi * eccTimesSinphi);
    return cosphi / denominator;
  }
 
  double TProjection::tCompute(const double phi, const double sinphi) const {
    if ((HALFPI) - fabs(phi) < DBL_EPSILON) return 0.0;
 
    double eccTimesSinphi = Eccentricity() * sinphi;
    double denominator  = pow((1.0 - eccTimesSinphi) /
                              (1.0 + eccTimesSinphi),
                              0.5 * Eccentricity());
    return tan(0.5 * (HALFPI - phi)) / denominator;
  }
 
  double TProjection::e4Compute() const {
    double onePlusEcc = 1.0 + Eccentricity();
    double oneMinusEcc = 1.0 - Eccentricity();
 
    return sqrt(pow(onePlusEcc, onePlusEcc) *
                pow(oneMinusEcc, oneMinusEcc));
  }
 
} //end namespace isis