Robinson.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#include "Robinson.h"
 
#include <cmath>
#include <cfloat>
 
#include "Constants.h"
#include "IException.h"
#include "Projection.h"
#include "Pvl.h"
#include "PvlGroup.h"
#include "PvlKeyword.h"
 
const double EPSILON = 1.0e-10;
 
using namespace std;
namespace Isis {
 
  Robinson::Robinson(Pvl &label, bool allowDefaults) :
      TProjection::TProjection(label) {
    try {
 
      // Initialize coefficients
      // PR  Add xtra element at beginning to mimic 1-based arrays to match Snyder's code for
      // easier debugging.
      m_pr << 0 << -.062 <<
                         0 <<
                     0.062 <<
                     0.124 <<
                     0.186 <<
                     0.248 <<
                     0.310 <<
                     0.372 <<
                     0.434 <<
                     0.4958 <<
                     0.5571 <<
                     0.6176 <<
                     0.6769 <<
                     0.7346 <<
                     0.7903 <<
                     0.8435 <<
                     0.8936 <<
                     0.9394 <<
                     0.9761 <<
                     1.;
      // XLR
      m_xlr << 0 << 0.9986 <<
                          1.0 <<
                       0.9986 <<
                       0.9954 <<
                         0.99 <<
                       0.9822 <<
                        0.973 << 
                         0.96 <<
                       0.9427 <<
                       0.9216 <<
                       0.8962 <<
                       0.8679 <<
                        0.835 << 
                       0.7986 <<
                       0.7597 <<
                       0.7186 <<
                       0.6732 <<
                       0.6213 <<
                       0.5722 <<
                       0.5322;
 
      // 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 ((allowDefaults) && (!mapGroup.hasKeyword("CenterLongitude"))) {
        double lon = (m_minimumLongitude + m_maximumLongitude) / 2.0;
        mapGroup += PvlKeyword("CenterLongitude", toString(lon));
      }
 
      // Get the center longitude
      m_centerLongitude = mapGroup["CenterLongitude"];
 
      // convert to radians, adjust for longitude direction
      m_centerLongitude *= DEG2RAD;
      if (m_longitudeDirection == PositiveWest) m_centerLongitude *= -1.0;
    }
    catch(IException &e) {
      QString message = "Invalid label group [Mapping]";
      throw IException(e, IException::Io, message, _FILEINFO_);
    }
  }
 
  Robinson::~Robinson() {
  }
 
  bool Robinson::operator== (const Projection &proj) {
    if (!Projection::operator==(proj)) return false;
    // dont do the below it is a recusive plunge
    //  if (Projection::operator!=(proj)) return false;
    Robinson *robin = (Robinson *) &proj;
    if (robin->m_centerLongitude != m_centerLongitude) return false;
    return true;
  }
 
  QString Robinson::Name() const {
    return "Robinson";
  }
 
  QString Robinson::Version() const {
    return "1.0";
  }
 
  bool Robinson::SetGround(const double lat, const double lon) {
    // Convert to radians
    m_latitude = lat;
    m_longitude = lon;
    double latRadians = lat * DEG2RAD;
    double lonRadians = lon * DEG2RAD;
    if (m_longitudeDirection == PositiveWest) lonRadians *= -1.0;
 
    // Compute the coordinate
    double deltaLon = (lonRadians - m_centerLongitude);
    double p2 = fabs(latRadians / 5.0 / DEG2RAD);
    long ip1 = (long) (p2 - EPSILON);
    if (ip1 > 17) {
      return false;
    }
 
    // Stirling's interpolation formula (using 2nd Diff.)
    //
    p2 -= (double) ip1;
    double x = 0.8487 * m_equatorialRadius * (m_xlr[ip1 + 2] + p2 * (m_xlr[ip1 + 3] -
                         m_xlr[ip1 + 1]) / 2.0 +
                         p2 * p2 * (m_xlr[ip1 + 3] - 2.0 * m_xlr[ip1 + 2] +
                         m_xlr[ip1 + 1])/2.0) * deltaLon;
 
    double y = 1.3523 * m_equatorialRadius * (m_pr[ip1 + 2] + p2 * (m_pr[ip1 + 3] -
                         m_pr[ip1 +1]) / 2.0 + p2 * p2 * (m_pr[ip1 + 3] -
                         2.0 * m_pr[ip1 + 2] + m_pr[ip1 + 1]) / 2.0);
    if (lat < 0) y *= -1.;
 
    SetComputedXY(x, y);
    m_good = true;
    return m_good;
  }
 
 
  bool Robinson::SetCoordinate(const double x, const double y) {
    m_latitude = 0;
    m_longitude = 0;
    // Save the coordinate
    SetXY(x, y);
 
    double yy = y / m_equatorialRadius / 1.3523;
    double phid = yy * 90.0;
    double p2 = fabs(phid / 5.0);
    long ip1 = (long) (p2 - EPSILON);
    if (ip1 == 0) ip1 = 1;
    if (ip1 > 17) {
      return false;
    }
 
    // Stirling's interpolation formula as used in forward transformation is 
    //   reversed for first estimation of LAT. from rectangular coordinates. LAT.
    //   is then adjusted by iteration until use of forward series provides correct 
    //   value of Y within tolerance.
    //
    double u, v, t, c;
    double y1;
 
    for (int i = 0;;) {
      u = m_pr[ip1 + 3] - m_pr[ip1 + 1];
      v = m_pr[ip1 + 3] - 2.0 * m_pr[ip1 + 2] + m_pr[ip1 + 1];
      t = 2.0 * (fabs(yy) - m_pr[ip1 + 2]) / u;
      c = v / u;
      p2 = t * (1.0 - c * t * (1.0 - 2.0 * c * t));
 
      if ((p2 >= 0.0) || (ip1 == 1)) {
        phid = (p2 + (double) ip1 ) * 5.0;
        if (y < 0) phid *= -1;
 
        do {
          p2 = fabs(phid / 5.0);
          ip1 = (long) (p2 - EPSILON);
          if (ip1 > 17) {
            return false;
          }
          p2 -= (double) ip1;
 
          y1 = 1.3523 * m_equatorialRadius * (m_pr[ip1 +2] + p2 *(m_pr[ip1 + 3] -
                        m_pr[ip1 +1]) / 2.0 + p2 * p2 * (m_pr[ip1 + 3] -
                        2.0 * m_pr[ip1 + 2] + m_pr[ip1 + 1])/2.0);
          if (y < 0) y1 *= -1.;
 
          phid -= (90.0 * (y1 - y) / m_equatorialRadius / 1.3523);
          i++;
          if (i > 75) {
            return false;
          }
        } while (fabs(y1 - y) > .00001);
        break;
      }
      else {
        ip1--;
        if (ip1 < 0) {
          return false;
        }
      }
    }
 
    // Compute latitude and make sure it is not above 90
    m_latitude  = phid;
 
    // Calculate  longitude. using final latitude. with transposed forward Stirling's 
    //   interpolation formula.
    m_longitude = m_centerLongitude + x / m_equatorialRadius / 0.8487 / (m_xlr[ip1 + 2] +
                          p2 * (m_xlr[ip1 + 3] - m_xlr[ip1 + 1]) / 2.0 +
                          p2 * p2 * (m_xlr[ip1 + 3] - 2.0 * m_xlr[ip1 + 2] + 
                          m_xlr[ip1 + 1]) / 2.0);
    m_longitude *= RAD2DEG;
    if (m_longitudeDirection == PositiveWest) m_longitude *= -1.0;
 
    // Our double precision is not good once we pass a certain magnitude of
    //   longitude. Prevent failures down the road by failing now.
    m_good = (fabs(m_longitude) < 1E10);
    return m_good;
  }
 
 
  bool Robinson::XYRange(double &minX, double &maxX,
                           double &minY, double &maxY) {
 
    // Check the corners of the lat/lon range
    XYRangeCheck(m_minimumLatitude, m_minimumLongitude);
    XYRangeCheck(m_maximumLatitude, m_minimumLongitude);
    XYRangeCheck(m_minimumLatitude, m_maximumLongitude);
    XYRangeCheck(m_maximumLatitude, m_maximumLongitude);
 
    // If the latitude crosses the equator check there
    if ((m_minimumLatitude < 0.0) && (m_maximumLatitude > 0.0)) {
      XYRangeCheck(0.0, m_minimumLongitude);
      XYRangeCheck(0.0, m_maximumLongitude);
    }
 
    // 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 Robinson::Mapping()  {
    PvlGroup mapping = TProjection::Mapping();
 
    mapping += m_mappingGrp["CenterLongitude"];
 
    return mapping;
  }
 
  PvlGroup Robinson::MappingLatitudes() {
    PvlGroup mapping = TProjection::MappingLatitudes();
 
    return mapping;
  }
 
  PvlGroup Robinson::MappingLongitudes() {
    PvlGroup mapping = TProjection::MappingLongitudes();
 
    mapping += m_mappingGrp["CenterLongitude"];
 
    return mapping;
  }
 
} // end namespace isis
 
extern "C" Isis::Projection *RobinsonPlugin(Isis::Pvl &lab,
    bool allowDefaults) {
  return new Isis::Robinson(lab, allowDefaults);
}