Sensor.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#include "Sensor.h"
 
#include <iomanip>
 
#include <QDebug>
#include <QString>
 
#include "MathUtils.h"
 
#include "Angle.h"
#include "Constants.h"
#include "CubeManager.h"
#include "Distance.h"
#include "EllipsoidShape.h"
#include "IException.h"
#include "IString.h"
#include "iTime.h"
#include "Latitude.h"
#include "Longitude.h"
#include "NaifStatus.h"
#include "Projection.h"
#include "ShapeModel.h"
#include "SpecialPixel.h"
#include "SurfacePoint.h"
#include "Target.h"
#include "UniqueIOCachingAlgorithm.h"
 
#define MAX(x,y) (((x) > (y)) ? (x) : (y))
 
using namespace std;
 
namespace Isis {
 
  Sensor::Sensor(Cube &cube) : Spice(cube) {
    m_newLookB = false;
  }
 
 
  Sensor::~Sensor() {
  }
 
 
  void Sensor::IgnoreElevationModel(bool ignore) {
    // if we have an elevation model and are not ignoring it,
    //   set p_hasElevationModel to true
    if (!ignore) {
      target()->restoreShape();
    }
    else {
      target()->setShapeEllipsoid();
    }
  }
 
 
  QList<QPointF> Sensor::PixelIfovOffsets() {
 
    QString message = "Pixel Ifov offsets not implemented for this camera.";
    throw IException(IException::Programmer, message, _FILEINFO_);
  }
 
 
  void Sensor::setTime(const iTime &time) {
    Spice::setTime(time);
    target()->shape()->clearSurfacePoint();
  }
 
 
  bool Sensor::SetLookDirection(const double v[3]) {
    // std::cout << "Sensor::SetLookDirection()\n";
    // The look vector must be in the camera coordinate system
 
    // copy v to LookC
    // lookC[0] = v[0];
    // lookC[1] = v[1];
    // lookC[2] = v[2]
    vector<double> lookC(v, v + 3);
 
    // Convert it to body-fixed
    const vector<double> &lookJ = instrumentRotation()->J2000Vector(lookC);
    const vector<double> &lookB = bodyRotation()->ReferenceVector(lookJ);
 
    // This memcpy does:
    // m_lookB[0] = lookB[0];
    // m_lookB[1] = lookB[1];
    // m_lookB[2] = lookB[2];
    memcpy(m_lookB, &lookB[0], sizeof(double) * 3);
    m_newLookB = true;
 
    // Don't try to intersect the sky
    if (target()->isSky()) {
      target()->shape()->setHasIntersection(false);
      return false;
    }
 
    // See if it intersects the planet
    const vector<double> &sB = bodyRotation()->ReferenceVector(
        instrumentPosition()->Coordinate());
 
    // double tolerance = resolution() / 100.0; return
    // target()->shape()->intersectSurface(sB, lookB, tolerance);
    return target()->shape()->intersectSurface(sB, lookB);
  }
 
 
  bool Sensor::HasSurfaceIntersection() const {
    return target()->shape()->hasIntersection();
  }
 
 
  void Sensor::Coordinate(double p[3]) const {
    ShapeModel *shape = target()->shape();
    p[0] = shape->surfaceIntersection()->GetX().kilometers();
    p[1] = shape->surfaceIntersection()->GetY().kilometers();
    p[2] = shape->surfaceIntersection()->GetZ().kilometers();
  }
 
 
  double Sensor::UniversalLatitude() const {
    return GetLatitude().degrees();
  }
 
 
  Latitude Sensor::GetLatitude() const {
    return target()->shape()->surfaceIntersection()->GetLatitude();
  }
 
 
  double Sensor::UniversalLongitude() const {
    return GetLongitude().degrees();
  }
 
 
  Longitude Sensor::GetLongitude() const {
    return target()->shape()->surfaceIntersection()->GetLongitude();
  }
 
 
  SurfacePoint Sensor::GetSurfacePoint() const {
    return *(target()->shape()->surfaceIntersection());
  }
 
 
  Distance Sensor::LocalRadius() const {
    //TODO: We probably need to be validating surface intersect point is not NULL
    // Here? or in ShapeModel? or what, man?
    return target()->shape()->surfaceIntersection()->GetLocalRadius();
  }
 
 
  Distance Sensor::LocalRadius(Latitude lat, Longitude lon) {
    return target()->shape()->localRadius(lat, lon);
  }
 
 
  Distance Sensor::LocalRadius(double lat, double lon) {
    return target()->shape()->localRadius(Latitude(lat, Angle::Degrees),
                                          Longitude(lon, Angle::Degrees));
  }
 
 
  double Sensor::PhaseAngle() const {
    double sensorPosBf[3];
    instrumentBodyFixedPosition(sensorPosBf);
    double groundPtBf[3];
    Coordinate(groundPtBf);
    return RAD2DEG * SensorUtilities::sepAngle(sensorPosBf, groundPtBf, m_uB);
  }
 
 
  double Sensor::EmissionAngle() const {
    return target()->shape()->emissionAngle(
        bodyRotation()->ReferenceVector(instrumentPosition()->Coordinate()));
  }
 
 
  double Sensor::IncidenceAngle() const {
    std::vector<double> sunB(m_uB, m_uB+3);
    return target()->shape()->incidenceAngle(sunB);
  }
 
 
  bool Sensor::SetUniversalGround(const double latitude,
                                  const double longitude,
                                  bool backCheck) {
 
    ShapeModel *shape = target()->shape();
    shape->clearSurfacePoint();
 
    // Can't intersect the sky
    if (target()->isSky()) {
      return false;
    }
 
    // Load the latitude/longitude
    Latitude lat(latitude, Angle::Degrees);
    Longitude lon(longitude, Angle::Degrees);
    // Local radius is deferred to (possible derived) shape model method
    shape->intersectSurface(lat, lon,
                            bodyRotation()->ReferenceVector(instrumentPosition()->Coordinate()),
                            backCheck);
 
    return SetGroundLocal(backCheck);
  }
 
  bool Sensor::SetUniversalGround(const double latitude,
                                  const double longitude,
                                  const double radius,
                                  bool backCheck) {
 
    ShapeModel *shape = target()->shape();
    shape->clearSurfacePoint();
 
   // Can't intersect the sky
    if (target()->isSky()) {
      return false;
    }
 
    Latitude lat(latitude, Angle::Degrees);
    Longitude lon(longitude, Angle::Degrees);
    Distance rad(radius, Distance::Meters);
 
    shape->intersectSurface(SurfacePoint(lat, lon, rad),
                            bodyRotation()->ReferenceVector(instrumentPosition()->Coordinate()),
                            backCheck);
 
    return SetGroundLocal(backCheck);
  }
 
 
  bool Sensor::SetGround(const SurfacePoint &surfacePt, bool backCheck) {
    //std::cout << "Sensor::SetGround()\n";
    ShapeModel *shape = target()->shape();
    shape->clearSurfacePoint();
 
    // Can't intersect the sky
    if (target()->isSky()) {
      return false;
    }
 
    shape->intersectSurface(surfacePt,
                            bodyRotation()->ReferenceVector(instrumentPosition()->Coordinate()),
                            backCheck);
 
    return SetGroundLocal(backCheck);
  }
 
 
  bool Sensor::SetGroundLocal(bool backCheck) {
    ShapeModel *shape = target()->shape();
    // With the 3 spherical value compute the x/y/z coordinate
    //latrec_c(m_radius, (m_longitude * PI / 180.0), (m_latitude * PI / 180.0), m_pB);
 
 
    if (!(shape->hasIntersection())) {
      return false;
    }
 
    // Make sure the point isn't on the backside of the body
 
    // This is static purely for performance reasons. A significant speedup
    // is achieved here.
    const vector<double> &sB =
        bodyRotation()->ReferenceVector(instrumentPosition()->Coordinate());
 
    m_lookB[0] = shape->surfaceIntersection()->GetX().kilometers() - sB[0];
    m_lookB[1] = shape->surfaceIntersection()->GetY().kilometers() - sB[1];
    m_lookB[2] = shape->surfaceIntersection()->GetZ().kilometers() - sB[2];
    m_newLookB = true;
 
    // See if the point is on the backside of the target. Note occlusion handling
    // now handles this in derived shape models that can support it. (KJB 2017-03-23)
    // This may be good if the computation of the look direction is more sophisticated.
    if (backCheck) {
      // Assume the intersection point is good in order to get the emission angle
      // shape->setHasIntersection(true);  //KJB there should be a formal intersection in ShapeModel
      std::vector<double> lookdir = lookDirectionBodyFixed();
      if ( !shape->isVisibleFrom(sB, lookdir) ) {
        shape->clearSurfacePoint();
        shape->setHasIntersection(false);
        return false;
      }
    }
 
    // return with success
    shape->setHasIntersection(true);
 
    return true;
  }
 
 
  void Sensor::LookDirection(double v[3]) const {
    vector<double> lookC = instrumentRotation()->ReferenceVector(lookDirectionJ2000());
    v[0] = lookC[0];
    v[1] = lookC[1];
    v[2] = lookC[2];
  }
 
  vector<double> Sensor::lookDirectionBodyFixed() const {
    vector<double> lookB(3);
    lookB[0] = m_lookB[0];
    lookB[1] = m_lookB[1];
    lookB[2] = m_lookB[2];
    return lookB;
  }
 
 
  vector<double> Sensor::lookDirectionJ2000() const {
    vector<double> lookJ = bodyRotation()->J2000Vector(lookDirectionBodyFixed());
    return lookJ;
  }
 
 
 
  double Sensor::RightAscension() {
    if (m_newLookB) {
      computeRaDec();
    }
    return m_ra;
  }
 
 
  double Sensor::Declination() {
    if (m_newLookB) {
      computeRaDec();
    }
    return m_dec;
  }
 
 
  void Sensor::computeRaDec() {
    m_newLookB = false;
    vector<double> lookB = {m_lookB[0], m_lookB[1], m_lookB[2]};
    vector<double> lookJ = bodyRotation()->J2000Vector(lookB);
 
    SensorUtilities::GroundPt3D sphericalPt = SensorUtilities::rectToSpherical(&lookJ[0]);
    m_ra = sphericalPt.lon;
    // Convert to [0, 2pi] domain
    if (m_ra < 0) {
      m_ra += 2 * PI;
    }
    m_ra *= RAD2DEG;
    m_dec = sphericalPt.lat * RAD2DEG;
  }
 
 
  bool Sensor::SetRightAscensionDeclination(const double ra, const double dec) {
    double raRad = ra * DEG2RAD;
    double decRad = dec * DEG2RAD;
 
    // Make the radius bigger, some multiple of the body radius -or- use sensor position at the reference point 
    SensorUtilities::GroundPt3D sphericalPt = {decRad, raRad, 1};
    SensorUtilities::Vec rectPt = SensorUtilities::sphericalToRect(sphericalPt);
 
    vector<double> lookC = instrumentRotation()->ReferenceVector(rectPt);
    return SetLookDirection(&lookC[0]);
  }
 
 
  void Sensor::SpacecraftSurfaceVector(double scSurfaceVector[3]) const {
    scSurfaceVector[0] = m_lookB[0];
    scSurfaceVector[1] = m_lookB[1];
    scSurfaceVector[2] = m_lookB[2];
  }
 
 
  double Sensor::SlantDistance() const {
    double sensorPos[3];
    instrumentBodyFixedPosition(sensorPos);
 
    SurfacePoint *surfPoint = target()->shape()->surfaceIntersection();
    SensorUtilities::Vec groundPt = {
          surfPoint->GetX().kilometers(),
          surfPoint->GetY().kilometers(),
          surfPoint->GetZ().kilometers()};
 
    return SensorUtilities::distance(sensorPos, groundPt);
  }
 
 
  double Sensor::LocalSolarTime() {
    double slat, slon;
    subSolarPoint(slat, slon);
 
    double lst = UniversalLongitude() - slon + 180.0;
    lst = lst / 15.0;  // 15 degress per hour
    if (lst < 0.0) lst += 24.0;
    if (lst > 24.0) lst -= 24.0;
    return lst;
  }
 
 
  double Sensor::SolarDistance() const {
    // Get the sun coord
    double sB[3];
    Spice::sunPosition(sB);
 
    SurfacePoint *intersection = target()->shape()->surfaceIntersection();
    SensorUtilities::Vec groundPt = {
          intersection->GetX().kilometers(),
          intersection->GetY().kilometers(),
          intersection->GetZ().kilometers()};
 
    return SensorUtilities::distance(groundPt, sB) / 149597870.691;
  }
 
 
  double Sensor::SpacecraftAltitude() {
    // Get the spacecraft coord
    double sensorPos[3];
    Spice::instrumentPosition(sensorPos);
 
    // Get subspacecraft point
    double lat, lon;
    subSpacecraftPoint(lat, lon);
 
    // Compute radius
    Distance rad = LocalRadius(lat, lon);
 
    // Take the difference
    return SensorUtilities::magnitude(sensorPos) - rad.kilometers();
  }
}