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

Last change on this file since 2227 was 2227, 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...

Frequency is stored as LSRK value although original frequency is
in an arbitrary frequency frame (e.g. TOPO). Conversion is done
in filler if necessary.
File size: 22.0 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  for ( uInt icon = 0 ; icon < numConfigDescId ; icon++ ) {
105    for ( unsigned int ifield = 0 ; ifield < numFieldId ; ifield++ ) {
106      //logsink_->postLocally( LogMessage("start configDescId "+String::toString(configDescIdList[icon])+" fieldId "+String::toString(fieldIdList[ifield]),LogOrigin(className_,funcName,WHERE)) ) ;
107
108      //Bool status = reader_->setMainRow( configDescIdList[icon], fieldIdList[ifield] ) ;
109      if ( !(reader_->setMainRow( configDescIdList[icon], fieldIdList[ifield] )) ) {
110        //logsink_->postLocally( LogMessage("skip configDescId "+String::toString(configDescIdList[icon])+" fieldId "+String::toString(fieldIdList[ifield]),LogOrigin(className_,funcName,WHERE)) ) ;
111        continue ;
112      }
113
114      // number of rows
115      uInt nrow = reader_->getNumMainRow() ;
116
117      //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)) ) ;
118     
119      // CYCLENO
120      unsigned int cycleno = 0 ;
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),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),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 = reader_->getTsys( idata ) ;
288          vector< vector<float> > tc = reader_->getTcal( idata ) ;
289          Matrix<casa::Float> spectra = toMatrix( sp, numPol, numChan ) ;
290          Vector<uChar> flagtra( numChan, 0 ) ;
291          Matrix<casa::Float> tsys = toMatrix( ts, numPol, numChan ) ;
292          Matrix<casa::Float> tcal = toMatrix( tc, numPol, numChan ) ;
293//           String caltime = "" ;
294//           if ( anyNE( tcal, (casa::Float)1.0 ) )
295//             caltime = toTcalTime( mjd ) ;
296          String caltime = toTcalTime( mjd ) ;
297
298          for ( unsigned int ipol = 0 ; ipol < numPol ; ipol++ ) {
299
300            // fill SCANNO, CYCLENO, IFNO, POLNO, and BEAMNO
301            setIndex( (uInt)scanno-1, (uInt)cycleno, ifno, ipol, beamno ) ;
302
303            // fill SPECTRA, FLAGTRA, TSYS
304            setSpectrum( spectra.row(ipol), flagtra, tsys.row(ipol) ) ;
305
306            // fill TCAL_ID and add TCAL row if necessary
307            setTcal2( caltime, tcal.row(ipol) ) ;
308
309            // fill OPACITY
310            setOpacity( opacity[ipol] ) ;
311           
312            // commit row
313            commitRow() ;
314          }
315
316          // increment CYCLENO
317          cycleno++ ;
318        }
319      }
320    }
321  }
322
323  return ;
324}
325
326void ASDMFiller::close()
327{
328  reader_->close() ;
329  reader_ = 0 ;
330
331  return ;
332}
333
334void ASDMFiller::fillHeader()
335{
336  STHeader hdr ;
337
338  reader_->fillHeader( hdr.nchan,
339                       hdr.npol,
340                       hdr.nif,
341                       hdr.nbeam,
342                       hdr.observer,
343                       hdr.project,
344                       hdr.obstype,
345                       hdr.antennaname,
346                       hdr.antennaposition,
347                       hdr.equinox,
348                       hdr.freqref,
349                       hdr.reffreq,
350                       hdr.bandwidth,
351                       hdr.utc,
352                       hdr.fluxunit,
353                       hdr.epoch,
354                       hdr.poltype ) ;
355
356  if ( hdr.freqref != "LSRK" ) {
357//     String freqref ;
358//     if (hdr.freqref == "TOPOCENT") {
359//       freqref = "TOPO";
360//     } else if (hdr.freqref == "GEOCENTR") {
361//       freqref = "GEO";
362//     } else if (hdr.freqref == "BARYCENT") {
363//       freqref = "BARY";
364//     } else if (hdr.freqref == "GALACTOC") {
365//       freqref = "GALACTO";
366//     } else if (hdr.freqref == "LOCALGRP") {
367//       freqref = "LGROUP";
368//     } else if (hdr.freqref == "CMBDIPOL") {
369//       freqref = "CMB";
370//     } else if (hdr.freqref == "SOURCE") {
371//       freqref = "REST";
372//     }
373    vector<double> sdir ;
374    string ref ;
375    reader_->getSourceDirection( sdir, ref ) ;
376    Vector<casa::Double> sourceDir( sdir ) ;
377    hdr.reffreq = toLSRK( hdr.reffreq, hdr.freqref, hdr.utc, hdr.antennaposition, sdir, String(ref) ) ;
378    hdr.freqref = "LSRK" ;
379  }
380
381  setHeader( hdr ) ;
382}
383
384String ASDMFiller::getIFKey( uInt ifno )
385{
386  return "IFNO"+String::toString( ifno ) ;
387}
388
389void ASDMFiller::getFrequencyRec( String key,
390                                       double &refpix,
391                                       double &refval,
392                                       double &incr )
393{
394  Record frec = ifrec_.asRecord( key ) ;
395  refpix = frec.asdouble( "REFPIX" ) ;
396  refval = frec.asdouble( "REFVAL" ) ;
397  incr = frec.asdouble( "INCREMENT" ) ;
398}
399
400void ASDMFiller::setFrequencyRec( String key,
401                                       double refpix,
402                                       double refval,
403                                       double incr )
404{
405  Record frec ;
406  frec.define( "REFPIX", refpix ) ;
407  frec.define( "REFVAL", refval ) ;
408  frec.define( "INCREMENT", incr ) ;
409  ifrec_.defineRecord( key, frec ) ;
410}
411
412Matrix<casa::Float> ASDMFiller::toMatrix( float *sp,
413                                         unsigned int npol,
414                                         unsigned int nchan )
415{
416  Matrix<casa::Float> mSp( npol, nchan ) ;
417  if ( npol <= 2 ) {
418    // 1 or 2 polarization case
419    for ( unsigned int ich = 0 ; ich < nchan ; ich++ ) {
420      for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ ) {
421        mSp(ipol,ich) = (casa::Float)(sp[npol*ich+ipol]) ;
422      }
423    }
424  }
425  else {
426    // 4 polarization case
427    for ( unsigned int ich = 0 ; ich < nchan ; ich++ ) {
428      mSp(0,ich) = (casa::Float)(sp[4*ich]) ;   // Re(XX)
429      mSp(1,ich) = (casa::Float)(sp[4*ich+4]) ; // Re(YY)
430      mSp(2,ich) = (casa::Float)(sp[4*ich+2]) ; // Re(XY)
431      mSp(3,ich) = (casa::Float)(sp[4*ich+3]) ; // Im(XY)
432    }
433  }
434  return mSp ;
435}
436
437Matrix<casa::Float> ASDMFiller::toMatrix( vector< vector<float> > &tsys,
438                                               unsigned int npol,
439                                               unsigned int nchan )
440{
441  unsigned int numRec = tsys.size() ;
442  unsigned int numChan = tsys[0].size() ;
443  Matrix<casa::Float> ret ;
444  if ( npol == numRec && nchan == numChan ) {
445    ret.resize( npol, nchan ) ;
446    for ( unsigned int ip = 0 ; ip < npol ; ip++ )
447      for ( unsigned int ic = 0 ; ic < nchan ; ic++ )
448        ret( ip, ic ) = (casa::Float)(tsys[ip][ic]) ;
449  }
450  else if ( npol == numRec && numChan == 1 ) {
451    ret.resize( npol, 1 ) ;
452    for ( unsigned int ip = 0 ; ip < npol ; ip++ )
453      ret( ip, 0 ) = (casa::Float)(tsys[0][0]) ;
454  }
455  else if ( numRec == 1 && 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[0][ic]) ;
460  }
461  else if ( numRec == 1 && 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 == 2 && npol == 4 && numChan == nchan ) {
467    // TODO: How to determine Tsys for XY?
468    //       at the moment Tsys[XY] = 0.5*(Tsys[X]+Tsys[Y])
469    ret.resize( npol, nchan ) ;
470    for ( unsigned int ic = 0 ; ic < nchan ; ic++ ) {
471      casa::Float tsysxy = (casa::Float)(0.5*(tsys[0][ic]+tsys[1][ic])) ;
472      ret( 0, ic ) = (casa::Float)(tsys[0][ic]) ;
473      ret( 1, ic ) = (casa::Float)(tsys[1][ic]) ;
474      ret( 2, ic ) = tsysxy ;
475      ret( 3, ic ) = tsysxy ;
476    }
477  }
478  else if ( numRec == 2 && npol == 4 && numChan == 1 ) {
479    // TODO: How to determine Tsys for XY?
480    //       at the moment Tsys[XY] = 0.5*(Tsys[X]+Tsys[Y])
481    ret.resize( npol, 1 ) ;
482    casa::Float tsysxy = (casa::Float)(0.5*(tsys[0][0]+tsys[1][0])) ;
483    ret( 0, 0 ) = (casa::Float)(tsys[0][0]) ;
484    ret( 1, 0 ) = (casa::Float)(tsys[1][0]) ;
485    ret( 2, 0 ) = tsysxy ;
486    ret( 3, 0 ) = tsysxy ;
487  }
488  else {
489    // I don't know how to handle ...
490    for ( unsigned int ip = 0 ; ip < npol ; ip++ )
491      for ( unsigned int ic = 0 ; ic < nchan ; ic++ )
492        ret( ip, ic ) = (casa::Float)(tsys[0][ic]) ;   
493  }
494  return ret ;
495}
496
497Vector<casa::Float> ASDMFiller::toVector( vector<float> &tau,
498                                               unsigned int npol )
499{
500  Vector<casa::Float> ret( npol ) ;
501
502  if ( npol == 4 ) {
503    ret[0] = (casa::Float)tau[0] ;
504    ret[1] = (casa::Float)tau[1] ;
505    ret[2] = 0.5 * ( ret[0] + ret[1] ) ;
506    ret[3] = ret[2] ;
507  }
508  else if ( npol == tau.size() ) {
509    for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ )
510      ret[ipol] = (casa::Float)tau[ipol] ;
511  }
512  else {
513    // I don't know how to handle...
514    for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ )
515      ret[ipol] = (casa::Float)tau[0] ;
516  }
517  return ret ;
518}
519
520String ASDMFiller::toTcalTime( casa::Double mjd )
521{
522  return MVTime( mjd ).string( MVTime::YMD ) ;
523}
524
525void ASDMFiller::toJ2000( Vector<casa::Double> &dir,
526                               double az,
527                               double el,
528                               casa::Double mjd,
529                               Vector<casa::Double> antpos )
530{
531  String funcName = "toJ2000" ;
532
533  Vector<casa::Double> azel( 2 ) ;
534  azel[0] = az ;
535  azel[1] = el ;
536//   MEpoch me( Quantity( mjd, "d" ), MEpoch::UTC ) ;
537//   Vector<Quantity> qantpos( 3 ) ;
538//   qantpos[0] = Quantity( antpos[0], "m" ) ;
539//   qantpos[1] = Quantity( antpos[1], "m" ) ;
540//   qantpos[2] = Quantity( antpos[2], "m" ) ;
541//   MPosition mp( MVPosition( qantpos ),
542//                 MPosition::ITRF ) ;
543// //   mp.print( os_.output() ) ;
544//   MeasFrame mf( me, mp ) ;
545//   MDirection::Convert toj2000( MDirection::AZELGEO,
546//                                MDirection::Ref( MDirection::J2000, mf ) ) ;
547//   dir = toj2000( azel ).getAngle( "rad" ).getValue() ;
548  dir = toJ2000( azel, "AZELGEO", mjd, antpos ) ;
549  //logsink_->postLocally( LogMessage("dir = "+String::toString(dir),LogOrigin(className_,funcName,WHERE)) ) ;
550}
551
552Vector<casa::Double> ASDMFiller::toJ2000( Vector<casa::Double> dir,
553                                          String dirref,
554                                          casa::Double mjd,
555                                          Vector<casa::Double> antpos )
556{
557  Vector<casa::Double> newd( dir ) ;
558  if ( dirref != "J2000" ) {
559    MEpoch me( Quantity( mjd, "d" ), MEpoch::UTC ) ;
560    Vector<Quantity> qantpos( 3 ) ;
561    qantpos[0] = Quantity( antpos[0], "m" ) ;
562    qantpos[1] = Quantity( antpos[1], "m" ) ;
563    qantpos[2] = Quantity( antpos[2], "m" ) ;
564    MPosition mp( MVPosition( qantpos ),
565                  MPosition::ITRF ) ;
566    //   mp.print( os_.output() ) ;
567    MeasFrame mf( me, mp ) ;
568    MDirection::Types dirtype ;
569    Bool b = MDirection::getType( dirtype, dirref ) ;
570    if ( b ) {
571      MDirection::Convert toj2000( dirtype,
572                                   MDirection::Ref( MDirection::J2000, mf ) ) ;
573      newd = toj2000( dir ).getAngle( "rad" ).getValue() ;
574    }
575  }
576  return newd ;
577}
578
579MFrequency::Types ASDMFiller::toFrameType( string &s )
580{
581  MFrequency::Types ftype = MFrequency::DEFAULT ;
582  if ( s == "LABREST" )
583    ftype = MFrequency::REST ;
584  else {
585    Bool b = MFrequency::getType( ftype, String(s) ) ;
586    if (!b)
587      ftype = MFrequency::DEFAULT ;
588  }
589  return ftype ;
590}
591
592casa::Double ASDMFiller::toLSRK( casa::Double freq,
593                                 String freqref,
594                                 casa::Double utc,
595                                 Vector<casa::Double> antpos,
596                                 Vector<casa::Double> dir,
597                                 String dirref )
598{
599  String funcName = "toLSRK" ;
600
601  //logsink_->postLocally( LogMessage("freqref = "+freqref,LogOrigin(className_,funcName,WHERE)) ) ;
602  casa::Double newf = freq ;
603  if ( freqref != "LSRK" ) {
604    MEpoch me( Quantum<casa::Double>( utc, Unit("d") ), MEpoch::UTC ) ;
605    Vector< Quantum<casa::Double> > antposQ( 3 ) ;
606    for ( int i = 0 ; i < 3 ; i++ )
607      antposQ[i] = Quantum<casa::Double>( antpos[i], Unit("m") ) ;
608    MPosition mp( antposQ, MPosition::ITRF ) ;
609    MDirection::Types dirtype ;
610    Bool b = MDirection::getType( dirtype, dirref ) ;
611    if ( !b )
612      dirtype = MDirection::J2000 ;
613    MDirection md( Quantum<casa::Double>( dir[0], Unit("rad") ),
614                   Quantum<casa::Double>( dir[1], Unit("rad") ),
615                   dirtype ) ;
616    MeasFrame mf( me, mp, md ) ;
617    MFrequency::Types freqtype ;
618    b = MFrequency::getType( freqtype, freqref ) ;
619    if ( !b )
620      freqtype = MFrequency::TOPO ;
621    MFrequency::Convert tolsr( freqtype,
622                               MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
623    newf = tolsr( Quantum<casa::Double>( freq, Unit("Hz") ) ).get( "Hz" ).getValue() ;
624    //logsink_->postLocally( LogMessage("freq = "+String::toString(freq)+", newf = "+String::toString(newf),LogOrigin(className_,funcName,WHERE)) ) ;
625  }
626  return newf ;
627}
Note: See TracBrowser for help on using the repository browser.