USGS

Isis 3.0 Application Source Code Reference

Home

MdisCalUtils.h

Go to the documentation of this file.
00001 #ifndef MdisCalUtils_h
00002 #define MdisCalUtils_h
00003 /**
00004  * @file
00005  * $Revision: 6392 $
00006  * $Date: 2015-10-09 14:29:19 -0700 (Fri, 09 Oct 2015) $
00007  *
00008  *   Unless noted otherwise, the portions of Isis written by the USGS are
00009  *   public domain. See individual third-party library and package descriptions
00010  *   for intellectual property information, user agreements, and related
00011  *   information.
00012  *
00013  *   Although Isis has been used by the USGS, no warranty, expressed or
00014  *   implied, is made by the USGS as to the accuracy and functioning of such
00015  *   software and related material nor shall the fact of distribution
00016  *   constitute any such warranty, and no responsibility is assumed by the
00017  *   USGS in connection therewith.
00018  *
00019  *   For additional information, launch
00020  *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html
00021  *   in a browser or see the Privacy & Disclaimers page on the Isis website,
00022  *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on
00023  *   http://www.usgs.gov/privacy.html.
00024  */
00025 #include <cmath>
00026 #include <string>
00027 #include <vector>
00028 
00029 #include "CSVReader.h"
00030 #include "FileName.h"
00031 #include "IString.h"
00032 #include "Spice.h"
00033 /**
00034  * @author ????-??-?? Kris Becker
00035  *
00036  * @internal
00037  *   @history 2015-08-31 Jeannie Backer - Changed method name from loadContaminationEvent() to
00038  *                           loadEmpiricalCorrection(). Brought code closer to coding standards.
00039  */
00040 namespace Isis {
00041   /**
00042    * @brief Helper function to convert values to doubles
00043    *
00044    * @param T Type of value to convert
00045    * @param value Value to convert
00046    *
00047    * @return double Converted value
00048    */
00049   template <typename T> double toDouble(const T &value) {
00050     return toDouble(QString(value).trimmed());
00051   }
00052 
00053   template <typename T> int toInteger(const T &value) {
00054     return toInt(QString(value).trimmed());
00055   }
00056 
00057   template <typename T> inline T MIN(const T &A, const T &B) {
00058     return (((A) < (B)) ? (A) : (B));
00059   }
00060 
00061   template <typename T> inline T MAX(const T &A, const T &B) {
00062     return (((A) > (B)) ? (A) : (B));
00063   }
00064 
00065   inline QString quote(const QString &value) {
00066     if (value.isEmpty()) return (value);
00067     if (value[0] == '"') return (value);
00068     return (QString('"' + value + '"'));
00069   }
00070 
00071   /**
00072    * @brief Load required NAIF kernels required for timing needs
00073    *
00074    * This method maintains the loading of kernels for MESSENGER timing and
00075    * planetary body ephemerides to support time and relative positions of planet
00076    * bodies.
00077    */
00078   static void loadNaifTiming() {
00079     static bool naifLoaded = false;
00080     if (!naifLoaded) {
00081 //  Load the NAIF kernels to determine timing data
00082       Isis::FileName leapseconds("$base/kernels/lsk/naif????.tls");
00083       leapseconds = leapseconds.highestVersion();
00084 
00085       Isis::FileName sclk("$messenger/kernels/sclk/messenger_????.tsc");
00086       sclk = sclk.highestVersion();
00087 
00088       Isis::FileName pck("$base/kernels/spk/de???.bsp");
00089       pck = pck.highestVersion();
00090 
00091 //  Load the kernels
00092       QString leapsecondsName(leapseconds.expanded());
00093       QString sclkName(sclk.expanded());
00094       QString pckName(pck.expanded());
00095       furnsh_c(leapsecondsName.toAscii().data());
00096       furnsh_c(sclkName.toAscii().data());
00097       furnsh_c(pckName.toAscii().data());
00098 
00099 //  Ensure it is loaded only once
00100       naifLoaded = true;
00101     }
00102     return;
00103   }
00104 
00105   /**
00106    * @brief Computes the distance from the Sun to the observed body
00107    *
00108    * This method requires the appropriate NAIK kernels to be loaded that
00109    * provides instrument time support, leap seconds and planet body ephemeris.
00110    *
00111    * @return double Distance in AU between Sun and observed body
00112    */
00113   static bool sunDistanceAU(const QString &scStartTime,
00114                             const QString &target,
00115                             double &sunDist) {
00116 
00117     //  Ensure NAIF kernels are loaded
00118     loadNaifTiming();
00119     sunDist = 1.0;
00120 
00121     //  Determine if the target is a valid NAIF target
00122     SpiceInt tcode;
00123     SpiceBoolean found;
00124     bodn2c_c(target.toAscii().data(), &tcode, &found);
00125     if (!found) return (false);
00126 
00127     //  Convert starttime to et
00128     double obsStartTime;
00129     scs2e_c(-236, scStartTime.toAscii().data(), &obsStartTime);
00130 
00131     //  Get the vector from target to sun and determine its length
00132     double sunv[3];
00133     double lt;
00134     spkpos_c(target.toAscii().data(), obsStartTime, "J2000", "LT+S", "sun",
00135                     sunv, &lt);
00136     double sunkm = vnorm_c(sunv);
00137 
00138     //  Return in AU units
00139     sunDist = sunkm / 1.49597870691E8;
00140     return (true);
00141   }
00142 
00143   std::vector<double> loadWACCSV(const QString &fname, int filter,
00144                                  int nvalues, bool header=true, int skip=0) {
00145     //  Open the CSV file
00146     FileName csvfile(fname);
00147     CSVReader csv(csvfile.expanded(), header, skip);
00148     for (int i = 0; i < csv.rows(); i++) {
00149       CSVReader::CSVAxis row = csv.getRow(i);
00150       if (toInteger(row[0]) == filter) {
00151         if ((row.dim1() - 1) < nvalues) {
00152           QString mess = "Number values (" + QString(row.dim1() - 1) +
00153                              ") in file " + fname +
00154                              " less than number requested (" +
00155                              QString(nvalues) + ")!";
00156           throw IException(IException::User, mess, _FILEINFO_);
00157         }
00158 
00159         std::vector<double> rsp;
00160         for (int c = 0; c < nvalues; c++) {
00161           rsp.push_back(toDouble(row[1+c]));
00162         }
00163         return (rsp);
00164       }
00165     }
00166 
00167     // If it reaches here, the filter was not found
00168     std::ostringstream mess;
00169     mess << "CSV Vector MDIS filter " << filter <<  ", not found in file "
00170          << fname << "!";
00171     throw IException(IException::User, mess.str(), _FILEINFO_);
00172   }
00173 
00174 
00175   std::vector<double> loadNACCSV(const QString &fname, int nvalues,
00176                                  bool header=true, int skip=0) {
00177     //  Open the CSV file
00178     FileName csvfile(fname);
00179     CSVReader csv(csvfile.expanded(), header, skip);
00180     CSVReader::CSVAxis row = csv.getRow(0);
00181     if (row.dim1() < nvalues) {
00182       QString mess = "Number values (" + QString(row.dim1()) +
00183                          ") in file " + fname + " less than number requested (" +
00184                          QString(nvalues) + ")!";
00185       throw IException(IException::User, mess, _FILEINFO_);
00186 
00187     }
00188     std::vector<double> rsp;
00189     for (int i = 0; i < nvalues; i++) {
00190       rsp.push_back(toDouble(row[i]));
00191     }
00192     return (rsp);
00193   }
00194 
00195 
00196   std::vector<double> loadResponsivity(bool isNAC, bool binned, int filter,
00197                                        QString &fname) {
00198 
00199     FileName resfile(fname);
00200     if (fname.isEmpty()) {
00201       QString camstr = (isNAC) ? "NAC" : "WAC";
00202       QString binstr = (binned)       ? "_BINNED" : "_NOTBIN";
00203       QString base   = "$messenger/calibration/RESPONSIVITY/";
00204       resfile = base + "MDIS" + camstr + binstr + "_RESP_?.TAB";
00205       resfile = resfile.highestVersion();
00206       fname = resfile.originalPath() + "/" + resfile.name();
00207     }
00208 
00209     // Unfortunately NAC has a slightly different format so must do it
00210     //  explicitly
00211     if (isNAC) {
00212       return (loadNACCSV(fname, 4, false, 0));
00213     }
00214     else {
00215       // Load the WAC parameters
00216       return (loadWACCSV(fname, filter, 4, false, 0));
00217     }
00218   }
00219 
00220 
00221   std::vector<double> loadSolarIrr(bool isNAC, bool binned, int filter,
00222                                    QString &fname)  {
00223 
00224     FileName solfile(fname);
00225     if (fname.isEmpty()) {
00226       QString camstr = (isNAC) ? "NAC" : "WAC";
00227       QString base   = "$messenger/calibration/SOLAR/";
00228       solfile = base + "MDIS" + camstr + "_SOLAR_?.TAB";
00229       solfile = solfile.highestVersion();
00230       fname = solfile.originalPath() + "/" + solfile.name();
00231     }
00232 
00233     if (isNAC) {
00234       return (loadNACCSV(fname, 3, false, 0));
00235     }
00236     else {
00237       return (loadWACCSV(fname, filter, 3, false, 0));
00238     }
00239   }
00240 
00241   double loadSmearComponent(bool isNAC, int filter, QString &fname) {
00242 
00243     FileName smearfile(fname);
00244     if (fname.isEmpty()) {
00245       QString camstr = (isNAC) ? "NAC" : "WAC";
00246       QString base   = "$messenger/calibration/smear/";
00247       smearfile = base + "MDIS" + camstr + "_FRAME_TRANSFER_??.TAB";
00248       smearfile = smearfile.highestVersion();
00249       fname = smearfile.originalPath() + "/" + smearfile.name();
00250     }
00251 
00252     std::vector<double> smear;
00253     if (isNAC) {
00254       smear = loadNACCSV(fname, 1, false, 0);
00255     }
00256     else {
00257       smear = loadWACCSV(fname, filter, 1, false, 0);
00258     }
00259     return (smear[0]);
00260   }
00261 
00262 /**
00263  * @brief Load and retrieve empirical correction factor 
00264  *  
00265  * This function determines the empirical correction factor for changes that 
00266  * that occured on the spacecraft after Mercury orbit insertion.  The 
00267  * affected dates are May 24, 2011 to January 3, 2012. 
00268  *  
00269  * The table of correction factors is expected to be stored in 
00270  * $messenger/calibration/events/event_table_ratioed_v?.txt.  However, the 
00271  * caller may provide a table that conforms to the expected format.  The 
00272  * expected format for the empirical correction file is a comma separated value
00273  * (CSV) table that contains 13 columns of data per row.  The first column is 
00274  * the UTC time during the event. The next 12 columns contain multiplicative 
00275  * correction factors for each WAC filter (NAC correction factors are not 
00276  * provided). These factors are expected to be around 1.0 (the default) as it 
00277  * is expected to directly scale DN values. 
00278  *  
00279  * Below is the expected mapping of column indexes to filter numbers as 
00280  * specfied in the BandBin/Number keyword from MDIS ISIS cube labels. Index is 
00281  * the column index from each row for a given filter, Number is the value of 
00282  * the BandBin/Number keyword from the label designating the filter number 
00283  * (corresponding to the filter parameter passed to this routine) and Letter is
00284  * the filter letter designation used in the last alpha numeric character in 
00285  * MDIS filenames: 
00286  *  
00287  *  Index   Number   Letter   Wavelength
00288  *    1       6        F         430 nm
00289  *    2       3        C         480 nm
00290  *    3       4        D         560 nm
00291  *    4       5        E         630 nm
00292  *    5       7        G         750 nm
00293  *    6      12        L         830 nm
00294  *    7      10        J         900 nm
00295  *    8       9        I        1000 nm
00296  *    9       1        A        Filter 1 (700 nm)
00297  *   10       2        B        Filter 2 (clear)
00298  *   11       8        H        Filter 8 (950 nm)
00299  *   12      11        K        Filter 11 (1010 nm)
00300  *  
00301  * The UTC dates in the first column are assumed to be strictly increasing in 
00302  * time.  The initial table (*_v2) contains dates that span the complete 
00303  * expected timeframe of the mission (launch at 2004-08-04T10:00:00.000000, 
00304  * termination at 2015-01-03T09:00:00.000000). 
00305  *  
00306  * The spacecraft clock time is provided as input (scStartTime) to this 
00307  * function.  This value is converted to ET (SCET) and used to determine the 
00308  * corresponding event time in the first column of the table.  The first table 
00309  * column time is represented in UTC time.  This time is converted to ET and 
00310  * then compared with the start time in ET. 
00311  *  
00312  * The algorithm searches linearly through the table essentially storing the 
00313  * time slot prior to the SCET and the next occuring one.  Ultimately, the 
00314  * factor returned by the algorithm is the one whose event time is closest to 
00315  * the SCET. 
00316  *  
00317  * The empirical correction model and algorithm was developed by
00318  * Mary Ruth Keller of JHA/APL.
00319  *
00320  * @author Kris Becker - 10/23/2012
00321  * 
00322  * @param scStartTime - Start time of the image in SCLK format
00323  * @param filter      - WAC filter number to return event correction factor for
00324  * @param ename       - Returns the name of the event table file if not 
00325  *                      provided by caller.  If non-empty string is passed by
00326  *                      caller, it is assumed to be a fully qualified filename
00327  *                      of the event table.
00328  * @param eDate       - Returns the UTC date entry of the selected correction 
00329  *                      event factor
00330  * 
00331  * @return double     - Event correction factor at the selected time to apply 
00332  *                      to WAC filter data.
00333  */
00334  double loadEmpiricalCorrection(const QString &scStartTime, const int filter, 
00335                                 QString &ename, QString &eDate) {
00336 
00337    //  This table maps the filter number extracted from BandBin/Number keyword
00338    //  to the columns (index) in the empirical correction table
00339    const int filterMap[12] = { 6, 3, 4, 5, 7, 12, 10, 9, 1, 2, 8, 11 };
00340 
00341    //  Find the WAC filter column index
00342    int ncols = sizeof(filterMap)/sizeof(filterMap[0]);
00343    int column = -1;
00344    for (int c = 0; c < ncols; c++) {
00345      if (filterMap[c] == filter) {
00346        column = c + 1;  // indexes start after 1st (time) column
00347        break;
00348      }
00349    }
00350 
00351    // Ensure we have a valid filter number
00352    if (column <= 0) {
00353      std::ostringstream mess;
00354       mess << "Invalid MDIS WAC filter number (" << filter 
00355            <<  " - range:1-12) for determining index into empirical correction table.";
00356       throw IException(IException::User, mess.str(), _FILEINFO_);     
00357    }
00358 
00359    //  File name not provided by caller.  Determine the event table name
00360     if ( ename.isEmpty() ) {
00361       FileName eventfile("$messenger/calibration/events/event_table_ratioed_v?.txt");
00362       eventfile = eventfile.highestVersion();
00363       ename = eventfile.originalPath() + "/" + eventfile.name();
00364     }
00365 
00366     //  Open/read the CSV empirical correction file
00367     FileName csvfile(ename);
00368     const bool header = false;  // No header in file
00369     const int skip = 0;         // No lines to skip to data
00370     const int nvalues = 13;     // Expected columns in table
00371     CSVReader csv(csvfile.expanded(), header, skip);
00372     if (csv.columns() < nvalues) {  // All rows should have same # columns 
00373       QString mess = "Number values (" + QString(csv.columns()) +
00374                          ") in file " + ename + " less than number requested (" +
00375                          QString(nvalues) + ")!";
00376       throw IException(IException::User, mess, _FILEINFO_);
00377     }
00378 
00379     // Ensure NAIF kernels are loaded for NAIF time computations
00380     loadNaifTiming();
00381 
00382     //  Convert s/c clock start time to et
00383     double obsStartTime;
00384     scs2e_c(-236, scStartTime.toAscii().data(), &obsStartTime);
00385 
00386     // Set initial conditions and loop through all rows in the event table
00387     double evalue = 1.0;
00388     eDate = "N/A";  // Will attain a valid time on guaranteed first pass
00389     double preEventTime = 0.0;
00390     for (int i = 0; i < csv.rows(); i++) {
00391       CSVReader::CSVAxis eRow = csv.getRow(i);
00392       QString utcTime = eRow[0];
00393       double eTime;
00394       utc2et_c(utcTime.toAscii().data(), &eTime);
00395 
00396       // If current time is greater than start time this is the post event case
00397       if (eTime > obsStartTime) {
00398         //  Get closest pre or post event correction factor
00399         if ( fabs(obsStartTime-preEventTime) > fabs(eTime-obsStartTime) ) {
00400           //  Post-event time closer to SCLK than Pre-event time
00401           eDate = utcTime;
00402           evalue = toDouble(eRow[column]);
00403         }
00404 
00405         break;  //  Terminate loop and return
00406       }
00407 
00408       // Record pre-event time slot - Sets return variables as well
00409       eDate = utcTime;
00410       preEventTime = eTime;
00411       evalue = toDouble(eRow[column]);
00412     }
00413  
00414     // Return the factor
00415     return (evalue);
00416   }
00417 
00418 
00419 };
00420 #endif