source: branches/casa-prerelease/pre-asap/external-alma/oldasdm2ASAP/OldASDMFiller.cc @ 2302

Last change on this file since 2302 was 2302, 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: test_importasdm_sd

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

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