source: branches/plotter2/external-alma/asdm2ASAP/ASDMFiller.cc@ 2900

Last change on this file since 2900 was 2754, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: Yes CSV-2532 (may be related to CSV-1908 and CSV-2161)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: In tool level, added parameter 'freq_tolsr' to

scantable constructor and function sd.splitant.

Test Programs: test_sdsave, test_importasdm_sd

Put in Release Notes: Yes

Module(s): Module Names change impacts.

Description: Describe your changes here...

In importing MS to Scantable, frequency frame information is
imported as is by default, i.e., base frame in Scantable is
TOPO for ALMA data, which is forcibly converted to LSRK with
wrong time and direction reference.

Some functions have a boolean parameter 'freq_tolsr' that controls
the above behavior. If freq_tolsr is False (default), frequency
is imported as is, while frequency is converted to LSRK (wrongly)
when it is True.


File size: 24.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/MFrequency.h>
10#include <measures/Measures/MCFrequency.h>
11#include <measures/Measures/MeasFrame.h>
12#include <measures/Measures/MeasConvert.h>
13#include <casa/Arrays/Vector.h>
14#include <casa/Arrays/Matrix.h>
15#include <casa/Quanta/MVTime.h>
16#include <casa/Logging/LogMessage.h>
17
18#include "ASDMFiller.h"
19
20using namespace std ;
21using namespace casa ;
22using namespace asap ;
23
24ASDMFiller::ASDMFiller( CountedPtr<Scantable> stable )
25 : FillerBase( stable ),
26 antennaId_( -1 ),
27 antennaName_( "" ),
28 className_("ASDMFiller"),
29 freqToLsr_(False)
30{
31 reader_ = new ASDMReader() ;
32}
33
34ASDMFiller::~ASDMFiller()
35{
36 // nothing to do?
37 logsink_ = 0 ;
38}
39
40void ASDMFiller::setLogger( CountedPtr<LogSinkInterface> &logsink )
41{
42 logsink_ = logsink ;
43 if ( !(reader_.null()) ) {
44 reader_->setLogger( logsink ) ;
45 }
46}
47
48bool ASDMFiller::open( const string &filename, const Record &rec )
49{
50 String funcName = "open" ;
51 bool status = reader_->open( filename, rec ) ;
52
53 antennaId_ = reader_->getAntennaId() ;
54 antennaName_ = reader_->getAntennaName() ;
55
56 if (rec.isDefined("asdm")) {
57 Record asdmrec = rec.asRecord("asdm") ;
58 if (asdmrec.isDefined("freq_tolsr")) {
59 freqToLsr_ = asdmrec.asBool("freq_tolsr");
60 }
61 }
62 logsink_->postLocally(LogMessage("freqToLsr_ = "+String(freqToLsr_?"True":"False"), LogOrigin(className_, funcName, WHERE)));
63
64 //logsink_->postLocally( LogMessage("antennaId_ = "+String::toString(antennaId_),LogOrigin(className_,funcName,WHERE)) ) ;
65 //logsink_->postLocally( LogMessage("antennaName_ = "+antennaName_,LogOrigin(className_,funcName,WHERE)) ) ;
66
67 return status ;
68}
69
70void ASDMFiller::fill()
71{
72 String funcName = "fill" ;
73
74 // header
75 fillHeader() ;
76 const STHeader hdr = table_->getHeader();
77
78 // set Frame for FREQUENCIES table
79 string sFreqFrame = reader_->getFrame() ;
80 //MFrequency::Types freqFrame = toFrameType( sFreqFrame ) ;
81 MFrequency::Types freqFrame = MFrequency::LSRK ;
82 table_->frequencies().setFrame( freqFrame, false ) ;
83 if ( freqToLsr_ ) {
84 table_->frequencies().setFrame( freqFrame, true ) ;
85 }
86 else {
87 string baseFrame = reader_->getFrame() ;
88 table_->frequencies().setFrame( baseFrame, true ) ;
89 }
90 //logsink_->postLocally( LogMessage("sFreqFrame = "+sFreqFrame,LogOrigin(className_,funcName,WHERE)) ) ;
91
92 Vector<casa::Double> antpos = table_->getHeader().antennaposition ;
93
94 // data selection
95 reader_->select() ;
96
97 // pick up valid configDescriptionId
98 Vector<uInt> configDescIdList = reader_->getConfigDescriptionIdList() ;
99 uInt numConfigDescId = configDescIdList.size() ;
100
101 //logsink_->postLocally( LogMessage("configDescIdList = "+String::toString(configDescIdList),LogOrigin(className_,funcName,WHERE)) ) ;
102
103 // get field list
104 Vector<uInt> fieldIdList = reader_->getFieldIdList() ;
105 uInt numFieldId = fieldIdList.size() ;
106
107 //logsink_->postLocally( LogMessage("fieldIdList = "+String::toString(fieldIdList),LogOrigin(className_,funcName,WHERE)) ) ;
108
109 // BEAMNO is always 0 since ALMA antenna is single beam
110 uInt beamno = 0 ;
111
112 // REFBEAMNO is -1
113 setReferenceBeam() ;
114
115 // fill FOCUS_ID and add FOCUS row if necessary
116 setFocus() ;
117
118 // CYCLENO
119 map< unsigned int, unsigned int > cycleno ;
120 map< unsigned int, unsigned int >::iterator citer ;
121
122 for ( unsigned int ifield = 0 ; ifield < numFieldId ; ifield++ ) {
123 for ( uInt icon = 0 ; icon < numConfigDescId ; icon++ ) {
124 //logsink_->postLocally( LogMessage("start configDescId "+String::toString(configDescIdList[icon])+" fieldId "+String::toString(fieldIdList[ifield]),LogOrigin(className_,funcName,WHERE)) ) ;
125
126 if ( !(reader_->setMainRow( configDescIdList[icon], fieldIdList[ifield] )) ) {
127 //logsink_->postLocally( LogMessage("skip configDescId "+String::toString(configDescIdList[icon])+" fieldId "+String::toString(fieldIdList[ifield]),LogOrigin(className_,funcName,WHERE)) ) ;
128 continue ;
129 }
130
131 // number of rows
132 uInt nrow = reader_->getNumMainRow() ;
133
134 //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)) ) ;
135
136 for ( uInt irow = 0 ; irow < nrow ; irow++ ) {
137
138 // set main row
139 if ( !(reader_->setMainRow( irow )) ) {
140 // skip row since the row doesn't have valid configDescId
141 //logsink_->postLocally( LogMessage("skip "+String::toString(irow)+" since ASDMReader::setMainrow() returns False",LogOrigin(className_,funcName,WHERE)) ) ;
142 continue ;
143 }
144
145 // scan and subscan
146 unsigned int scanno = reader_->getScanNoOfCurrentRow() ;
147 //logsink_->postLocally( LogMessage("scanno = "+String::toString(scanno),LogOrigin(className_,funcName,WHERE)) ) ;
148 citer = cycleno.find( scanno ) ;
149 if ( citer == cycleno.end() )
150 cycleno[scanno] = 0 ;
151
152 // set data
153 if ( !(reader_->setData()) ) {
154 // skip row since reader failed to retrieve data
155 //logsink_->postLocally( LogMessage("skip "+String::toString(irow)+" since ASDMReader::setData() returns False",LogOrigin(className_,funcName,WHERE)) ) ;
156 continue ;
157 }
158
159 unsigned int numData = reader_->getNumData() ;
160 double refpix = 0.0 ;
161 double refval = 0.0 ;
162 double incr = 0.0 ;
163 string freqref = "" ;
164
165 //logsink_->postLocally( LogMessage("numData = "+String::toString(numData),LogOrigin(className_,funcName,WHERE)) ) ;
166 for ( unsigned int idata = 0 ; idata < numData ; idata++ ) {
167 // prepare to extract binary data
168 //logsink_->postLocally( LogMessage("prepare data...",LogOrigin(className_,funcName,WHERE)) ) ;
169 reader_->prepareData( idata ) ;
170
171 // subscan number
172 unsigned int subscanno = reader_->getSubscanNo() ;
173 //logsink_->postLocally( LogMessage("subscanno = "+String::toString(subscanno),LogOrigin(className_,funcName,WHERE)) ) ;
174
175 // IFNO
176 uInt ifno = reader_->getIFNo() ;
177 //logsink_->postLocally( LogMessage("ifno = "+String::toString(ifno),LogOrigin(className_,funcName,WHERE)) ) ;
178 // source spec
179 int srctype = reader_->getSrcType( scanno, subscanno ) ;
180 //logsink_->postLocally( LogMessage("srctype = "+String::toString(srctype),LogOrigin(className_,funcName,WHERE)) ) ;
181 string srcname ;
182 string fieldname ;
183 vector<double> srcDirection ;
184 vector<double> srcProperMotion ;
185 double sysVel ;
186 vector<double> rf ;
187 reader_->getSourceProperty( srcname,
188 fieldname,
189 srcDirection,
190 srcProperMotion,
191 sysVel,
192 rf ) ;
193 //logsink_->postLocally( LogMessage("srcname = "+String::toString(srcname),LogOrigin(className_,funcName,WHERE)) ) ;
194
195 // fill MOLECULE_ID and add MOLECULES row if necessary
196 Vector<casa::Double> restFreqs( rf.size() ) ;
197 for ( uInt i = 0 ; i < rf.size() ; i++ )
198 restFreqs[i] = (casa::Double)(rf[i]) ;
199 setMolecule( restFreqs ) ;
200
201 // time and interval
202 casa::Double mjd = (casa::Double)(reader_->getTime()) ;
203 casa::Double interval = (casa::Double)(reader_->getInterval()) ;
204 //logsink_->postLocally( LogMessage("mjd = "+String::toString(mjd),LogOrigin(className_,funcName,WHERE)) ) ;
205
206 // fill TIME and INTERVAL
207 setTime( mjd, interval ) ;
208
209 // fill SRCNAME, SRCTYPE, FIELDNAME, SRCDIRECTION, SRCPROPERMOTION, and SRCVELOCITY
210 Vector<casa::Double> srcDir( 2 ) ;
211 srcDir[0] = (casa::Double)(srcDirection[0]) ;
212 srcDir[1] = (casa::Double)(srcDirection[1]) ;
213 Vector<casa::Double> srcPM( 2 ) ;
214 srcPM[0] = (casa::Double)(srcProperMotion[0]) ;
215 srcPM[1] = (casa::Double)(srcProperMotion[1]) ;
216 setSource( srcname, srctype, fieldname, srcDir, srcPM, (casa::Double)sysVel ) ;
217
218 // fill FLAGROW
219 unsigned int flagrow = reader_->getFlagRow() ;
220 setFlagrow( (uInt)flagrow ) ;
221 //logsink_->postLocally( LogMessage("flagrow = "+String::toString(flagrow),LogOrigin(className_,funcName,WHERE)) ) ;
222 // fill WEATHER_ID and add WEATHER row if necessary
223 float temperature ;
224 float pressure ;
225 float humidity ;
226 float windspeed ;
227 float windaz ;
228 reader_->getWeatherInfo( temperature,
229 pressure,
230 humidity,
231 windspeed,
232 windaz ) ;
233 setWeather2( (casa::Float)temperature,
234 (casa::Float)pressure,
235 (casa::Float)humidity,
236 (casa::Float)windspeed,
237 (casa::Float)windaz ) ;
238 //logsink_->postLocally( LogMessage("temperature = "+String::toString(temperature),LogOrigin(className_,funcName,WHERE)) ) ;
239 // fill AZIMUTH, ELEVATION, DIRECTION and SCANRATE
240 vector<double> dir ;
241 double az ;
242 double el ;
243 vector<double> srate ;
244 reader_->getPointingInfo( dir,
245 az,
246 el,
247 srate ) ;
248 Vector<casa::Double> scanRate( 2, 0.0 ) ;
249 Vector<casa::Double> direction( 2, 0.0 ) ;
250 if ( srate.size() > 0 ) {
251 scanRate[0] = (casa::Double)(srate[0]) ;
252 scanRate[1] = (casa::Double)(srate[1]) ;
253 }
254 setScanRate( scanRate ) ;
255 if ( dir.size() > 0 ) {
256 direction[0] = (casa::Double)(dir[0]) ;
257 direction[1] = (casa::Double)(dir[1]) ;
258 }
259 else {
260 toJ2000( direction, az, el, mjd, antpos ) ;
261 }
262 //logsink_->postLocally( LogMessage("direction = "+String::toString(direction),LogOrigin(className_,funcName,WHERE)) ) ;
263 setDirection( direction, (casa::Float)az, (casa::Float)el ) ;
264
265 // REFPIX, REFVAL, INCREMENT
266 String ifkey = getIFKey( ifno ) ;
267 if ( ifrec_.isDefined( ifkey ) ) {
268 getFrequencyRec( ifkey, refpix, refval, incr ) ;
269 }
270 else {
271 reader_->getFrequency( refpix, refval, incr, freqref ) ;
272 refval = (double)toLSRK( casa::Double(refval),
273 String(freqref),
274 hdr.utc,
275 antpos,
276 //direction,
277 srcDir,
278 "J2000" ) ;
279 setFrequencyRec( ifkey, refpix, refval, incr ) ;
280 }
281
282 // fill FREQ_ID and add FREQUENCIES row if necessary
283 setFrequency( (casa::Double)refpix, (casa::Double)refval, (casa::Double)incr ) ;
284
285 // loop on polarization
286 vector<unsigned int> dataShape = reader_->getDataShape() ;
287// ostringstream oss ;
288// for ( unsigned int i = 0 ; i < dataShape.size() ; i++ ) {
289// if ( i == 0 )
290// oss << "dataShape=[" << dataShape[i] << ", " ;
291// else if ( i == dataShape.size()-1 )
292// oss << dataShape[i] << "]" ;
293// else
294// oss << dataShape[i] << ", " ;
295// }
296// logsink_->postLocally( LogMessage(oss.str(),LogOrigin(className_,funcName,WHERE)) ) ;
297
298 unsigned int numPol = dataShape[0] ;
299 unsigned int numChan = dataShape[1] ;
300
301 //logsink_->postLocally( LogMessage("numPol = "+String::toString(numPol),LogOrigin(className_,funcName,WHERE)) ) ;
302
303 // OPACITY
304 vector<float> tau = reader_->getOpacity() ;
305 Vector<casa::Float> opacity = toVector( tau, numPol ) ;
306
307 // SPECTRA, FLAGTRA, TSYS, TCAL
308 float *sp = reader_->getSpectrum() ;
309 //logsink_->postLocally( LogMessage("sp[0] = "+String::toString(sp[0]),LogOrigin(className_,funcName,WHERE)) ) ;
310 vector< vector<float> > ts ;
311 vector< vector<float> > tc ;
312 reader_->getTcalAndTsys( tc, ts ) ;
313 Matrix<casa::Float> spectra = toMatrix( sp, numPol, numChan ) ;
314 //logsink_->postLocally( LogMessage("spectra(0,0) = "+String::toString(spectra(0,0)),LogOrigin(className_,funcName,WHERE)) ) ;
315 Vector<uChar> flagtra( numChan, 0 ) ;
316 Matrix<casa::Float> tsys = toMatrix( ts, numPol, numChan ) ;
317 Matrix<casa::Float> tcal = toMatrix( tc, numPol, numChan ) ;
318 //logsink_->postLocally( LogMessage("tsys(0,0) = "+String::toString(tsys(0,0)),LogOrigin(className_,funcName,WHERE)) ) ;
319// String caltime = "" ;
320// if ( anyNE( tcal, (casa::Float)1.0 ) )
321// caltime = toTcalTime( mjd ) ;
322 String caltime = toTcalTime( mjd ) ;
323
324 for ( unsigned int ipol = 0 ; ipol < numPol ; ipol++ ) {
325
326 // fill SCANNO, CYCLENO, IFNO, POLNO, and BEAMNO
327 setIndex( (uInt)scanno-1, (uInt)cycleno[scanno], ifno, ipol, beamno ) ;
328
329 // fill SPECTRA, FLAGTRA, TSYS
330 setSpectrum( spectra.row(ipol), flagtra, tsys.row(ipol) ) ;
331
332 // fill TCAL_ID and add TCAL row if necessary
333 setTcal2( caltime, tcal.row(ipol) ) ;
334
335 // fill OPACITY
336 setOpacity( opacity[ipol] ) ;
337
338 // commit row
339 commitRow() ;
340 }
341
342 // increment CYCLENO
343 cycleno[scanno]++ ;
344 }
345 }
346 }
347 }
348
349 return ;
350}
351
352void ASDMFiller::close()
353{
354 reader_->close() ;
355 reader_ = 0 ;
356
357 return ;
358}
359
360void ASDMFiller::fillHeader()
361{
362 STHeader hdr ;
363
364 reader_->fillHeader( hdr.nchan,
365 hdr.npol,
366 hdr.nif,
367 hdr.nbeam,
368 hdr.observer,
369 hdr.project,
370 hdr.obstype,
371 hdr.antennaname,
372 hdr.antennaposition,
373 hdr.equinox,
374 hdr.freqref,
375 hdr.reffreq,
376 hdr.bandwidth,
377 hdr.utc,
378 hdr.fluxunit,
379 hdr.epoch,
380 hdr.poltype ) ;
381
382 if ( hdr.freqref != "LSRK" ) {
383// String freqref ;
384// if (hdr.freqref == "TOPOCENT") {
385// freqref = "TOPO";
386// } else if (hdr.freqref == "GEOCENTR") {
387// freqref = "GEO";
388// } else if (hdr.freqref == "BARYCENT") {
389// freqref = "BARY";
390// } else if (hdr.freqref == "GALACTOC") {
391// freqref = "GALACTO";
392// } else if (hdr.freqref == "LOCALGRP") {
393// freqref = "LGROUP";
394// } else if (hdr.freqref == "CMBDIPOL") {
395// freqref = "CMB";
396// } else if (hdr.freqref == "SOURCE") {
397// freqref = "REST";
398// }
399 vector<double> sdir ;
400 string ref ;
401 reader_->getSourceDirection( sdir, ref ) ;
402 Vector<casa::Double> sourceDir( sdir ) ;
403 hdr.reffreq = toLSRK( hdr.reffreq, hdr.freqref, hdr.utc, hdr.antennaposition, sdir, String(ref) ) ;
404 hdr.freqref = "LSRK" ;
405 }
406
407 setHeader( hdr ) ;
408}
409
410String ASDMFiller::getIFKey( uInt ifno )
411{
412 return "IFNO"+String::toString( ifno ) ;
413}
414
415void ASDMFiller::getFrequencyRec( String key,
416 double &refpix,
417 double &refval,
418 double &incr )
419{
420 Record frec = ifrec_.asRecord( key ) ;
421 refpix = frec.asdouble( "REFPIX" ) ;
422 refval = frec.asdouble( "REFVAL" ) ;
423 incr = frec.asdouble( "INCREMENT" ) ;
424}
425
426void ASDMFiller::setFrequencyRec( String key,
427 double refpix,
428 double refval,
429 double incr )
430{
431 Record frec ;
432 frec.define( "REFPIX", refpix ) ;
433 frec.define( "REFVAL", refval ) ;
434 frec.define( "INCREMENT", incr ) ;
435 ifrec_.defineRecord( key, frec ) ;
436}
437
438Matrix<casa::Float> ASDMFiller::toMatrix( float *sp,
439 unsigned int npol,
440 unsigned int nchan )
441{
442 Matrix<casa::Float> mSp( npol, nchan ) ;
443 if ( npol <= 2 ) {
444 // 1 or 2 polarization case
445 for ( unsigned int ich = 0 ; ich < nchan ; ich++ ) {
446 for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ ) {
447 mSp(ipol,ich) = (casa::Float)(sp[npol*ich+ipol]) ;
448 }
449 }
450 }
451 else {
452 // 4 polarization case
453 for ( unsigned int ich = 0 ; ich < nchan ; ich++ ) {
454 mSp(0,ich) = (casa::Float)(sp[4*ich]) ; // Re(XX)
455 mSp(1,ich) = (casa::Float)(sp[4*ich+4]) ; // Re(YY)
456 mSp(2,ich) = (casa::Float)(sp[4*ich+2]) ; // Re(XY)
457 mSp(3,ich) = (casa::Float)(sp[4*ich+3]) ; // Im(XY)
458 }
459 }
460 return mSp ;
461}
462
463Matrix<casa::Float> ASDMFiller::toMatrix( vector< vector<float> > &tsys,
464 unsigned int npol,
465 unsigned int nchan )
466{
467 unsigned int numRec = tsys.size() ;
468 unsigned int numChan = tsys[0].size() ;
469 Matrix<casa::Float> ret ;
470 if ( npol == numRec && nchan == numChan ) {
471 ret.resize( npol, nchan ) ;
472 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
473 for ( unsigned int ic = 0 ; ic < nchan ; ic++ )
474 ret( ip, ic ) = (casa::Float)(tsys[ip][ic]) ;
475 }
476 else if ( npol == numRec && numChan == 1 ) {
477 ret.resize( npol, 1 ) ;
478 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
479 ret( ip, 0 ) = (casa::Float)(tsys[0][0]) ;
480 }
481 else if ( numRec == 1 && nchan == numChan ) {
482 ret.resize( npol, nchan ) ;
483 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
484 for ( unsigned int ic = 0 ; ic < nchan ; ic++ )
485 ret( ip, ic ) = (casa::Float)(tsys[0][ic]) ;
486 }
487 else if ( numRec == 1 && numChan == 1 ) {
488 ret.resize( npol, 1 ) ;
489 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
490 ret( ip, 0 ) = (casa::Float)(tsys[0][0]) ;
491 }
492 else if ( numRec == 2 && npol == 4 && numChan == nchan ) {
493 // TODO: How to determine Tsys for XY?
494 // at the moment Tsys[XY] = 0.5*(Tsys[X]+Tsys[Y])
495 ret.resize( npol, nchan ) ;
496 for ( unsigned int ic = 0 ; ic < nchan ; ic++ ) {
497 casa::Float tsysxy = (casa::Float)(0.5*(tsys[0][ic]+tsys[1][ic])) ;
498 ret( 0, ic ) = (casa::Float)(tsys[0][ic]) ;
499 ret( 1, ic ) = (casa::Float)(tsys[1][ic]) ;
500 ret( 2, ic ) = tsysxy ;
501 ret( 3, ic ) = tsysxy ;
502 }
503 }
504 else if ( numRec == 2 && npol == 4 && numChan == 1 ) {
505 // TODO: How to determine Tsys for XY?
506 // at the moment Tsys[XY] = 0.5*(Tsys[X]+Tsys[Y])
507 ret.resize( npol, 1 ) ;
508 casa::Float tsysxy = (casa::Float)(0.5*(tsys[0][0]+tsys[1][0])) ;
509 ret( 0, 0 ) = (casa::Float)(tsys[0][0]) ;
510 ret( 1, 0 ) = (casa::Float)(tsys[1][0]) ;
511 ret( 2, 0 ) = tsysxy ;
512 ret( 3, 0 ) = tsysxy ;
513 }
514 else {
515 // I don't know how to handle ...
516 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
517 for ( unsigned int ic = 0 ; ic < nchan ; ic++ )
518 ret( ip, ic ) = (casa::Float)(tsys[0][ic]) ;
519 }
520 return ret ;
521}
522
523Vector<casa::Float> ASDMFiller::toVector( vector<float> &tau,
524 unsigned int npol )
525{
526 String funcName = "toVector" ;
527
528 Vector<casa::Float> ret( npol ) ;
529 //logsink_->postLocally( LogMessage("tau0="+String::toString(tau[0]),LogOrigin(className_,funcName,WHERE)) ) ;
530 if ( npol == 4 ) {
531 ret[0] = (casa::Float)tau[0] ;
532 ret[1] = (casa::Float)tau[1] ;
533 ret[2] = 0.5 * ( ret[0] + ret[1] ) ;
534 ret[3] = ret[2] ;
535 }
536 else if ( npol == tau.size() ) {
537 for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ )
538 ret[ipol] = (casa::Float)tau[ipol] ;
539 }
540 else {
541 // I don't know how to handle...
542 for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ )
543 ret[ipol] = (casa::Float)tau[0] ;
544 }
545 //logsink_->postLocally( LogMessage("tau="+String::toString(ret),LogOrigin(className_,funcName,WHERE)) ) ;
546 return ret ;
547}
548
549String ASDMFiller::toTcalTime( casa::Double mjd )
550{
551 return MVTime( mjd ).string( MVTime::YMD ) ;
552}
553
554void ASDMFiller::toJ2000( Vector<casa::Double> &dir,
555 double az,
556 double el,
557 casa::Double mjd,
558 Vector<casa::Double> antpos )
559{
560 String funcName = "toJ2000" ;
561
562 Vector<casa::Double> azel( 2 ) ;
563 azel[0] = az ;
564 azel[1] = el ;
565// MEpoch me( Quantity( mjd, "d" ), MEpoch::UTC ) ;
566// Vector<Quantity> qantpos( 3 ) ;
567// qantpos[0] = Quantity( antpos[0], "m" ) ;
568// qantpos[1] = Quantity( antpos[1], "m" ) ;
569// qantpos[2] = Quantity( antpos[2], "m" ) ;
570// MPosition mp( MVPosition( qantpos ),
571// MPosition::ITRF ) ;
572// // mp.print( os_.output() ) ;
573// MeasFrame mf( me, mp ) ;
574// MDirection::Convert toj2000( MDirection::AZELGEO,
575// MDirection::Ref( MDirection::J2000, mf ) ) ;
576// dir = toj2000( azel ).getAngle( "rad" ).getValue() ;
577 dir = toJ2000( azel, "AZELGEO", mjd, antpos ) ;
578 //logsink_->postLocally( LogMessage("dir = "+String::toString(dir),LogOrigin(className_,funcName,WHERE)) ) ;
579}
580
581Vector<casa::Double> ASDMFiller::toJ2000( Vector<casa::Double> dir,
582 String dirref,
583 casa::Double mjd,
584 Vector<casa::Double> antpos )
585{
586 Vector<casa::Double> newd( dir ) ;
587 if ( dirref != "J2000" ) {
588 MEpoch me( Quantity( mjd, "d" ), MEpoch::UTC ) ;
589 Vector<Quantity> qantpos( 3 ) ;
590 qantpos[0] = Quantity( antpos[0], "m" ) ;
591 qantpos[1] = Quantity( antpos[1], "m" ) ;
592 qantpos[2] = Quantity( antpos[2], "m" ) ;
593 MPosition mp( MVPosition( qantpos ),
594 MPosition::ITRF ) ;
595 // mp.print( os_.output() ) ;
596 MeasFrame mf( me, mp ) ;
597 MDirection::Types dirtype ;
598 Bool b = MDirection::getType( dirtype, dirref ) ;
599 if ( b ) {
600 MDirection::Convert toj2000( dirtype,
601 MDirection::Ref( MDirection::J2000, mf ) ) ;
602 newd = toj2000( dir ).getAngle( "rad" ).getValue() ;
603 }
604 }
605 return newd ;
606}
607
608MFrequency::Types ASDMFiller::toFrameType( string &s )
609{
610 MFrequency::Types ftype = MFrequency::DEFAULT ;
611 if ( s == "LABREST" )
612 ftype = MFrequency::REST ;
613 else {
614 Bool b = MFrequency::getType( ftype, String(s) ) ;
615 if (!b)
616 ftype = MFrequency::DEFAULT ;
617 }
618 return ftype ;
619}
620
621casa::Double ASDMFiller::toLSRK( casa::Double freq,
622 String freqref,
623 casa::Double utc,
624 Vector<casa::Double> antpos,
625 Vector<casa::Double> dir,
626 String dirref )
627{
628 String funcName = "toLSRK" ;
629
630 //logsink_->postLocally( LogMessage("freqref = "+freqref,LogOrigin(className_,funcName,WHERE)) ) ;
631 casa::Double newf = freq ;
632 if ( freqToLsr_ && freqref != "LSRK" ) {
633 MEpoch me( Quantum<casa::Double>( utc, Unit("d") ), MEpoch::UTC ) ;
634 Vector< Quantum<casa::Double> > antposQ( 3 ) ;
635 for ( int i = 0 ; i < 3 ; i++ )
636 antposQ[i] = Quantum<casa::Double>( antpos[i], Unit("m") ) ;
637 MPosition mp( antposQ, MPosition::ITRF ) ;
638 MDirection::Types dirtype ;
639 Bool b = MDirection::getType( dirtype, dirref ) ;
640 if ( !b )
641 dirtype = MDirection::J2000 ;
642 MDirection md( Quantum<casa::Double>( dir[0], Unit("rad") ),
643 Quantum<casa::Double>( dir[1], Unit("rad") ),
644 dirtype ) ;
645 MeasFrame mf( me, mp, md ) ;
646 MFrequency::Types freqtype ;
647 b = MFrequency::getType( freqtype, freqref ) ;
648 if ( !b )
649 freqtype = MFrequency::TOPO ;
650 MFrequency::Convert tolsr( freqtype,
651 MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
652 newf = tolsr( Quantum<casa::Double>( freq, Unit("Hz") ) ).get( "Hz" ).getValue() ;
653 //logsink_->postLocally( LogMessage("freq = "+String::toString(freq)+", newf = "+String::toString(newf),LogOrigin(className_,funcName,WHERE)) ) ;
654 }
655 return newf ;
656}
Note: See TracBrowser for help on using the repository browser.