Changeset 354 for trunk/src/SDMath.cc

Timestamp:

02/03/05 15:43:18 (19 years ago)

Author:

kil064

Message:

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

File:

• trunk/src/SDMath.cc (modified) (10 diffs)

trunk/src/SDMath.cc

@@ -51 +51 @@
 #include <casa/Quanta/Unit.h>
 #include <casa/Quanta/MVEpoch.h>
-#include <casa/Quanta/QC.h>
 #include <casa/Quanta/MVTime.h>
 #include <casa/Utilities/Assert.h>
@@ -79 +78 @@
 #include "MathUtils.h"
 #include "SDDefs.h"
+#include "SDAttr.h"
 #include "SDContainer.h"
 #include "SDMemTable.h"
@@ -1061 +1061 @@
 
 
-SDMemTable* SDMath::convertFlux (const SDMemTable& in, Float a, Float eta, Bool doAll) const
+SDMemTable* SDMath::convertFlux (const SDMemTable& in, Float D, Float etaAp, 
+                                 Float JyPerK, Bool doAll) const
 // 
-// As it is, this function could be implemented with 'simpleOperate'
-// However, I anticipate that eventually we will look the conversion
-// values up in a Table and apply them in a frequency dependent way,
-// so I have implemented it fully here
+// etaAp = aperture efficiency
+// D     = geometric diameter (m)
+// JyPerK 
 //
 {
@@ -1072 +1072 @@
   SDMemTable* pTabOut = new SDMemTable(in, True);
 
-// FInd out how to convert values into Jy and K (e.g. units might be mJy or mK)
+// Find out how to convert values into Jy and K (e.g. units might be mJy or mK)
 // Also automatically find out what we are converting to according to the
 // flux unit
@@ -1081 +1081 @@
 //
   Bool toKelvin = True;
-  Double inFac = 1.0;
+  Double cFac = 1.0;    
   if (fluxUnit==JY) {
      cerr << "Converting to K" << endl;
@@ -1087 +1087 @@
      Quantum<Double> t(1.0,fluxUnit);
      Quantum<Double> t2 = t.get(JY);
-     inFac = (t2 / t).getValue();
+     cFac = (t2 / t).getValue();               // value to Jy
 //
      toKelvin = True;
@@ -1096 +1096 @@
      Quantum<Double> t(1.0,fluxUnit);
      Quantum<Double> t2 = t.get(K);
-     inFac = (t2 / t).getValue();
+     cFac = (t2 / t).getValue();              // value to K
 //
      toKelvin = False;
@@ -1105 +1105 @@
   pTabOut->putSDHeader(sh);
 
-// Compute conversion factor. 'a' and 'eta' are really frequency, time and  
-// telescope dependent and should be looked// up in a table
-
-  Float factor = 2.0 * inFac * 1.0e-7 * 1.0e26 * 
-                 QC::k.getValue(Unit(String("erg/K"))) / a / eta;
-  if (toKelvin) {
-    factor = 1.0 / factor;
-  }
-  cerr << "Applying conversion factor = " << factor << endl;
-
-// Generate correction vector.  Apply same factor regardless
-// of beam/pol/IF.  This will need to change somewhen.
-
-  Vector<Float> factors(in.nRow(), factor);
-
-// Correct
-
-  correctFromVector (pTabOut, in, doAll, factors);
+// Make sure input values are converted to either Jy or K first...
+
+  Float factor = cFac;
+
+// Select method
+
+  if (JyPerK>0.0) {
+     factor *= JyPerK;
+     if (toKelvin) factor = 1.0 / JyPerK;
+//
+     cerr << "Applying supplied conversion factor = " << factor << endl;
+     Vector<Float> factors(in.nRow(), factor);
+     correctFromVector (pTabOut, in, doAll, factors);
+  } else if (etaAp>0.0) {
+     factor *= SDAttr::findJyPerKFac (etaAp, D);
+     if (toKelvin) {
+        factor = 1.0 / factor;
+     }
+//
+     cerr << "Applying supplied conversion factor = " << factor << endl;
+     Vector<Float> factors(in.nRow(), factor);
+     correctFromVector (pTabOut, in, doAll, factors);
+  } else {
+
+// OK now we must deal with automatic look up of values.
+// We must also deal with the fact that the factors need
+// to be computed per IF and may be different and may
+// change per integration.
+
+     cerr << "Looking up conversion factors" << endl;
+     convertBrightnessUnits (pTabOut, in, toKelvin, cFac, doAll);
+  }
 //
   return pTabOut;
 }
+
+
+
 
 
@@ -1253 +1270 @@
 // 'private' functions
 
+void SDMath::convertBrightnessUnits (SDMemTable* pTabOut, const SDMemTable& in, 
+                                     Bool toKelvin, Float cFac, Bool doAll) const
+{
+
+// Get header
+
+   SDHeader sh = in.getSDHeader();
+   const uInt nChan = sh.nchan;
+
+// Get instrument
+
+   Bool throwIt = True;
+   Instrument inst = SDMemTable::convertInstrument (sh.antennaname, throwIt);
+
+// Get Diameter (m)
+
+   SDAttr sdAtt;
+
+// Get epoch of first row
+
+   MEpoch dateObs = in.getEpoch(0);
+
+// Generate a Vector of correction factors. One per FreqID
+
+   SDFrequencyTable sdft = in.getSDFreqTable();
+   Vector<uInt> freqIDs;
+//
+   Vector<Float> freqs(sdft.length());
+   for (uInt i=0; i<sdft.length(); i++) {
+      freqs(i) = (nChan/2 - sdft.referencePixel(i))*sdft.increment(i) + sdft.referenceValue(i);
+   }
+//
+   Vector<Float> JyPerK = sdAtt.JyPerK(inst, dateObs, freqs);
+   Vector<Float> factors = cFac * JyPerK;
+   if (toKelvin) factors = Float(1.0) / factors;
+
+// For operations only on specified cursor location
+
+   IPosition start, end;
+   getCursorLocation(start, end, in);
+   const uInt ifAxis = in.getIF();
+
+// Iteration axes
+
+   IPosition axes(asap::nAxes-1,0);
+   for (uInt i=0,j=0; i<asap::nAxes; i++) {
+      if (i!=asap::IFAxis) {
+         axes(j++) = i;
+      }
+   }
+
+// Loop over rows and apply correction factor
+
+   Float factor = 1.0;  
+   const uInt axis = asap::ChanAxis;
+   for (uInt i=0; i < in.nRow(); ++i) {
+
+// Get data
+
+      MaskedArray<Float> dataIn(in.rowAsMaskedArray(i));
+      Array<Float>& values = dataIn.getRWArray();
+
+// Get SDCOntainer
+
+      SDContainer sc = in.getSDContainer(i);
+
+// Get FreqIDs
+
+      freqIDs = sc.getFreqMap();
+
+// Now the conversion factor depends only upon frequency
+// So we need to iterate through by IF only giving 
+// us BEAM/POL/CHAN cubes
+
+      if (doAll) {
+         ArrayIterator<Float> itIn(values, axes);
+         uInt ax = 0;
+         while (!itIn.pastEnd()) {
+           itIn.array() *= factors(freqIDs(ax));         // Writes back to dataIn
+           itIn.next();
+         }
+      } else {
+         MaskedArray<Float> dataIn2 = dataIn(start,end);  // reference
+         dataIn2 *= factors(freqIDs(ifAxis));
+      }
+
+// Write out
+
+      putDataInSDC(sc, dataIn.getArray(), dataIn.getMask());
+//
+      pTabOut->putSDContainer(sc);
+   }
+}
+
+
 
 SDMemTable* SDMath::frequencyAlign (const SDMemTable& in,
@@ -1800 +1912 @@
    }
 }
-
-