Spice.cpp

Go to the documentation of this file.

#include "Spice.h"

#include <cfloat>
#include <iomanip>

#include <QDebug>
#include <QVector>

#include "Constants.h"
#include "Distance.h"
#include "EllipsoidShape.h"
#include "EndianSwapper.h"
#include "FileName.h"
#include "IException.h"
#include "IString.h"
#include "iTime.h"
#include "Longitude.h"
#include "LightTimeCorrectionState.h"
#include "NaifStatus.h"
#include "ShapeModel.h"
#include "SpacecraftPosition.h"
#include "Target.h"

#include "getSpkAbCorrState.hpp"

using namespace std;

namespace Isis {
  // TODO: DOCUMENT EVERYTHING
  Spice::Spice(Cube &cube) {
    Pvl &lab = *cube.label();
    PvlGroup kernels = lab.findGroup("Kernels", Pvl::Traverse);
    bool hasTables = (kernels["TargetPosition"][0] == "Table");

    init(lab, !hasTables);
  }

  Spice::Spice(Cube &cube, bool noTables) {
    init(*cube.label(), noTables);
  }

  void Spice::init(Pvl &lab, bool noTables) {
    NaifStatus::CheckErrors();

    // Initialize members
    m_solarLongitude = new Longitude;
    m_et = NULL;
    m_kernels = new QVector<QString>;
    m_target = new Target(this, lab);

    m_startTime = new iTime;
    m_endTime = new iTime;
    m_cacheSize = new SpiceDouble;
    *m_cacheSize = 0;

    m_startTimePadding = new SpiceDouble;
    *m_startTimePadding = 0;
    m_endTimePadding = new SpiceDouble;
    *m_endTimePadding = 0;

    m_instrumentPosition = NULL;
    m_instrumentRotation = NULL;
    m_sunPosition = NULL;
    m_bodyRotation = NULL;

    m_allowDownsizing = false;

    m_spkCode = new SpiceInt;
    m_ckCode = new SpiceInt;
    m_ikCode = new SpiceInt;
    m_sclkCode = new SpiceInt;
    m_spkBodyCode = new SpiceInt;
    m_bodyFrameCode = new SpiceInt;

    m_naifKeywords = new PvlObject("NaifKeywords");

    // m_sky = false;

    // Get the kernel group and load main kernels
    PvlGroup kernels = lab.findGroup("Kernels", Pvl::Traverse);

    // Get the time padding first
    if (kernels.hasKeyword("StartPadding")) {
      *m_startTimePadding = toDouble(kernels["StartPadding"][0]);
    }
    else {
      *m_startTimePadding = 0.0;
    }

    if (kernels.hasKeyword("EndPadding")) {
      *m_endTimePadding  = toDouble(kernels["EndPadding"][0]);
    }
    else {
      *m_endTimePadding = 0.0;
    }

    m_usingNaif = !lab.hasObject("NaifKeywords") || noTables;

    //  Modified  to load planetary ephemeris SPKs before s/c SPKs since some
    //  missions (e.g., MESSENGER) may augment the s/c SPK with new planet
    //  ephemerides. (2008-02-27 (KJB))
    if (m_usingNaif) {
      if (noTables) {
        load(kernels["TargetPosition"], noTables);
        load(kernels["InstrumentPosition"], noTables);
        load(kernels["InstrumentPointing"], noTables);
      }

      if (kernels.hasKeyword("Frame")) {
        load(kernels["Frame"], noTables);
      }

      load(kernels["TargetAttitudeShape"], noTables);
      if (kernels.hasKeyword("Instrument")) {
        load(kernels["Instrument"], noTables);
      }
      // Always load after instrument
      if (kernels.hasKeyword("InstrumentAddendum")) {
        load(kernels["InstrumentAddendum"], noTables);
      }
      load(kernels["LeapSecond"], noTables);
      if ( kernels.hasKeyword("SpacecraftClock")) {
        load(kernels["SpacecraftClock"], noTables);
      }

      // Modified to load extra kernels last to allow overriding default values
      // (2010-04-07) (DAC)
      if (kernels.hasKeyword("Extra")) {
        load(kernels["Extra"], noTables);
      }

      // If Target is Saturn and ShapeModel is RingPlane, load the extra rings pck file
      //  which changes the prime meridian values to report longitudes with respect to
      // the ascending node of the ringplane.
      if (m_target->name().toUpper() == "SATURN" && m_target->shape()->name().toUpper() == "PLANE") {
        PvlKeyword ringPck = PvlKeyword("RingPCK","$cassini/kernels/pck/saturnRings_v001.tpc");
        load(ringPck, noTables);
      }
    }
    else {
      *m_naifKeywords = lab.findObject("NaifKeywords");
    }

    // Get NAIF ik, spk, sclk, and ck codes
    //
    //    Use ikcode to get parameters from instrument kernel such as focal
    //    length, distortions, focal plane maps, etc
    //
    //    Use spkcode to get spacecraft position from spk file
    //
    //    Use sclkcode to transform times from et to tics
    //
    //    Use ckcode to transform between frames
    //
    //    Use bodycode to obtain radii and attitude (pole position/omega0)
    //
    //    Use spkbodycode to read body position from spk

    QString trykey = "NaifIkCode";
    if (kernels.hasKeyword("NaifFrameCode")) trykey = "NaifFrameCode";
    *m_ikCode = toInt(kernels[trykey][0]);

    *m_spkCode  = *m_ikCode / 1000;
    *m_sclkCode = *m_spkCode;
    *m_ckCode   = *m_ikCode;

    if (!m_target->isSky()) {
      QString radiiKey = "BODY" + toString(m_target->naifBodyCode()) + "_RADII";
      std::vector<Distance> radii(3,Distance());
      radii[0] = Distance(getDouble(radiiKey, 0), Distance::Kilometers);
      radii[1] = Distance(getDouble(radiiKey, 1), Distance::Kilometers);
      radii[2] = Distance(getDouble(radiiKey, 2), Distance::Kilometers);
      // m_target doesn't have the getDouble method so Spice gets the radii for it
      m_target->setRadii(radii);
    }
    *m_spkBodyCode = m_target->naifBodyCode();

    // Override them if they exist in the labels
    if (kernels.hasKeyword("NaifSpkCode")) {
      *m_spkCode = (int) kernels["NaifSpkCode"];
    }

    if (kernels.hasKeyword("NaifCkCode")) {
      *m_ckCode = (int) kernels["NaifCkCode"];
    }

    if (kernels.hasKeyword("NaifSclkCode")) {
      *m_sclkCode = (int) kernels["NaifSclkCode"];
    }

    if (!m_target->isSky()) {
      if (kernels.hasKeyword("NaifSpkBodyCode")) {
        *m_spkBodyCode = (int) kernels["NaifSpkBodyCode"];
      }
    }

    if (m_target->isSky()) {
      // Create the identity rotation for sky targets
      // Everything in bodyfixed will really be J2000
      m_bodyRotation = new SpiceRotation(1);
    }
    else {
      // JAA - Modified to stored and look for the frame body code in the
      // cube labels
      SpiceInt frameCode;
      if ((m_usingNaif) || (!m_naifKeywords->hasKeyword("BODY_FRAME_CODE"))) {
        char frameName[32];
        SpiceBoolean found;
        cidfrm_c(*m_spkBodyCode, sizeof(frameName), &frameCode, frameName,
                 &found);

        if (!found) {
          QString naifTarget = "IAU_" + QString(m_target->name()).toUpper();
          namfrm_c(naifTarget.toAscii().data(), &frameCode);
          if (frameCode == 0) {
            string msg = "Can not find NAIF BODY_FRAME_CODE for target [" 
                         + IString(m_target->name()) + "]";
            throw IException(IException::Io, msg, _FILEINFO_);
          }
        }

        QVariant result = (int)frameCode;
        storeValue("BODY_FRAME_CODE",0,SpiceIntType,result);
      }
      else {
        frameCode = getInteger("BODY_FRAME_CODE",0);
      }

      m_bodyRotation = new SpiceRotation(frameCode);
      *m_bodyFrameCode = frameCode;
    }

    m_instrumentRotation = new SpiceRotation(*m_ckCode);

    //  Set up for observer/target and light time correction to between s/c 
    // and target body.
    LightTimeCorrectionState ltState(*m_ikCode, this);
    ltState.checkSpkKernelsForAberrationCorrection();

    vector<Distance> radius = m_target->radii();
    Distance targetRadius((radius[0] + radius[2])/2.0);
    m_instrumentPosition = new SpacecraftPosition(*m_spkCode, *m_spkBodyCode,
                                                  ltState, targetRadius);

    m_sunPosition = new SpicePosition(10, m_target->naifBodyCode());

    // Check to see if we have nadir pointing that needs to be computed &
    // See if we have table blobs to load
    if (kernels["TargetPosition"][0].toUpper() == "TABLE") {
      Table t("SunPosition", lab.fileName(), lab);
      m_sunPosition->LoadCache(t);

      Table t2("BodyRotation", lab.fileName(), lab);
      m_bodyRotation->LoadCache(t2);
      if (t2.Label().hasKeyword("SolarLongitude")) {
        *m_solarLongitude = Longitude(t2.Label()["SolarLongitude"],
            Angle::Degrees);
      }
      else {
        solarLongitude();
      }
    }

    //  We can't assume InstrumentPointing & InstrumentPosition exist, old
    //  files may be around with the old keywords, SpacecraftPointing &
    //  SpacecraftPosition.  The old keywords were in existance before the
    //  Table option, so we don't need to check for Table under the old
    //  keywords.

    if (kernels["InstrumentPointing"].size() == 0) {
      throw IException(IException::Unknown,
                       "No camera pointing available",
                       _FILEINFO_);
    }

    //  2009-03-18  Tracie Sucharski - Removed test for old keywords, any files
    // with the old keywords should be re-run through spiceinit.
    if (kernels["InstrumentPointing"][0].toUpper() == "NADIR") {
      if (m_instrumentRotation) {
        delete m_instrumentRotation;
        m_instrumentRotation = NULL;
      }

      m_instrumentRotation = new SpiceRotation(*m_ikCode, *m_spkBodyCode);
    }
    else if (kernels["InstrumentPointing"][0].toUpper() == "TABLE") {
      Table t("InstrumentPointing", lab.fileName(), lab);
      m_instrumentRotation->LoadCache(t);
    }

    if (kernels["InstrumentPosition"].size() == 0) {
      throw IException(IException::Unknown,
                       "No instrument position available",
                       _FILEINFO_);
    }

    if (kernels["InstrumentPosition"][0].toUpper() == "TABLE") {
      Table t("InstrumentPosition", lab.fileName(), lab);
      m_instrumentPosition->LoadCache(t);
    }

    NaifStatus::CheckErrors();
  }


  void Spice::load(PvlKeyword &key, bool noTables) {
    NaifStatus::CheckErrors();

    for (int i = 0; i < key.size(); i++) {
      if (key[i] == "") continue;
      if (IString(key[i]).UpCase() == "NULL") break;
      if (IString(key[i]).UpCase() == "NADIR") break;
      if (IString(key[i]).UpCase() == "TABLE" && !noTables) break;
      if (IString(key[i]).UpCase() == "TABLE" && noTables) continue;
      FileName file(key[i]);
      if (!file.fileExists()) {
        QString msg = "Spice file does not exist [" + file.expanded() + "]";
        throw IException(IException::Io, msg, _FILEINFO_);
      }
      QString fileName(file.expanded());
      furnsh_c(fileName.toAscii().data());
      m_kernels->push_back((QString)key[i]);
    }

    NaifStatus::CheckErrors();
  }

  Spice::~Spice() {
    NaifStatus::CheckErrors();

    if (m_solarLongitude != NULL) {
      delete m_solarLongitude;
      m_solarLongitude = NULL;
    }

    if (m_et != NULL) {
      delete m_et;
      m_et = NULL;
    }

    if (m_startTime != NULL) {
      delete m_startTime;
      m_startTime = NULL;
    }

    if (m_endTime != NULL) {
      delete m_endTime;
      m_endTime = NULL;
    }

    if (m_cacheSize != NULL) {
      delete m_cacheSize;
      m_cacheSize = NULL;
    }

    if (m_startTimePadding != NULL) {
      delete m_startTimePadding;
      m_startTimePadding = NULL;
    }

    if (m_endTimePadding != NULL) {
      delete m_endTimePadding;
      m_endTimePadding = NULL;
    }

    if (m_instrumentPosition != NULL) {
      delete m_instrumentPosition;
      m_instrumentPosition = NULL;
    }

    if (m_instrumentRotation != NULL) {
      delete m_instrumentRotation;
      m_instrumentRotation = NULL;
    }

    if (m_sunPosition != NULL) {
      delete m_sunPosition;
      m_sunPosition = NULL;
    }

    if (m_bodyRotation != NULL) {
      delete m_bodyRotation;
      m_bodyRotation = NULL;
    }

    if (m_spkCode != NULL) {
      delete m_spkCode;
      m_spkCode = NULL;
    }

    if (m_ckCode != NULL) {
      delete m_ckCode;
      m_ckCode = NULL;
    }

    if (m_ikCode != NULL) {
      delete m_ikCode;
      m_ikCode = NULL;
    }

    if (m_sclkCode != NULL) {
      delete m_sclkCode;
      m_sclkCode = NULL;
    }

    if (m_spkBodyCode != NULL) {
      delete m_spkBodyCode;
      m_spkBodyCode = NULL;
    }

    if (m_bodyFrameCode != NULL) {
      delete m_bodyFrameCode;
      m_bodyFrameCode = NULL;
    }

    if (m_target != NULL) {
      delete m_target;
      m_target = NULL;
    }

    // Unload the kernels (TODO: Can this be done faster)
    for (int i = 0; m_kernels && i < m_kernels->size(); i++) {
      FileName file(m_kernels->at(i));
      QString fileName(file.expanded());
      unload_c(fileName.toAscii().data());
    }

    if (m_kernels != NULL) {
      delete m_kernels;
      m_kernels = NULL;
    }

    NaifStatus::CheckErrors();
  }

  void Spice::createCache(iTime startTime, iTime endTime,
      int cacheSize, double tol) {
    NaifStatus::CheckErrors();

    // Check for errors
    if (cacheSize <= 0) {
      string msg = "Argument cacheSize must be greater than zero";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    if (startTime > endTime) {
      string msg = "Argument startTime must be less than or equal to endTime";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    if (*m_cacheSize > 0) {
      string msg = "A cache has already been created";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    if (cacheSize == 1 && (*m_startTimePadding != 0 || *m_endTimePadding != 0)) {
      string msg = "This instrument does not support time padding";
      throw IException(IException::User, msg, _FILEINFO_);
    }

    string abcorr;
    if (getSpkAbCorrState(abcorr) ) {
      instrumentPosition()->SetAberrationCorrection("NONE");
    }

    iTime avgTime((startTime.Et() + endTime.Et()) / 2.0);
    computeSolarLongitude(avgTime);

    // Cache everything
    if (!m_bodyRotation->IsCached()) {
      int bodyRotationCacheSize = cacheSize;
      if (cacheSize > 2) bodyRotationCacheSize = 2;
      m_bodyRotation->LoadCache(
          startTime.Et() - *m_startTimePadding,
          endTime.Et() + *m_endTimePadding,
          bodyRotationCacheSize);
    }

    if (m_instrumentRotation->GetSource() < SpiceRotation::Memcache) {
      if (cacheSize > 3) m_instrumentRotation->MinimizeCache(SpiceRotation::Yes);
      m_instrumentRotation->LoadCache(
          startTime.Et() - *m_startTimePadding,
          endTime.Et() + *m_endTimePadding,
          cacheSize);
    }

    if (m_instrumentPosition->GetSource() < SpicePosition::Memcache) {
      m_instrumentPosition->LoadCache(
          startTime.Et() - *m_startTimePadding,
          endTime.Et() + *m_endTimePadding,
          cacheSize);
      if (cacheSize > 3) m_instrumentPosition->Memcache2HermiteCache(tol);
    }

    if (!m_sunPosition->IsCached()) {
      int sunPositionCacheSize = cacheSize;
      if (cacheSize > 2) sunPositionCacheSize = 2;
      m_sunPosition->LoadCache(
          startTime.Et() - *m_startTimePadding,
          endTime.Et() + *m_endTimePadding,
          sunPositionCacheSize);
    }

    // Save the time and cache size
    *m_startTime = startTime;
    *m_endTime = endTime;
    *m_cacheSize = cacheSize;
    m_et = NULL;

    // Unload the kernels (TODO: Can this be done faster)
    for (int i = 0; i < m_kernels->size(); i++) {
      FileName file(m_kernels->at(i));
      QString fileName(file.expanded());
      unload_c(fileName.toAscii().data());
    }

    m_kernels->clear();

    NaifStatus::CheckErrors();
  }


  iTime Spice::cacheStartTime() const {
    if (m_startTime) {
      return *m_startTime;
    }

    return iTime();
  }

  iTime Spice::cacheEndTime() const {
    if (m_endTime) {
      return *m_endTime;
    }

    return iTime();
  }

  void Spice::setTime(const iTime &et) {

    if (m_et == NULL) {
      m_et = new iTime();
      // Before the Spice is cached, but after the camera aberration correction
      // is set, check to see if the instrument position kernel was created
      // by spkwriter.  If so turn off aberration corrections because the camera
      // set aberration corrections are included in the spk already.
      string abcorr;
      if (*m_cacheSize == 0) {
        if (m_startTime->Et() == 0.0 && m_endTime->Et() == 0.0
            && getSpkAbCorrState(abcorr)) {
          instrumentPosition()->SetAberrationCorrection("NONE");
        }
      }
    }

    *m_et = et;

    m_bodyRotation->SetEphemerisTime(et.Et());
    m_instrumentRotation->SetEphemerisTime(et.Et());
    m_instrumentPosition->SetEphemerisTime(et.Et());
    m_sunPosition->SetEphemerisTime(et.Et());

    std::vector<double> uB = m_bodyRotation->ReferenceVector(m_sunPosition->Coordinate());
    m_uB[0] = uB[0];
    m_uB[1] = uB[1];
    m_uB[2] = uB[2];

    computeSolarLongitude(*m_et);
  }

  void Spice::instrumentPosition(double p[3]) const {
    instrumentBodyFixedPosition(p);
  }

  void Spice::instrumentBodyFixedPosition(double p[3]) const {
    if (m_et == NULL) {
      std::string msg = "You must call SetTime first";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    std::vector<double> sB = m_bodyRotation->ReferenceVector(m_instrumentPosition->Coordinate());
    p[0] = sB[0];
    p[1] = sB[1];
    p[2] = sB[2];
  }

  void Spice::instrumentBodyFixedVelocity(double v[3]) const {
    if (m_et == NULL) {
      std::string msg = "You must call SetTime first";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    std::vector<double> state;
    state.push_back(m_instrumentPosition->Coordinate()[0]);
    state.push_back(m_instrumentPosition->Coordinate()[1]);
    state.push_back(m_instrumentPosition->Coordinate()[2]);
    state.push_back(m_instrumentPosition->Velocity()[0]);
    state.push_back(m_instrumentPosition->Velocity()[1]);
    state.push_back(m_instrumentPosition->Velocity()[2]);

    std::vector<double> vB = m_bodyRotation->ReferenceVector(state);
    v[0] = vB[3];
    v[1] = vB[4];
    v[2] = vB[5];
  }

  iTime Spice::time() const {
    return *m_et;
  }

  void Spice::sunPosition(double p[3]) const {
    if (m_et == NULL) {
      std::string msg = "You must call SetTime first";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    p[0] = m_uB[0];
    p[1] = m_uB[1];
    p[2] = m_uB[2];
  }

  double Spice::targetCenterDistance() const {
    std::vector<double> sB = m_bodyRotation->ReferenceVector(m_instrumentPosition->Coordinate());
    return sqrt(pow(sB[0], 2) + pow(sB[1], 2) + pow(sB[2], 2));
  }

  void Spice::radii(Distance r[3]) const {
    for (int i = 0; i < 3; i++) 
       r[i] =m_target->radii()[i];
  }

  SpiceInt Spice::naifBodyCode() const {
    return (int) m_target->naifBodyCode();
  }

  SpiceInt Spice::naifSpkCode() const {
    return *m_spkCode;
  }

  SpiceInt Spice::naifCkCode() const {
    return *m_ckCode;
  }

  SpiceInt Spice::naifIkCode() const {
    return *m_ikCode;
  }

  SpiceInt Spice::naifSclkCode() const {
    return *m_sclkCode;
  }

  SpiceInt Spice::naifBodyFrameCode() const {
    return *m_bodyFrameCode;
  }


  PvlObject Spice::getStoredNaifKeywords() const {
    return *m_naifKeywords;
  }

  double Spice::resolution() {
    return 1.;
  };


  SpiceInt Spice::getInteger(const QString &key, int index) {
    return readValue(key, SpiceIntType, index).toInt();
  }

  SpiceDouble Spice::getDouble(const QString &key, int index) {
    return readValue(key, SpiceDoubleType, index).toDouble();
  }


  iTime Spice::getClockTime(QString clockValue, int sclkCode) {
    if (sclkCode == -1) {
      sclkCode = naifSclkCode();
    }

    iTime result;

    QString key = "CLOCK_ET_" + toString(sclkCode) + "_" + clockValue;
    QVariant storedClockTime = getStoredResult(key, SpiceDoubleType);

    if (storedClockTime.isNull()) {
      SpiceDouble timeOutput;
      scs2e_c(sclkCode, clockValue.toAscii().data(), &timeOutput);
      storedClockTime = timeOutput;
      storeResult(key, SpiceDoubleType, timeOutput);
    }

    result = storedClockTime.toDouble();

    return result;
  }


  QVariant Spice::readValue(QString key, SpiceValueType type, int index) {
    QVariant result;

    if (m_usingNaif) {
      NaifStatus::CheckErrors();

      // This is the success status of the naif call
      SpiceBoolean found = false;

      // Naif tells us how many values were read, but we always just read one.
      //   Use this variable to make naif happy.
      SpiceInt numValuesRead;

      if (type == SpiceDoubleType) {
        SpiceDouble kernelValue;
        gdpool_c(key.toAscii().data(), (SpiceInt)index, 1,
                 &numValuesRead, &kernelValue, &found);

        if (found)
          result = kernelValue;
      }
      else if (type == SpiceStringType) {
        char kernelValue[512];
        gcpool_c(key.toAscii().data(), (SpiceInt)index, 1, sizeof(kernelValue),
                 &numValuesRead, kernelValue, &found);

        if (found)
          result = kernelValue;
      }
      else if (type == SpiceIntType) {
        SpiceInt kernelValue;
        gipool_c(key.toAscii().data(), (SpiceInt)index, 1, &numValuesRead,
                 &kernelValue, &found);

        if (found)
          result = (int)kernelValue;
      }

      if (!found) {
        QString msg = "Can not find [" + key + "] in text kernels";
        throw IException(IException::Io, msg, _FILEINFO_);
      }

      storeValue(key, index, type, result);
    }
    else {
      // Read from PvlObject that is our naif keywords
      result = readStoredValue(key, type, index);

      if (result.isNull()) {
        IString msg = "The camera is requesting spice data [" + key + "] that "
            "was not attached, please re-run spiceinit";
        throw IException(IException::Unknown, msg, _FILEINFO_);
      }
    }

    return result;
  }


  void Spice::storeResult(QString name, SpiceValueType type, QVariant value) {
    if (type == SpiceDoubleType) {
      EndianSwapper swapper("LSB");

      double doubleVal = value.toDouble();
      doubleVal = swapper.Double(&doubleVal);
      QByteArray byteCode((char *) &doubleVal, sizeof(double));
      value = byteCode;
      type = SpiceByteCodeType;
    }

    storeValue(name + "_COMPUTED", 0, type, value);
  }


  QVariant Spice::getStoredResult(QString name, SpiceValueType type) {
    bool wasDouble = false;

    if (type == SpiceDoubleType) {
      wasDouble = true;
      type = SpiceByteCodeType;
    }

    QVariant stored = readStoredValue(name + "_COMPUTED", type, 0);

    if (wasDouble && !stored.isNull()) {
      EndianSwapper swapper("LSB");
      double doubleVal = swapper.Double((void *)QByteArray::fromHex(
          stored.toByteArray()).data());
      stored = doubleVal;
    }

    return stored;
  }


  void Spice::storeValue(QString key, int index, SpiceValueType type,
                         QVariant value) {
    if (!m_naifKeywords->hasKeyword(key)) {
      m_naifKeywords->addKeyword(PvlKeyword(key));
    }

    PvlKeyword &storedKey = m_naifKeywords->findKeyword(key);

    while(index >= storedKey.size()) {
      storedKey.addValue("");
    }

    if (type == SpiceByteCodeType) {
      storedKey[index] = QString(value.toByteArray().toHex().data());
    }
    else if (type == SpiceStringType) {
      storedKey[index] = value.toString();
    }
    else if (type == SpiceDoubleType) {
      storedKey[index] = toString(value.toDouble());
    }
    else if (type == SpiceIntType) {
      storedKey[index] = toString(value.toInt());
    }
    else {
      IString msg = "Unable to store variant in labels for key [" + key + "]";
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
  }


  QVariant Spice::readStoredValue(QString key, SpiceValueType type,
                                  int index) {
    // Read from PvlObject that is our naif keywords
    QVariant result;

    if (m_naifKeywords->hasKeyword(key) && !m_usingNaif) {
      PvlKeyword &storedKeyword = m_naifKeywords->findKeyword(key);

      try {
        if (type == SpiceDoubleType) {
          result = toDouble(storedKeyword[index]);
        }
        else if (type == SpiceStringType) {
          result = storedKeyword[index];
        }
        else if (type == SpiceByteCodeType || SpiceStringType) {
          result = storedKeyword[index].toAscii();
        }
        else if (type == SpiceIntType) {
          result = toInt(storedKeyword[index]);
        }
      }
      catch(IException &) {
      }
    }

    return result;
  }

  QString Spice::getString(const QString &key, int index) {
    return readValue(key, SpiceStringType, index).toString();
  }


  void Spice::subSpacecraftPoint(double &lat, double &lon) {
    NaifStatus::CheckErrors();

    if (m_et == NULL) {
      std::string msg = "You must call SetTime first";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    SpiceDouble usB[3], dist;
    std::vector<double> vsB = m_bodyRotation->ReferenceVector(m_instrumentPosition->Coordinate());
    SpiceDouble sB[3];
    sB[0] = vsB[0];
    sB[1] = vsB[1];
    sB[2] = vsB[2];
    unorm_c(sB, usB, &dist);

    std::vector<Distance> radii = target()->radii();
    SpiceDouble a = radii[0].kilometers();
    SpiceDouble b = radii[1].kilometers();
    SpiceDouble c = radii[2].kilometers();

    SpiceDouble originB[3];
    originB[0] = originB[1] = originB[2] = 0.0;

    SpiceBoolean found;
    SpiceDouble subB[3];
    surfpt_c(originB, usB, a, b, c, subB, &found);

    SpiceDouble mylon, mylat;
    reclat_c(subB, &a, &mylon, &mylat);
    lat = mylat * 180.0 / PI;
    lon = mylon * 180.0 / PI;
    if (lon < 0.0) lon += 360.0;

    NaifStatus::CheckErrors();
  }

  void Spice::subSolarPoint(double &lat, double &lon) {
    NaifStatus::CheckErrors();

    if (m_et == NULL) {
      std::string msg = "You must call SetTime first";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    SpiceDouble uuB[3], dist;
    unorm_c(m_uB, uuB, &dist);

    std::vector<Distance> radii = target()->radii();
    SpiceDouble a = radii[0].kilometers();
    SpiceDouble b = radii[1].kilometers();
    SpiceDouble c = radii[2].kilometers();

    SpiceDouble originB[3];
    originB[0] = originB[1] = originB[2] = 0.0;

    SpiceBoolean found;
    SpiceDouble subB[3];
    surfpt_c(originB, uuB, a, b, c, subB, &found);

    SpiceDouble mylon, mylat;
    reclat_c(subB, &a, &mylon, &mylat);
    lat = mylat * 180.0 / PI;
    lon = mylon * 180.0 / PI;
    if (lon < 0.0) lon += 360.0;

    NaifStatus::CheckErrors();
  }


  Target *Spice::target() const {
    return m_target;
  }


  QString Spice::targetName() const {
    return m_target->name();
  }


  void Spice::computeSolarLongitude(iTime et) {
    NaifStatus::CheckErrors();

    if (m_target->isSky()) {
      *m_solarLongitude = Longitude();
      return;
    }

    if (m_bodyRotation->IsCached()) return;

    double tipm[3][3], npole[3];
    char frameName[32];
    SpiceInt frameCode;
    SpiceBoolean found;

    cidfrm_c(*m_spkBodyCode, sizeof(frameName), &frameCode, frameName, &found);

    if (found) {
      pxform_c("J2000", frameName, et.Et(), tipm);
    }
    else {
      tipbod_c("J2000", *m_spkBodyCode, et.Et(), tipm);
    }

    for (int i = 0; i < 3; i++) {
      npole[i] = tipm[2][i];
    }

    double state[6], lt;
    spkez_c(*m_spkBodyCode, et.Et(), "J2000", "NONE", 10, state, &lt);

    double uavel[3];
    ucrss_c(state, &state[3], uavel);

    double x[3], y[3], z[3];
    vequ_c(uavel, z);
    ucrss_c(npole, z, x);
    ucrss_c(z, x, y);

    double trans[3][3];
    for (int i = 0; i < 3; i++) {
      trans[0][i] = x[i];
      trans[1][i] = y[i];
      trans[2][i] = z[i];
    }

    spkez_c(10, et.Et(), "J2000", "LT+S", *m_spkBodyCode, state, &lt);

    double pos[3];
    mxv_c(trans, state, pos);

    double radius, ls, lat;
    reclat_c(pos, &radius, &ls, &lat);

    *m_solarLongitude = Longitude(ls, Angle::Radians).force360Domain();

    NaifStatus::CheckErrors();
  }


  Longitude Spice::solarLongitude() {
    if (m_et) {
      computeSolarLongitude(*m_et);
      return *m_solarLongitude;
    }

    return Longitude();
  }


  bool Spice::hasKernels(Pvl &lab) {

    // Get the kernel group and check main kernels
    PvlGroup kernels = lab.findGroup("Kernels", Pvl::Traverse);
    std::vector<string> keywords;
    keywords.push_back("TargetPosition");

    if (kernels.hasKeyword("SpacecraftPosition")) {
      keywords.push_back("SpacecraftPosition");
    }
    else {
      keywords.push_back("InstrumentPosition");
    }

    if (kernels.hasKeyword("SpacecraftPointing")) {
      keywords.push_back("SpacecraftPointing");
    }
    else {
      keywords.push_back("InstrumentPointing");
    }

    if (kernels.hasKeyword("Frame")) {
      keywords.push_back("Frame");
    }

    if (kernels.hasKeyword("Extra")) {
      keywords.push_back("Extra");
    }

    PvlKeyword key;
    for (int ikey = 0; ikey < (int) keywords.size(); ikey++) {
      key = kernels[ikey];

      for (int i = 0; i < key.size(); i++) {
        if (key[i] == "") return false;
        if (IString(key[i]).UpCase() == "NULL") return false;
        if (IString(key[i]).UpCase() == "NADIR") return false;
        if (IString(key[i]).UpCase() == "TABLE") return false;
      }
    }
    return true;
  }

  SpicePosition *Spice::sunPosition() const {
    return m_sunPosition;
  }

  SpicePosition *Spice::instrumentPosition() const {
    return m_instrumentPosition;
  }

  SpiceRotation *Spice::bodyRotation() const {
    return m_bodyRotation;
  }

  SpiceRotation *Spice::instrumentRotation() const {
    return m_instrumentRotation;
  }

}