Changeset 354 for trunk/src


Ignore:
Timestamp:
02/03/05 15:43:18 (20 years ago)
Author:
kil064
Message:

rework convertFLux stuff using new SDAttr and looking up values
where possible

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/SDMath.cc

    r330 r354  
    5151#include <casa/Quanta/Unit.h>
    5252#include <casa/Quanta/MVEpoch.h>
    53 #include <casa/Quanta/QC.h>
    5453#include <casa/Quanta/MVTime.h>
    5554#include <casa/Utilities/Assert.h>
    … …  
    7978#include "MathUtils.h"
    8079#include "SDDefs.h"
     80#include "SDAttr.h"
    8181#include "SDContainer.h"
    8282#include "SDMemTable.h"
    … …  
    10611061
    10621062
    1063 SDMemTable* SDMath::convertFlux (const SDMemTable& in, Float a, Float eta, Bool doAll) const
     1063SDMemTable* SDMath::convertFlux (const SDMemTable& in, Float D, Float etaAp,
     1064                                 Float JyPerK, Bool doAll) const
    10641065//
    1065 // As it is, this function could be implemented with 'simpleOperate'
    1066 // However, I anticipate that eventually we will look the conversion
    1067 // values up in a Table and apply them in a frequency dependent way,
    1068 // so I have implemented it fully here
     1066// etaAp = aperture efficiency
     1067// D     = geometric diameter (m)
     1068// JyPerK
    10691069//
    10701070{
    … …  
    10721072  SDMemTable* pTabOut = new SDMemTable(in, True);
    10731073
    1074 // FInd out how to convert values into Jy and K (e.g. units might be mJy or mK)
     1074// Find out how to convert values into Jy and K (e.g. units might be mJy or mK)
    10751075// Also automatically find out what we are converting to according to the
    10761076// flux unit
    … …  
    10811081//
    10821082  Bool toKelvin = True;
    1083   Double inFac = 1.0;
     1083  Double cFac = 1.0;   
    10841084  if (fluxUnit==JY) {
    10851085     cerr << "Converting to K" << endl;
    … …  
    10871087     Quantum<Double> t(1.0,fluxUnit);
    10881088     Quantum<Double> t2 = t.get(JY);
    1089      inFac = (t2 / t).getValue();
     1089     cFac = (t2 / t).getValue();               // value to Jy
    10901090//
    10911091     toKelvin = True;
    … …  
    10961096     Quantum<Double> t(1.0,fluxUnit);
    10971097     Quantum<Double> t2 = t.get(K);
    1098      inFac = (t2 / t).getValue();
     1098     cFac = (t2 / t).getValue();              // value to K
    10991099//
    11001100     toKelvin = False;
    … …  
    11051105  pTabOut->putSDHeader(sh);
    11061106
    1107 // Compute conversion factor. 'a' and 'eta' are really frequency, time and 
    1108 // telescope dependent and should be looked// up in a table
    1109 
    1110   Float factor = 2.0 * inFac * 1.0e-7 * 1.0e26 *
    1111                  QC::k.getValue(Unit(String("erg/K"))) / a / eta;
    1112   if (toKelvin) {
    1113     factor = 1.0 / factor;
    1114   }
    1115   cerr << "Applying conversion factor = " << factor << endl;
    1116 
    1117 // Generate correction vector.  Apply same factor regardless
    1118 // of beam/pol/IF.  This will need to change somewhen.
    1119 
    1120   Vector<Float> factors(in.nRow(), factor);
    1121 
    1122 // Correct
    1123 
    1124   correctFromVector (pTabOut, in, doAll, factors);
     1107// Make sure input values are converted to either Jy or K first...
     1108
     1109  Float factor = cFac;
     1110
     1111// Select method
     1112
     1113  if (JyPerK>0.0) {
     1114     factor *= JyPerK;
     1115     if (toKelvin) factor = 1.0 / JyPerK;
     1116//
     1117     cerr << "Applying supplied conversion factor = " << factor << endl;
     1118     Vector<Float> factors(in.nRow(), factor);
     1119     correctFromVector (pTabOut, in, doAll, factors);
     1120  } else if (etaAp>0.0) {
     1121     factor *= SDAttr::findJyPerKFac (etaAp, D);
     1122     if (toKelvin) {
     1123        factor = 1.0 / factor;
     1124     }
     1125//
     1126     cerr << "Applying supplied conversion factor = " << factor << endl;
     1127     Vector<Float> factors(in.nRow(), factor);
     1128     correctFromVector (pTabOut, in, doAll, factors);
     1129  } else {
     1130
     1131// OK now we must deal with automatic look up of values.
     1132// We must also deal with the fact that the factors need
     1133// to be computed per IF and may be different and may
     1134// change per integration.
     1135
     1136     cerr << "Looking up conversion factors" << endl;
     1137     convertBrightnessUnits (pTabOut, in, toKelvin, cFac, doAll);
     1138  }
    11251139//
    11261140  return pTabOut;
    11271141}
     1142
     1143
     1144
    11281145
    11291146
    … …  
    12531270// 'private' functions
    12541271
     1272void SDMath::convertBrightnessUnits (SDMemTable* pTabOut, const SDMemTable& in,
     1273                                     Bool toKelvin, Float cFac, Bool doAll) const
     1274{
     1275
     1276// Get header
     1277
     1278   SDHeader sh = in.getSDHeader();
     1279   const uInt nChan = sh.nchan;
     1280
     1281// Get instrument
     1282
     1283   Bool throwIt = True;
     1284   Instrument inst = SDMemTable::convertInstrument (sh.antennaname, throwIt);
     1285
     1286// Get Diameter (m)
     1287
     1288   SDAttr sdAtt;
     1289
     1290// Get epoch of first row
     1291
     1292   MEpoch dateObs = in.getEpoch(0);
     1293
     1294// Generate a Vector of correction factors. One per FreqID
     1295
     1296   SDFrequencyTable sdft = in.getSDFreqTable();
     1297   Vector<uInt> freqIDs;
     1298//
     1299   Vector<Float> freqs(sdft.length());
     1300   for (uInt i=0; i<sdft.length(); i++) {
     1301      freqs(i) = (nChan/2 - sdft.referencePixel(i))*sdft.increment(i) + sdft.referenceValue(i);
     1302   }
     1303//
     1304   Vector<Float> JyPerK = sdAtt.JyPerK(inst, dateObs, freqs);
     1305   Vector<Float> factors = cFac * JyPerK;
     1306   if (toKelvin) factors = Float(1.0) / factors;
     1307
     1308// For operations only on specified cursor location
     1309
     1310   IPosition start, end;
     1311   getCursorLocation(start, end, in);
     1312   const uInt ifAxis = in.getIF();
     1313
     1314// Iteration axes
     1315
     1316   IPosition axes(asap::nAxes-1,0);
     1317   for (uInt i=0,j=0; i<asap::nAxes; i++) {
     1318      if (i!=asap::IFAxis) {
     1319         axes(j++) = i;
     1320      }
     1321   }
     1322
     1323// Loop over rows and apply correction factor
     1324
     1325   Float factor = 1.0; 
     1326   const uInt axis = asap::ChanAxis;
     1327   for (uInt i=0; i < in.nRow(); ++i) {
     1328
     1329// Get data
     1330
     1331      MaskedArray<Float> dataIn(in.rowAsMaskedArray(i));
     1332      Array<Float>& values = dataIn.getRWArray();
     1333
     1334// Get SDCOntainer
     1335
     1336      SDContainer sc = in.getSDContainer(i);
     1337
     1338// Get FreqIDs
     1339
     1340      freqIDs = sc.getFreqMap();
     1341
     1342// Now the conversion factor depends only upon frequency
     1343// So we need to iterate through by IF only giving
     1344// us BEAM/POL/CHAN cubes
     1345
     1346      if (doAll) {
     1347         ArrayIterator<Float> itIn(values, axes);
     1348         uInt ax = 0;
     1349         while (!itIn.pastEnd()) {
     1350           itIn.array() *= factors(freqIDs(ax));         // Writes back to dataIn
     1351           itIn.next();
     1352         }
     1353      } else {
     1354         MaskedArray<Float> dataIn2 = dataIn(start,end);  // reference
     1355         dataIn2 *= factors(freqIDs(ifAxis));
     1356      }
     1357
     1358// Write out
     1359
     1360      putDataInSDC(sc, dataIn.getArray(), dataIn.getMask());
     1361//
     1362      pTabOut->putSDContainer(sc);
     1363   }
     1364}
     1365
     1366
    12551367
    12561368SDMemTable* SDMath::frequencyAlign (const SDMemTable& in,
    … …  
    18001912   }
    18011913}
    1802 
    1803 

  • trunk/src/SDMath.h

    r330 r354  
    102102
    103103// Flux conversion between Jansky and Kelvin
    104    SDMemTable* convertFlux (const SDMemTable& in, casa::Float area,
    105                             casa::Float eta, casa::Bool doAll) const;
     104   SDMemTable* convertFlux (const SDMemTable& in, casa::Float D, casa::Float etaAp,
     105                            casa::Float JyPerK, casa::Bool doAll) const;
    106106
    107107// Gain-elevation correction
    … …  
    142142                   casa::Bool useMask, WeightType wtType) const;
    143143
     144// Work out conversion factor for converting Jy<->K per IF per row and apply
     145   void convertBrightnessUnits (SDMemTable* pTabOut, const SDMemTable& in,
     146                                casa::Bool toKelvin, casa::Float sFac, casa::Bool doAll) const;
     147
     148// Convert weight string to enum value
     149
     150   void convertWeightString (WeightType& wt, const casa::String& weightStr) const;
     151
     152// Convert interpolation type string
     153//   void convertInterpString(casa::Int& type, const casa::String& interp) const;
     154   void convertInterpString(casa::InterpolateArray1D<casa::Double,casa::Float>::InterpolationMethod& method, 
     155                             const casa::String& interp) const;
     156
     157// Correct data from an ascii Table
     158   void correctFromAsciiTable(SDMemTable* pTabOut, const SDMemTable& in,
     159                              const casa::String& fileName,
     160                              const casa::String& col0, const casa::String& col1,
     161                              const casa::String& methodStr, casa::Bool doAll,
     162                              const casa::Vector<casa::Float>& xOut) const;
     163
     164// Correct data from a Table
     165   void correctFromTable(SDMemTable* pTabOut, const SDMemTable& in, const casa::Table& tTable,
     166                         const casa::String& col0, const casa::String& col1,
     167                         const casa::String& methodStr, casa::Bool doAll,
     168                         const casa::Vector<casa::Float>& xOut) const;
     169
     170// Correct data from a Vector
     171   void correctFromVector (SDMemTable* pTabOut, const SDMemTable& in,
     172                           casa::Bool doAll, const casa::Vector<casa::Float>& factor) const;
     173
     174// Convert time String to Epoch
     175   casa::MEpoch epochFromString (const casa::String& str, casa::MEpoch::Types timeRef) const;
     176
    144177// Function to fill Scan Container when averaging in time
    145178
    … …  
    151184                const casa::Vector<casa::uInt>& freqID) const;
    152185
    153 // Put the data and mask into the SDContainer
    154    void putDataInSDC (SDContainer& sc, const casa::Array<casa::Float>& data,
    155                       const casa::Array<casa::Bool>& mask) const;
    156 
    157 // Function to normalize data when averaging in time
    158 
    159   void normalize (casa::MaskedArray<casa::Float>& data,
    160                   const casa::Array<casa::Float>& sumSq,
    161                   const casa::Array<casa::Float>& nPts,
    162                   WeightType wtType, casa::Int axis, casa::Int nAxes) const;
    163 
    164 // Function to get the current cursor location
    165    void getCursorLocation (casa::IPosition& start, casa::IPosition& end,
    166                            const SDMemTable& in) const;
    167 
    168 // Convert weight string to enum value
    169 
    170    void convertWeightString (WeightType& wt, const casa::String& weightStr) const;
    171 
    172 // Convert interpolation type string
    173 //   void convertInterpString(casa::Int& type, const casa::String& interp) const;
    174    void convertInterpString(casa::InterpolateArray1D<casa::Double,casa::Float>::InterpolationMethod& method, 
    175                              const casa::String& interp) const;
    176 
    177 // Correct data from an ascii Table
    178    void correctFromAsciiTable(SDMemTable* pTabOut, const SDMemTable& in,
    179                               const casa::String& fileName,
    180                               const casa::String& col0, const casa::String& col1,
    181                               const casa::String& methodStr, casa::Bool doAll,
    182                               const casa::Vector<casa::Float>& xOut) const;
    183 
    184 // Correct data from a Table
    185    void correctFromTable(SDMemTable* pTabOut, const SDMemTable& in, const casa::Table& tTable,
    186                          const casa::String& col0, const casa::String& col1,
    187                          const casa::String& methodStr, casa::Bool doAll,
    188                          const casa::Vector<casa::Float>& xOut) const;
    189 
    190 // Correct data from a Vector
    191    void correctFromVector (SDMemTable* pTabOut, const SDMemTable& in,
    192                            casa::Bool doAll, const casa::Vector<casa::Float>& factor) const;
    193 
    194 // Read ascii file into a Table
    195 
    196    casa::Table readAsciiFile (const casa::String& fileName) const;
     186// Format EPoch
     187   casa::String formatEpoch(const casa::MEpoch& epoch)  const;
     188
     189// Align in Frequency
     190   SDMemTable* frequencyAlign (const SDMemTable& in,
     191                              casa::MFrequency::Types system,
     192                              const casa::String& timeRef,
     193                              const casa::String& method) const;
    197194
    198195// Generate frequency aligners
    … …  
    213210                               const casa::ROArrayColumn<casa::uInt>& fqIDCol) const;
    214211
    215 // Align in Frequency
    216    SDMemTable* frequencyAlign (const SDMemTable& in,
    217                               casa::MFrequency::Types system,
    218                               const casa::String& timeRef,
    219                               const casa::String& method) const;
    220 
    221 // Convert time String to Epoch
    222    casa::MEpoch epochFromString (const casa::String& str, casa::MEpoch::Types timeRef) const;
    223 
    224 // Format EPoch
    225    casa::String formatEpoch(const casa::MEpoch& epoch)  const;
     212// Function to get the current cursor location
     213   void getCursorLocation (casa::IPosition& start, casa::IPosition& end,
     214                           const SDMemTable& in) const;
     215
     216// Function to normalize data when averaging in time
     217
     218  void normalize (casa::MaskedArray<casa::Float>& data,
     219                  const casa::Array<casa::Float>& sumSq,
     220                  const casa::Array<casa::Float>& nPts,
     221                  WeightType wtType, casa::Int axis, casa::Int nAxes) const;
     222
     223// Put the data and mask into the SDContainer
     224   void putDataInSDC (SDContainer& sc, const casa::Array<casa::Float>& data,
     225                      const casa::Array<casa::Bool>& mask) const;
     226
     227// Read ascii file into a Table
     228
     229   casa::Table readAsciiFile (const casa::String& fileName) const;
    226230};
    227231
Note: See TracChangeset for help on using the changeset viewer.