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

Last change on this file since 2407 was 2407, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: Yes/No?

What Interface Changed: Please list interface changes

Test Programs: test_importasdm_sd

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

some clean up

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