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

Last change on this file since 3106 was 3106, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/No?

Interface Changes: Yes/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...


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