Sensor.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#include "Sensor.h"
 
#include <iomanip>
 
#include <QDebug>
#include <QString>
 
#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 {
    std::vector<double> sunB(m_uB, m_uB+3);
    return target()->shape()->phaseAngle(
                               bodyRotation()->ReferenceVector(instrumentPosition()->Coordinate()), sunB);
  }
 
 
  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(3);
    lookB[0] = m_lookB[0];
    lookB[1] = m_lookB[1];
    lookB[2] = m_lookB[2];
    vector<double> lookJ = bodyRotation()->J2000Vector(lookB);
 
    SpiceDouble range;
    recrad_c((SpiceDouble *)&lookJ[0], &range, &m_ra, &m_dec);
    m_ra *= 180.0 / PI;
    m_dec *= 180.0 / PI;
  }
 
 
  bool Sensor::SetRightAscensionDeclination(const double ra, const double dec) {
    vector<double> lookJ(3);
    radrec_c(1.0, ra * PI / 180.0, dec * PI / 180.0, (SpiceDouble *)&lookJ[0]);
 
    vector<double> lookC = instrumentRotation()->ReferenceVector(lookJ);
    return SetLookDirection((double *)&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 {
    SpiceDouble psB[3], upsB[3];
    SpiceDouble dist;
 
    std::vector<double> sB = bodyRotation()->ReferenceVector(instrumentPosition()->Coordinate());
 
    SpiceDouble pB[3];
    ShapeModel *shape = target()->shape();
    pB[0] = shape->surfaceIntersection()->GetX().kilometers();
    pB[1] = shape->surfaceIntersection()->GetY().kilometers();
    pB[2] = shape->surfaceIntersection()->GetZ().kilometers();
 
    vsub_c(pB, (SpiceDouble *) &sB[0], psB);
    unorm_c(psB, upsB, &dist);
    return dist;
  }
 
 
  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);
 
    // Calc the change
    ShapeModel *shape = target()->shape();
    double xChange = sB[0] - shape->surfaceIntersection()->GetX().kilometers();
    double yChange = sB[1] - shape->surfaceIntersection()->GetY().kilometers();
    double zChange = sB[2] - shape->surfaceIntersection()->GetZ().kilometers();
 
    // Calc the distance and convert to AU
    double dist = sqrt(xChange*xChange + yChange*yChange + zChange*zChange);
    dist /= 149597870.691;
    return dist;
  }
 
 
  double Sensor::SpacecraftAltitude() {
    // Get the spacecraft coord
    double spB[3];
    Spice::instrumentPosition(spB);
 
    // Get subspacecraft point
    double lat, lon;
    subSpacecraftPoint(lat, lon);
    double rlat = lat * PI / 180.0;
    double rlon = lon * PI / 180.0;
 
    // Compute radius
    Distance rad = LocalRadius(lat, lon);
 
    // Now with the 3 spherical value compute the x/y/z coordinate
    double ssB[3];
    latrec_c(rad.kilometers(), rlon, rlat, ssB);
 
    // Calc the change
    double xChange = spB[0] - ssB[0];
    double yChange = spB[1] - ssB[1];
    double zChange = spB[2] - ssB[2];
 
    // Calc the distance
    double dist = sqrt(xChange*xChange + yChange*yChange + zChange*zChange);
    return dist;
  }
 
 
//  Distance Sensor::DemRadius(const SurfacePoint &pt) {
//    return DemRadius(pt.GetLatitude(), pt.GetLongitude());
//  }
 
 
//  Distance Sensor::DemRadius(const Latitude &lat, const Longitude &lon) {
//    if (!m_hasElevationModel) return Distance();
//    //if (!lat.Valid() || !lon.Valid()) return Distance();
//    m_demProj->SetUniversalGround(lat.degrees(), lon.degrees());
//    if (!m_demProj->IsGood()) {
//      return Distance();
//    }
 
//    m_portal->SetPosition(m_demProj->WorldX(), m_demProj->WorldY(), 1);
 
//    m_demCube->read(*m_portal);
 
//    const double &radius = m_interp->Interpolate(m_demProj->WorldX(),
//                                                 m_demProj->WorldY(),
//                                                 m_portal->DoubleBuffer());
 
//    return Distance(radius, Distance::Meters);
//      Distance fred;
//      return fred;
//  }
 
 
//  double Sensor::DemRadius(double lat, double lon) {
//    if (!m_demProj->SetUniversalGround(lat, lon)) {
//      return Isis::Null;
//    }
 
//    m_portal->SetPosition(m_demProj->WorldX(), m_demProj->WorldY(), 1);
//    m_demCube->read(*m_portal);
 
//    double radius = m_interp->Interpolate(m_demProj->WorldX(),
//                                          m_demProj->WorldY(),
//                                          m_portal->DoubleBuffer());
//    if (Isis::IsSpecial(radius)) {
//      return Isis::Null;
//    }
 
//    return radius / 1000.0;
//      double fred;
//      return fred;
//  }
//   bool HasElevationModel() {
//     return m_hasElevationModel;
//   };
 
}