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.