UpturnedEllipsoidTransverseAzimuthal.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#include "UpturnedEllipsoidTransverseAzimuthal.h"
 
#include <cmath>
#include <cfloat>
#include <iostream>
#include <iomanip>
 
#include <QDebug>
 
#include "Constants.h"
#include "IException.h"
#include "IString.h"
#include "TProjection.h"
#include "Pvl.h"
#include "PvlGroup.h"
#include "PvlKeyword.h"
#include "SpecialPixel.h"
 
using namespace std;
namespace Isis {
  UpturnedEllipsoidTransverseAzimuthal::UpturnedEllipsoidTransverseAzimuthal(Pvl &label,
                                               bool allowDefaults) :
      TProjection::TProjection(label) {
    try {
      // This algorithm can be found in the professional paper 
      // Cartographic Projections for Small Bodies of the Solar System
      // by Maria E Fleis, Kira B Shingareva,
      // Michael M Borisov, Michael V Alexandrovich, and Philip Stooke
      
      // Try to read the mapping group
      PvlGroup &mapGroup = label.findGroup("Mapping", Pvl::Traverse);
 
      // Compute and write the default center longitude if allowed and
      // necessary
      double centerLongitude = 0.0;
      if (!mapGroup.hasKeyword("CenterLongitude")) {
        if (allowDefaults) {
          // centerLongitude still 0.0 here
          mapGroup += PvlKeyword("CenterLongitude", toString(centerLongitude), "Degrees");
        }
        else {
          QString message = "Cannot project using upturned ellipsoid Transverse Azimuthal";
          message += " without [CenterLongitude] value.  Keyword does not exist";
          message += " in labels and defaults are not allowed.";
          throw IException(IException::Unknown, message, _FILEINFO_);
        }
      }
      else {
        centerLongitude = mapGroup["CenterLongitude"];
      }
 
      if (MinimumLongitude() < centerLongitude - 90.0) {
        QString message = "MinimumLongitude ["  
                          + toString(MinimumLongitude()) 
                          + "] is invalid. Must be within -90 degrees of the CenterLongitude ["
                          + toString(centerLongitude) + "].";
        throw IException(IException::Unknown, message, _FILEINFO_);
      }
      if (MaximumLongitude() > centerLongitude + 90.0) {
        QString message = "MaximumLongitude ["  
                          + toString(MaximumLongitude()) 
                          + "] is invalid. Must be within +90 degrees of the CenterLongitude ["
                          + toString(centerLongitude) + "].";
        throw IException(IException::Unknown, message, _FILEINFO_);
      }
 
      // Get the center longitude  & latitude
      // initialize member variables
      init(centerLongitude);
    }
    catch(IException &e) {
      QString message = "Invalid label group [Mapping]";
      throw IException(e, IException::Unknown, message, _FILEINFO_);
    }
  }
  
 
  UpturnedEllipsoidTransverseAzimuthal::~UpturnedEllipsoidTransverseAzimuthal() {
  }
 
 
  bool UpturnedEllipsoidTransverseAzimuthal::operator== (const Projection &proj) {
    // don't do the below it is a recursive plunge
    //  if (Projection::operator!=(proj)) return false;
 
    // all other data members calculated using m_lambda0, m_a and m_b.
    UpturnedEllipsoidTransverseAzimuthal *tcyl = (UpturnedEllipsoidTransverseAzimuthal *) &proj;
    if ((tcyl->m_lambda0 != m_lambda0) ||
        (tcyl->m_a != m_a) ||
        (tcyl->m_b != m_b)) return false;
 
    return true;
  }
 
 
  QString UpturnedEllipsoidTransverseAzimuthal::Name() const {
    return "UpturnedEllipsoidUpturnedEllipsoidTransverseAzimuthal";
  }
 
 
  QString  UpturnedEllipsoidTransverseAzimuthal::Version() const {
    return "1.0";
  }
 
 
  void UpturnedEllipsoidTransverseAzimuthal::init(double centerLongitude) {
 
    // adjust for positive east longitude direction and convert to radians 
    if (IsPositiveEast()) {
      m_lambda0 = centerLongitude * DEG2RAD;
    }
    else {
      m_lambda0 = ToPositiveEast(centerLongitude, 360) * DEG2RAD;
    }
    if (Has180Domain()) {
      m_lambda0 = To180Domain(m_lambda0);
    }
    else {
      m_lambda0 = To360Domain(m_lambda0);
    }
 
    // Initialize miscellaneous protected data elements
    m_good = false;
 
    m_minimumX = DBL_MAX;
    m_maximumX = -DBL_MAX;
    m_minimumY = DBL_MAX;
    m_maximumY = -DBL_MAX;
 
    // Descriptions of the variables a, e
    double axis1 = EquatorialRadius();
    double axis2 = PolarRadius();
    m_a = qMax(axis1, axis2); // major axis. generally, this is the larger of the equatorial axes
    m_b = qMin(axis1, axis2); // minor axis. generally, this is the polar axis
    m_e = Eccentricity();  // e = sqrt(1 - (b/a)^2) 
                           // must be 0 <= e < 1
    if (qFuzzyCompare(0.0, m_e)) {
      m_e = 0.0;
    }
 
    // to reduce calculations, calculate some scalar constants:
    double t0 = m_b / m_a; // equals sqrt( 1 - m_e * m_e );
    m_t  = t0 * t0; // equals 1 - m_e * m_e ;
    m_t1 = m_e * t0; // equals e/sqrt( 1 - m_e * m_e );
    double k1 = 2 * m_a * exp(m_t1 * atan(m_t1));
    // k = radius of the Equator of transverse graticule on Azimuthal projection
    //     under condition of no distortion in the center of projection
    m_k = k1 * m_t;
  }
 
 
  bool UpturnedEllipsoidTransverseAzimuthal::SetGround(const double lat, 
                                                       const double lon) {
    // when lat/lon is Null we can't set ground
    if (lat == Null || lon == Null) {
      m_good = false;
      return m_good;
    }
 
    // convert given lat to radians, planetocentric
    // phiNorm = The planetocentric latitude, in normal aspect.
    double phiNorm;
    // when lat just barely beyond pole, set equal to pole
    if (qFuzzyCompare(90.0, qAbs(lat)) && qAbs(lat) > 90.0) {
      phiNorm = copysign(HALFPI, lat); // returns sign(lat) * PI/2
    }
    // if well beyond pole, it's bad
    else if (qAbs(lat) > 90.0) {
      m_good = false;
      return m_good;
    }
    else if (IsPlanetocentric()) {
      phiNorm = lat * DEG2RAD;
    }
    else {
      // equations expect ocentric latitudes.
      phiNorm = ToPlanetocentric(lat) * DEG2RAD;
    }
 
    // convert given lon to positive east, then subtract center lon to get 
    // lambdaNorm = longitude east of the center of projection 
    double lambdaNorm;
    if (IsPositiveEast()) {
      lambdaNorm = lon * DEG2RAD - m_lambda0;
    }
    else {
      lambdaNorm = ToPositiveEast(lon, 360) * DEG2RAD - m_lambda0;
    }
 
    // Compute the rectangular x,y coordinate using the appropriate set of equations below
    double x,y;
 
    // NOTE: z = the angular distance from the center of projection.
    double cosz = cos(phiNorm) * cos(lambdaNorm);
 
    // First, we take care of the edge cases where
    // 1. a rounding error causes cosz to be outside of the range of cosine
    // 2. z == 0 (this could be handled by the next set of equations,
    //            but taking care of it here reduces calculations.)
    if (cosz >= 1) {
      // This is the origin,
      // at lat = equator, lon = centerlon 
      x = 0.0;
      y = 0.0;
    }
 
    // This condition is used for the set of equations that are
    // written to exclude the singularities where z == 0 or z == PI
    // (i.e. the given longitude is equal to the center longitude
    // or 180 degrees from it).
    // Use these equations for 0.5 < cosz < 1.0, i.e. 0 < z < PI/3
    else if (cosz > 0.5) { //  && cosz < 1 is implied by bounds on cosine
      // use pythagorean identity to get sine
      double sinz = sqrt( 1 - cosz*cosz ); // due to restrictions above sinz != 0
      // phi = the latitude at "upturned" ellipsoid of revolution
      // NOTE: Since cos(z) > 0.5 or cos(z) < -0.5,
      //       there is no risk of zero denominator when dividing by cosz
      double phi = HALFPI - atan2( sinz, m_t * cosz );
      double sinPhi = sin(phi);
      // rhoOverTanZ = rho/tanz
      //       where rho = the radius of latitude circle (transverse graticule)
      // NOTE: We know sin(phi) = -1 only if
      //       phi = PI/2 - arctan(angle) = -PI/2 or 3PI/2
      //       but the range of arctangent is (-PI/2, PI/2),
      //       so there is no risk of zero denominator dividing by (1+sin(phi))
      double rhoOverTanZ = m_k * sinPhi / ( ( 1 + sinPhi )
                                            * m_t
                                            * exp( m_t1 * atan( m_t1 * sinPhi ) ) );
      x = rhoOverTanZ * tan(lambdaNorm);
      y = rhoOverTanZ * (sin(phiNorm) / cosz ); // Note: cos(z) > 0.5, so no dividing by zero
    }
 
    // This condition is used for the set of equations that are
    // written to exclude the singularity where z == PI/2.
    // Use these equations for -1 <= cosz <= 0.5, i.e. PI/3 <= z < PI.
    else { 
      // For this case, we will restrict z near multiples of PI using a
      // tolerance value in order to avoid singularities at 0 and PI. We
      // define zmin and zmax by:
      // 
      // zmin = minimal value of angular distance from the center of projection
      //      = 0 + angular tolerance
      // zmax = maximal value of angular distance from the center of projection
      //        in the neighborhood of point opposite to the center of projection
      //      = PI - angular tolerance
      double tolerance = 0.0016;//???
      double coszmax = cos(PI - tolerance);
      double lambdaModulus = fmod(lambdaNorm, TWOPI);
      if (cosz < coszmax) {
        cosz = coszmax; // prevents cosz from equalling -1
        // The following prevent lambdaNorm from being too close to +/-PI
        
        // if the given longitude is within tolerance of -PI and less than -PI,
        // set it to -PI - tolerance)
        //if (qFuzzyCompare(lambdaNorm, -PI) && lambdaNorm <= -PI) {
        if (-PI - tolerance < lambdaModulus && lambdaModulus <= -PI) {
          lambdaNorm = -PI - tolerance;
        }
        // if the given longitude is within tolerance of -PI and greater than -PI, 
        // set it to -zmax (i.e. -PI + tolerance)
        // 
        //if (qFuzzyCompare(lambdaNorm, -PI) && lambdaNorm > -PI) {
        else if (-PI < lambdaModulus && lambdaModulus <= -PI + tolerance) {
          lambdaNorm = -PI + tolerance;
        }
        // if the given longitude is within tolerance of PI and less than PI,
        // set it to zmax (i.e. PI - tolerance)
        // 
        //if (qFuzzyCompare(lambdaNorm, PI) && lambdaNorm <= PI) {
        else if (PI - tolerance < lambdaModulus && lambdaModulus <= PI) {
          lambdaNorm = PI - tolerance;
        }
        // if the given longitude is within tolerance of PI and greater than PI,
        // set it to PI + tolerance
        // 
        //if (qFuzzyCompare(lambdaNorm, PI) && lambdaNorm > PI) {
        else if (PI < lambdaModulus && lambdaModulus < PI + tolerance) {
          lambdaNorm = PI + tolerance;
        }
      }
      // use pythagorean identity to get sine
      double sinz = sqrt( 1 - cosz*cosz ); // due to restrictions above, we know sinz != 0
 
      // phi = latitude at "upturned" ellipsoid of revolution
      // NOTE: On the interval, PI/3 <= z < PI,
      //       we know 0 < sin(z) <= 1,
      //       so there is no risk of zero denom when dividing by sinz
      double phi = atan2( m_t * cosz, sinz );
      double sinPhi = sin(phi);
      // Note: We know sin(phi) = -1 only if
      //       phi = arctan(angle) = -PI/2
      //       but the range of arctangent is the open interval (-PI/2, PI/2),
      //       so there is no risk of zero denominator dividing by (1+sin(phi))
      double rhoOverSinZ = m_k * cos(phi) / ( ( 1 + sinPhi )
                                              * sinz
                                              * exp( m_t1 * atan(m_t1 * sinPhi ) ) );
      x = rhoOverSinZ * cos(phiNorm) * sin(lambdaNorm);
      y = rhoOverSinZ * sin(phiNorm);
    }
 
    // set member data and return status
    m_good = true;
    m_latitude = lat;
    m_longitude = lon; 
    SetComputedXY(x, y); // sets m_x = x and m_y = y and handles rotation
                         // Note: if x or y are Null, sets m_good back to false
    return m_good;
  }
 
 
  bool UpturnedEllipsoidTransverseAzimuthal::SetCoordinate(const double x, 
                                                           const double y) {
    if (x == Null || y == Null) {
      m_good = false;
      return m_good;
    }
    // Save the coordinate
    SetXY(x, y);
 
    if (qFuzzyCompare(x + 1.0,  1.0) && qFuzzyCompare(y + 1.0,  1.0)) {
      // origin
      m_latitude = 0.0;
      m_longitude = m_lambda0 * RAD2DEG;
    }
    else {
 
      // We are given the following forward equations:
      // 
      //     x = (rho/sin(z)) * cos(phiNorm) * sin(lambdaNorm)
      //     y = (rho/sin(z)) * sin(phiNorm)
      //     rho/sin(z) = ( k*cos(phi) / [(1+sin(phi))*sin(z)*e^{t1*arctan(t1*sin(phi)}] )
      //     phi = arctan( (1-e^2)*cos(z) / sin(z) )
      // 
      // So, after simplifying, we can verify x^2 + y^2 = rho^2
      // Thus, we have two solutions for rho:
      // 
      //     rho = +/- sqrt(x^2 + y^2)
      //     rho = ( k*cos(phi) / [(1+sin(phi))*e^{t1*arctan(t1*sin(phi)}] )
      // 
      // Now, we can use Newton's root finding method for the function
      //     f(phi) = g(phi) - h(phi) on [-pi/2, pi/2]
      // where
      //     g(phi) = k*cos(phi) / [(1+sin(phi))*e^{t1*arctan(t1*sin(phi)}]
      // and
      //     h(phi) = sqrt(x^2 + y^2).
      // 
      // After simplifying, the derivative with respect to phi is
      //     f'(phi) = -k * (1 + t1^2) /
      //               [(1 + sin(phi)) * e^{t1 * arctan(t1 * sin(phi)} * (1 + t1^2 * sin^2(phi))]
 
      double phi0, fphi0, fprimephi0, phi1;
      bool converged = false;
      int iterations = 0;
      double tolerance = 10e-10;
      phi0 = 0.0; // start at the equator???
      while (!converged && iterations < 1000) {
        fphi0 = m_k * cos(phi0) / ((1 + sin(phi0)) * exp(m_t1 * atan( m_t1 * sin(phi0) )))
                - sqrt(x*x + y*y);
 
        fprimephi0 = -m_k * (1 + m_t1 * m_t1)
                     / ((1 + sin(phi0)) 
                        * exp(m_t1 * atan( m_t1 * sin(phi0) )) 
                        * (1 + m_t1 * m_t1 * sin(phi0) * sin(phi0)));
        phi1 = phi0 - fphi0 / fprimephi0;
        // if phi wrapped at the poles, make sure phi is on [-pi/2, pi/2]
        if (qAbs(phi1) > HALFPI) {
          double phiDegrees = To180Domain(phi1*RAD2DEG);
          if (phiDegrees > 90.0) {
            phiDegrees -= 90.0;
          }
          if (phiDegrees < -180.0) {
            phiDegrees += 90.0;
          }
          phi1 = phiDegrees * DEG2RAD;
        }
        if (qAbs(phi0 - phi1) < tolerance) {
          converged = true;
        }
        else {
          phi0 = phi1;
          iterations++;
        }
      }
 
      if (!converged) {
        m_good = false;
        return m_good;
      }
 
      // Now we have phi, the latitude at "upturned" ellipsoid of revolution.
      double phi = phi1;
 
      // Now use the forward equation for phi to solve for z, the angular distance from the center
      // of projection:
      // 
      //     phi = arctan( (1-e^2) * cos(z) / sin(z))
      // 
      double z = atan2((1 - m_e * m_e), tan(phi)); // see handling of cylind???
 
      // Get phiNorm,the planetocentric latitude, in normal aspect. The range of arcsine is
      // [-pi/2, pi/2] so we are guaranteed to get an angle between the poles. Use the forward
      // equation:
      // 
      //     y = (rho / sinz ) * sin(phiNorm)
      //double rho = ( m_k * cos(phi) / ((1+sin(phi))*exp(m_e, m_t1*arctan(m_t1*sin(phi))) );
      double rho = sqrt(x*x + y*y);
      double phiNorm = asin( y * sin(z) / rho );
 
      // Get lambdaNorm, the longitude east of lambda0, using the forward equation
      //     cos(z) = cos(phiNorm) * cos(lambdaNorm)
      double lambdaNorm;
 
      // make sure we are in the correct quadrant...
      double cosLambdaNorm = cos(z) / cos(phiNorm);    
      if (cosLambdaNorm > 1.0) {
        lambdaNorm = 0.0;
      }
      else if (cosLambdaNorm < -1.0) {
        lambdaNorm = PI; // +/- PI doesn't matter here... same answer either way
      }
      else if (x >= 0 && y >= 0) {
        lambdaNorm = acos(cosLambdaNorm);
      }
      else if (x < 0 && y >= 0) {// fix conditional???
        lambdaNorm = -acos(cosLambdaNorm);
      }
      else if (y < 0 && x >= 0) {// fix conditional???
        lambdaNorm = acos(cosLambdaNorm);
      }
      else { // y < 0 && x < 0
        lambdaNorm = -acos(cosLambdaNorm);
      }
 
      // calculations give positive east longitude
      m_longitude = (lambdaNorm + m_lambda0) * RAD2DEG;
      m_latitude = phiNorm * RAD2DEG;
    }
 
    // Cleanup the latitude
    if (IsPlanetocentric()) {
      m_latitude = ToPlanetocentric(m_latitude);
    }
 
    // Cleanup the longitude
    if (IsPositiveWest()) {
      m_longitude = ToPositiveWest(m_longitude, m_longitudeDomain);
    }
 
    if (Has180Domain()) {
      m_longitude = To180Domain(m_longitude);
    }
    else {
      // Do this because longitudeDirection could cause (-360,0)
      m_longitude = To360Domain(m_longitude);
    }
 
    m_good = true;
    return m_good;
  }
 
 
  bool UpturnedEllipsoidTransverseAzimuthal::XYRange(double &minX, double &maxX, 
                                      double &minY, double &maxY) {
 
    // First check combinations of the lat/lon range boundaries
    XYRangeCheck(m_minimumLatitude, m_minimumLongitude);
    XYRangeCheck(m_maximumLatitude, m_minimumLongitude);
    XYRangeCheck(m_minimumLatitude, m_maximumLongitude);
    XYRangeCheck(m_maximumLatitude, m_maximumLongitude);
 
    double centerLongitude = m_lambda0 * RAD2DEG;
 
    bool centerLongitudeInRange =    TProjection::inLongitudeRange(centerLongitude);
    bool centerLongitude90InRange =  TProjection::inLongitudeRange(centerLongitude + 90.0);
    bool centerLongitude180InRange = TProjection::inLongitudeRange(centerLongitude + 180.0);
    bool centerLongitude270InRange = TProjection::inLongitudeRange(centerLongitude + 270.0);
 
    if (centerLongitudeInRange) {
      XYRangeCheck(m_minimumLatitude, centerLongitude);
      XYRangeCheck(m_maximumLatitude, centerLongitude);
    }
    if (centerLongitude90InRange) {
      XYRangeCheck(m_minimumLatitude, centerLongitude + 90.0);
      XYRangeCheck(m_maximumLatitude, centerLongitude + 90.0);
    }
    if (centerLongitude180InRange) {
      XYRangeCheck(m_minimumLatitude, centerLongitude + 180.0);
      XYRangeCheck(m_maximumLatitude, centerLongitude + 180.0);
    }
    if (centerLongitude270InRange) {
      XYRangeCheck(m_minimumLatitude, centerLongitude + 270.0);
      XYRangeCheck(m_maximumLatitude, centerLongitude + 270.0);
    }
 
    if (TProjection::inLatitudeRange(0.0)) {
      XYRangeCheck(0.0, m_minimumLongitude);
      XYRangeCheck(0.0, m_maximumLongitude);
      if (centerLongitudeInRange) {
        XYRangeCheck(0.0, centerLongitude);
      }
      if (centerLongitude90InRange) {
        XYRangeCheck(0.0, centerLongitude + 90.0);
      }
      if (centerLongitude180InRange) {
        XYRangeCheck(0.0, centerLongitude + 180.0);
      }
      if (centerLongitude270InRange) {
        XYRangeCheck(0.0, centerLongitude + 270.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 UpturnedEllipsoidTransverseAzimuthal::Mapping() {
    PvlGroup mapping = TProjection::Mapping();
    mapping += PvlKeyword("CenterLongitude", toString(m_lambda0 * RAD2DEG));
    return mapping;
  }
 
 
  PvlGroup UpturnedEllipsoidTransverseAzimuthal::MappingLatitudes() {
    return TProjection::MappingLatitudes();
  }
 
 
  PvlGroup UpturnedEllipsoidTransverseAzimuthal::MappingLongitudes() {
    PvlGroup mapping = TProjection::MappingLongitudes();
    mapping += PvlKeyword("CenterLongitude", toString(m_lambda0 * RAD2DEG));
    return mapping;
  }
 
} // end namespace isis
 
extern "C" Isis::Projection *UpturnedEllipsoidTransverseAzimuthalPlugin(Isis::Pvl &lab,
                                                                          bool allowDefaults) {
  return new Isis::UpturnedEllipsoidTransverseAzimuthal(lab, allowDefaults);
}