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

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

Use getSourceDirection( uInt ), which always returns J2000 direction,
instead of getSourceDirection( uInt, vector<double>, string ).
File size: 22.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{
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),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  String funcName = "toVector" ;
501
502  Vector<casa::Float> ret( npol ) ;
503  //logsink_->postLocally( LogMessage("tau0="+String::toString(tau[0]),LogOrigin(className_,funcName,WHERE)) ) ;
504  if ( npol == 4 ) {
505    ret[0] = (casa::Float)tau[0] ;
506    ret[1] = (casa::Float)tau[1] ;
507    ret[2] = 0.5 * ( ret[0] + ret[1] ) ;
508    ret[3] = ret[2] ;
509  }
510  else if ( npol == tau.size() ) {
511    for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ )
512      ret[ipol] = (casa::Float)tau[ipol] ;
513  }
514  else {
515    // I don't know how to handle...
516    for ( unsigned int ipol = 0 ; ipol < npol ; ipol++ )
517      ret[ipol] = (casa::Float)tau[0] ;
518  }
519  //logsink_->postLocally( LogMessage("tau="+String::toString(ret),LogOrigin(className_,funcName,WHERE)) ) ;
520  return ret ;
521}
522
523String ASDMFiller::toTcalTime( casa::Double mjd )
524{
525  return MVTime( mjd ).string( MVTime::YMD ) ;
526}
527
528void ASDMFiller::toJ2000( Vector<casa::Double> &dir,
529                               double az,
530                               double el,
531                               casa::Double mjd,
532                               Vector<casa::Double> antpos )
533{
534  String funcName = "toJ2000" ;
535
536  Vector<casa::Double> azel( 2 ) ;
537  azel[0] = az ;
538  azel[1] = el ;
539//   MEpoch me( Quantity( mjd, "d" ), MEpoch::UTC ) ;
540//   Vector<Quantity> qantpos( 3 ) ;
541//   qantpos[0] = Quantity( antpos[0], "m" ) ;
542//   qantpos[1] = Quantity( antpos[1], "m" ) ;
543//   qantpos[2] = Quantity( antpos[2], "m" ) ;
544//   MPosition mp( MVPosition( qantpos ),
545//                 MPosition::ITRF ) ;
546// //   mp.print( os_.output() ) ;
547//   MeasFrame mf( me, mp ) ;
548//   MDirection::Convert toj2000( MDirection::AZELGEO,
549//                                MDirection::Ref( MDirection::J2000, mf ) ) ;
550//   dir = toj2000( azel ).getAngle( "rad" ).getValue() ;
551  dir = toJ2000( azel, "AZELGEO", mjd, antpos ) ;
552  //logsink_->postLocally( LogMessage("dir = "+String::toString(dir),LogOrigin(className_,funcName,WHERE)) ) ;
553}
554
555Vector<casa::Double> ASDMFiller::toJ2000( Vector<casa::Double> dir,
556                                          String dirref,
557                                          casa::Double mjd,
558                                          Vector<casa::Double> antpos )
559{
560  Vector<casa::Double> newd( dir ) ;
561  if ( dirref != "J2000" ) {
562    MEpoch me( Quantity( mjd, "d" ), MEpoch::UTC ) ;
563    Vector<Quantity> qantpos( 3 ) ;
564    qantpos[0] = Quantity( antpos[0], "m" ) ;
565    qantpos[1] = Quantity( antpos[1], "m" ) ;
566    qantpos[2] = Quantity( antpos[2], "m" ) ;
567    MPosition mp( MVPosition( qantpos ),
568                  MPosition::ITRF ) ;
569    //   mp.print( os_.output() ) ;
570    MeasFrame mf( me, mp ) ;
571    MDirection::Types dirtype ;
572    Bool b = MDirection::getType( dirtype, dirref ) ;
573    if ( b ) {
574      MDirection::Convert toj2000( dirtype,
575                                   MDirection::Ref( MDirection::J2000, mf ) ) ;
576      newd = toj2000( dir ).getAngle( "rad" ).getValue() ;
577    }
578  }
579  return newd ;
580}
581
582MFrequency::Types ASDMFiller::toFrameType( string &s )
583{
584  MFrequency::Types ftype = MFrequency::DEFAULT ;
585  if ( s == "LABREST" )
586    ftype = MFrequency::REST ;
587  else {
588    Bool b = MFrequency::getType( ftype, String(s) ) ;
589    if (!b)
590      ftype = MFrequency::DEFAULT ;
591  }
592  return ftype ;
593}
594
595casa::Double ASDMFiller::toLSRK( casa::Double freq,
596                                 String freqref,
597                                 casa::Double utc,
598                                 Vector<casa::Double> antpos,
599                                 Vector<casa::Double> dir,
600                                 String dirref )
601{
602  String funcName = "toLSRK" ;
603
604  //logsink_->postLocally( LogMessage("freqref = "+freqref,LogOrigin(className_,funcName,WHERE)) ) ;
605  casa::Double newf = freq ;
606  if ( freqref != "LSRK" ) {
607    MEpoch me( Quantum<casa::Double>( utc, Unit("d") ), MEpoch::UTC ) ;
608    Vector< Quantum<casa::Double> > antposQ( 3 ) ;
609    for ( int i = 0 ; i < 3 ; i++ )
610      antposQ[i] = Quantum<casa::Double>( antpos[i], Unit("m") ) ;
611    MPosition mp( antposQ, MPosition::ITRF ) ;
612    MDirection::Types dirtype ;
613    Bool b = MDirection::getType( dirtype, dirref ) ;
614    if ( !b )
615      dirtype = MDirection::J2000 ;
616    MDirection md( Quantum<casa::Double>( dir[0], Unit("rad") ),
617                   Quantum<casa::Double>( dir[1], Unit("rad") ),
618                   dirtype ) ;
619    MeasFrame mf( me, mp, md ) ;
620    MFrequency::Types freqtype ;
621    b = MFrequency::getType( freqtype, freqref ) ;
622    if ( !b )
623      freqtype = MFrequency::TOPO ;
624    MFrequency::Convert tolsr( freqtype,
625                               MFrequency::Ref( MFrequency::LSRK, mf ) ) ;
626    newf = tolsr( Quantum<casa::Double>( freq, Unit("Hz") ) ).get( "Hz" ).getValue() ;
627    //logsink_->postLocally( LogMessage("freq = "+String::toString(freq)+", newf = "+String::toString(newf),LogOrigin(className_,funcName,WHERE)) ) ;
628  }
629  return newf ;
630}
Note: See TracBrowser for help on using the repository browser.