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

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