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