source: trunk/external-alma/asdm2ASAP/ASDMFiller.cc@ 2208

Last change on this file since 2208 was 2208, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: Yes CAS-1913

Ready for Test: No

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Update of asdm2ASAP and related classes.

  • replaced LogIO with StreamLogSink
  • all log messages are written to logger via StreamLogSink
  • no output to stdout/stderr
  • commented out junk log
  • added time_sampling selection
  • supported wvr_corrected_data='both'
  • added logfile specification
  • added corr_mode selection


File size: 17.2 KB
Line 
1#include<iostream>
2
3#include <STHeader.h>
4
5#include <measures/Measures/MEpoch.h>
6#include <measures/Measures/MPosition.h>
7#include <measures/Measures/MDirection.h>
8#include <measures/Measures/MCDirection.h>
9#include <measures/Measures/MeasFrame.h>
10#include <measures/Measures/MeasConvert.h>
11#include <casa/Arrays/Vector.h>
12#include <casa/Arrays/Matrix.h>
13#include <casa/Quanta/MVTime.h>
14#include <casa/Logging/LogMessage.h>
15
16#include "ASDMFiller.h"
17
18using namespace std ;
19using namespace casa ;
20using namespace asap ;
21
22ASDMFiller::ASDMFiller( CountedPtr<Scantable> stable )
23 : FillerBase( stable ),
24 antennaId_( -1 ),
25 antennaName_( "" ),
26 className_("ASDMFiller")
27{
28 reader_ = new ASDMReader() ;
29}
30
31ASDMFiller::~ASDMFiller()
32{
33 // nothing to do?
34 logsink_ = 0 ;
35}
36
37void ASDMFiller::setLogger( CountedPtr<LogSinkInterface> &logsink )
38{
39 logsink_ = logsink ;
40 if ( !(reader_.null()) ) {
41 reader_->setLogger( logsink ) ;
42 }
43}
44
45bool ASDMFiller::open( const string &filename, const Record &rec )
46{
47 String funcName = "open" ;
48 bool status = reader_->open( filename, rec ) ;
49
50 antennaId_ = reader_->getAntennaId() ;
51 antennaName_ = reader_->getAntennaName() ;
52
53 //logsink->postLocally( LogMessage("antennaId_ = "+String::toString(antennaId_),LogOrigin(className_,funcName,WHERE)) ) ;
54 //logsink->postLocally( LogMessage("antennaName_ = "+antennaName_,LogOrigin(className_,funcName,WHERE)) ) ;
55
56 return status ;
57}
58
59void ASDMFiller::fill()
60{
61 String funcName = "fill" ;
62
63 // header
64 fillHeader() ;
65
66 Vector<casa::Double> antpos = table_->getHeader().antennaposition ;
67
68 //STHeader hdr = table_->getHeader() ;
69
70 // data selection
71 reader_->select() ;
72
73 // pick up valid configDescriptionId
74 Vector<uInt> configDescIdList = reader_->getConfigDescriptionIdList() ;
75 uInt numConfigDescId = configDescIdList.size() ;
76
77 //logsink->postLocally( LogMessage("configDescIdList = "+String::toString(configDescIdList),LogOrigin(className_,funcName,WHERE)) ) ;
78
79 // get field list
80 Vector<uInt> fieldIdList = reader_->getFieldIdList() ;
81 uInt numFieldId = fieldIdList.size() ;
82
83 //logsink->postLocally( LogMessage("fieldIdList = "+String::toString(fieldIdList),LogOrigin(className_,funcName,WHERE)) ) ;
84
85 // BEAMNO is always 0 since ALMA antenna is single beam
86 uInt beamno = 0 ;
87
88 // REFBEAMNO is -1
89 setReferenceBeam() ;
90
91 // fill FOCUS_ID and add FOCUS row if necessary
92 setFocus() ;
93
94 for ( uInt icon = 0 ; icon < numConfigDescId ; icon++ ) {
95 for ( unsigned int ifield = 0 ; ifield < numFieldId ; ifield++ ) {
96 //logsink_->postLocally( LogMessage("start configDescId "+String::toString(configDescIdList[icon])+" fieldId "+String::toString(fieldIdList[ifield]),LogOrigin(className_,funcName,WHERE)) ) ;
97
98 //Bool status = reader_->setMainRow( configDescIdList[icon], fieldIdList[ifield] ) ;
99 if ( !(reader_->setMainRow( configDescIdList[icon], fieldIdList[ifield] )) ) {
100 //logsink_->postLocally( LogMessage("skip configDescId "+String::toString(configDescIdList[icon])+" fieldId "+String::toString(fieldIdList[ifield]),LogOrigin(className_,funcName,WHERE)) ) ;
101 continue ;
102 }
103
104 // number of rows
105 uInt nrow = reader_->getNumMainRow() ;
106
107 //logsink_->postLocally( LogMessage("There are "+String::toString(nrow)+" rows in Main table corresponding to configDescId "+String::toString(configDescIdList[icon]+" fieldId "+String::toString(fieldIdList[ifield]),LogOrigin(className_,funcName,WHERE)) ) ;
108
109 // CYCLENO
110 unsigned int cycleno = 0 ;
111
112 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
113
114 // set main row
115 if ( !(reader_->setMainRow( irow )) ) {
116 // skip row since the row doesn't have valid configDescId
117 //logsink_->postLocally( LogMessage("skip "+String::toString(irow),LogOrigin(className_,funcName,WHERE)) ) ;
118 continue ;
119 }
120
121 // scan and subscan
122 unsigned int scanno = reader_->getScanNo() ;
123 //uInt subscanno = reader_->getSubscanNo() ;
124
125 // set data
126 if ( !(reader_->setData()) ) {
127 // skip row since reader failed to retrieve data
128 //logsink_->postLocally( LogMessage("skip "+String::toString(irow),LogOrigin(className_,funcName,WHRER)) ) ;
129 continue ;
130 }
131
132 unsigned int numData = reader_->getNumData() ;
133 double refpix = 0.0 ;
134 double refval = 0.0 ;
135 double incr = 0.0 ;
136
137 for ( unsigned int idata = 0 ; idata < numData ; idata++ ) {
138
139 // subscan number
140 unsigned int subscanno = reader_->getSubscanNo( idata ) ;
141
142 // IFNO
143 uInt ifno = reader_->getIFNo( idata ) ;
144
145
146 // REFPIX, REFVAL, INCREMENT
147 String ifkey = getIFKey( ifno ) ;
148 if ( ifrec_.isDefined( ifkey ) ) {
149 getFrequencyRec( ifkey, refpix, refval, incr ) ;
150 }
151 else {
152 reader_->getFrequency( idata, refpix, refval, incr ) ;
153 setFrequencyRec( ifkey, refpix, refval, incr ) ;
154 }
155
156 // fill FREQ_ID and add FREQUENCIES row if necessary
157 setFrequency( (casa::Double)refpix, (casa::Double)refval, (casa::Double)incr ) ;
158
159
160 // rest frequency
161 vector<double> rf = reader_->getRestFrequency( idata ) ;
162
163 // fill MOLECULE_ID and add MOLECULES row if necessary
164 Vector<casa::Double> restFreqs( rf.size() ) ;
165 for ( uInt i = 0 ; i < rf.size() ; i++ )
166 restFreqs[i] = (casa::Double)(rf[i]) ;
167 setMolecule( restFreqs ) ;
168
169
170 // time and interval
171 casa::Double mjd = (casa::Double)(reader_->getTime( idata )) ;
172 casa::Double interval = (casa::Double)(reader_->getInterval( idata )) ;
173
174 // fill TIME and INTERVAL
175 setTime( mjd, interval ) ;
176
177
178 // source spec
179 string srcname = reader_->getSourceName( idata ) ;
180 string fieldname = reader_->getFieldName( idata ) ;
181 int srctype = reader_->getSrcType( scanno, subscanno ) ;
182 vector<double> srcDirection = reader_->getSourceDirection( idata ) ;
183 vector<double> srcProperMotion = reader_->getSourceProperMotion( idata ) ;
184 double sysVel = reader_->getSysVel( idata ) ;
185
186 // fill SRCNAME, SRCTYPE, FIELDNAME, SRCDIRECTION, SRCPROPERMOTION, and SRCVELOCITY
187 Vector<casa::Double> srcDir( 2 ) ;
188 srcDir[0] = (casa::Double)(srcDirection[0]) ;
189 srcDir[1] = (casa::Double)(srcDirection[1]) ;
190 Vector<casa::Double> srcPM( 2 ) ;
191 srcPM[0] = (casa::Double)(srcProperMotion[0]) ;
192 srcPM[1] = (casa::Double)(srcProperMotion[1]) ;
193 setSource( srcname, srctype, fieldname, srcDir, srcPM, (casa::Double)sysVel ) ;
194
195 // fill FLAGROW
196 unsigned int flagrow = reader_->getFlagRow( idata ) ;
197 setFlagrow( (uInt)flagrow ) ;
198
199 // fill WEATHER_ID and add WEATHER row if necessary
200 float temperature ;
201 float pressure ;
202 float humidity ;
203 float windspeed ;
204 float windaz ;
205 reader_->getWeatherInfo( idata,
206 temperature,
207 pressure,
208 humidity,
209 windspeed,
210 windaz ) ;
211 setWeather2( (casa::Float)temperature,
212 (casa::Float)pressure,
213 (casa::Float)humidity,
214 (casa::Float)windspeed,
215 (casa::Float)windaz ) ;
216
217 // fill AZIMUTH, ELEVATION, DIRECTION and SCANRATE
218 vector<double> dir ;
219 double az ;
220 double el ;
221 vector<double> srate ;
222 reader_->getPointingInfo( idata,
223 dir,
224 az,
225 el,
226 srate ) ;
227 Vector<casa::Double> scanRate( 2, 0.0 ) ;
228 Vector<casa::Double> direction( 2, 0.0 ) ;
229 if ( srate.size() > 0 ) {
230 scanRate[0] = (casa::Double)(srate[0]) ;
231 scanRate[1] = (casa::Double)(srate[1]) ;
232 }
233 setScanRate( scanRate ) ;
234 if ( dir.size() > 0 ) {
235 direction[0] = (casa::Double)(dir[0]) ;
236 direction[1] = (casa::Double)(dir[1]) ;
237 }
238 else {
239 toJ2000( direction, az, el, mjd, antpos ) ;
240 }
241 //logsink_->postLocally( LogMessage("direction = "+String::toString(direction),LogOrigin(className_,funcName,WHERE)) ) ;
242 setDirection( direction, (casa::Float)az, (casa::Float)el ) ;
243
244 // loop on polarization
245 vector<unsigned int> dataShape = reader_->getDataShape( idata ) ;
246// ostringstream oss ;
247// for ( unsigned int i = 0 ; i < dataShape.size() ; i++ ) {
248// if ( i == 0 )
249// oss << "dataShape=[" << dataShape[i] << ", " ;
250// else if ( i == dataShape.size()-1 )
251// oss << dataShape[i] << "]" ;
252// else
253// oss << dataShape[i] << ", " ;
254// }
255// logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
256
257 //int numPol = reader_->getNumPol( idata ) ;
258 unsigned int numPol = dataShape[0] ;
259 unsigned int numChan = dataShape[1] ;
260
261 //logsink_->postLocally( LogMessage("numPol = "+String::toString(numPol),LogOrigin(className_,funcName,WHERE)) ) ;
262
263 // OPACITY
264 vector<float> tau = reader_->getOpacity( idata ) ;
265 Vector<casa::Float> opacity = toVector( tau, numPol ) ;
266
267 // SPECTRA, FLAGTRA, TSYS, TCAL
268 float *sp = reader_->getSpectrum( idata ) ;
269 vector< vector<float> > ts = reader_->getTsys( idata ) ;
270 vector< vector<float> > tc = reader_->getTcal( idata ) ;
271 Matrix<casa::Float> spectra = toMatrix( sp, numPol, numChan ) ;
272 Vector<uChar> flagtra( numChan, 0 ) ;
273 Matrix<casa::Float> tsys = toMatrix( ts, numPol, numChan ) ;
274 Matrix<casa::Float> tcal = toMatrix( tc, numPol, numChan ) ;
275// String caltime = "" ;
276// if ( anyNE( tcal, (casa::Float)1.0 ) )
277// caltime = toTcalTime( mjd ) ;
278 String caltime = toTcalTime( mjd ) ;
279
280 for ( unsigned int ipol = 0 ; ipol < numPol ; ipol++ ) {
281
282 // fill SCANNO, CYCLENO, IFNO, POLNO, and BEAMNO
283 setIndex( (uInt)scanno-1, (uInt)cycleno, ifno, ipol, beamno ) ;
284
285 // fill SPECTRA, FLAGTRA, TSYS
286 setSpectrum( spectra.row(ipol), flagtra, tsys.row(ipol) ) ;
287
288 // fill TCAL_ID and add TCAL row if necessary
289 setTcal2( caltime, tcal.row(ipol) ) ;
290
291 // fill OPACITY
292 setOpacity( opacity[ipol] ) ;
293
294 // commit row
295 commitRow() ;
296 }
297
298 // increment CYCLENO
299 cycleno++ ;
300 }
301 }
302 }
303 }
304
305 return ;
306}
307
308void ASDMFiller::close()
309{
310 reader_->close() ;
311 reader_ = 0 ;
312
313 return ;
314}
315
316void ASDMFiller::fillHeader()
317{
318 STHeader hdr ;
319
320 reader_->fillHeader( hdr.nchan,
321 hdr.npol,
322 hdr.nif,
323 hdr.nbeam,
324 hdr.observer,
325 hdr.project,
326 hdr.obstype,
327 hdr.antennaname,
328 hdr.antennaposition,
329 hdr.equinox,
330 hdr.freqref,
331 hdr.reffreq,
332 hdr.bandwidth,
333 hdr.utc,
334 hdr.fluxunit,
335 hdr.epoch,
336 hdr.poltype ) ;
337
338 setHeader( hdr ) ;
339}
340
341String ASDMFiller::getIFKey( uInt ifno )
342{
343 return "IFNO"+String::toString( ifno ) ;
344}
345
346void ASDMFiller::getFrequencyRec( String key,
347 double &refpix,
348 double &refval,
349 double &incr )
350{
351 Record frec = ifrec_.asRecord( key ) ;
352 refpix = frec.asdouble( "REFPIX" ) ;
353 refval = frec.asdouble( "REFVAL" ) ;
354 incr = frec.asdouble( "INCREMENT" ) ;
355}
356
357void ASDMFiller::setFrequencyRec( String key,
358 double refpix,
359 double refval,
360 double incr )
361{
362 Record frec ;
363 frec.define( "REFPIX", refpix ) ;
364 frec.define( "REFVAL", refval ) ;
365 frec.define( "INCREMENT", incr ) ;
366 ifrec_.defineRecord( key, frec ) ;
367}
368
369Matrix<casa::Float> ASDMFiller::toMatrix( float *sp,
370 unsigned int npol,
371 unsigned int nchan )
372{
373 Matrix<casa::Float> mSp( npol, nchan ) ;
374 if ( npol <= 2 ) {
375 // 1 or 2 polarization case
376 for ( unsigned int ich = 0 ; ich < nchan ; ich++ ) {
377 for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ ) {
378 mSp(ipol,ich) = (casa::Float)(sp[npol*ich+ipol]) ;
379 }
380 }
381 }
382 else {
383 // 4 polarization case
384 for ( unsigned int ich = 0 ; ich < nchan ; ich++ ) {
385 mSp(0,ich) = (casa::Float)(sp[4*ich]) ; // Re(XX)
386 mSp(1,ich) = (casa::Float)(sp[4*ich+4]) ; // Re(YY)
387 mSp(2,ich) = (casa::Float)(sp[4*ich+2]) ; // Re(XY)
388 mSp(3,ich) = (casa::Float)(sp[4*ich+3]) ; // Im(XY)
389 }
390 }
391 return mSp ;
392}
393
394Matrix<casa::Float> ASDMFiller::toMatrix( vector< vector<float> > &tsys,
395 unsigned int npol,
396 unsigned int nchan )
397{
398 unsigned int numRec = tsys.size() ;
399 unsigned int numChan = tsys[0].size() ;
400 Matrix<casa::Float> ret ;
401 if ( npol == numRec && nchan == numChan ) {
402 ret.resize( npol, nchan ) ;
403 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
404 for ( unsigned int ic = 0 ; ic < nchan ; ic++ )
405 ret( ip, ic ) = (casa::Float)(tsys[ip][ic]) ;
406 }
407 else if ( npol == numRec && numChan == 1 ) {
408 ret.resize( npol, 1 ) ;
409 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
410 ret( ip, 0 ) = (casa::Float)(tsys[0][0]) ;
411 }
412 else if ( numRec == 1 && nchan == numChan ) {
413 ret.resize( npol, nchan ) ;
414 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
415 for ( unsigned int ic = 0 ; ic < nchan ; ic++ )
416 ret( ip, ic ) = (casa::Float)(tsys[0][ic]) ;
417 }
418 else if ( numRec == 1 && numChan == 1 ) {
419 ret.resize( npol, 1 ) ;
420 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
421 ret( ip, 0 ) = (casa::Float)(tsys[0][0]) ;
422 }
423 else if ( numRec == 2 && npol == 4 && numChan == nchan ) {
424 // TODO: How to determine Tsys for XY?
425 // at the moment Tsys[XY] = 0.5*(Tsys[X]+Tsys[Y])
426 ret.resize( npol, nchan ) ;
427 for ( unsigned int ic = 0 ; ic < nchan ; ic++ ) {
428 casa::Float tsysxy = (casa::Float)(0.5*(tsys[0][ic]+tsys[1][ic])) ;
429 ret( 0, ic ) = (casa::Float)(tsys[0][ic]) ;
430 ret( 1, ic ) = (casa::Float)(tsys[1][ic]) ;
431 ret( 2, ic ) = tsysxy ;
432 ret( 3, ic ) = tsysxy ;
433 }
434 }
435 else if ( numRec == 2 && npol == 4 && numChan == 1 ) {
436 // TODO: How to determine Tsys for XY?
437 // at the moment Tsys[XY] = 0.5*(Tsys[X]+Tsys[Y])
438 ret.resize( npol, 1 ) ;
439 casa::Float tsysxy = (casa::Float)(0.5*(tsys[0][0]+tsys[1][0])) ;
440 ret( 0, 0 ) = (casa::Float)(tsys[0][0]) ;
441 ret( 1, 0 ) = (casa::Float)(tsys[1][0]) ;
442 ret( 2, 0 ) = tsysxy ;
443 ret( 3, 0 ) = tsysxy ;
444 }
445 else {
446 // I don't know how to handle ...
447 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
448 for ( unsigned int ic = 0 ; ic < nchan ; ic++ )
449 ret( ip, ic ) = (casa::Float)(tsys[0][ic]) ;
450 }
451 return ret ;
452}
453
454Vector<casa::Float> ASDMFiller::toVector( vector<float> &tau,
455 unsigned int npol )
456{
457 Vector<casa::Float> ret( npol ) ;
458
459 if ( npol == 4 ) {
460 ret[0] = (casa::Float)tau[0] ;
461 ret[1] = (casa::Float)tau[1] ;
462 ret[2] = 0.5 * ( ret[0] + ret[1] ) ;
463 ret[3] = ret[2] ;
464 }
465 else if ( npol == tau.size() ) {
466 for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ )
467 ret[ipol] = (casa::Float)tau[ipol] ;
468 }
469 else {
470 // I don't know how to handle...
471 for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ )
472 ret[ipol] = (casa::Float)tau[0] ;
473 }
474 return ret ;
475}
476
477String ASDMFiller::toTcalTime( casa::Double mjd )
478{
479 return MVTime( mjd ).string( MVTime::YMD ) ;
480}
481
482void ASDMFiller::toJ2000( Vector<casa::Double> &dir,
483 double az,
484 double el,
485 casa::Double mjd,
486 Vector<casa::Double> antpos )
487{
488 String funcName = "toJ2000" ;
489
490 Vector<casa::Double> azel( 2 ) ;
491 azel[0] = az ;
492 azel[1] = el ;
493 MEpoch me( Quantity( mjd, "d" ), MEpoch::UTC ) ;
494 Vector<Quantity> qantpos( 3 ) ;
495 qantpos[0] = Quantity( antpos[0], "m" ) ;
496 qantpos[1] = Quantity( antpos[1], "m" ) ;
497 qantpos[2] = Quantity( antpos[2], "m" ) ;
498 MPosition mp( MVPosition( qantpos ),
499 MPosition::ITRF ) ;
500// mp.print( os_.output() ) ;
501 MeasFrame mf( me, mp ) ;
502 MDirection::Convert toj2000( MDirection::AZELGEO,
503 MDirection::Ref( MDirection::J2000, mf ) ) ;
504 dir = toj2000( azel ).getAngle( "rad" ).getValue() ;
505 //logsink->postLocally( LogMessage("dir = "+String::toString(dir),LogOrigin(className_,funcName,WHERE)) ) ;
506}
507
Note: See TracBrowser for help on using the repository browser.