Changeset 822
- Timestamp:
- 02/17/06 10:05:16 (18 years ago)
- Location:
- trunk/src
- Files:
-
- 2 edited
-
STWriter.cpp (modified) (8 diffs)
-
STWriter.h (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/STWriter.cpp
r820 r822 1 1 //#--------------------------------------------------------------------------- 2 //# S DWriter.cc: ASAP class to write out single dish spectra.2 //# STWriter.cc: ASAP class to write out single dish spectra. 3 3 //#--------------------------------------------------------------------------- 4 4 //# Copyright (C) 2004 … … 49 49 #include "SDContainer.h" 50 50 #include "SDMemTable.h" 51 #include "S DWriter.h"51 #include "STWriter.h" 52 52 #include "SDFITSImageWriter.h" 53 53 #include "SDAsciiWriter.h" … … 57 57 using namespace asap; 58 58 59 //--------------------------------------------------------- S DWriter::SDWriter59 //--------------------------------------------------------- STWriter::STWriter 60 60 61 61 // Default constructor. 62 62 63 S DWriter::SDWriter(const std::string &format)63 STWriter::STWriter(const std::string &format) 64 64 { 65 65 cFormat = format; … … 80 80 } 81 81 82 //-------------------------------------------------------- S DWriter::~SDWriter82 //-------------------------------------------------------- STWriter::~STWriter 83 83 84 84 // Destructor. 85 85 86 S DWriter::~SDWriter()86 STWriter::~STWriter() 87 87 { 88 88 if (cWriter) { … … 91 91 } 92 92 93 //-------------------------------------------------------- S DWriter::setFormat93 //-------------------------------------------------------- STWriter::setFormat 94 94 95 95 // Reset the output format. 96 96 97 Int S DWriter::setFormat(const std::string &format)97 Int STWriter::setFormat(const std::string &format) 98 98 { 99 99 if (format != cFormat) { … … 118 118 } 119 119 120 //------------------------------------------------------------ S DWriter::write120 //------------------------------------------------------------ STWriter::write 121 121 122 122 // Write an SDMemTable to file in the desired format, closing the file when 123 123 // finished. 124 124 125 Int S DWriter::write(const CountedPtr<SDMemTable> in,125 Int STWriter::write(const CountedPtr<SDMemTable> in, 126 126 const std::string &filename, Bool toStokes) 127 127 { … … 180 180 } 181 181 182 Int scanNo = -1; 183 Int cycleNo; 184 Double mjd0 = 0.0; 185 // 182 Vector<Double> srcPM(2, 0.0); 183 Double srcVel = 0.0; 184 186 185 Array<Float> spectra, tSys, stokes; 187 186 Array<uChar> flags; 188 187 Bool doLinear = True; 189 // 188 189 String fieldName, srcName, tcalTime; 190 Vector<Float> calFctr, sigma, tcal, tsys; 191 Vector<Double> direction(2), scanRate(2), srcDir(2), srcPM(2,0.0); 192 Matrix<Float> spectra; 193 Matrix<uChar> flagtra; 194 Complex xCalFctr; 190 195 Int count = 0; 191 for (Int iRow = 0; iRow < in->nRow(); iRow++) { 192 // Extract the next integration from the table. 193 SDContainer sd = in->getSDContainer(iRow); 194 if (sd.scanid != scanNo) { 195 scanNo = sd.scanid; 196 mjd0 = sd.timestamp; 197 cycleNo = 1; 198 } else if (fabs(sd.timestamp-mjd0) > sd.interval) { 199 cycleNo++; 196 Int scanno = 1; 197 // use spearate iterators to ensure renumbering of all numbers 198 TableIterator scanit(table, "SCANNO"); 199 while (!scanit.pastEnd() ) { 200 Table stable = scanit.table(); 201 TableIterator beamit(table, "BEAMNO"); 202 Int beamno = 1; 203 while (!beamit.pastEnd() ) { 204 Table btable = beamit.table(); 205 TableIterator ifit(btable, "IFNO"); 206 Int ifno = 1; 207 while (!ifit.pastEnd() ) { 208 Table itable = ifit.table(); 209 TableIterator cycit(itable, "CYCLENO"); 210 Int cycno = 1; 211 while (!cycit.pastEnd() ) { 212 Table ctable = cycit.table(); 213 TableRow row(ctable); 214 // use the first row to fill in all the "metadata" 215 const TableRecord& rec = row.record(0); 216 ROArrayColumn<Float> specCol(ctable, "SPECTRA"); 217 uInt nchan = specCol(0).nelements(); 218 Double cdelt,crval,crpix, restfreq; 219 String tmp,tmp2; 220 in->frequencies().getEntry(crpix,crval,cdelt, rec.asuInt("FREQ_ID")); 221 in->molecules().getEntry(restfreq,tmp,tmp2,rec.asuInt("RESTFREQ_ID")); 222 Double pixel = Double(nchan/2); 223 Double refFreqNew = (pixel-crpix)*cdelt + crval; 224 // ok, now we have nrows for the n polarizations in this table 225 Matrix<Float> specs; 226 Matrix<uChar> flags; 227 Vector<Complex> xpol; 228 polConversion(specs, flags, xpol, ctable); 229 // dummy data 230 uInt nPol = specs.ncolumns(); 231 Matrix<Float> baseLin(npol,2, 0.0f); 232 Matrix<Float> baseSub(npol,9, 0.0f); 233 Complex xCalFctr; 234 Vector<Double> scanRate(2, 0.0); 235 Vector<Float> sigma(npol, 0.0f); 236 Vector<Float> calFctr(npol, 0.0f); 237 238 239 if (status = cWriter->write(scano, cycno, rec.asDouble("TIME"), 240 rec.asDouble("INTERVAL"), 241 rec.asString("FIELDNAME"), 242 rec.asString("SRCNAME"), 243 direction ,// 244 srcPM, srcVel, // not in scantable yet 245 ifno, 246 refFreqNew, nchan*abs(cdelt), cdelt, 247 restfreq, 248 tcal,// 249 tcaltime,// 250 rec.asFloat("AZIMUTH"), 251 rec.asFloat("ELEVATION"), 252 rec.asFloat("PARANGLE"), 253 focusAxi, focusTan, focusRot,// 254 temperature,// 255 pressure, humidity, windSpeed, windAz,// 256 rec.asInt("REFBEAM"), beamno, 257 direction,// 258 scanRate,// not in scantable 259 tSys, // 260 sigma, calFctr,// not in scantable 261 baseLin, baseSub,// not in scantable 262 specs, flags, 263 xCalFctr,// 264 xpol) 265 ) { 266 cerr << "Error writing output file." << endl; 267 return 1; 268 } 269 270 ++cycno; 271 ++cycit; 272 } 273 ++ifno; 274 ++ifit; 275 } 276 ++beamno; 277 ++beamit; 200 278 } 201 202 // Get FreqID vector 203 freqIDCol.get(iRow, freqIDs); 204 205 // Get Stokes array (all beams/IFs/Stokes). Its not in SDContainer 206 if (toStokes) { 207 stokes = in->getStokesSpectrum(iRow,-1,-1); 208 } 209 210 IPosition start(asap::nAxes,0); 211 IPosition end(stokes.shape()-1); 212 213 // Write it out beam by beam. 214 for (Int iBeam = 0; iBeam < hdr.nbeam; iBeam++) { 215 start(asap::BeamAxis) = iBeam; 216 end(asap::BeamAxis) = iBeam; 217 218 // Write it out IF by IF. 219 for (Int iIF = 0; iIF < hdr.nif; iIF++) { 220 start(asap::IFAxis) = iIF; 221 end(asap::IFAxis) = iIF; 222 uInt freqID = freqIDs(iIF); 223 224 // None of these are stored in SDMemTable by SDReader. 225 //String fieldName = ""; 226 //Vector<Double> srcDir(2, 0.0); 227 Vector<Double> srcPM(2, 0.0); 228 Double srcVel = 0.0; 229 230 // The writer will assume refPix = nChan/2 (0-rel). So recompute 231 // the frequency at this location. 232 233 Double cdelt = sdft.increment(freqID); 234 Double crval = sdft.referenceValue(freqID); 235 Double crpix = sdft.referencePixel(freqID); 236 Double pixel = nChan/2; 237 Double refFreqNew = (pixel-crpix)*cdelt + crval; 238 // 239 //Vector<Float> tcal(2, 0.0f); 240 //String tcalTime = ""; 241 //Float azimuth = 0.0f; 242 //Float elevation = 0.0f; 243 //Float parAngle = 0.0f; 244 Float focusAxi = 0.0f; 245 Float focusTan = 0.0f; 246 Float focusRot = 0.0f; 247 Float temperature = 0.0f; 248 Float pressure = 0.0f; 249 Float humidity = 0.0f; 250 Float windSpeed = 0.0f; 251 Float windAz = 0.0f; 252 //Int refBeam = 0; 253 //Vector<Double> direction(2, 0.0); 254 Vector<Double> scanRate(2, 0.0); 255 Vector<Float> sigma(nPol, 0.0f); 256 Vector<Float> calFctr(nPol, 0.0f); 257 Matrix<Float> baseLin(nPol,2, 0.0f); 258 Matrix<Float> baseSub(nPol,9, 0.0f); 259 Complex xCalFctr; 260 Vector<Complex> xPol; 261 262 // Get data. 263 264 if (toStokes) { 265 266 // Big pain in the ass because 267 // 1) Stokes vector may have less polarizations than input 268 // 2) Stokes column virtual and not in SDContainer 269 // 3) Tsys and FLags already selected for IF and Beam 270 // 3) Have to flip order 271 272 Array<Float> t = sd.getTsys(iBeam,iIF); 273 tSys = SDPolUtil::computeStokesTSysForWriter (t, doLinear); 274 Array<uChar> t2 = sd.getFlags(iBeam,iIF); 275 flags = SDPolUtil::computeStokesFlagsForWriter (t2, doLinear); 276 spectra = SDPolUtil::extractStokesForWriter (stokes,start,end); 277 } else { 278 tSys = sd.getTsys(iBeam,iIF); 279 flags = sd.getFlags(iBeam,iIF); 280 spectra = sd.getSpectrum(iBeam,iIF); 281 } 282 // 283 if (status = cWriter->write(sd.scanid+1, cycleNo, sd.timestamp, 284 sd.interval, sd.fieldname, sd.sourcename, 285 sd.getDirection(iBeam), 286 srcPM, srcVel, iIF+1, refFreqNew, 287 nChan*abs(cdelt), cdelt, restFreq, sd.tcal, 288 sd.tcaltime, sd.azimuth, sd.elevation, 289 sd.parangle, 290 focusAxi, focusTan, focusRot, temperature, 291 pressure, humidity, windSpeed, windAz, 292 sd.refbeam, iBeam+1, 293 sd.getDirection(iBeam), 294 scanRate, 295 tSys, sigma, calFctr, 296 baseLin, baseSub, 297 spectra, flags, 298 xCalFctr, xPol)) { 299 cerr << "Error writing output file." << endl; 300 return 1; 301 } 302 303 count++; 304 } 305 } 279 ++scanno; 280 ++scanit; 306 281 } 307 282 ostringstream oss; 308 oss << "S DWriter: wrote " << count << " rows to " << filename << endl;283 oss << "STWriter: wrote " << count << " rows to " << filename << endl; 309 284 pushLog(String(oss)); 310 285 cWriter->close(); … … 313 288 } 314 289 315 290 void STWriter::polConversion( Matrix< Float >& specs, Matrix< Float >& flags, 291 Vector< Complex > & xpol, const Table & tab ) 292 { 293 TableRow row(tab); 294 String poltype = tab.keywordSet().asString("POLTYPE"); 295 if ( poltype != "linear") 296 String msg = "poltype = " + poltype + " not yet supported in output."; 297 throw(AipsError("msg")); 298 // use the first row to fill in all the "metadata" 299 const TableRecord& rec = row.record(0); 300 ROArrayColumn<Float> specCol(ctable, "SPECTRA"); 301 ROArrayColumn<uChar> flagCol(ctable, "FLAGTRA"); 302 uInt nchan = specCol(0).nelements(); 303 uInt ncols = ( ctable.nrow()==1 ? 1: 2 ); 304 specs.resize(nchan, ncols); 305 flags.resize(nchan, ncols); 306 // the linears 307 for (uInt i=0; i<ncols; ++i) { 308 specs.column(i) = specCol(i); 309 flags.column(i) = flagCol(i); 310 } 311 // now the complex if exists 312 Vector<Complex> xpol; 313 Bool hasxpol = False; 314 xpol.resize(); 315 if ( ctable.nrow() == 4 ) { 316 hasxpol = True; 317 xpol.resize(nchan); 318 Vector<Float> reals, imags; 319 reals = specCol(2); imags = specCol(3); 320 for (uInt k=0; k < nchan; ++k) { 321 ` xpol[k] = Complex(reals[k], imags[k]); 322 } 323 } 324 } -
trunk/src/STWriter.h
r820 r822 1 1 //#--------------------------------------------------------------------------- 2 //# S DWriter.h: ASAP class to write out single dish spectra.2 //# STWriter.h: ASAP class to write out single dish spectra. 3 3 //#--------------------------------------------------------------------------- 4 4 //# Copyright (C) 2004 … … 29 29 //# $Id$ 30 30 //#--------------------------------------------------------------------------- 31 #ifndef S DWRITER_H32 #define S DWRITER_H31 #ifndef STWRITER_H 32 #define STWRITER_H 33 33 34 34 #include <string> … … 46 46 namespace asap { 47 47 48 class S DWriter : public SDLog {48 class STWriter : public SDLog { 49 49 public: 50 S DWriter(const string &format = "SDFITS");51 ~S DWriter();50 STWriter(const string &format = "SDFITS"); 51 ~STWriter(); 52 52 53 53 // Format can be "SDFITS", "FITS", "MS2" or "ASCII" … … 58 58 59 59 private: 60 void polConversion( casa::Matrix<casa::Float>& spec, 61 casa::Matrix<casa::Float>& flag, 62 casa::Vector<casa::Complex>& xpol, 63 const Table& tab); 60 64 std::string cFormat; 61 65 PKSwriter *cWriter;
Note: See TracChangeset
for help on using the changeset viewer.
