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

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

New Development: Yes

JIRA Issue: Yes CAS-1913

Ready for Test: No

Interface Changes: Yes/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...

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