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

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

New Development: No

JIRA Issue: Yes CAS-1913

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

Support latest ASDM definition in asdm2ASAP.
Now asdm2ASAP is linked to libalma_v3.so to be able to read latest ASDM.

The importasdm task invokes asdm2ASAP when singledish is True.

Older version of asdm2ASAP is renamec as oldasdm2ASAP.
The oldasdm2ASAP and related classes are built by linking to libalma.so.
It is installed but not used by any tasks.


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