diff --git a/include/usgscsm/UsgsAstroLsSensorModel.h b/include/usgscsm/UsgsAstroLsSensorModel.h
index 1bb18239582e79e5f08ad5b749dd69abf1c862df..78bb3633fe9ac28c8e037fc9a3dba63486b05341 100644
--- a/include/usgscsm/UsgsAstroLsSensorModel.h
+++ b/include/usgscsm/UsgsAstroLsSensorModel.h
@@ -32,6 +32,9 @@
 #include <CorrelationModel.h>
 #include "Distortion.h"
 
+#include "spdlog/spdlog.h"
+#include "spdlog/sinks/basic_file_sink.h"
+
 class UsgsAstroLsSensorModel : public csm::RasterGM, virtual public csm::SettableEllipsoid
 {
 public:
@@ -115,6 +118,10 @@ public:
    std::vector<double> m_covariance;
    int          m_imageFlipFlag;
 
+   // Define logging pointer and file content
+   std::string m_logFile;
+   std::shared_ptr<spdlog::logger> m_logger;
+
    // Hardcoded
    static const std::string      _SENSOR_MODEL_NAME; // state date element 0
 
@@ -683,6 +690,9 @@ public:
    //  reference by the given index.
    //<
 
+   virtual std::shared_ptr<spdlog::logger> getLogger();
+   virtual void setLogger(std::shared_ptr<spdlog::logger> logger);
+
 
    //---
    // Uncertainty Propagation
diff --git a/src/UsgsAstroFrameSensorModel.cpp b/src/UsgsAstroFrameSensorModel.cpp
index 86562d2d1e269af57eb57e09730a92acda977e75..d48276e5a5cc633b7257c4a1d2fc4a0504d48c6d 100644
--- a/src/UsgsAstroFrameSensorModel.cpp
+++ b/src/UsgsAstroFrameSensorModel.cpp
@@ -97,7 +97,7 @@ csm::ImageCoord UsgsAstroFrameSensorModel::groundToImage(const csm::EcefCoord &g
                               double desiredPrecision,
                               double *achievedPrecision,
                               csm::WarningList *warnings) const {
-  MESSAGE_LOG(this->m_logger, "Computeing groundToImage(No adjustments) for {}, {}, {}, with desired precision {}",
+  MESSAGE_LOG(this->m_logger, "Computing groundToImage(No adjustments) for {}, {}, {}, with desired precision {}",
               groundPt.x, groundPt.y, groundPt.z, desiredPrecision);
 
   return groundToImage(groundPt, m_noAdjustments, desiredPrecision, achievedPrecision, warnings);
@@ -121,7 +121,7 @@ csm::ImageCoord UsgsAstroFrameSensorModel::groundToImage(
     double*                    achieved_precision,
     csm::WarningList*          warnings ) const {
 
-  MESSAGE_LOG(this->m_logger, "Computeing groundToImage for {}, {}, {}, with desired precision {}",
+  MESSAGE_LOG(this->m_logger, "Computing groundToImage for {}, {}, {}, with desired precision {}",
               groundPt.x, groundPt.y, groundPt.z, desired_precision);
 
   double x = groundPt.x;
@@ -171,6 +171,7 @@ csm::ImageCoordCovar UsgsAstroFrameSensorModel::groundToImage(const csm::EcefCoo
                                    csm::WarningList *warnings) const {
     MESSAGE_LOG(this->m_logger, "Computeing groundToImage(Covar) for {}, {}, {}, with desired precision {}",
                 groundPt.x, groundPt.y, groundPt.z, desiredPrecision);
+
     csm::EcefCoord gp;
     gp.x = groundPt.x;
     gp.y = groundPt.y;
diff --git a/src/UsgsAstroLsSensorModel.cpp b/src/UsgsAstroLsSensorModel.cpp
index 282843f22717c78e701cd3a253d6c7a48cbd7548..44986ffe0f2735d4d13be77730feb0adc579d357 100644
--- a/src/UsgsAstroLsSensorModel.cpp
+++ b/src/UsgsAstroLsSensorModel.cpp
@@ -28,7 +28,11 @@
 #include <sstream>
 #include <Error.h>
 #include <json/json.hpp>
+
+#define MESSAGE_LOG(logger, ...) if (logger) { logger->info(__VA_ARGS__); }
+
 using json = nlohmann::json;
+using namespace std;
 
 const std::string UsgsAstroLsSensorModel::_SENSOR_MODEL_NAME
          = "USGS_ASTRO_LINE_SCANNER_SENSOR_MODEL";
@@ -121,9 +125,13 @@ const csm::param::Type
    csm::param::FIXED
 };
 
-
+//***************************************************************************
+// UsgsAstroLineScannerSensorModel::replaceModelState
+//***************************************************************************
 void UsgsAstroLsSensorModel::replaceModelState(const std::string &stateString )
 {
+   MESSAGE_LOG(m_logger, "Replacing model state")
+
    reset();
    auto j = json::parse(stateString);
    int num_params    = NUM_PARAMETERS;
@@ -134,6 +142,16 @@ void UsgsAstroLsSensorModel::replaceModelState(const std::string &stateString )
    m_nLines = j["m_nLines"];
    m_nSamples = j["m_nSamples"];
    m_platformFlag = j["m_platformFlag"];
+   MESSAGE_LOG(m_logger, "m_imageIdentifier: {} "
+                         "m_sensorName: {} "
+                         "m_nLines: {} "
+                         "m_nSamples: {} "
+                         "m_platformFlag: {} ",
+                         j["m_imageIdentifier"].dump(),
+                         j["m_sensorName"].dump(),
+                         j["m_nLines"].dump(),
+                         j["m_nSamples"].dump(),
+                         j["m_platformFlag"].dump())
 
    m_intTimeLines = j["m_intTimeLines"].get<std::vector<double>>();
    m_intTimeStartTimes = j["m_intTimeStartTimes"].get<std::vector<double>>();
@@ -144,38 +162,91 @@ void UsgsAstroLsSensorModel::replaceModelState(const std::string &stateString )
    m_detectorSampleSumming = j["m_detectorSampleSumming"];
    m_startingSample = j["m_startingDetectorSample"];
    m_ikCode = j["m_ikCode"];
-
+   MESSAGE_LOG(m_logger, "m_startingEphemerisTime: {} "
+                         "m_centerEphemerisTime: {} "
+                         "m_detectorSampleSumming: {} "
+                         "m_startingSample: {} "
+                         "m_ikCode: {} ",
+                         j["m_startingEphemerisTime"].dump(),
+                         j["m_centerEphemerisTime"].dump(),
+                         j["m_detectorSampleSumming"].dump(),
+                         j["m_startingSample"].dump(), j["m_ikCode"].dump())
 
    m_focalLength = j["m_focalLength"];
    m_zDirection = j["m_zDirection"];
    m_distortionType = (DistortionType)j["m_distortionType"].get<int>();
    m_opticalDistCoeffs = j["m_opticalDistCoeffs"].get<std::vector<double>>();
+   MESSAGE_LOG(m_logger, "m_focalLength: {} "
+                         "m_zDirection: {} "
+                         "m_distortionType: {} ",
+                         j["m_focalLength"].dump(), j["m_zDirection"].dump(),
+                         j["m_distortionType"].dump())
+
    for (int i = 0; i < 3; i++) {
      m_iTransS[i] = j["m_iTransS"][i];
      m_iTransL[i] = j["m_iTransL"][i];
    }
+
    m_detectorSampleOrigin = j["m_detectorSampleOrigin"];
    m_detectorLineOrigin = j["m_detectorLineOrigin"];
    m_majorAxis = j["m_majorAxis"];
    m_minorAxis = j["m_minorAxis"];
+   MESSAGE_LOG(m_logger, "m_detectorSampleOrigin: {} "
+                         "m_detectorLineOrigin: {} "
+                         "m_majorAxis: {} "
+                         "m_minorAxis: {} ",
+                         j["m_detectorSampleOrigin"].dump(),
+                         j["m_detectorLineOrigin"].dump(),
+                         j["m_majorAxis"].dump(), j["m_minorAxis"].dump())
+
    m_platformIdentifier = j["m_platformIdentifier"];
    m_sensorIdentifier = j["m_sensorIdentifier"];
+   MESSAGE_LOG(m_logger, "m_platformIdentifier: {} "
+                         "m_sensorIdentifier: {} ",
+                         j["m_platformIdentifier"].dump(),
+                         j["m_sensorIdentifier"].dump())
+
    m_minElevation = j["m_minElevation"];
    m_maxElevation = j["m_maxElevation"];
+   MESSAGE_LOG(m_logger, "m_minElevation: {} "
+                         "m_maxElevation: {} ",
+                         j["m_minElevation"].dump(),
+                         j["m_maxElevation"].dump())
+
    m_dtEphem = j["m_dtEphem"];
    m_t0Ephem = j["m_t0Ephem"];
    m_dtQuat = j["m_dtQuat"];
    m_t0Quat = j["m_t0Quat"];
    m_numPositions = j["m_numPositions"];
+   MESSAGE_LOG(m_logger, "m_dtEphem: {} "
+                         "m_t0Ephem: {} "
+                         "m_dtQuat: {} "
+                         "m_t0Quat: {} ",
+                         j["m_dtEphem"].dump(), j["m_t0Ephem"].dump(),
+                         j["m_dtQuat"].dump(), j["m_t0Quat"].dump())
 
    m_numQuaternions = j["m_numQuaternions"];
    m_referencePointXyz.x = j["m_referencePointXyz"][0];
    m_referencePointXyz.y = j["m_referencePointXyz"][1];
    m_referencePointXyz.z = j["m_referencePointXyz"][2];
+   MESSAGE_LOG(m_logger, "m_numQuaternions: {} "
+                         "m_referencePointX: {} "
+                         "m_referencePointY: {} "
+                         "m_referencePointZ: {} ",
+                         j["m_numQuaternions"].dump(), j["m_referencePointXyz"][0].dump(),
+                         j["m_referencePointXyz"][1].dump(),
+                         j["m_referencePointXyz"][2].dump())
+
    m_gsd = j["m_gsd"];
    m_flyingHeight = j["m_flyingHeight"];
    m_halfSwath = j["m_halfSwath"];
    m_halfTime = j["m_halfTime"];
+   MESSAGE_LOG(m_logger, "m_gsd: {} "
+                         "m_flyingHeight: {} "
+                         "m_halfSwath: {} "
+                         "m_halfTime: {} ",
+                         j["m_gsd"].dump(), j["m_flyingHeight"].dump(),
+                         j["m_halfSwath"].dump(), j["m_halfTime"].dump())
    // Vector = is overloaded so explicit get with type required.
 
    m_positions = j["m_positions"].get<std::vector<double>>();
@@ -184,6 +255,17 @@ void UsgsAstroLsSensorModel::replaceModelState(const std::string &stateString )
    m_currentParameterValue = j["m_currentParameterValue"].get<std::vector<double>>();
    m_covariance = j["m_covariance"].get<std::vector<double>>();
 
+   m_logFile = j["m_logFile"].get<std::string>();
+   if (m_logFile.empty()) {
+     m_logger.reset();
+   }
+   else {
+     m_logger = spdlog::get(m_logFile);
+     if (!m_logger) {
+       m_logger = spdlog::basic_logger_mt(m_logFile, m_logFile);
+     }
+   }
+
    for (int i = 0; i < num_params; i++) {
      for (int k = 0; k < NUM_PARAM_TYPES; k++) {
        if (j["m_parameterType"][i] == PARAM_STRING_ALL[k]) {
@@ -209,6 +291,9 @@ void UsgsAstroLsSensorModel::replaceModelState(const std::string &stateString )
    }
 }
 
+//***************************************************************************
+// UsgsAstroLineScannerSensorModel::getModelNameFromModelState
+//***************************************************************************
 std::string UsgsAstroLsSensorModel::getModelNameFromModelState(
    const std::string& model_state)
 {
@@ -236,7 +321,12 @@ std::string UsgsAstroLsSensorModel::getModelNameFromModelState(
    return model_name;
 }
 
+//***************************************************************************
+// UsgsAstroLineScannerSensorModel::getModelState
+//***************************************************************************
 std::string UsgsAstroLsSensorModel::getModelState() const {
+      MESSAGE_LOG(m_logger, "Running getModelState")
+
       json state;
       state["m_modelName"] = _SENSOR_MODEL_NAME;
       state["m_startingDetectorSample"] = m_startingSample;
@@ -245,55 +335,124 @@ std::string UsgsAstroLsSensorModel::getModelState() const {
       state["m_nLines"] = m_nLines;
       state["m_nSamples"] = m_nSamples;
       state["m_platformFlag"] = m_platformFlag;
+      MESSAGE_LOG(m_logger, "m_imageIdentifier: {} "
+                                  "m_sensorName: {} "
+                                  "m_nLines: {} "
+                                  "m_nSamples: {} "
+                                  "m_platformFlag: {} ",
+                                  m_imageIdentifier, m_sensorName,
+                                  m_nLines, m_nSamples, m_platformFlag)
+
       state["m_intTimeLines"] = m_intTimeLines;
       state["m_intTimeStartTimes"] = m_intTimeStartTimes;
       state["m_intTimes"] = m_intTimes;
       state["m_startingEphemerisTime"] = m_startingEphemerisTime;
       state["m_centerEphemerisTime"] = m_centerEphemerisTime;
+      MESSAGE_LOG(m_logger, "m_startingEphemerisTime: {} "
+                                  "m_centerEphemerisTime: {} ",
+                                  m_startingEphemerisTime,
+                                  m_centerEphemerisTime)
+
       state["m_detectorSampleSumming"] = m_detectorSampleSumming;
       state["m_startingSample"] = m_startingSample;
       state["m_ikCode"] = m_ikCode;
+      MESSAGE_LOG(m_logger, "m_detectorSampleSumming: {} "
+                                  "m_startingSample: {} "
+                                  "m_ikCode: {} ",
+                                  m_detectorSampleSumming, m_startingSample,
+                                  m_ikCode)
+
       state["m_focalLength"] = m_focalLength;
       state["m_zDirection"] = m_zDirection;
       state["m_distortionType"] = m_distortionType;
       state["m_opticalDistCoeffs"] = m_opticalDistCoeffs;
+      MESSAGE_LOG(m_logger, "m_focalLength: {} "
+                                  "m_zDirection: {} "
+                                  "m_distortionType (0-Radial, 1-Transverse): {} ",
+                                  m_focalLength, m_zDirection,
+                                  m_distortionType)
+
       state["m_iTransS"] = std::vector<double>(m_iTransS, m_iTransS+3);
       state["m_iTransL"] = std::vector<double>(m_iTransL, m_iTransL+3);
+
       state["m_detectorSampleOrigin"] = m_detectorSampleOrigin;
       state["m_detectorLineOrigin"] = m_detectorLineOrigin;
       state["m_majorAxis"] = m_majorAxis;
       state["m_minorAxis"] = m_minorAxis;
+      MESSAGE_LOG(m_logger, "m_detectorSampleOrigin: {} "
+                                  "m_detectorLineOrigin: {} "
+                                  "m_majorAxis: {} "
+                                  "m_minorAxis: {} ",
+                                  m_detectorSampleOrigin, m_detectorLineOrigin,
+                                  m_majorAxis, m_minorAxis)
+
       state["m_platformIdentifier"] = m_platformIdentifier;
       state["m_sensorIdentifier"] = m_sensorIdentifier;
       state["m_minElevation"] = m_minElevation;
       state["m_maxElevation"] = m_maxElevation;
+      MESSAGE_LOG(m_logger, "m_platformIdentifier: {} "
+                                  "m_sensorIdentifier: {} "
+                                  "m_minElevation: {} "
+                                  "m_maxElevation: {} ",
+                                  m_platformIdentifier, m_sensorIdentifier,
+                                  m_minElevation, m_maxElevation)
+
       state["m_dtEphem"] = m_dtEphem;
       state["m_t0Ephem"] = m_t0Ephem;
       state["m_dtQuat"] = m_dtQuat;
       state["m_t0Quat"] = m_t0Quat;
+      MESSAGE_LOG(m_logger, "m_dtEphem: {} "
+                                  "m_t0Ephem: {} "
+                                  "m_dtQuat: {} "
+                                  "m_t0Quat: {} ",
+                                  m_dtEphem, m_t0Ephem,
+                                  m_dtQuat, m_t0Quat)
+
       state["m_numPositions"] = m_numPositions;
       state["m_numQuaternions"] = m_numQuaternions;
       state["m_positions"] = m_positions;
       state["m_velocities"] = m_velocities;
       state["m_quaternions"] = m_quaternions;
+      MESSAGE_LOG(m_logger, "m_numPositions: {} "
+                            "m_numQuaternions: {} ",
+                             m_numPositions, m_numQuaternions)
+
       state["m_currentParameterValue"] = m_currentParameterValue;
       state["m_parameterType"] = m_parameterType;
+
       state["m_gsd"] = m_gsd;
       state["m_flyingHeight"] = m_flyingHeight;
       state["m_halfSwath"] = m_halfSwath;
       state["m_halfTime"] = m_halfTime;
       state["m_covariance"] = m_covariance;
+      MESSAGE_LOG(m_logger, "m_gsd: {} "
+                            "m_flyingHeight: {} "
+                            "m_halfSwath: {} "
+                            "m_halfTime: {} ",
+                            m_gsd, m_flyingHeight,
+                            m_halfSwath, m_halfTime)
 
       state["m_referencePointXyz"] = json();
       state["m_referencePointXyz"][0] = m_referencePointXyz.x;
       state["m_referencePointXyz"][1] = m_referencePointXyz.y;
       state["m_referencePointXyz"][2] = m_referencePointXyz.z;
+      MESSAGE_LOG(m_logger, "m_referencePointXyz: {} "
+                            "m_referencePointXyz: {} "
+                            "m_referencePointXyz: {} ",
+                             m_referencePointXyz.x, m_referencePointXyz.y,
+                             m_referencePointXyz.z)
+
+      state["m_logFile"] = m_logFile;
 
       return state.dump();
  }
 
+ //***************************************************************************
+ // UsgsAstroLineScannerSensorModel::reset
+ //***************************************************************************
 void UsgsAstroLsSensorModel::reset()
 {
+  MESSAGE_LOG(m_logger, "Running reset()")
   _linear = false; // default until a linear model is made
   _u0    = 0.0;
   _du_dx = 0.0;
@@ -359,8 +518,10 @@ void UsgsAstroLsSensorModel::reset()
   m_halfTime = 10.0;
 
   m_covariance = std::vector<double>(NUM_PARAMETERS * NUM_PARAMETERS,0.0); // 52
-}
 
+  m_logFile = "";
+  m_logger.reset();
+}
 
 //*****************************************************************************
 // UsgsAstroLsSensorModel Constructor
@@ -385,13 +546,20 @@ void UsgsAstroLsSensorModel::updateState()
 {
    // If sensor model is being created for the first time
    // This routine will set some parameters not found in the ISD.
-
+   MESSAGE_LOG(m_logger, "Updating State")
    // Reference point (image center)
    double lineCtr = m_nLines / 2.0;
    double sampCtr = m_nSamples / 2.0;
    csm::ImageCoord ip(lineCtr, sampCtr);
+   MESSAGE_LOG(m_logger, "updateState: center image coordinate set to {} {}",
+                               lineCtr, sampCtr)
+
    double refHeight = 0;
    m_referencePointXyz = imageToGround(ip, refHeight);
+   MESSAGE_LOG(m_logger, "updateState: reference point (x, y, z) {} {} {}",
+                                m_referencePointXyz.x, m_referencePointXyz.y,
+                                m_referencePointXyz.z)
+
    // Compute ground sample distance
    ip.line += 1;
    ip.samp += 1;
@@ -400,6 +568,9 @@ void UsgsAstroLsSensorModel::updateState()
    double dy = delta.y - m_referencePointXyz.y;
    double dz = delta.z - m_referencePointXyz.z;
    m_gsd = sqrt((dx*dx + dy*dy + dz*dz) / 2.0);
+   MESSAGE_LOG(m_logger, "updateState: ground sample distance set to {} "
+                              "based on dx {} dy {} dz {}",
+                               m_gsd, dx, dy, dz)
 
    // Compute flying height
    csm::EcefCoord sensorPos = getSensorPosition(0.0);
@@ -407,14 +578,21 @@ void UsgsAstroLsSensorModel::updateState()
    dy = sensorPos.y - m_referencePointXyz.y;
    dz = sensorPos.z - m_referencePointXyz.z;
    m_flyingHeight = sqrt(dx*dx + dy*dy + dz*dz);
+   MESSAGE_LOG(m_logger, "updateState: flight height set to {}"
+                              "based on dx {} dy {} dz {}",
+                               m_flyingHeight, dx, dy, dz)
 
    // Compute half swath
    m_halfSwath = m_gsd * m_nSamples / 2.0;
+   MESSAGE_LOG(m_logger, "updateState: half swath set to {}",
+                               m_halfSwath)
 
    // Compute half time duration
    double fullImageTime = m_intTimeStartTimes.back() - m_intTimeStartTimes.front()
                           + m_intTimes.back() * (m_nLines - m_intTimeLines.back());
    m_halfTime = fullImageTime / 2.0;
+   MESSAGE_LOG(m_logger, "updateState: half time duration set to {}",
+                               m_halfTime)
 
    // Parameter covariance, hardcoded accuracy values
    int num_params = NUM_PARAMETERS;
@@ -444,6 +622,10 @@ csm::ImageCoord UsgsAstroLsSensorModel::groundToImage(
    double*               achieved_precision,
    csm::WarningList*     warnings) const
 {
+
+  MESSAGE_LOG(m_logger, "Computing groundToImage(No adjustments) for {}, {}, {}, with desired precision {}",
+              ground_pt.x, ground_pt.y, ground_pt.z, desired_precision);
+
    // The public interface invokes the private interface with no adjustments.
    return groundToImage(
       ground_pt, _no_adjustment,
@@ -460,6 +642,8 @@ csm::ImageCoord UsgsAstroLsSensorModel::groundToImage(
    double*               achieved_precision,
    csm::WarningList*     warnings) const
 {
+   MESSAGE_LOG(m_logger, "Computing groundToImage (with adjustments) for {}, {}, {}, with desired precision {}",
+              ground_pt.x, ground_pt.y, ground_pt.z, desired_precision);
    // Search for the line, sample coordinate that viewed a given ground point.
    // This method uses an iterative secant method to search for the image
    // line.
@@ -587,6 +771,8 @@ csm::ImageCoordCovar UsgsAstroLsSensorModel::groundToImage(
    double* achieved_precision,
    csm::WarningList* warnings) const
 {
+  MESSAGE_LOG(m_logger, "Computing groundToImage(Covar) for {}, {}, {}, with desired precision {}",
+                groundPt.x, groundPt.y, groundPt.z, desired_precision);
    // Ground to image with error propagation
    // Compute corresponding image point
    csm::EcefCoord gp;
@@ -658,6 +844,9 @@ csm::EcefCoord UsgsAstroLsSensorModel::imageToGround(
    double* achieved_precision,
    csm::WarningList* warnings) const
 {
+   MESSAGE_LOG(m_logger, "Computing imageToGround for {}, {}, {}, with desired precision {}",
+               image_pt.line, image_pt.samp, height, desired_precision);
+
    double xc, yc, zc;
    double vx, vy, vz;
    double xl, yl, zl;
@@ -683,14 +872,25 @@ csm::EcefCoord UsgsAstroLsSensorModel::imageToGround(
             "UsgsAstroLsSensorModel::imageToGround()"));
    }
 
+/*
+    MESSAGE_LOG(m_logger, "imageToGround for {} {} {} achieved precision {}",
+                                image_pt.line, image_pt.samp, height, achieved_precision)
+*/
+
    return csm::EcefCoord(x, y, z);
 }
 
-
+//***************************************************************************
+// UsgsAstroLineScannerSensorModel::determineSensorCovarianceInImageSpace
+//***************************************************************************
 void UsgsAstroLsSensorModel::determineSensorCovarianceInImageSpace(
    csm::EcefCoord &gp,
    double         sensor_cov[4] ) const
 {
+   MESSAGE_LOG(m_logger, "Calculating determineSensorCovarianceInImageSpace for {} {} {}",
+                              gp.x, gp.y, gp.z)
+
+
    int i, j, totalAdjParams;
    totalAdjParams = getNumParameters();
 
@@ -725,6 +925,8 @@ csm::EcefCoordCovar UsgsAstroLsSensorModel::imageToGround(
    double* achieved_precision,
    csm::WarningList* warnings) const
 {
+  MESSAGE_LOG(m_logger, "Calculating imageToGround (with error propagation) for {}, {}, {} with height varinace {} and desired precision {}",
+                              image_pt.line, image_pt.samp, height, heightVariance, desired_precision)
    // Image to ground with error propagation
    // Use numerical partials
 
@@ -810,6 +1012,10 @@ csm::EcefLocus UsgsAstroLsSensorModel::imageToProximateImagingLocus(
    double* achieved_precision,
    csm::WarningList* warnings) const
 {
+
+  MESSAGE_LOG(m_logger, "Computing imageToProximateImagingLocus (ground {}, {}, {}) for image point {}, {} with desired precision {}",
+                              ground_pt.x, ground_pt.y, ground_pt.z, image_pt.line, image_pt.samp, desired_precision);
+
    // Object ray unit direction near given ground location
    const double DELTA_GROUND = m_gsd;
 
@@ -871,6 +1077,10 @@ csm::EcefLocus UsgsAstroLsSensorModel::imageToRemoteImagingLocus(
    double* achieved_precision,
    csm::WarningList* warnings) const
 {
+
+  MESSAGE_LOG(m_logger, "Calculating imageToRemoteImagingLocus for point {}, {} with desired precision {}",
+                              image_pt.line, image_pt.samp, desired_precision)
+
   double vx, vy, vz;
   csm::EcefLocus locus;
   losToEcf(
@@ -895,6 +1105,10 @@ csm::EcefLocus UsgsAstroLsSensorModel::imageToRemoteImagingLocus(
 std::vector<double> UsgsAstroLsSensorModel::computeGroundPartials(
    const csm::EcefCoord& ground_pt) const
 {
+
+   MESSAGE_LOG(m_logger, "Computing computeGroundPartials for point {}, {}, {}",
+                                ground_pt.x, ground_pt.y, ground_pt.z)
+
    double GND_DELTA = m_gsd;
    // Partial of line, sample wrt X, Y, Z
    double x = ground_pt.x;
@@ -928,6 +1142,10 @@ csm::RasterGM::SensorPartials UsgsAstroLsSensorModel::computeSensorPartials(
    double* achieved_precision,
    csm::WarningList* warnings) const
 {
+
+  MESSAGE_LOG(m_logger, "Calculating computeSensorPartials for ground point {}, {}, {} with desired precision {}",
+                              ground_pt.x, ground_pt.y, ground_pt.z, desired_precision)
+
    // Compute image coordinate first
    csm::ImageCoord img_pt = groundToImage(
       ground_pt, desired_precision, achieved_precision);
@@ -948,6 +1166,10 @@ csm::RasterGM::SensorPartials UsgsAstroLsSensorModel::computeSensorPartials(
    double*                achieved_precision,
    csm::WarningList*      warnings) const
 {
+
+  MESSAGE_LOG(m_logger, "Calculating computeSensorPartials (with image points {}, {}) for ground point {}, {}, {} with desired precision {}",
+                              image_pt.line, image_pt.samp, ground_pt.x, ground_pt.y, ground_pt.z, desired_precision)
+
    // Compute numerical partials ls wrt specific parameter
 
    const double DELTA = m_gsd;
@@ -974,6 +1196,8 @@ UsgsAstroLsSensorModel::computeAllSensorPartials(
    double*               achieved_precision,
    csm::WarningList*     warnings) const
 {
+  MESSAGE_LOG(m_logger, "Computing computeAllSensorPartials for ground point {}, {}, {} with desired precision {}",
+                              ground_pt.x, ground_pt.y, ground_pt.z, desired_precision)
    csm::ImageCoord image_pt = groundToImage(
       ground_pt, desired_precision, achieved_precision, warnings);
 
@@ -993,6 +1217,10 @@ UsgsAstroLsSensorModel::computeAllSensorPartials(
    double*                achieved_precision,
    csm::WarningList*      warnings) const
 {
+
+   MESSAGE_LOG(m_logger, "Computing computeAllSensorPartials for image {} {} and ground {}, {}, {} with desired precision {}",
+                                image_pt.line, image_pt.samp, ground_pt.x, ground_pt.y, ground_pt.z, desired_precision)
+
    std::vector<int> indices = getParameterSetIndices(pSet);
    size_t num = indices.size();
    std::vector<csm::RasterGM::SensorPartials> partials;
@@ -1013,7 +1241,12 @@ double UsgsAstroLsSensorModel::getParameterCovariance(
    int index1,
    int index2) const
 {
+
    int index = UsgsAstroLsSensorModel::NUM_PARAMETERS * index1 + index2;
+
+   MESSAGE_LOG(m_logger, "getParameterCovariance for {} {} is {}",
+                                index1, index2, m_covariance[index])
+
    return m_covariance[index];
 }
 
@@ -1026,6 +1259,10 @@ void UsgsAstroLsSensorModel::setParameterCovariance(
    double covariance)
 {
    int index = UsgsAstroLsSensorModel::NUM_PARAMETERS * index1 + index2;
+
+   MESSAGE_LOG(m_logger, "setParameterCovariance for {} {} is {}",
+                                index1, index2, m_covariance[index])
+
    m_covariance[index] = covariance;
 }
 
@@ -1058,6 +1295,7 @@ std::string UsgsAstroLsSensorModel::getReferenceDateAndTime() const
 double UsgsAstroLsSensorModel::getImageTime(
    const csm::ImageCoord& image_pt) const
 {
+
    // Flip image taken backwards
    double line1 = image_pt.line;
 
@@ -1080,6 +1318,9 @@ double UsgsAstroLsSensorModel::getImageTime(
    double time = m_intTimeStartTimes[referenceIndex]
       + m_intTimes[referenceIndex] * (lineUSGSFull - m_intTimeLines[referenceIndex]);
 
+  MESSAGE_LOG(m_logger, "getImageTime for image line {} is {}",
+                                image_pt.line, time)
+
    return time;
 
 }
@@ -1090,6 +1331,9 @@ double UsgsAstroLsSensorModel::getImageTime(
 csm::EcefCoord UsgsAstroLsSensorModel::getSensorPosition(
    const csm::ImageCoord& imagePt) const
 {
+   MESSAGE_LOG(m_logger, "getSensorPosition at line {}",
+                                imagePt.line)
+
    return getSensorPosition(getImageTime(imagePt));
 }
 
@@ -1102,6 +1346,9 @@ csm::EcefCoord UsgsAstroLsSensorModel::getSensorPosition(double time) const
    double x, y, z, vx, vy, vz;
    getAdjSensorPosVel(time, _no_adjustment, x, y, z, vx, vy, vz);
 
+   MESSAGE_LOG(m_logger, "getSensorPosition at {}",
+                                time)
+
    return csm::EcefCoord(x, y, z);
 }
 
@@ -1111,6 +1358,8 @@ csm::EcefCoord UsgsAstroLsSensorModel::getSensorPosition(double time) const
 csm::EcefVector UsgsAstroLsSensorModel::getSensorVelocity(
    const csm::ImageCoord& imagePt) const
 {
+  MESSAGE_LOG(m_logger, "getSensorVelocity at {}",
+                               imagePt.line)
    return getSensorVelocity(getImageTime(imagePt));
 }
 
@@ -1122,6 +1371,9 @@ csm::EcefVector UsgsAstroLsSensorModel::getSensorVelocity(double time) const
    double x, y, z, vx, vy, vz;
    getAdjSensorPosVel(time, _no_adjustment, x, y, z, vx, vy, vz);
 
+   MESSAGE_LOG(m_logger, "getSensorVelocity at {}",
+                                time)
+
    return csm::EcefVector(vx, vy, vz);
 }
 
@@ -1379,6 +1631,10 @@ UsgsAstroLsSensorModel::getValidImageRange() const
 csm::EcefVector UsgsAstroLsSensorModel::getIlluminationDirection(
    const csm::EcefCoord& groundPt) const
 {
+  MESSAGE_LOG(m_logger, "Accessing illimination direction of ground point"
+              "{} {} {}."
+              "Illimination direction is not supported, throwing exception",
+              groundPt.x, groundPt.y, groundPt.z);
    throw csm::Error(
       csm::Error::UNSUPPORTED_FUNCTION,
       "Unsupported function",
@@ -1400,9 +1656,11 @@ int UsgsAstroLsSensorModel::getNumGeometricCorrectionSwitches() const
 //***************************************************************************
 // UsgsAstroLsSensorModel::getGeometricCorrectionName
 //***************************************************************************
-
 std::string UsgsAstroLsSensorModel::getGeometricCorrectionName(int index) const
 {
+  MESSAGE_LOG(m_logger, "Accessing name of geometric correction switch {}. "
+              "Geometric correction switches are not supported, throwing exception",
+              index);
    // Since there are no geometric corrections, all indices are out of range
    throw csm::Error(
       csm::Error::INDEX_OUT_OF_RANGE,
@@ -1419,6 +1677,10 @@ void UsgsAstroLsSensorModel::setGeometricCorrectionSwitch(
    csm::param::Type pType)
 
 {
+   MESSAGE_LOG(m_logger, "Setting geometric correction switch {} to {} "
+               "with parameter type {}. "
+               "Geometric correction switches are not supported, throwing exception",
+               index, value, pType);
    // Since there are no geometric corrections, all indices are out of range
    throw csm::Error(
       csm::Error::INDEX_OUT_OF_RANGE,
@@ -1431,6 +1693,9 @@ void UsgsAstroLsSensorModel::setGeometricCorrectionSwitch(
 //***************************************************************************
 bool UsgsAstroLsSensorModel::getGeometricCorrectionSwitch(int index) const
 {
+   MESSAGE_LOG(m_logger, "Accessing value of geometric correction switch {}. "
+               "Geometric correction switches are not supported, throwing exception",
+               index);
    // Since there are no geometric corrections, all indices are out of range
    throw csm::Error(
       csm::Error::INDEX_OUT_OF_RANGE,
@@ -1455,6 +1720,9 @@ std::vector<double> UsgsAstroLsSensorModel::getCrossCovarianceMatrix(
    return std::vector<double>(num_rows * num_cols, 0.0);
 }
 
+//***************************************************************************
+// UsgsAstroLineScannerSensorModel::getCorrelationModel
+//***************************************************************************
 const csm::CorrelationModel&
 UsgsAstroLsSensorModel::getCorrelationModel() const
 {
@@ -1506,6 +1774,8 @@ bool UsgsAstroLsSensorModel::isParameterShareable(int index) const
 csm::SharingCriteria UsgsAstroLsSensorModel::getParameterSharingCriteria(
    int index) const
 {
+  MESSAGE_LOG(m_logger, "Checking sharing criteria for parameter {}. "
+              "Sharing is not supported, throwing exception", index);
    // Parameter sharing is not supported for this sensor,
    // all indices are out of range
    throw csm::Error(
@@ -1538,7 +1808,9 @@ csm::Version UsgsAstroLsSensorModel::getVersion() const
    return csm::Version(1, 0, 0);
 }
 
-
+//***************************************************************************
+// UsgsAstroLineScannerSensorModel::getEllipsoid
+//***************************************************************************
 csm::Ellipsoid UsgsAstroLsSensorModel::getEllipsoid() const
 {
    return csm::Ellipsoid(m_majorAxis, m_minorAxis);
@@ -1551,8 +1823,9 @@ void UsgsAstroLsSensorModel::setEllipsoid(
    m_minorAxis = ellipsoid.getSemiMinorRadius();
 }
 
-
-
+//***************************************************************************
+// UsgsAstroLineScannerSensorModel::getValue
+//***************************************************************************
 double UsgsAstroLsSensorModel::getValue(
    int   index,
    const std::vector<double> &adjustments) const
@@ -1564,7 +1837,6 @@ double UsgsAstroLsSensorModel::getValue(
 //***************************************************************************
 // Functions pulled out of losToEcf and computeViewingPixel
 // **************************************************************************
-
 void UsgsAstroLsSensorModel::getQuaternions(const double& time, double q[4]) const{
   int nOrder = 8;
   if (m_platformFlag == 0)
@@ -1572,16 +1844,24 @@ void UsgsAstroLsSensorModel::getQuaternions(const double& time, double q[4]) con
   int nOrderQuat = nOrder;
   if (m_numQuaternions < 6 && nOrder == 8)
      nOrderQuat = 4;
+
+  MESSAGE_LOG(m_logger, "Calculating getQuaternions for time {} with {}"
+                              "order lagrange",
+                              time, nOrder)
   lagrangeInterp(
      m_numQuaternions/4, &m_quaternions[0], m_t0Quat, m_dtQuat, time, 4, nOrderQuat, q);
 }
 
-
+//***************************************************************************
+// UsgsAstroLineScannerSensorModel::calculateAttitudeCorrection
+//***************************************************************************
 void UsgsAstroLsSensorModel::calculateAttitudeCorrection(
    const double& time,
    const std::vector<double>& adj,
    double attCorr[9]) const
 {
+  MESSAGE_LOG(m_logger, "Computing calculateAttitudeCorrection (with adjustment)"
+                              "for time {}", time)
   double aTime = time - m_t0Quat;
   double euler[3];
   double nTime = aTime / m_halfTime;
@@ -1592,6 +1872,8 @@ void UsgsAstroLsSensorModel::calculateAttitudeCorrection(
     (getValue(7, adj) + getValue(10, adj)* nTime + getValue(13, adj)* nTime2) / m_flyingHeight;
   euler[2] =
     (getValue(8, adj) + getValue(11, adj)* nTime + getValue(14, adj)* nTime2) / m_halfSwath;
+  MESSAGE_LOG(m_logger, "calculateAttitudeCorrection: euler {} {} {}",
+                              euler[0], euler[1], euler[2])
 
   calculateRotationMatrixFromEuler(euler, attCorr);
 }
@@ -1617,6 +1899,9 @@ void UsgsAstroLsSensorModel::losToEcf(
    //# private_func_description
    // Computes image ray (look vector) in ecf coordinate system.
    // Compute adjusted sensor position and velocity
+   MESSAGE_LOG(m_logger, "Computing losToEcf (with adjustments) for"
+                               "line {} sample {}",
+                                line, sample)
 
    double time = getImageTime(csm::ImageCoord(line, sample));
    getAdjSensorPosVel(time, adj, xc, yc, zc, vx, vy, vz);
@@ -1636,7 +1921,7 @@ void UsgsAstroLsSensorModel::losToEcf(
          m_iTransS, m_iTransL,
          distortedFocalPlaneX, distortedFocalPlaneY);
 
-   // Remove lens distortion
+   // Remove lens
    double undistortedFocalPlaneX, undistortedFocalPlaneY;
    removeDistortion(distortedFocalPlaneX, distortedFocalPlaneY,
                     undistortedFocalPlaneX, undistortedFocalPlaneY,
@@ -1661,7 +1946,9 @@ void UsgsAstroLsSensorModel::losToEcf(
    correctedCameraLook[2] = attCorr[6] * cameraLook[0]
                           + attCorr[7] * cameraLook[1]
                           + attCorr[8] * cameraLook[2];
-
+   MESSAGE_LOG(m_logger, "losToEcf: corrected camera look vector {} {} {}",
+                               correctedCameraLook[0], correctedCameraLook[1],
+                               correctedCameraLook[2])
 // Rotate the look vector into the body fixed frame from the camera reference frame by applying the rotation matrix from the sensor quaternions
    double quaternions[4];
    getQuaternions(time, quaternions);
@@ -1677,6 +1964,9 @@ void UsgsAstroLsSensorModel::losToEcf(
    bodyLookZ = cameraToBody[6] * correctedCameraLook[0]
              + cameraToBody[7] * correctedCameraLook[1]
              + cameraToBody[8] * correctedCameraLook[2];
+   MESSAGE_LOG(m_logger, "losToEcf: body look vector {} {} {}",
+                                bodyLookX, bodyLookY, bodyLookZ)
+
 }
 
 
@@ -1694,6 +1984,9 @@ void UsgsAstroLsSensorModel::lightAberrationCorr(
    double&       dyl,
    double&       dzl) const
 {
+   MESSAGE_LOG(m_logger, "Computing lightAberrationCorr for camera velocity"
+                               "{} {} {} and image ray {} {} {}",
+                                vx, vy, vz, xl, yl, zl)
    //# func_description
    //  Computes light aberration correction vector
 
@@ -1712,6 +2005,8 @@ void UsgsAstroLsSensorModel::lightAberrationCorr(
       dxl = 0.0;
       dyl = 0.0;
       dzl = 0.0;
+      MESSAGE_LOG(m_logger, "lightAberrationCorr: image ray is parallel"
+                                  "to velocity vector")
    }
 
    // Compute the angle between the corrected image ray and spacecraft
@@ -1730,6 +2025,8 @@ void UsgsAstroLsSensorModel::lightAberrationCorr(
    dxl = cfac * vx;
    dyl = cfac * vy;
    dzl = cfac * vz;
+   MESSAGE_LOG(m_logger, "lightAberrationCorr: light of sight correction"
+                               "{} {} {}", dxl, dyl, dzl)
 }
 
 //***************************************************************************
@@ -1743,6 +2040,9 @@ void UsgsAstroLsSensorModel::computeElevation(
    double&       achieved_precision,
    const double& desired_precision) const
 {
+   MESSAGE_LOG(m_logger, "Calculating computeElevation for {} {} {}"
+                               "with desired precision {}",
+                              x, y, z, desired_precision)
    // Compute elevation given xyz
    // Requires semi-major-axis and eccentricity-square
    const int MKTR = 10;
@@ -1770,6 +2070,7 @@ void UsgsAstroLsSensorModel::computeElevation(
          tanPhi = zz / d;
          ktr++;
       } while (MKTR > ktr && fabs(h - hPrev) > desired_precision);
+      MESSAGE_LOG(m_logger, "computeElevation: point is near equator")
    }
 
    // Suited for points near the poles
@@ -1788,10 +2089,14 @@ void UsgsAstroLsSensorModel::computeElevation(
          cotPhi = dd / z;
          ktr++;
       } while (MKTR > ktr && fabs(h - hPrev) > desired_precision);
+      MESSAGE_LOG(m_logger, "computeElevation: point is near poles")
    }
 
    height = h;
    achieved_precision = fabs(h - hPrev);
+   MESSAGE_LOG(m_logger, "computeElevation: height {} with achieved"
+                               "precision of {}",
+                                height, achieved_precision)
 }
 
 //***************************************************************************
@@ -1811,6 +2116,10 @@ void UsgsAstroLsSensorModel::losEllipsoidIntersect(
    double&       achieved_precision,
    const double& desired_precision) const
 {
+   MESSAGE_LOG(m_logger, "Computing losEllipsoidIntersect for camera position"
+                               "{} {} {} looking {} {} {} with desired precision"
+                               "{}",
+                                xc, yc, zc, xl, yl, zl, desired_precision)
    // Helper function which computes the intersection of the image ray
    // with the ellipsoid.  All vectors are in earth-centered-fixed
    // coordinate system with origin at the center of the earth.
@@ -1860,6 +2169,9 @@ void UsgsAstroLsSensorModel::losEllipsoidIntersect(
    slope = -1;
 
    achieved_precision = fabs(height - h);
+   MESSAGE_LOG(m_logger, "losEllipsoidIntersect: found intersect at {} {} {}"
+                                "with achieved precision of {}",
+                                x, y, z, achieved_precision)
 }
 
 //***************************************************************************
@@ -1879,6 +2191,9 @@ void UsgsAstroLsSensorModel::losPlaneIntersect(
                           //         0(X), 1(Y), or 2(Z) fixed
                           // output: 0(X), 1(Y), or 2(Z) fixed
 {
+  MESSAGE_LOG(m_logger, "Calculating losPlaneIntersect for camera position"
+                              "{} {} {} and image ray {} {} {}",
+                              xc, yc, zc, xl, yl, zl)
    //# func_description
    //  Computes 2 of the 3 coordinates of a ground point given the 3rd
    //  coordinate.  The 3rd coordinate that is held fixed corresponds
@@ -1901,7 +2216,8 @@ void UsgsAstroLsSensorModel::losPlaneIntersect(
          mode = 2;
       }
    }
-
+   MESSAGE_LOG(m_logger, "losPlaneIntersect: largest/fixed image ray component"
+                                "{} (1-x, 2-y, 3-z)", mode)
    // X is the fixed or largest component
 
    if (0 == mode)
@@ -1925,6 +2241,7 @@ void UsgsAstroLsSensorModel::losPlaneIntersect(
       x = xc + (z - zc) * xl / zl;
       y = yc + (z - zc) * yl / zl;
    }
+   MESSAGE_LOG(m_logger, "ground coordinate {} {} {}", x, y, z)
 }
 
 //***************************************************************************
@@ -1940,6 +2257,7 @@ void UsgsAstroLsSensorModel::imageToPlane(
    double&       z,
    int&          mode) const
 {
+  MESSAGE_LOG(m_logger, "Computing imageToPlane")
    //# func_description
    //  Computes ground coordinates by intersecting image ray with
    //  a plane perpendicular to the coordinate axis with the largest
@@ -1971,6 +2289,9 @@ void UsgsAstroLsSensorModel::getAdjSensorPosVel(
    double&       vy,
    double&       vz) const
 {
+  MESSAGE_LOG(m_logger, "Calculating getAdjSensorPosVel at time {}",
+                              time)
+
    // Sensor position and velocity (4th or 8th order Lagrange).
    int nOrder = 8;
    if (m_platformFlag == 0)
@@ -1981,8 +2302,10 @@ void UsgsAstroLsSensorModel::getAdjSensorPosVel(
    double sensVelNom[3];
    lagrangeInterp(m_numPositions/3, &m_velocities[0], m_t0Ephem, m_dtEphem,
       time, 3, nOrder, sensVelNom);
-   // Compute rotation matrix from ICR to ECF
+   MESSAGE_LOG(m_logger, "getAdjSensorPosVel: using {} order Lagrange",
+                                nOrder)
 
+   // Compute rotation matrix from ICR to ECF
    double radialUnitVec[3];
    double radMag = sqrt(sensPosNom[0] * sensPosNom[0] +
       sensPosNom[1] * sensPosNom[1] +
@@ -2018,8 +2341,8 @@ void UsgsAstroLsSensorModel::getAdjSensorPosVel(
    ecfFromIcr[6] = inTrackUnitVec[2];
    ecfFromIcr[7] = crossTrackUnitVec[2];
    ecfFromIcr[8] = radialUnitVec[2];
-   // Apply position and velocity corrections
 
+   // Apply position and velocity corrections
    double aTime = time - m_t0Ephem;
    double dvi = getValue(3, adj) / m_halfTime;
    double dvc = getValue(4, adj) / m_halfTime;
@@ -2039,6 +2362,10 @@ void UsgsAstroLsSensorModel::getAdjSensorPosVel(
       + ecfFromIcr[3] * di + ecfFromIcr[4] * dc + ecfFromIcr[5] * dr;
    zc = sensPosNom[2]
       + ecfFromIcr[6] * di + ecfFromIcr[7] * dc + ecfFromIcr[8] * dr;
+
+   MESSAGE_LOG(m_logger, "getAdjSensorPosVel: postition {} {} {}"
+                                "and velocity {} {} {}",
+                              xc, yc, zc, vx, vy, vz)
 }
 
 
@@ -2051,6 +2378,10 @@ csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel(
    const std::vector<double>& adj,
    const double& desiredPrecision) const
 {
+  MESSAGE_LOG(m_logger, "Computing computeViewingPixel (with adjusments)"
+                              "for ground point {} {} {} at time {} ",
+                              groundPoint.x, groundPoint.y, groundPoint.z, time)
+
   // Helper function to compute the CCD pixel that views a ground point based
   // on the exterior orientation at a given time.
 
@@ -2062,6 +2393,8 @@ csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel(
    double bodyLookX = groundPoint.x - xc;
    double bodyLookY = groundPoint.y - yc;
    double bodyLookZ = groundPoint.z - zc;
+   MESSAGE_LOG(m_logger, "computeViewingPixel: look vector {} {} {}",
+                                bodyLookX, bodyLookY, bodyLookZ)
 
    // Rotate the look vector into the camera reference frame
    double quaternions[4];
@@ -2079,6 +2412,9 @@ csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel(
    double cameraLookZ = bodyToCamera[2] * bodyLookX
                       + bodyToCamera[5] * bodyLookY
                       + bodyToCamera[8] * bodyLookZ;
+   MESSAGE_LOG(m_logger, "computeViewingPixel: look vector (camrea ref frame)"
+                               "{} {} {}",
+                                cameraLookX, cameraLookY, cameraLookZ)
 
    // Invert the attitude correction
    double attCorr[9];
@@ -2094,12 +2430,18 @@ csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel(
    double adjustedLookZ = attCorr[2] * cameraLookX
                         + attCorr[5] * cameraLookY
                         + attCorr[8] * cameraLookZ;
+   MESSAGE_LOG(m_logger, "computeViewingPixel: adjusted look vector"
+                               "{} {} {}",
+                                adjustedLookX, adjustedLookY, adjustedLookZ)
 
    // Convert to focal plane coordinate
    double lookScale = m_focalLength / adjustedLookZ;
    double focalX = adjustedLookX * lookScale;
    double focalY = adjustedLookY * lookScale;
    double distortedFocalX, distortedFocalY;
+   MESSAGE_LOG(m_logger, "computeViewingPixel: focal plane coordinates"
+                                "x:{} y:{} scale:{}",
+                                focalX, focalY, lookScale)
 
    // Invert distortion
    applyDistortion(focalX, focalY, distortedFocalX, distortedFocalY,
@@ -2117,6 +2459,8 @@ csm::ImageCoord UsgsAstroLsSensorModel::computeViewingPixel(
    double line = detectorLine + m_detectorLineOrigin;
    double sample = (detectorSample + m_detectorSampleOrigin - m_startingSample)
                  / m_detectorSampleSumming;
+   MESSAGE_LOG(m_logger, "computeViewingPixel: image line sample {} {}",
+                                line, sample)
 
    return csm::ImageCoord(line, sample);
 }
@@ -2143,11 +2487,15 @@ void UsgsAstroLsSensorModel::computeLinearApproximation(
 
       if (ip.samp < 0.0)     ip.samp = 0.0;
       if (ip.samp > numCols) ip.samp = numCols;
+      MESSAGE_LOG(m_logger, "Computing computeLinearApproximation"
+                                  "with linear approximation")
    }
    else
    {
       ip.line = m_nLines / 2.0;
       ip.samp = m_nSamples / 2.0;
+      MESSAGE_LOG(m_logger, "Computing computeLinearApproximation"
+                                  "nonlinear approx line/2 and sample/2")
    }
 }
 
@@ -2157,124 +2505,136 @@ void UsgsAstroLsSensorModel::computeLinearApproximation(
 //***************************************************************************
 void UsgsAstroLsSensorModel::setLinearApproximation()
 {
-   double u_factors[4] = { 0.0, 0.0, 1.0, 1.0 };
-   double v_factors[4] = { 0.0, 1.0, 0.0, 1.0 };
-
-   csm::EcefCoord refPt = getReferencePoint();
-
-   double height, aPrec;
-   double desired_precision = 0.01;
-   computeElevation(refPt.x, refPt.y, refPt.z, height, aPrec, desired_precision);
-   if (isnan(height))
-   {
-      _linear = false;
-      return;
-   }
-
-
-   double numRows = m_nLines;
-   double numCols = m_nSamples;
-
-   csm::ImageCoord imagePt;
-   csm::EcefCoord  gp[8];
-
-   int i;
-   for (i = 0; i < 4; i++)
-   {
-      imagePt.line = u_factors[i] * numRows;
-      imagePt.samp = v_factors[i] * numCols;
-      gp[i] = imageToGround(imagePt, height);
-   }
+  MESSAGE_LOG(m_logger, "Calculating setLinearApproximation")
+  double u_factors[4] = { 0.0, 0.0, 1.0, 1.0 };
+  double v_factors[4] = { 0.0, 1.0, 0.0, 1.0 };
+
+  csm::EcefCoord refPt = getReferencePoint();
+
+  double height, aPrec;
+  double desired_precision = 0.01;
+  computeElevation(refPt.x, refPt.y, refPt.z, height, aPrec, desired_precision);
+  if (isnan(height))
+  {
+    MESSAGE_LOG(m_logger, "setLinearApproximation: computeElevation of"
+                                "reference point {} {} {} with desired precision"
+                                "{} returned nan height; nonlinear",
+                                refPt.x, refPt.y, refPt.z, desired_precision)
+    _linear = false;
+    return;
+  }
+  MESSAGE_LOG(m_logger, "setLinearApproximation: computeElevation of"
+                              "reference point {} {} {} with desired precision"
+                              "{} returned {} height",
+                              refPt.x, refPt.y, refPt.z, desired_precision, height)
+
+  double numRows = m_nLines;
+  double numCols = m_nSamples;
+
+  csm::ImageCoord imagePt;
+  csm::EcefCoord  gp[8];
+
+  int i;
+  for (i = 0; i < 4; i++)
+  {
+    imagePt.line = u_factors[i] * numRows;
+    imagePt.samp = v_factors[i] * numCols;
+    gp[i] = imageToGround(imagePt, height);
+  }
 
-   double delz = 100.0;
-   height += delz;
-   for (i = 0; i < 4; i++)
-   {
-      imagePt.line = u_factors[i] * numRows;
-      imagePt.samp = v_factors[i] * numCols;
-      gp[i + 4] = imageToGround(imagePt, height);
-   }
+  double delz = 100.0;
+  height += delz;
+  for (i = 0; i < 4; i++)
+  {
+    imagePt.line = u_factors[i] * numRows;
+    imagePt.samp = v_factors[i] * numCols;
+    gp[i + 4] = imageToGround(imagePt, height);
+  }
 
-   csm::EcefCoord d_du;
-   d_du.x = (
-      (gp[2].x + gp[3].x + gp[6].x + gp[7].x) -
-      (gp[0].x + gp[1].x + gp[4].x + gp[5].x)) / numRows / 4.0;
-   d_du.y = (
-      (gp[2].y + gp[3].y + gp[6].y + gp[7].y) -
-      (gp[0].y + gp[1].y + gp[4].y + gp[5].y)) / numRows / 4.0;
-   d_du.z = (
-      (gp[2].z + gp[3].z + gp[6].z + gp[7].z) -
-      (gp[0].z + gp[1].z + gp[4].z + gp[5].z)) / numRows / 4.0;
-
-   csm::EcefCoord d_dv;
-   d_dv.x = (
-      (gp[1].x + gp[3].x + gp[5].x + gp[7].x) -
-      (gp[0].x + gp[2].x + gp[4].x + gp[6].x)) / numCols / 4.0;
-   d_dv.y = (
-      (gp[1].y + gp[3].y + gp[5].y + gp[7].y) -
-      (gp[0].y + gp[2].y + gp[4].y + gp[6].y)) / numCols / 4.0;
-   d_dv.z = (
-      (gp[1].z + gp[3].z + gp[5].z + gp[7].z) -
-      (gp[0].z + gp[2].z + gp[4].z + gp[6].z)) / numCols / 4.0;
-
-   double mat3x3[9];
-
-   mat3x3[0] = d_du.x;
-   mat3x3[1] = d_dv.x;
-   mat3x3[2] = 1.0;
-   mat3x3[3] = d_du.y;
-   mat3x3[4] = d_dv.y;
-   mat3x3[5] = 1.0;
-   mat3x3[6] = d_du.z;
-   mat3x3[7] = d_dv.z;
-   mat3x3[8] = 1.0;
-
-   double denom = determinant3x3(mat3x3);
-
-   if (fabs(denom) < 1.0e-8) // can not get derivatives this way
-   {
-      _linear = false;
-      return;
-   }
+  csm::EcefCoord d_du;
+  d_du.x = (
+    (gp[2].x + gp[3].x + gp[6].x + gp[7].x) -
+    (gp[0].x + gp[1].x + gp[4].x + gp[5].x)) / numRows / 4.0;
+  d_du.y = (
+    (gp[2].y + gp[3].y + gp[6].y + gp[7].y) -
+    (gp[0].y + gp[1].y + gp[4].y + gp[5].y)) / numRows / 4.0;
+  d_du.z = (
+    (gp[2].z + gp[3].z + gp[6].z + gp[7].z) -
+    (gp[0].z + gp[1].z + gp[4].z + gp[5].z)) / numRows / 4.0;
+
+  csm::EcefCoord d_dv;
+  d_dv.x = (
+    (gp[1].x + gp[3].x + gp[5].x + gp[7].x) -
+    (gp[0].x + gp[2].x + gp[4].x + gp[6].x)) / numCols / 4.0;
+  d_dv.y = (
+    (gp[1].y + gp[3].y + gp[5].y + gp[7].y) -
+    (gp[0].y + gp[2].y + gp[4].y + gp[6].y)) / numCols / 4.0;
+  d_dv.z = (
+    (gp[1].z + gp[3].z + gp[5].z + gp[7].z) -
+    (gp[0].z + gp[2].z + gp[4].z + gp[6].z)) / numCols / 4.0;
+
+  double mat3x3[9];
+
+  mat3x3[0] = d_du.x;
+  mat3x3[1] = d_dv.x;
+  mat3x3[2] = 1.0;
+  mat3x3[3] = d_du.y;
+  mat3x3[4] = d_dv.y;
+  mat3x3[5] = 1.0;
+  mat3x3[6] = d_du.z;
+  mat3x3[7] = d_dv.z;
+  mat3x3[8] = 1.0;
+
+  double denom = determinant3x3(mat3x3);
+
+  if (fabs(denom) < 1.0e-8) // can not get derivatives this way
+  {
+    MESSAGE_LOG(m_logger, "setLinearApproximation: determinant3x3 of"
+                                "matrix of partials is {}; nonlinear",
+                                denom)
+    _linear = false;
+    return;
+  }
 
-   mat3x3[0] = 1.0;
-   mat3x3[3] = 0.0;
-   mat3x3[6] = 0.0;
-   _du_dx = determinant3x3(mat3x3) / denom;
+  mat3x3[0] = 1.0;
+  mat3x3[3] = 0.0;
+  mat3x3[6] = 0.0;
+  _du_dx = determinant3x3(mat3x3) / denom;
 
-   mat3x3[0] = 0.0;
-   mat3x3[3] = 1.0;
-   mat3x3[6] = 0.0;
-   _du_dy = determinant3x3(mat3x3) / denom;
+  mat3x3[0] = 0.0;
+  mat3x3[3] = 1.0;
+  mat3x3[6] = 0.0;
+  _du_dy = determinant3x3(mat3x3) / denom;
 
-   mat3x3[0] = 0.0;
-   mat3x3[3] = 0.0;
-   mat3x3[6] = 1.0;
-   _du_dz = determinant3x3(mat3x3) / denom;
+  mat3x3[0] = 0.0;
+  mat3x3[3] = 0.0;
+  mat3x3[6] = 1.0;
+  _du_dz = determinant3x3(mat3x3) / denom;
 
-   mat3x3[0] = d_du.x;
-   mat3x3[3] = d_du.y;
-   mat3x3[6] = d_du.z;
+  mat3x3[0] = d_du.x;
+  mat3x3[3] = d_du.y;
+  mat3x3[6] = d_du.z;
 
-   mat3x3[1] = 1.0;
-   mat3x3[4] = 0.0;
-   mat3x3[7] = 0.0;
-   _dv_dx = determinant3x3(mat3x3) / denom;
+  mat3x3[1] = 1.0;
+  mat3x3[4] = 0.0;
+  mat3x3[7] = 0.0;
+  _dv_dx = determinant3x3(mat3x3) / denom;
 
-   mat3x3[1] = 0.0;
-   mat3x3[4] = 1.0;
-   mat3x3[7] = 0.0;
-   _dv_dy = determinant3x3(mat3x3) / denom;
+  mat3x3[1] = 0.0;
+  mat3x3[4] = 1.0;
+  mat3x3[7] = 0.0;
+  _dv_dy = determinant3x3(mat3x3) / denom;
 
-   mat3x3[1] = 0.0;
-   mat3x3[4] = 0.0;
-   mat3x3[7] = 1.0;
-   _dv_dz = determinant3x3(mat3x3) / denom;
+  mat3x3[1] = 0.0;
+  mat3x3[4] = 0.0;
+  mat3x3[7] = 1.0;
+  _dv_dz = determinant3x3(mat3x3) / denom;
 
-   _u0 = -gp[0].x * _du_dx - gp[0].y * _du_dy - gp[0].z * _du_dz;
-   _v0 = -gp[0].x * _dv_dx - gp[0].y * _dv_dy - gp[0].z * _dv_dz;
+  _u0 = -gp[0].x * _du_dx - gp[0].y * _du_dy - gp[0].z * _du_dz;
+  _v0 = -gp[0].x * _dv_dx - gp[0].y * _dv_dy - gp[0].z * _dv_dz;
 
-   _linear = true;
+  _linear = true;
+  MESSAGE_LOG(m_logger, "Completed setLinearApproximation")
 }
 
 //***************************************************************************
@@ -2282,17 +2642,19 @@ void UsgsAstroLsSensorModel::setLinearApproximation()
 //***************************************************************************
 double UsgsAstroLsSensorModel::determinant3x3(double mat[9]) const
 {
-   return
-      mat[0] * (mat[4] * mat[8] - mat[7] * mat[5]) -
-      mat[1] * (mat[3] * mat[8] - mat[6] * mat[5]) +
-      mat[2] * (mat[3] * mat[7] - mat[6] * mat[4]);
+  return
+    mat[0] * (mat[4] * mat[8] - mat[7] * mat[5]) -
+    mat[1] * (mat[3] * mat[8] - mat[6] * mat[5]) +
+    mat[2] * (mat[3] * mat[7] - mat[6] * mat[4]);
 }
 
 
-
-
+//***************************************************************************
+// UsgsAstroLineScannerSensorModel::constructStateFromIsd
+//***************************************************************************
 std::string UsgsAstroLsSensorModel::constructStateFromIsd(const std::string imageSupportData, csm::WarningList *warnings) const
 {
+  MESSAGE_LOG(m_logger, "Constructing state from Isd")
   // Instantiate UsgsAstroLineScanner sensor model
   json isd = json::parse(imageSupportData);
   json state = {};
@@ -2305,47 +2667,100 @@ std::string UsgsAstroLsSensorModel::constructStateFromIsd(const std::string imag
   state["m_imageIdentifier"] = getImageId(isd, parsingWarnings);
   state["m_sensorName"] = getSensorName(isd, parsingWarnings);
   state["m_platformName"] = getPlatformName(isd, parsingWarnings);
+  MESSAGE_LOG(m_logger, "m_modelName: {} "
+                        "m_imageIdentifier: {} "
+                        "m_sensorName: {} "
+                        "m_platformName: {} ",
+                        state["m_modelName"].dump(),
+                        state["m_imageIdentifier"].dump(),
+                        state["m_sensorName"].dump(),
+                        state["m_platformName"].dump())
 
   state["m_focalLength"] = getFocalLength(isd, parsingWarnings);
+  MESSAGE_LOG(m_logger, "m_focalLength: {} ", state["m_focalLength"].dump())
 
   state["m_nLines"] = getTotalLines(isd, parsingWarnings);
   state["m_nSamples"] = getTotalSamples(isd, parsingWarnings);
+  MESSAGE_LOG(m_logger, "m_nLines: {} "
+                        "m_nSamples: {} ",
+                        state["m_nLines"].dump(), state["m_nSamples"].dump())
 
   state["m_iTransS"] = getFocal2PixelSamples(isd, parsingWarnings);
   state["m_iTransL"] = getFocal2PixelLines(isd, parsingWarnings);
+  MESSAGE_LOG(m_logger, "m_iTransS: {} "
+                        "m_iTransL: {} ",
+                        state["m_iTransS"].dump(), state["m_iTransL"].dump())
 
   state["m_platformFlag"] = 1;
   state["m_ikCode"] = 0;
   state["m_zDirection"] = 1;
+  MESSAGE_LOG(m_logger, "m_platformFlag: {} "
+                        "m_ikCode: {} "
+                        "m_zDirection: {} ",
+                        state["m_platformFlag"].dump(), state["m_ikCode"].dump(),
+                        state["m_zDirection"].dump())
 
   state["m_distortionType"] = getDistortionModel(isd, parsingWarnings);
   state["m_opticalDistCoeffs"] = getDistortionCoeffs(isd, parsingWarnings);
+  MESSAGE_LOG(m_logger, "m_distortionType: {} "
+                        "m_opticalDistCoeffs: {} ",
+                        state["m_distortionType"].dump(),
+                        state["m_opticalDistCoeffs"].dump())
 
   // Zero computed state values
   state["m_referencePointXyz"] = std::vector<double>(3, 0.0);
+  MESSAGE_LOG(m_logger, "m_referencePointXyz: {} ",
+                        state["m_referencePointXyz"].dump())
 
   // leave these be for now.
   state["m_gsd"] = 1.0;
   state["m_flyingHeight"] = 1000.0;
   state["m_halfSwath"] = 1000.0;
   state["m_halfTime"] = 10.0;
+  MESSAGE_LOG(m_logger, "m_gsd: {} "
+                        "m_flyingHeight: {} "
+                        "m_halfSwath: {} "
+                        "m_halfTime: {} ",
+                        state["m_gsd"].dump(), state["m_flyingHeight"].dump(),
+                        state["m_halfSwath"].dump(), state["m_halfTime"].dump())
 
   state["m_centerEphemerisTime"] = getCenterTime(isd, parsingWarnings);
   state["m_startingEphemerisTime"] = getStartingTime(isd, parsingWarnings);
+  MESSAGE_LOG(m_logger, "m_centerEphemerisTime: {} "
+                        "m_startingEphemerisTime: {} ",
+                        state["m_centerEphemerisTime"].dump(),
+                        state["m_startingEphemerisTime"].dump())
 
   state["m_intTimeLines"] = getIntegrationStartLines(isd, parsingWarnings);
   state["m_intTimeStartTimes"] = getIntegrationStartTimes(isd, parsingWarnings);
   state["m_intTimes"] = getIntegrationTimes(isd, parsingWarnings);
+  MESSAGE_LOG(m_logger, "m_intTimeLines: {} "
+                        "m_intTimeStartTimes: {} "
+                        "m_intTimes: {} ",
+                        state["m_intTimeLines"].dump(),
+                        state["m_intTimeStartTimes"].dump(),
+                        state["m_intTimes"].dump())
 
   state["m_detectorSampleSumming"] = getSampleSumming(isd, parsingWarnings);
   state["m_startingDetectorSample"] = getDetectorStartingSample(isd, parsingWarnings);
   state["m_startingDetectorLine"] = getDetectorStartingLine(isd, parsingWarnings);
   state["m_detectorSampleOrigin"] = getDetectorCenterSample(isd, parsingWarnings);
   state["m_detectorLineOrigin"] = getDetectorCenterLine(isd, parsingWarnings);
+  MESSAGE_LOG(m_logger, "m_detectorSampleSumming: {} "
+                        "m_startingDetectorSample: {} "
+                        "m_startingDetectorLine: {} "
+                        "m_detectorSampleOrigin: {} "
+                        "m_detectorLineOrigin: {} ",
+                        state["m_detectorSampleSumming"].dump(),
+                        state["m_startingDetectorSample"].dump(),
+                        state["m_startingDetectorLine"].dump(),
+                        state["m_detectorSampleOrigin"].dump(),
+                        state["m_detectorLineOrigin"].dump())
 
   // These are exlusive to LineScanners, leave them here for now.
   try {
     state["m_dtEphem"] = isd.at("dt_ephemeris");
+    MESSAGE_LOG(m_logger, "m_dtEphem: {} ", state["m_dtEphem"].dump())
   }
   catch(...) {
     parsingWarnings->push_back(
@@ -2353,10 +2768,12 @@ std::string UsgsAstroLsSensorModel::constructStateFromIsd(const std::string imag
         csm::Warning::DATA_NOT_AVAILABLE,
         "dt_ephemeris not in ISD",
         "UsgsAstroFrameSensorModel::constructStateFromIsd()"));
+    MESSAGE_LOG(m_logger, "m_dtEphem not in ISD")
   }
 
   try {
     state["m_t0Ephem"] = isd.at("t0_ephemeris");
+    MESSAGE_LOG(m_logger, "t0_ephemeris: {}", state["m_t0Ephem"].dump())
   }
   catch(...) {
     parsingWarnings->push_back(
@@ -2364,10 +2781,12 @@ std::string UsgsAstroLsSensorModel::constructStateFromIsd(const std::string imag
         csm::Warning::DATA_NOT_AVAILABLE,
         "t0_ephemeris not in ISD",
         "UsgsAstroFrameSensorModel::constructStateFromIsd()"));
+    MESSAGE_LOG(m_logger, "t0_ephemeris not in ISD")
   }
 
   try{
     state["m_dtQuat"] =  isd.at("dt_quaternion");
+    MESSAGE_LOG(m_logger, "dt_quaternion: {}", state["m_dtQuat"].dump())
   }
   catch(...) {
     parsingWarnings->push_back(
@@ -2375,10 +2794,12 @@ std::string UsgsAstroLsSensorModel::constructStateFromIsd(const std::string imag
         csm::Warning::DATA_NOT_AVAILABLE,
         "dt_quaternion not in ISD",
         "UsgsAstroFrameSensorModel::constructStateFromIsd()"));
+    MESSAGE_LOG(m_logger, "dt_quaternion not in ISD")
   }
 
   try{
     state["m_t0Quat"] =  isd.at("t0_quaternion");
+    MESSAGE_LOG(m_logger, "m_t0Quat: {}", state["m_t0Quat"].dump())
   }
   catch(...) {
     parsingWarnings->push_back(
@@ -2386,32 +2807,55 @@ std::string UsgsAstroLsSensorModel::constructStateFromIsd(const std::string imag
         csm::Warning::DATA_NOT_AVAILABLE,
         "t0_quaternion not in ISD",
         "UsgsAstroFrameSensorModel::constructStateFromIsd()"));
+    MESSAGE_LOG(m_logger, "t0_quaternion not in ISD")
   }
 
   std::vector<double> positions = getSensorPositions(isd, parsingWarnings);
   state["m_positions"] = positions;
   state["m_numPositions"] = positions.size();
+  MESSAGE_LOG(m_logger, "m_positions: {}"
+                        "m_numPositions: {}",
+                        state["m_positions"].dump(),
+                        state["m_numPositions"].dump())
 
   state["m_velocities"] = getSensorVelocities(isd, parsingWarnings);
+  MESSAGE_LOG(m_logger, "m_velocities: {}",
+                        state["m_velocities"].dump())
 
   std::vector<double> quaternions = getSensorOrientations(isd, parsingWarnings);
   state["m_quaternions"] = quaternions;
   state["m_numQuaternions"] = quaternions.size();
+  MESSAGE_LOG(m_logger, "m_quaternions: {}"
+                        "m_numQuaternions: {}",
+                        state["m_quaternions"].dump(),
+                        state["m_numQuaternions"].dump())
 
   state["m_currentParameterValue"] = std::vector<double>(NUM_PARAMETERS, 0.0);
+  MESSAGE_LOG(m_logger, "m_currentParameterValue: {}",
+                        state["m_currentParameterValue"].dump())
 
   // get radii
   state["m_minorAxis"] = getSemiMinorRadius(isd, parsingWarnings);
   state["m_majorAxis"] = getSemiMajorRadius(isd, parsingWarnings);
+  MESSAGE_LOG(m_logger, "m_minorAxis: {}"
+                        "m_majorAxis: {}",
+                        state["m_minorAxis"].dump(), state["m_majorAxis"].dump())
 
   // set identifiers
   state["m_platformIdentifier"] = getPlatformName(isd, parsingWarnings);
   state["m_sensorIdentifier"] = getSensorName(isd, parsingWarnings);
+  MESSAGE_LOG(m_logger, "m_platformIdentifier: {}"
+                        "m_sensorIdentifier: {}",
+                        state["m_platformIdentifier"].dump(),
+                        state["m_sensorIdentifier"].dump())
 
   // get reference_height
   state["m_minElevation"] = getMinHeight(isd, parsingWarnings);
   state["m_maxElevation"] = getMaxHeight(isd, parsingWarnings);
-
+  MESSAGE_LOG(m_logger, "m_minElevation: {}"
+                        "m_maxElevation: {}",
+                        state["m_minElevation"].dump(),
+                        state["m_maxElevation"].dump())
 
   // Default to identity covariance
   state["m_covariance"] =
@@ -2420,8 +2864,24 @@ std::string UsgsAstroLsSensorModel::constructStateFromIsd(const std::string imag
    state["m_covariance"][i * NUM_PARAMETERS + i] = 1.0;
   }
 
+  // Get the optional logging file
+  state["m_logFile"] = getLogFile(isd);
 
    // The state data will still be updated when a sensor model is created since
    // some state data is notin the ISD and requires a SM to compute them.
    return state.dump();
 }
+
+
+//***************************************************************************
+// UsgsAstroLineScannerSensorModel::getLogger
+//***************************************************************************
+std::shared_ptr<spdlog::logger> UsgsAstroLsSensorModel::getLogger() {
+  // MESSAGE_LOG(m_logger, "Getting log pointer, logger is {}",
+  //                             m_logger)
+  return m_logger;
+}
+
+void UsgsAstroLsSensorModel::setLogger(std::shared_ptr<spdlog::logger> logger) {
+  m_logger = logger;
+}