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

Last change on this file since 2710 was 2710, checked in by Takeshi Nakazato, 11 years ago

New Development: No

JIRA Issue: No

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

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