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

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_importasdm_sd

Put in Release Notes: es/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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