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

Last change on this file since 2754 was 2754, checked in by Takeshi Nakazato, 11 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.