LambertAzimuthalEqualArea.cpp

Go to the documentation of this file.

#include "LambertAzimuthalEqualArea.h"

#include <cmath>
#include <cfloat>

#include <iostream>
#include <iomanip>

#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 {
  LambertAzimuthalEqualArea::LambertAzimuthalEqualArea(
      Pvl &label, 
      bool allowDefaults) :
      TProjection::TProjection(label) {
    try {
      // This algorithm can be found in the USGS Professional Paper 1395
      // Map Projections--A Working Manual by John P. Snyder
      
      // Try to read the mapping group
      PvlGroup &mapGroup = label.findGroup("Mapping", Pvl::Traverse);

      // Compute and write the default center longitude if allowed and
      // necessary
      if (!mapGroup.hasKeyword("CenterLongitude")) {
        if (allowDefaults) {
          double centerLon = (MinimumLongitude() + MaximumLongitude()) / 2.0;
          mapGroup += PvlKeyword("CenterLongitude", toString(centerLon), "Degrees");
        }
        else {
          QString message = "Cannot project using Lambert Azimuthal equal-area";
          message += " without [CenterLongitude] value.  Keyword does not exist";
          message += " in labels and defaults are not allowed.";
          throw IException(IException::Unknown, message, _FILEINFO_);
        }
      }

      // Compute and write the default center latitude if allowed and
      // necessary
      if (!mapGroup.hasKeyword("CenterLatitude")) {
        if (allowDefaults) {
          double centerLat = (MinimumLatitude() + MaximumLatitude()) / 2.0;
          mapGroup += PvlKeyword("CenterLatitude", toString(centerLat), "Degrees");
        }
        else {
          QString message = "Cannot project using Lambert Azimuthal equal-area";
          message += " without [CenterLatitude] value.  Keyword does not exist";
          message += " in labels and defaults are not allowed.";
          throw IException(IException::Unknown, message, _FILEINFO_);
        }
      }

      // Get the center longitude  & latitude
      double centerLongitude = mapGroup["CenterLongitude"];
      double centerLatitude = mapGroup["CenterLatitude"];

      if (fabs(MinimumLongitude() - centerLongitude) > 360.0 ||
          fabs(MaximumLongitude() - centerLongitude) > 360.0) {
        IString message = "[MinimumLongitude,MaximumLongitude] range of [";
        message += IString(MinimumLongitude())+","+IString(MaximumLongitude());
        message += "] is invalid.  No more than 360 degrees from the "
                   "CenterLongitude [" + IString(centerLongitude) + "] is allowed.";
        throw IException(IException::Unknown, message, _FILEINFO_);
      }

      if (MaximumLongitude() - MinimumLongitude() > 360.0) {
        IString message = "[MinimumLongitude,MaximumLongitude] range of ["  
                        + IString(MinimumLongitude()) + "," 
                        + IString(MaximumLongitude()) + "] is invalid. "
                        "No more than 360 degree range width is allowed.";
        throw IException(IException::Unknown, message, _FILEINFO_);
      }

      // initialize member variables
      init(centerLatitude, centerLongitude);
    }
    catch(IException &e) {
      IString message = "Invalid label group [Mapping]";
      throw IException(e, IException::Unknown, message, _FILEINFO_);
    }
  }
  
  LambertAzimuthalEqualArea::~LambertAzimuthalEqualArea() {
  }

  bool LambertAzimuthalEqualArea::operator== (const Projection &proj) {
    if (!Projection::operator==(proj)) return false;
    // don't do the below it is a recursive plunge
    //  if (Projection::operator!=(proj)) return false;
    LambertAzimuthalEqualArea *lama = (LambertAzimuthalEqualArea *) &proj;
    if ((lama->m_phi1 != m_phi1) ||
        (lama->m_lambda0 != m_lambda0) ||
        (lama->m_a != m_a) ||
        (lama->m_e != m_e)) return false;
    // all other data members calculated using m_lambda0, m_phi1, m_a and m_e.
    return true;
  }

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

  QString  LambertAzimuthalEqualArea::Version() const {
    return "1.0";
  }

  double LambertAzimuthalEqualArea::TrueScaleLatitude() const {
    //Snyder pg. 
    // no distortion at center of projection (centerLatitude, centerLongitude)
    return m_phi1 * 180.0 / PI;
  }

  
  void LambertAzimuthalEqualArea::init(double centerLatitude, 
                                       double centerLongitude) {
    // Initialize miscellaneous protected data elements
    m_good = false;

    m_minimumX = DBL_MAX;
    m_maximumX = -DBL_MAX;
    m_minimumY = DBL_MAX;
    m_maximumY = -DBL_MAX;

    // Test to make sure center longitude is valid
    if (fabs(centerLongitude) > 360.0) {
      IString message = "CenterLongitude [" + IString(centerLongitude);
      message += "] is outside the range of [-360, 360]";
      throw IException(IException::Unknown, message, _FILEINFO_);
    } 
 
    // Test to make sure center lat is valid
    if (fabs(centerLatitude) > 90.0) {
      IString message = "CenterLatitude [" + IString(centerLatitude);
      message += "] is outside the range of [-90, 90]";
      throw IException(IException::Unknown, message, _FILEINFO_);
    }

    // We need to convert to planetographic since phi is geographic latitude 
    // (see phi and phi1 - Snyder pg ix)
    if (IsPlanetocentric()) {
      centerLatitude = ToPlanetographic(centerLatitude);
    }

    // adjust for positive east longitude direction 
    // (see lambda and lambda0 - Snyder pg ix)
    if (m_longitudeDirection == PositiveWest) {
      centerLongitude = -1.0*centerLongitude;
    }

    // Descriptions of the variables a, e, lambda0, and phi1 
    // can be found in Snyder text on pgs viii-ix, pg 187
    m_a = EquatorialRadius();
    m_e = Eccentricity();  // e = sqrt(1 - (PolarRadius/EqRadius)^2) 
                           // must be 0 <= e < 1

    // Snyder variables used for spherical projection
    m_lambda0 = centerLongitude*PI / 180.0; // convert to radians
    m_phi1 = centerLatitude*PI / 180.0;     // convert to radians
    m_sinPhi1 = sin(m_phi1);
    m_cosPhi1 = cos(m_phi1);

    // flags for determinining which formulas to use
    m_spherical = false;
    m_northPolarAspect = false;
    m_southPolarAspect = false;
    m_equatorialAspect = false;

    if (qFuzzyCompare(0.0, m_e)) {
      m_e = 0.0;
      m_spherical = true;
    }
    if (qFuzzyCompare(HALFPI, m_phi1)) {
      m_phi1 = HALFPI;
      m_northPolarAspect = true;
    }
    if (qFuzzyCompare(-HALFPI, m_phi1)) {
      m_phi1 = -HALFPI;
      m_southPolarAspect = true;
    }
    if (qFuzzyCompare(0.0, m_phi1)) {
      m_phi1 = 0.0;
      m_equatorialAspect = true;
    }

    // Snyder ellipsoid variables will not be used
    m_qp = Null;
    m_q1 = Null;
    m_m1 = Null;
    m_beta1 = Null;
    m_sinBeta1 = Null;
    m_cosBeta1 = Null;
    m_Rq = Null;
    m_D = Null;

    // other Snyder variables
    m_h = Null;
    m_k = Null;

    // if eccentricity = 0, we are projecting a sphere.
    if (!m_spherical) {
      // Calculate Snyder ellipsoid variables
      initEllipsoid();
    }

    // Check if the antipodal point is in the lat/lon ranges.
    // 
    // We can only allow this if we have a polar projection.
    // Otherwise, we cannot SetGround() for the antipodal point since it is
    // projected as a circle, not a single point with (x,y) value.
    // 
    // The antipodal point is defined by the coordinates
    //     (-centerLat, centerLon-180) or (-centerLat, centerLon+180)
    if (!m_northPolarAspect && !m_southPolarAspect) {
      if (-centerLatitude >= MinimumLatitude() 
          && -centerLatitude <= MaximumLatitude()
          && (  (MinimumLongitude() <= centerLongitude - 180 &&
                 MaximumLongitude() >= centerLongitude - 180)
             || (MinimumLongitude() <= centerLongitude + 180 &&
                 MaximumLongitude() >= centerLongitude + 180)  ) ) {
        IString message = "Invalid latitude and/or longitude range. ";
        message += "Non-polar Lambert Azimuthal equal-area "
                   "projections can not project the antipodal point on "
                   "the body.";
        throw IException(IException::Unknown, message, _FILEINFO_);
      }
    }
  }

  void LambertAzimuthalEqualArea::initEllipsoid() {
    m_spherical = false;

    // Decided not to use qCompute method here to reduce chance of roundoff 
    // error since q for polar simplifies 
    // m_qp = qCompute(1.0); //q for polar, sin(phi) = sin(pi/2) = 1.0
    m_qp = 1 - (1 - m_e * m_e) / (2 * m_e) * log( (1 - m_e) / (1 + m_e) );
           // Snyder eqn (3-12) pg 187, use phi = pi/2 (polar)
           // Note: We know that qp is well defined since 0 < e < 1, thus
           //       1+e != 0 (no 0 denominator) and (1-e)/(1+e) > 0 (log domain)
    if (!m_northPolarAspect && !m_southPolarAspect) {
      // we only need to compute these variables for oblique projections

      m_q1 = qCompute(m_sinPhi1); // Snyder eqn (3-12) pg 187, use phi = phi1
      m_m1 = mCompute(m_sinPhi1, m_cosPhi1); // Snyder eqn (14-15) pg 187, 
             //use phi = phi1, thus m1 = cosPhi1/sqrt(1-eSinPhi1*eSinPhi1)
      m_beta1 = asin(m_q1 / m_qp); // Snyder eqn (3-11) pg 187, use q = q1
            // Note: We can show mathematically that beta1 is well defined
            //       First, since 0 < eccentricity < 1, we have (1-e)/(1+e) < 1
            //       Thus, ln((1-e)/(1+e)) < 0 and so qp > 1 (no 0 denominator)
            //       Next, verify the domain for arcsine, -1 <= q1/qp <= 1 
            //       Since not polar, |phi1| < 90, and since ellipsoid 0 < e < 1
            //       Then |sinPhi1| < 1, so e|sinPhi1| < e
            //       So, 1 < (1+e|sinPhi1|)/(1-e|sinPhi1|) < (1+e)/(1-e).
            //       0 < ln[(1+e|sinPhi1|)/(1-e|sinPhi1|)] < ln[(1+e)/(1-e)]
            //       Note that (1-e^2)/(2e) > 0 and recall -ln(x) = ln(1/x) 
            // So |q1| = |sinPhi1 - (1-e^2)/(2e)ln[(1-eSinPhi1)/(1+eSinPhi1)]|
            //         = |sinPhi1 + (1-e^2)/(2e)ln[(1+eSinPhi1)/(1-eSinPhi1)]|
            //        <= |sinPhi1| + |(1-e^2)/(2e)ln[(1+eSinPhi1)/(1-eSinPhi1)]|
            //         = |sinPhi1| +(1-e^2)/(2e) |ln[(1+eSinPhi1)/(1-eSinPhi1)]|
            //         <   1   + (1-e^2)/(2e) |ln[(1+eSinPhi1)/(1-eSinPhi1)]|
            //        <=   1   + (1-e^2)/(2e) ln[(1+e|SinPhi1|)/(1-e|SinPhi1|)]
            //         <   1   + (1-e^2)/(2e) ln[(1+e)/(1-e)]
            //         = 1 - (1-e^2)/(2e)ln[(1-e)/(1+e)]
            //         = qp
      m_sinBeta1 = sin(m_beta1);
      m_cosBeta1 = cos(m_beta1);
      m_Rq = m_a * sqrt(m_qp / 2); // Snyder eqn (3-13) pg 16, 187
            // Note: We know Rq is well defined since qp > 1 (see m_beta1)
      m_D  = m_a * m_m1 / (m_Rq * m_cosBeta1); // Snyder eqn (24-20) pg 187
            // Note: We know D is well-defined since
            //       cosBeta1 = 0 implies beta1 = +/-pi/2, 
            //       which in turn implies q1 = +/-qp
            //       This should only happen if polar aspect, 
            //       Furthermore, Rq = 0 is impossible since a > 0 and q > 1
    }
  }

  bool LambertAzimuthalEqualArea::SetGround(const double lat, 
                                            const double lon) {
    // when lat/lon is Null or lat is "beyond" the poles, we can't set ground
    if (fabs(lat) - 90.0 > DBL_EPSILON || lat == Null || lon == Null) {
      if (!qFuzzyCompare(90.0, fabs(lat))) {
        m_good = false;
        return m_good;
      }
    }
 
    m_longitude = lon; 
    m_latitude = lat;

    // assign input values to Snyder's phi and lambda variables 
    // convert to radians, positive east longitude, planetographic
    double phi     = lat * PI / 180.0;
    double lambda  = lon * PI / 180.0;

    // when lat just barely beyond pole, set equal to pole
    if (lat > 90.0 && qFuzzyCompare(90.0, lat)) {
      phi = HALFPI;
    }
    if (lat < -90.0 && qFuzzyCompare(-90.0, lat)) {
      phi = -HALFPI;
    }

    if (m_longitudeDirection == PositiveWest) lambda *= -1.0;
    if (IsPlanetocentric()) phi = ToPlanetographic(phi);

    if (!m_spherical) return setGroundEllipsoid(phi, lambda);

    // calculate the following once to reduce computations 
    double sinPhi   = sin(phi);
    double cosPhi   = cos(phi);
    double sinLambdaDiff = sin(lambda - m_lambda0);
    double cosLambdaDiff = cos(lambda - m_lambda0);

    // Use the variables defined above to compute the x,y coordinate
    double x,y;

    // In this case, Spherical Radius = Equitorial Radius, 
    // so Snyder's R = Snyder's a variable pg ix
    double R = m_a;
    if (m_northPolarAspect ) {
      double sinQuarterPiMinusHalfPhi = sin(PI/4-phi/2); 
      x = 2*R*sinQuarterPiMinusHalfPhi*sinLambdaDiff; //Snyder eqn (24-3) pg 186
      y = -2*R*sinQuarterPiMinusHalfPhi*cosLambdaDiff;//Snyder eqn (24-4) pg 186
      m_h = cos(PI/4-phi/2); // Snyder eqn (24-6) pg 186
      m_k = 1/m_h; // Snyder eqn (24-5) and (24-6) pg 186 
                   // - recall sec(theta) = 1/cos(theta)
    }
    else if (m_southPolarAspect ) {
      double cosQuarterPiMinusHalfPhi = cos(PI/4-phi/2);
      x = 2*R*cosQuarterPiMinusHalfPhi*sinLambdaDiff; //Snyder eqn (24-8) pg 186
      y = 2*R*cosQuarterPiMinusHalfPhi*cosLambdaDiff; //Snyder eqn (24-9) pg 186
      m_h = sin(PI/4-phi/2); // Snyder eqn (24-11) pg 186
      m_k = 1/m_h; // Snyder eqn (24-10) and (24-11) pg 186
    }
    else { // spherical oblique aspect (this includes equatorial)

      // Check if (phi, lambda) is the antipodal point (point diametrically
      // opposite the center point of projection)
      if (qFuzzyCompare(-m_phi1, phi) // phi = -phi1
          && fabs(fmod(lambda-m_lambda0, PI)) < DBL_EPSILON //lambda-lambda0=k*PI
          && fabs(fmod(lambda-m_lambda0, 2*PI)) > DBL_EPSILON) {  // k is odd
        // We don't allow this for oblique aspect projections.  This point,
        // at the opposite side of the planet, is projected as a circle of
        // radius 2R - we can't return a unique (x,y) value.
        m_good = false;
        return m_good;
      }

      // First, be sure that the denominator will not be zero
      double trigTerms;
      if (m_equatorialAspect) {
        // If m_phi1 == 0, the general case below cancels alegbraically
        // to the following equations, however, the simplified equations
        // are used here to reduce chance of roundoff errors
        // -- Snyder eq (24-14) pg 186
        trigTerms = cosPhi*cosLambdaDiff;
      }
      else {
        // general case for oblique projections -- Snyder eq (24-2) pg 185
        double sinProduct = m_sinPhi1 * sinPhi;
        double cosProduct = m_cosPhi1 * cosPhi * cosLambdaDiff;
        trigTerms = sinProduct + cosProduct;
      }
      // make sure when we add 1 to trigTerms that we are not near zero
      if (qFuzzyCompare(-1.0, trigTerms)) {
          m_good = false;
          return m_good;
      }
      double denominator = 1 + trigTerms;
      double kprime = sqrt(2/denominator); // Snyder eqn (24-2) or (24-14) pg 185-186
      x = R*kprime*cosPhi*sinLambdaDiff; // Snyder eqn (22-4) pg 185
      if (m_equatorialAspect) {
        // If m_phi1 == 0, the general case below cancels alegbraically
        // to the following, however, use the simplified equations for y
        // to reduce the likelihood of roundoff errors
        // -- Snyder eq (24-13) pg 186
        y = R*kprime*sinPhi;
      }
      else { // spherical general oblique aspect -- Snyder eqn (22-5) pg 185
        y = R*kprime*(m_cosPhi1*sinPhi - m_sinPhi1*cosPhi*cosLambdaDiff); 
      }
      // these are the relative scale factors for oblique aspect
      m_k = kprime;
      m_h = 1/m_k;
    }

    SetComputedXY(x, y); // sets m_x = x and m_y = y and handles rotation
    m_good = true;
    return m_good;
  }

  bool LambertAzimuthalEqualArea::setGroundEllipsoid(double phi, 
                                                     double lambda) {
    // calculate the following once to reduce computations 
    double sinPhi   = sin(phi);
    double cosPhi   = cos(phi);
    double sinLambdaDiff = sin(lambda - m_lambda0);
    double cosLambdaDiff = cos(lambda - m_lambda0);

    // Use the variables passed in to compute the x,y coordinate
    double x,y;

    // projecting an ellipsoid
    double q  = qCompute(sinPhi); //q for given lat
    if (qFuzzyCompare(PI/2, phi)) {
      q = m_qp; // reduce roundoff error
    }
    if (qFuzzyCompare(-PI/2, phi)) {
      q = -m_qp; // reduce roundoff error
    }
    // q = (1-e^2)*(sinphi/(1-e^2sinphi^2))-1/(2e)*log((1-esinphi)/(1+esinphi)))

    double m  = mCompute(sinPhi, cosPhi); // m = cosPhi/sqrt(1-eSinPhi*eSinPhi)
    if (m_northPolarAspect) {
      double rho = m_a* sqrt(m_qp - q); // Snyder eqn (24-23) pg 188
          // Note: Rho is well defined (qp >= q) as shown in initEllipsoid()
      x = rho * sinLambdaDiff; // Snyder eqn (21-30) pg 188
      y = -1.0 * rho * cosLambdaDiff; // Snyder eqn (21-31) pg 188
      if (qFuzzyCompare(m_phi1, phi)) {
        m_k = 1; // Since True scale at clat/clon - see Snyder page 190, table 29
      }
      else{
        m_k = rho / (m_a * m); // Snyder eqn (21-32) pg 188
      }
      m_h = 1 / m_k; // Snyder paragraph after eqn (24-23) pg 188

    }
    else if (m_southPolarAspect) {
      double rho = m_a* sqrt(m_qp + q); // Snyder eqn (24-25) pg 188
          // Note: Rho is well defined (qp >= -q) as shown in initEllipsoid()
      x = rho * sinLambdaDiff; // Snyder eqn (21-30) pg 188
      y = rho * cosLambdaDiff; // Snyder eqn (24-24) pg 188
      if (qFuzzyCompare(m_phi1, phi)) {
        m_k = 1; // Since True scale at clat/clon - see Snyder page 190, table 29
      }
      else{
        m_k = rho / (m_a * m); // Snyder eqn (21-32) pg 188
      }
      m_h = 1 / m_k;
    }
    else { // ellipsoidal oblique aspect
      // following comment is impossible -- see initEllipsoid(),  |q1| > |qp|
      // if (fabs(q) > fabs(m_qp)) {
      //  m_good = false;
      //  return m_good;
      // }
      double beta  = asin( q/m_qp); // Snyder eqn (3-11) pg 187
          // Note: beta is well defined- see initEllipsoid() m_beta1

      // calculate the following once to reduce computations when defining 
      // Snyder's variable B
      double sinBeta  = sin(beta);
      double cosBeta  = cos(beta);

      if (m_equatorialAspect) {
        // If phi1 == 0, the general case below cancels alegbraically to 
        // the following equations, however, the simplified equations are 
        // used here to reduce chance of roundoff errors
        double trigTerm = cosBeta*cosLambdaDiff;
        if (qFuzzyCompare(-1, trigTerm)) { // avoid divide by zero & possible roundoff errors
          // This happens when trying to project the antipodal point.
          // lambda = lambda0 +/- 180 and beta = 0 (that is, phi = 0)
          // Thus denominator = 1 + (1)(-1) = 0
          // We don't allow this for equatorial aspect projections.  This point,
          // at the opposite side of the planet, is projected as a circle of
          // radius 2R - we can't return a unique (x,y) value.
          m_good = false;
          return m_good;
        }
        double denominator = 1 + trigTerm;
        x = m_a * cosBeta * sinLambdaDiff * sqrt(2/denominator); 
            // Snyder eqn (24-21) pg 187
        y = (m_Rq*m_Rq/m_a) * sinBeta * sqrt(2/denominator); 
            // Snyder eqn (24-22) pg 187
            // Note: Clearly 1 + cosBeta*cosLambdaDiff >= 0
            // 1 + cosBeta*cosLambdaDiff = 0 implies 2 cases
            // 1) cosBeta = 1 and cosLambdaDiff = -1
            //    So, beta = 0 and lambdaDiff = 180
            //    Then, q = 0, which implies phi = 0 = phi1
            //    this only happens when projecting opposite side of equator
            // 2) cosBeta = -1 and cosLambdaDiff = 1
            //    So, beta = 180 impossible since beta is defined by asin 
            //    (range of asin is -90 to 90)
            //    Then, q = 0, which implies phi = 0 = phi1
            //    this only happens when projecting opposite side of equator
            // Therefore x and y are well defined
      }
      else { // ellipsoidal General oblique aspect
        double sinProduct = m_sinBeta1 * sinBeta;
        double cosProduct = m_cosBeta1 * cosBeta * cosLambdaDiff;
        double trigTerms = sinProduct + cosProduct;
        if (qFuzzyCompare(-1, trigTerms)) { // avoid divide by zero & possible roundoff errors
          // This happens when trying to project the antipodal point.
          // lambda = lambda0 + 180 and beta = -beta1 (that is, phi = -phi1)
          // Then, by odd symmetry of sine function, sinBeta = -sinBeta1,
          // by even symmetry of cosine, cosBeta = cosBeta1
          // and cosLambdaDiff = cos(180) = -1
          // Thus denominator = 1 - sin^2(beta) - cos^2(beta) = 1-1 = 0
          // We don't allow this for oblique aspect projections.  This point,
          // at the opposite side of the planet, is projected as a circle of
          // radius 2R - we can't return a unique (x,y) value.
          m_good = false;
          return m_good;
        }
        double denominator = 1 + trigTerms;
        double kprime = sqrt(2/denominator); // Snyder eqn (24-2) page 185
        double B = m_Rq * kprime; // Snyder eqns (24-19) pg 187

        x = B * m_D * cosBeta*sinLambdaDiff; // Snyder eqn (24-17) pg 187
        y = (B/m_D) * (m_cosBeta1*sinBeta - m_sinBeta1*cosBeta*cosLambdaDiff); 
            // Snyder eqn (24-18) pg 187
            // Note: y is well defined since D = 0 implies a = 0 or m1 = 0
            //       we know a > 0,
            //       so m1 = 0 would imply cosPhi1 = 0, thus phi1 = +/- 90
            //       but this calculation is not used for polar aspect,
            //       so y doesn't have zero denominator
      }
      // There are no ellipsoidal values for the scale factors h and k unless 
      // using the polar aspects, Snyder pg 26.
    }
    SetComputedXY(x, y); // sets m_x = x and m_y = y and handles rotation

    m_good = true;
    return m_good;
  }

  bool LambertAzimuthalEqualArea::SetCoordinate(const double x, 
                                                const double y) {
    if (x == Null || y == Null) {
      m_good = false;
      return m_good;
    }
    // Save the coordinate
    SetXY(x, y);

    // 2012-06-11
    // Uncomment following line if ellipsoid projection is implemented
    if (!m_spherical) return setCoordinateEllipsoid(x,y);

    double phi;    // latitude value to be calculated
    double lambda; // longitude value to be calculated

    // Spherical Radius = Equatorial Radius, so
    // Snyder's R = Snyder's a variable pg ix
    double R = m_a;

    double rho = sqrt(x*x+y*y); // Snyder eqn (20-18) pg 187
    if (rho < DBL_EPSILON) { // this implies (x,y) = (0,0), 
                             // so project to center lat/lon values
      phi = m_phi1;
      lambda = m_lambda0;
    }
    else {
      if (fabs(rho/(2*R)) > 1 + DBL_EPSILON) {// given (x,y) point is off planet 
                                             // projection if distance from
                                             // origin (rho) is greater than
                                             // 2*radius of sphere 
          m_good = false;
          return m_good;
      }
      else if (fabs(rho/(2*R)) > 1) { // if ratio is slightly larger than 1 
                                     // due to rounding error, then fix so we 
                                     // can take the arcsin
        rho = 2*R;
      }
      double c = 2*asin(rho/(2*R)); // c is angular distance, 
                                    //Snyder eqn (24-16) pg 187
      double sinC = sin(c);
      double cosC = cos(c);

      // verify the following is in the domain of the arcsin function
      if (fabs(cosC*m_sinPhi1+y*sinC*m_cosPhi1/rho) > 1) {
        m_good = false;
        return m_good;
      }
      phi = asin(cosC*m_sinPhi1+y*sinC*m_cosPhi1/rho);//Snyder eq (20-14) pg 186
      if (m_northPolarAspect ) { // spherical north polar aspect
        lambda = m_lambda0 + atan2(x,-y); // Snyder eqn (20-16) pg 187, 
                                          // need to use atan2 to correct for
                                          // quadrant - see pg 150
      }
      else if (m_southPolarAspect ) { // spherical north polar aspect
        lambda = m_lambda0 + atan2(x,y); // Snyder eqn (20-17) pg 187, 150
      }
      else { // spherical oblique aspect
        double denominator = rho*m_cosPhi1*cosC - y*m_sinPhi1*sinC;
        // atan2 well defined for denominator = 0 
        // if numerator is also 0, returns 0
        lambda = m_lambda0 + atan2(x * sinC, denominator);
                 // Snyder eqn (20-15) pg 186
      }
    }

    m_latitude = phi * 180.0 / PI;
    m_longitude = lambda * 180.0 / PI;

    // Cleanup the longitude
    if (m_longitudeDirection == PositiveWest) {
      m_longitude *= -1.0;
    }

    if (m_longitudeDomain == 180) {
      m_longitude = To180Domain(m_longitude);
    }
    else {
      // Do this because longitudeDirection could cause (-360,0)
      m_longitude = To360Domain(m_longitude);
    }

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

    m_good = true;
    return m_good;
  }

  bool LambertAzimuthalEqualArea::setCoordinateEllipsoid(const double x, 
                                                         const double y) { 
    double phi;    // latitude value to be calculated
    double lambda; // longitude value to be calculated

    // ellipsoidal projection
    double q;
    // for polar aspects
    double rho = sqrt(x*x+y*y); // Snyder eqn (20-18) pg 190 
    if (m_northPolarAspect ) { // ellipsoidal north polar aspect
      q = m_qp - (rho*rho/(m_a*m_a)); // Snyder eqn (24-31) pg 190
      lambda = m_lambda0 + atan2(x,-y); // Snyder eqn (20-16) pgs 190, 150
    }
    else if (m_southPolarAspect ) { // ellipsoidal south polar aspect
      q = -1.0*(m_qp - (rho*rho/(m_a*m_a))); // Snyder eqn (24-31) pg 190
      lambda = m_lambda0 + atan2(x,y); // Snyder eqn (20-17) pg 190, 150
    }
    else {// ellipsoidal oblique aspect 
      double xD = x/m_D; // m_D = 0 implies polar aspect 
                         // (see setGroundEllipsoid() general oblique case
      double Dy = m_D*y;
      double rho = sqrt((xD*xD)+(Dy*Dy)); // Snyder eqn (24-28) pg 189
      if (fabs(rho) > fabs(2*m_Rq)) {
        m_good = false;
        return m_good;
      }
      double Ce = 2*asin(rho/(2*m_Rq)); // Snyder eqn (24-29) pg 189 
                                        // TYPO IN TEXT, parentheses required 
                                        // around denominator (see usage pg 335)
           // Note: Ce well defined since Rq = 0 is impossible (a > 0 and q > 1)
      double sinCe = sin(Ce);
      double cosCe = cos(Ce);

      if (rho < DBL_EPSILON) { // pg 189, first line appears to be typo in text, 
                               // below is our interpretation.
        // if rho = 0, then x=0 and (D=0 or y=0), so it
        // makes sense that lambda = lambda0 in this case.
        q = m_qp * m_sinBeta1;
        lambda = m_lambda0;
      }
      else {
        q = m_qp * (cosCe * m_sinBeta1 + m_D * y * sinCe * m_cosBeta1 / rho);
            // Snyder eqn (24-27) pg 188
        double numerator = x * sinCe;
        double denominator = m_D * rho * m_cosBeta1 * cosCe 
                             - m_D * m_D * y * m_sinBeta1 * sinCe;
        // atan2 well defined for denominator = 0 
        // if numerator is also 0, returns 0
        lambda = m_lambda0 + atan2(numerator, denominator);
                 // Snyder eqn (24-26) pg 188
      }
    }

    if (qFuzzyCompare(fabs(q), fabs(m_qp))) { 
      phi = copysign(HALFPI,q);// Snyder eqn (14-20) pg 189, 
                                        //(see definition of qp on pg 187)
    }
    else {
      if (fabs(q) > 2) {
        // verify that q/2 is in the domain of arcsin
        m_good = false;
        return m_good;
      }
      phi = asin(q/2); // Snyder pg 189 above (14-20) describes iteration 
                       // process and initial value for phi
      double phiOld = phi;
      bool converge = false;
      double tolerance = 1e-10; // tolerance same as Projection::phi2Compute()
                               // should we attempt to relate tolerance to pixel resolution ???
      int maximumIterations = 100;
      int i = 0;
      while (!converge) {
        i++;
        // calculate values to reduce calculations:
        double sinPhi = sin(phi);
        double eSinPhi = m_e*sinPhi;
        double cosPhi = cos(phi);
        double squareESinPhi = eSinPhi*eSinPhi;
        double oneMinusSquareESinPhi = 1 - squareESinPhi;
        // find next iteration of phi
        phi += oneMinusSquareESinPhi*oneMinusSquareESinPhi / (2 * cosPhi)
               * (q / (1 - m_e * m_e) 
                  - sinPhi / (oneMinusSquareESinPhi) 
                  + log( (1 - eSinPhi) / (1 + eSinPhi) ) / (2 * m_e));
               //eqn (3-16) pg 188
        if (fabs(phiOld - phi) < tolerance) {
          converge = true;
        }
        if (i > maximumIterations) {
          m_good = false;
          return m_good;
        }
        phiOld = phi;
      }
    }
    m_latitude = phi * 180.0 / PI;
    m_longitude = lambda * 180.0 / PI;

    // Cleanup the longitude
    if (m_longitudeDirection == PositiveWest) {
      m_longitude *= -1.0;
    }

    if (m_longitudeDomain == 180) {
      m_longitude = To180Domain(m_longitude);
    }
    else {
      // Do this because longitudeDirection could cause (-360,0)
      m_longitude = To360Domain(m_longitude);
    }

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

    m_good = true;
    return m_good;
  }

  bool LambertAzimuthalEqualArea::XYRange(double &minX, double &maxX, 
                                          double &minY, double &maxY) {
    // Global Polar projection
    if ((m_northPolarAspect && qFuzzyCompare(-90.0, MinimumLatitude()))
        || (m_southPolarAspect && qFuzzyCompare(90.0, MaximumLatitude()))) {
      // for polar aspect projections, the antipodal point is the opposite pole.
      // if it is included in lat range, the bounding circle will exist, no
      // matter the longitude range.
      double eRad = EquatorialRadius();
      double pRad = PolarRadius();
      double maxCoordValue = 0;

      if (m_spherical) {
        maxCoordValue = 2*EquatorialRadius();
      }
      else {
        maxCoordValue = sqrt(2*eRad*eRad 
                             + pRad*pRad*log((1+m_e)/(1-m_e))/m_e);
      }
      m_minimumX = -maxCoordValue;
      m_maximumX =  maxCoordValue;
      m_minimumY = -maxCoordValue;
      m_maximumY =  maxCoordValue;
    } 
    // Polar projection, not including antipodal point
    else if (m_northPolarAspect || m_southPolarAspect) {
      return xyRangeLambertAzimuthalPolar(minX, maxX, minY, maxY);
    }
    // Oblique projection
    else { // oblique cases (includes equatorial aspect)
      if (xyRangeOblique(minX, maxX, minY, maxY)) {
        // Make sure that the calculations didn't go beyond the accepatable x,y
        // values.  As far as we know, the x, y values should not exceed
        // 2*LocalRadius(-phi)
        double maxCoordValue = 2*LocalRadius(-m_phi1*180/PI);
        if (m_minimumX < -maxCoordValue) m_minimumX = -maxCoordValue;
        if (m_maximumX >  maxCoordValue) m_maximumX =  maxCoordValue;
        if (m_minimumX < -maxCoordValue) m_minimumX = -maxCoordValue;
        if (m_maximumY >  maxCoordValue) m_maximumY =  maxCoordValue;
      }
      else {
        return false;
      }
    }

    // If we haven't already returned, we have a Global Polar projection or an
    // Oblique projection that did not fail xyRangeOblique
    
    // 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;
  }
  
  bool LambertAzimuthalEqualArea::xyRangeLambertAzimuthalPolar(double &minX, 
                                                               double &maxX, 
                                                               double &minY, 
                                                               double &maxY) {
    // Test the four combinations of min/max lat/lon
    XYRangeCheck(MinimumLatitude(), MinimumLongitude());
    XYRangeCheck(MinimumLatitude(), MaximumLongitude());
    XYRangeCheck(MaximumLatitude(), MinimumLongitude());
    XYRangeCheck(MaximumLatitude(), MaximumLongitude());

    double centerLonDeg = m_lambda0 * 180.0 / PI;
    if (m_longitudeDirection == PositiveWest) centerLonDeg = centerLonDeg * -1.0;

    // Since this projection is Polar aspect, 4 longitudes will lie
    // directly on the horizontal and vertical axes of the
    // projection.
    // test combinations with each of these longitudes:
    // 
    // down  (negative vertical axis) - center long for north polar
    // up    (positive vertical axis) - center long for south polar
    // left  (negative horizontal axis)
    // right (positive horizontal axis)
    // 
    for (double lon = centerLonDeg; lon < centerLonDeg + 360; lon += 90) {
      checkLongitude(lon);
    }

    // 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 LambertAzimuthalEqualArea::checkLongitude(double longitude) {

    double centerLatDeg = m_phi1 * 180.0 / PI;

    double innerLatitude, outerLatitude;
    if (m_northPolarAspect) {
      innerLatitude = MaximumLatitude();
      outerLatitude = MinimumLatitude();
    }
    else if (m_southPolarAspect) {
      innerLatitude = MinimumLatitude();
      outerLatitude = MaximumLatitude();
    }
    else {
      IString message = "checkLongitude() should only be called for a "
                    "polar aspect projection. CenterLatitude is [";
      message = message + IString(centerLatDeg) + "] degrees.";
      throw IException(IException::Programmer, message, _FILEINFO_);
    }

    // Check whether longitude is in the min/max longitude range
    bool thisLongitudeInMinMaxRange = false;
    if (MaximumLongitude() - MinimumLongitude() == 360) {
      thisLongitudeInMinMaxRange = true;
    }
    double adjustedLon =    To360Domain(longitude);
    double adjustedMinLon = To360Domain(MinimumLongitude());
    double adjustedMaxLon = To360Domain(MaximumLongitude()); 

    if (adjustedMinLon > adjustedMaxLon) {
      if (adjustedLon > adjustedMinLon) {
        adjustedLon -= 360;
      }
      adjustedMinLon -= 360;
    }
    // if lon value for this axis is between min lon and max lon
    if (adjustedMinLon <= adjustedLon 
        && adjustedLon <= adjustedMaxLon) { 
      thisLongitudeInMinMaxRange = true;
    }
    else {
      thisLongitudeInMinMaxRange = false;
    }

    if (thisLongitudeInMinMaxRange) {
      // any shape that includes this longitude
      XYRangeCheck(outerLatitude, longitude);
    }
    else {
      // determine which boundary value (minlon or maxlon)
      // is closer to given longitude 
      double distToMinLon  = fabs(adjustedMinLon - adjustedLon);
      double distToMaxLon = fabs(adjustedMaxLon - adjustedLon);

      if (distToMinLon >= 180 ) {
        distToMinLon = fabs(360 - distToMinLon);
      }
      if (distToMaxLon >= 180 ) {
        distToMaxLon = fabs(360 - distToMaxLon);
      }

      double nearestBoundary = 0;
      if (distToMinLon < distToMaxLon) {
        nearestBoundary = MinimumLongitude();
      }
      else {
        nearestBoundary = MaximumLongitude();
      }

      // shapes that come within 90 degrees of the given longitude
      if (distToMinLon <= 90 + DBL_EPSILON 
          || distToMaxLon <= 90 + DBL_EPSILON) {
        XYRangeCheck(outerLatitude, nearestBoundary);
      }
      // shapes more than 90 degrees from the longitude that include the origin
      else if (qFuzzyCompare(MaximumLatitude(), centerLatDeg)) {
        XYRangeCheck(centerLatDeg,longitude);
      }
      // shapes more than 90 degrees from the longitude without the origin
      else {
        XYRangeCheck(innerLatitude, nearestBoundary);
      }
    }
    return;
  }

  PvlGroup LambertAzimuthalEqualArea::Mapping() {
    PvlGroup mapping = TProjection::Mapping();

    mapping += m_mappingGrp["CenterLatitude"];
    mapping += m_mappingGrp["CenterLongitude"];

    return mapping;
  }

  PvlGroup LambertAzimuthalEqualArea::MappingLatitudes() {
    PvlGroup mapping = TProjection::MappingLatitudes();

    mapping += m_mappingGrp["CenterLatitude"];

    return mapping;
  }

  PvlGroup LambertAzimuthalEqualArea::MappingLongitudes() {
    PvlGroup mapping = TProjection::MappingLongitudes();

    mapping += m_mappingGrp["CenterLongitude"];

    return mapping;
  }

  // There are no ellipsoidal values for the scale factors h and k unless 
  // using the polar aspects, Snyder pg 26.

  // These methods are not currently used in other Projection classes. However,
  // if similar methods are implemented for other projections, we should create
  // a virtual prototype in the parent Projection. Currently since they do not
  // exist in Projection, we can not call them if ProjectionFactory was used to
  // create the object.
  double LambertAzimuthalEqualArea::relativeScaleFactorLongitude() const {
    validateRelativeScaleFactor();
    return m_h;
  };

  double LambertAzimuthalEqualArea::relativeScaleFactorLatitude() const {
    validateRelativeScaleFactor();
    return m_k;
  }

  void LambertAzimuthalEqualArea::validateRelativeScaleFactor() const {
    if (!m_good) {
      IString message = "Projection failed or ground and coordinates have not "
                        "been set.  Relative scale factor can not be computed.";
      throw IException(IException::Unknown, message, _FILEINFO_);
    }
    if (!m_spherical && !(m_northPolarAspect || m_southPolarAspect)) {
      IString message = "For ellipsidal bodies, relative scale factor can only "
                        "be computed for polar aspect projections.";
      throw IException(IException::Unknown, message, _FILEINFO_);
    }
    if (m_northPolarAspect && qFuzzyCompare(-90.0, m_latitude)) {
        IString message = "Relative scale factor can not be computed for north "
                          "polar aspect projection when ground is set to "
                          "latitude -90.";
      throw IException(IException::Unknown, message, _FILEINFO_);
    }
    if (m_southPolarAspect && qFuzzyCompare(90.0, m_latitude)) {
        IString message = "Relative scale factor can not be computed for south "
                          "polar aspect projection when ground is set to "
                          "latitude 90.";
      throw IException(IException::Unknown, message, _FILEINFO_);
    }
    if (m_k == Null || m_h == Null) {
      IString message = "Relative scale factor can not be computed.";
      throw IException(IException::Unknown, message, _FILEINFO_);
    }
  }

} // end namespace isis

extern "C" Isis::Projection *LambertAzimuthalEqualAreaPlugin(Isis::Pvl &lab,
    bool allowDefaults) {
  return new Isis::LambertAzimuthalEqualArea(lab, allowDefaults);
}