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

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

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...

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