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

Last change on this file since 3123 was 3106, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/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...


Check-in asap modifications from Jim regarding casacore namespace conversion.
File size: 24.6 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<casacore::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<casacore::Double> restFreqs( rf.size() ) ;
197 for ( uInt i = 0 ; i < rf.size() ; i++ )
198 restFreqs[i] = (casacore::Double)(rf[i]) ;
199 setMolecule( restFreqs ) ;
200
201 // time and interval
202 casacore::Double mjd = (casacore::Double)(reader_->getTime()) ;
203 casacore::Double interval = (casacore::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<casacore::Double> srcDir( 2 ) ;
211 srcDir[0] = (casacore::Double)(srcDirection[0]) ;
212 srcDir[1] = (casacore::Double)(srcDirection[1]) ;
213 Vector<casacore::Double> srcPM( 2 ) ;
214 srcPM[0] = (casacore::Double)(srcProperMotion[0]) ;
215 srcPM[1] = (casacore::Double)(srcProperMotion[1]) ;
216 setSource( srcname, srctype, fieldname, srcDir, srcPM, (casacore::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( (casacore::Float)temperature,
234 (casacore::Float)pressure,
235 (casacore::Float)humidity,
236 (casacore::Float)windspeed,
237 (casacore::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<casacore::Double> scanRate( 2, 0.0 ) ;
249 Vector<casacore::Double> direction( 2, 0.0 ) ;
250 if ( srate.size() > 0 ) {
251 scanRate[0] = (casacore::Double)(srate[0]) ;
252 scanRate[1] = (casacore::Double)(srate[1]) ;
253 }
254 setScanRate( scanRate ) ;
255 if ( dir.size() > 0 ) {
256 direction[0] = (casacore::Double)(dir[0]) ;
257 direction[1] = (casacore::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, (casacore::Float)az, (casacore::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( casacore::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( (casacore::Double)refpix, (casacore::Double)refval, (casacore::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<casacore::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<casacore::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<casacore::Float> tsys = toMatrix( ts, numPol, numChan ) ;
317 Matrix<casacore::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, (casacore::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 // CAS-5841: SCANNO should be consistent with ASDM scanNumber
328 setIndex( (uInt)scanno, (uInt)cycleno[scanno], ifno, ipol, beamno ) ;
329
330 // fill SPECTRA, FLAGTRA, TSYS
331 setSpectrum( spectra.row(ipol), flagtra, tsys.row(ipol) ) ;
332
333 // fill TCAL_ID and add TCAL row if necessary
334 setTcal2( caltime, tcal.row(ipol) ) ;
335
336 // fill OPACITY
337 setOpacity( opacity[ipol] ) ;
338
339 // commit row
340 commitRow() ;
341 }
342
343 // increment CYCLENO
344 cycleno[scanno]++ ;
345 }
346 }
347 }
348 }
349
350 return ;
351}
352
353void ASDMFiller::close()
354{
355 reader_->close() ;
356 reader_ = 0 ;
357
358 return ;
359}
360
361void ASDMFiller::fillHeader()
362{
363 STHeader hdr ;
364
365 reader_->fillHeader( hdr.nchan,
366 hdr.npol,
367 hdr.nif,
368 hdr.nbeam,
369 hdr.observer,
370 hdr.project,
371 hdr.obstype,
372 hdr.antennaname,
373 hdr.antennaposition,
374 hdr.equinox,
375 hdr.freqref,
376 hdr.reffreq,
377 hdr.bandwidth,
378 hdr.utc,
379 hdr.fluxunit,
380 hdr.epoch,
381 hdr.poltype ) ;
382
383 if ( hdr.freqref != "LSRK" ) {
384// String freqref ;
385// if (hdr.freqref == "TOPOCENT") {
386// freqref = "TOPO";
387// } else if (hdr.freqref == "GEOCENTR") {
388// freqref = "GEO";
389// } else if (hdr.freqref == "BARYCENT") {
390// freqref = "BARY";
391// } else if (hdr.freqref == "GALACTOC") {
392// freqref = "GALACTO";
393// } else if (hdr.freqref == "LOCALGRP") {
394// freqref = "LGROUP";
395// } else if (hdr.freqref == "CMBDIPOL") {
396// freqref = "CMB";
397// } else if (hdr.freqref == "SOURCE") {
398// freqref = "REST";
399// }
400 vector<double> sdir ;
401 string ref ;
402 reader_->getSourceDirection( sdir, ref ) ;
403 Vector<casacore::Double> sourceDir( sdir ) ;
404 hdr.reffreq = toLSRK( hdr.reffreq, hdr.freqref, hdr.utc, hdr.antennaposition, sdir, String(ref) ) ;
405 hdr.freqref = "LSRK" ;
406 }
407
408 setHeader( hdr ) ;
409}
410
411String ASDMFiller::getIFKey( uInt ifno )
412{
413 return "IFNO"+String::toString( ifno ) ;
414}
415
416void ASDMFiller::getFrequencyRec( String key,
417 double &refpix,
418 double &refval,
419 double &incr )
420{
421 Record frec = ifrec_.asRecord( key ) ;
422 refpix = frec.asdouble( "REFPIX" ) ;
423 refval = frec.asdouble( "REFVAL" ) ;
424 incr = frec.asdouble( "INCREMENT" ) ;
425}
426
427void ASDMFiller::setFrequencyRec( String key,
428 double refpix,
429 double refval,
430 double incr )
431{
432 Record frec ;
433 frec.define( "REFPIX", refpix ) ;
434 frec.define( "REFVAL", refval ) ;
435 frec.define( "INCREMENT", incr ) ;
436 ifrec_.defineRecord( key, frec ) ;
437}
438
439Matrix<casacore::Float> ASDMFiller::toMatrix( float *sp,
440 unsigned int npol,
441 unsigned int nchan )
442{
443 Matrix<casacore::Float> mSp( npol, nchan ) ;
444 if ( npol <= 2 ) {
445 // 1 or 2 polarization case
446 for ( unsigned int ich = 0 ; ich < nchan ; ich++ ) {
447 for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ ) {
448 mSp(ipol,ich) = (casacore::Float)(sp[npol*ich+ipol]) ;
449 }
450 }
451 }
452 else {
453 // 4 polarization case
454 for ( unsigned int ich = 0 ; ich < nchan ; ich++ ) {
455 mSp(0,ich) = (casacore::Float)(sp[4*ich]) ; // Re(XX)
456 mSp(1,ich) = (casacore::Float)(sp[4*ich+4]) ; // Re(YY)
457 mSp(2,ich) = (casacore::Float)(sp[4*ich+2]) ; // Re(XY)
458 mSp(3,ich) = (casacore::Float)(sp[4*ich+3]) ; // Im(XY)
459 }
460 }
461 return mSp ;
462}
463
464Matrix<casacore::Float> ASDMFiller::toMatrix( vector< vector<float> > &tsys,
465 unsigned int npol,
466 unsigned int nchan )
467{
468 unsigned int numRec = tsys.size() ;
469 unsigned int numChan = tsys[0].size() ;
470 Matrix<casacore::Float> ret ;
471 if ( npol == numRec && nchan == numChan ) {
472 ret.resize( npol, nchan ) ;
473 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
474 for ( unsigned int ic = 0 ; ic < nchan ; ic++ )
475 ret( ip, ic ) = (casacore::Float)(tsys[ip][ic]) ;
476 }
477 else if ( npol == numRec && numChan == 1 ) {
478 ret.resize( npol, 1 ) ;
479 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
480 ret( ip, 0 ) = (casacore::Float)(tsys[0][0]) ;
481 }
482 else if ( numRec == 1 && nchan == numChan ) {
483 ret.resize( npol, nchan ) ;
484 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
485 for ( unsigned int ic = 0 ; ic < nchan ; ic++ )
486 ret( ip, ic ) = (casacore::Float)(tsys[0][ic]) ;
487 }
488 else if ( numRec == 1 && numChan == 1 ) {
489 ret.resize( npol, 1 ) ;
490 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
491 ret( ip, 0 ) = (casacore::Float)(tsys[0][0]) ;
492 }
493 else if ( numRec == 2 && npol == 4 && numChan == nchan ) {
494 // TODO: How to determine Tsys for XY?
495 // at the moment Tsys[XY] = 0.5*(Tsys[X]+Tsys[Y])
496 ret.resize( npol, nchan ) ;
497 for ( unsigned int ic = 0 ; ic < nchan ; ic++ ) {
498 casacore::Float tsysxy = (casacore::Float)(0.5*(tsys[0][ic]+tsys[1][ic])) ;
499 ret( 0, ic ) = (casacore::Float)(tsys[0][ic]) ;
500 ret( 1, ic ) = (casacore::Float)(tsys[1][ic]) ;
501 ret( 2, ic ) = tsysxy ;
502 ret( 3, ic ) = tsysxy ;
503 }
504 }
505 else if ( numRec == 2 && npol == 4 && numChan == 1 ) {
506 // TODO: How to determine Tsys for XY?
507 // at the moment Tsys[XY] = 0.5*(Tsys[X]+Tsys[Y])
508 ret.resize( npol, 1 ) ;
509 casacore::Float tsysxy = (casacore::Float)(0.5*(tsys[0][0]+tsys[1][0])) ;
510 ret( 0, 0 ) = (casacore::Float)(tsys[0][0]) ;
511 ret( 1, 0 ) = (casacore::Float)(tsys[1][0]) ;
512 ret( 2, 0 ) = tsysxy ;
513 ret( 3, 0 ) = tsysxy ;
514 }
515 else {
516 // I don't know how to handle ...
517 for ( unsigned int ip = 0 ; ip < npol ; ip++ )
518 for ( unsigned int ic = 0 ; ic < nchan ; ic++ )
519 ret( ip, ic ) = (casacore::Float)(tsys[0][ic]) ;
520 }
521 return ret ;
522}
523
524Vector<casacore::Float> ASDMFiller::toVector( vector<float> &tau,
525 unsigned int npol )
526{
527 String funcName = "toVector" ;
528
529 Vector<casacore::Float> ret( npol ) ;
530 //logsink_->postLocally( LogMessage("tau0="+String::toString(tau[0]),LogOrigin(className_,funcName,WHERE)) ) ;
531 if ( npol == 4 ) {
532 ret[0] = (casacore::Float)tau[0] ;
533 ret[1] = (casacore::Float)tau[1] ;
534 ret[2] = 0.5 * ( ret[0] + ret[1] ) ;
535 ret[3] = ret[2] ;
536 }
537 else if ( npol == tau.size() ) {
538 for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ )
539 ret[ipol] = (casacore::Float)tau[ipol] ;
540 }
541 else {
542 // I don't know how to handle...
543 for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ )
544 ret[ipol] = (casacore::Float)tau[0] ;
545 }
546 //logsink_->postLocally( LogMessage("tau="+String::toString(ret),LogOrigin(className_,funcName,WHERE)) ) ;
547 return ret ;
548}
549
550String ASDMFiller::toTcalTime( casacore::Double mjd )
551{
552 return MVTime( mjd ).string( MVTime::YMD ) ;
553}
554
555void ASDMFiller::toJ2000( Vector<casacore::Double> &dir,
556 double az,
557 double el,
558 casacore::Double mjd,
559 Vector<casacore::Double> antpos )
560{
561 String funcName = "toJ2000" ;
562
563 Vector<casacore::Double> azel( 2 ) ;
564 azel[0] = az ;
565 azel[1] = el ;
566// MEpoch me( Quantity( mjd, "d" ), MEpoch::UTC ) ;
567// Vector<Quantity> qantpos( 3 ) ;
568// qantpos[0] = Quantity( antpos[0], "m" ) ;
569// qantpos[1] = Quantity( antpos[1], "m" ) ;
570// qantpos[2] = Quantity( antpos[2], "m" ) ;
571// MPosition mp( MVPosition( qantpos ),
572// MPosition::ITRF ) ;
573// // mp.print( os_.output() ) ;
574// MeasFrame mf( me, mp ) ;
575// MDirection::Convert toj2000( MDirection::AZELGEO,
576// MDirection::Ref( MDirection::J2000, mf ) ) ;
577// dir = toj2000( azel ).getAngle( "rad" ).getValue() ;
578 dir = toJ2000( azel, "AZELGEO", mjd, antpos ) ;
579 //logsink_->postLocally( LogMessage("dir = "+String::toString(dir),LogOrigin(className_,funcName,WHERE)) ) ;
580}
581
582Vector<casacore::Double> ASDMFiller::toJ2000( Vector<casacore::Double> dir,
583 String dirref,
584 casacore::Double mjd,
585 Vector<casacore::Double> antpos )
586{
587 Vector<casacore::Double> newd( dir ) ;
588 if ( dirref != "J2000" ) {
589 MEpoch me( Quantity( mjd, "d" ), MEpoch::UTC ) ;
590 Vector<Quantity> qantpos( 3 ) ;
591 qantpos[0] = Quantity( antpos[0], "m" ) ;
592 qantpos[1] = Quantity( antpos[1], "m" ) ;
593 qantpos[2] = Quantity( antpos[2], "m" ) ;
594 MPosition mp( MVPosition( qantpos ),
595 MPosition::ITRF ) ;
596 // mp.print( os_.output() ) ;
597 MeasFrame mf( me, mp ) ;
598 MDirection::Types dirtype ;
599 Bool b = MDirection::getType( dirtype, dirref ) ;
600 if ( b ) {
601 MDirection::Convert toj2000( dirtype,
602 MDirection::Ref( MDirection::J2000, mf ) ) ;
603 newd = toj2000( dir ).getAngle( "rad" ).getValue() ;
604 }
605 }
606 return newd ;
607}
608
609MFrequency::Types ASDMFiller::toFrameType( string &s )
610{
611 MFrequency::Types ftype = MFrequency::DEFAULT ;
612 if ( s == "LABREST" )
613 ftype = MFrequency::REST ;
614 else {
615 Bool b = MFrequency::getType( ftype, String(s) ) ;
616 if (!b)
617 ftype = MFrequency::DEFAULT ;
618 }
619 return ftype ;
620}
621
622casacore::Double ASDMFiller::toLSRK( casacore::Double freq,
623 String freqref,
624 casacore::Double utc,
625 Vector<casacore::Double> antpos,
626 Vector<casacore::Double> dir,
627 String dirref )
628{
629 String funcName = "toLSRK" ;
630
631 //logsink_->postLocally( LogMessage("freqref = "+freqref,LogOrigin(className_,funcName,WHERE)) ) ;
632 casacore::Double newf = freq ;
633 if ( freqToLsr_ && freqref != "LSRK" ) {
634 MEpoch me( Quantum<casacore::Double>( utc, Unit("d") ), MEpoch::UTC ) ;
635 Vector< Quantum<casacore::Double> > antposQ( 3 ) ;
636 for ( int i = 0 ; i < 3 ; i++ )
637 antposQ[i] = Quantum<casacore::Double>( antpos[i], Unit("m") ) ;
638 MPosition mp( antposQ, MPosition::ITRF ) ;
639 MDirection::Types dirtype ;
640 Bool b = MDirection::getType( dirtype, dirref ) ;
641 if ( !b )
642 dirtype = MDirection::J2000 ;
643 MDirection md( Quantum<casacore::Double>( dir[0], Unit("rad") ),
644 Quantum<casacore::Double>( dir[1], Unit("rad") ),
645 dirtype ) ;
646 MeasFrame mf( me, mp, md ) ;
647 MFrequency::Types freqtype ;
648 b = MFrequency::getType( freqtype, freqref ) ;
649 if ( !b )
650 freqtype = MFrequency::TOPO ;
651 MFrequency::Convert tolsr( freqtype,
652 MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
653 newf = tolsr( Quantum<casacore::Double>( freq, Unit("Hz") ) ).get( "Hz" ).getValue() ;
654 //logsink_->postLocally( LogMessage("freq = "+String::toString(freq)+", newf = "+String::toString(newf),LogOrigin(className_,funcName,WHERE)) ) ;
655 }
656 return newf ;
657}
Note: See TracBrowser for help on using the repository browser.