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

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_importasdm_sd

Put in Release Notes: es/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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