source: trunk/src/STFiller.cpp @ 2289

Last change on this file since 2289 was 2289, checked in by ShinnosukeKawakami, 13 years ago

merged parallel branch to trunk
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.0 KB
Line 
1//
2// C++ Implementation: STFiller
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006-2007
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12#include <casa/iostream.h>
13#include <casa/iomanip.h>
14
15#include <casa/Exceptions.h>
16#include <casa/OS/Path.h>
17#include <casa/OS/File.h>
18#include <casa/Quanta/Unit.h>
19#include <casa/Arrays/ArrayMath.h>
20#include <casa/Arrays/ArrayLogical.h>
21#include <casa/Utilities/Regex.h>
22
23#include <casa/Containers/RecordField.h>
24
25#include <tables/Tables/TableRow.h>
26
27#include <measures/Measures/MDirection.h>
28#include <measures/Measures/MeasConvert.h>
29
30#include <atnf/PKSIO/PKSrecord.h>
31#include <atnf/PKSIO/PKSreader.h>
32#ifdef HAS_ALMA
33 #include <casa/System/ProgressMeter.h>
34#endif
35#include <casa/System/ProgressMeter.h>
36#include <atnf/PKSIO/NROReader.h>
37#include <casa/Logging/LogIO.h>
38
39#include <time.h>
40
41
42#include "STDefs.h"
43#include "STAttr.h"
44
45#include "STFiller.h"
46#include "STHeader.h"
47
48using namespace casa;
49
50namespace asap {
51
52STFiller::STFiller() :
53  reader_(0),
54  header_(0),
55  table_(0),
56  refRx_(".*(e|w|_R)$"),
57  nreader_(0)
58{
59}
60
61STFiller::STFiller( CountedPtr< Scantable > stbl ) :
62  reader_(0),
63  header_(0),
64  table_(stbl),
65  refRx_(".*(e|w|_R)$"),
66  nreader_(0)
67{
68}
69
70STFiller::STFiller(const std::string& filename, int whichIF, int whichBeam ) :
71  reader_(0),
72  header_(0),
73  table_(0),
74  refRx_(".*(e|w|_R)$"),
75  nreader_(0)
76{
77  open(filename, "", whichIF, whichBeam);
78}
79
80STFiller::~STFiller()
81{
82  close();
83}
84
85void STFiller::open( const std::string& filename, const std::string& antenna, int whichIF, int whichBeam, casa::Bool getPt )
86{
87  if (table_.null())  {
88    table_ = new Scantable();
89  }
90  if (reader_)  { delete reader_; reader_ = 0; }
91  Bool haveBase, haveSpectra;
92
93  String inName(filename);
94  Path path(inName);
95  inName = path.expandedName();
96
97  File file(inName);
98  if ( !file.exists() ) {
99     throw(AipsError("File does not exist"));
100  }
101  filename_ = inName;
102
103  // Create reader and fill in values for arguments
104  String format;
105  Vector<Bool> beams, ifs;
106  Vector<uInt> nchans,npols;
107
108  //
109  // if isNRO_ is true, try NROReader
110  //
111  // 2008/11/11 Takeshi Nakazato
112  isNRO_ = fileCheck() ;
113  if ( isNRO_ ) {
114    if ( (nreader_ = getNROReader( inName, format )) == 0 ) {
115      throw(AipsError("Creation of NROReader failed")) ;
116    }
117    else {
118      openNRO( whichIF, whichBeam ) ;
119      return ;
120    }
121  }
122  //
123
124  if ( (reader_ = getPKSreader(inName, antenna, 0, 0, format, beams, ifs,
125                              nchans, npols, haveXPol_,haveBase, haveSpectra
126                              )) == 0 )  {
127    throw(AipsError("Creation of PKSreader failed"));
128  }
129  if (!haveSpectra) {
130    delete reader_;
131    reader_ = 0;
132    throw(AipsError("No spectral data in file."));
133    return;
134  }
135  nBeam_ = beams.nelements();
136  nIF_ = ifs.nelements();
137  // Get basic parameters.
138  if ( anyEQ(haveXPol_, True) ) {
139    pushLog("Cross polarization present");
140    for (uInt i=0; i< npols.nelements();++i) {
141      if (npols[i] < 3) npols[i] += 2;// Convert Complex -> 2 Floats
142    }
143  }
144  if (header_) delete header_;
145  header_ = new STHeader();
146  header_->nchan = max(nchans);
147  header_->npol = max(npols);
148  header_->nbeam = nBeam_;
149
150  Int status = reader_->getHeader(header_->observer, header_->project,
151                                  header_->antennaname, header_->antennaposition,
152                                  header_->obstype,
153                                  header_->fluxunit,
154                                  header_->equinox,
155                                  header_->freqref,
156                                  header_->utc, header_->reffreq,
157                                  header_->bandwidth);
158
159  if (status) {
160    delete reader_;
161    reader_ = 0;
162    delete header_;
163    header_ = 0;
164    throw(AipsError("Failed to get header."));
165  }
166  if ((header_->obstype).matches("*SW*")) {
167    // need robust way here - probably read ahead of next timestamp
168    pushLog("Header indicates frequency switched observation.\n"
169               "setting # of IFs = 1 ");
170    nIF_ = 1;
171    header_->obstype = String("fswitch");
172  }
173  // Determine Telescope and set brightness unit
174
175
176  Bool throwIt = False;
177  Instrument inst = STAttr::convertInstrument(header_->antennaname, throwIt);
178
179  if (inst==ATMOPRA || inst==TIDBINBILLA) {
180    header_->fluxunit = "K";
181  } else {
182    // downcase for use with Quanta
183    if (header_->fluxunit == "JY") {
184      header_->fluxunit = "Jy";
185    }
186  }
187  STAttr stattr;
188  header_->poltype = stattr.feedPolType(inst);
189  header_->nif = nIF_;
190  header_->epoch = "UTC";
191  // *** header_->frequnit = "Hz"
192  // Apply selection criteria.
193  Vector<Int> ref;
194  ifOffset_ = 0;
195  if (whichIF>=0) {
196    if (whichIF>=0 && whichIF<nIF_) {
197      ifs = False;
198      ifs(whichIF) = True;
199      header_->nif = 1;
200      nIF_ = 1;
201      ifOffset_ = whichIF;
202    } else {
203      delete reader_;
204      reader_ = 0;
205      delete header_;
206      header_ = 0;
207      throw(AipsError("Illegal IF selection"));
208    }
209  }
210  beamOffset_ = 0;
211  if (whichBeam>=0) {
212    if (whichBeam>=0 && whichBeam<nBeam_) {
213      beams = False;
214      beams(whichBeam) = True;
215      header_->nbeam = 1;
216      nBeam_ = 1;
217      beamOffset_ = whichBeam;
218    } else {
219      delete reader_;
220      reader_ = 0;
221      delete header_;
222      header_ = 0;
223      throw(AipsError("Illegal Beam selection"));
224    }
225  }
226  Vector<Int> start(nIF_, 1);
227  Vector<Int> end(nIF_, 0);
228  reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False, getPt);
229  table_->setHeader(*header_);
230  //For MS, add the location of POINTING of the input MS so one get
231  //pointing data from there, if necessary.
232  //Also find nrow in MS
233  nInDataRow = 0;
234  if (format == "MS2") {
235    Path datapath(inName);
236    String ptTabPath = datapath.absoluteName();
237    Table inMS(ptTabPath);
238    nInDataRow = inMS.nrow();
239    ptTabPath.append("/POINTING");
240    table_->table().rwKeywordSet().define("POINTING", ptTabPath);
241    if ((header_->antennaname).matches("GBT")) {
242      String GOTabPath = datapath.absoluteName();
243      GOTabPath.append("/GBT_GO");
244      table_->table().rwKeywordSet().define("GBT_GO", GOTabPath);
245    }
246  }
247  String freqFrame = header_->freqref;
248  //translate frequency reference frame back to
249  //MS style (as PKSMS2reader converts the original frame
250  //in FITS standard style)
251  if (freqFrame == "TOPOCENT") {
252    freqFrame = "TOPO";
253  } else if (freqFrame == "GEOCENER") {
254    freqFrame = "GEO";
255  } else if (freqFrame == "BARYCENT") {
256    freqFrame = "BARY";
257  } else if (freqFrame == "GALACTOC") {
258    freqFrame = "GALACTO";
259  } else if (freqFrame == "LOCALGRP") {
260    freqFrame = "LGROUP";
261  } else if (freqFrame == "CMBDIPOL") {
262    freqFrame = "CMB";
263  } else if (freqFrame == "SOURCE") {
264    freqFrame = "REST";
265  }
266  // set both "FRAME" and "BASEFRAME"
267  table_->frequencies().setFrame(freqFrame, false);
268  table_->frequencies().setFrame(freqFrame,true);
269  //table_->focus().setParallactify(true);
270}
271
272void STFiller::close( )
273{
274  delete reader_;reader_=0;
275  delete nreader_;nreader_=0;
276  delete header_;header_=0;
277  table_ = 0;
278}
279
280int asap::STFiller::read( )
281{
282  int status = 0;
283
284  //
285  // for NRO data
286  //
287  // 2008/11/12 Takeshi Nakazato
288  if ( isNRO_ ) {
289    status = readNRO() ;
290    return status ;
291  }
292  //
293
294/**
295  Int    beamNo, IFno, refBeam, scanNo, cycleNo;
296  Float  azimuth, elevation, focusAxi, focusRot, focusTan,
297    humidity, parAngle, pressure, temperature, windAz, windSpeed;
298  Double bandwidth, freqInc, interval, mjd, refFreq, srcVel;
299  String          fieldName, srcName, tcalTime, obsType;
300  Vector<Float>   calFctr, sigma, tcal, tsys;
301  Matrix<Float>   baseLin, baseSub;
302  Vector<Double>  direction(2), scanRate(2), srcDir(2), srcPM(2), restFreq(1);
303  Matrix<Float>   spectra;
304  Matrix<uChar>   flagtra;
305  Complex         xCalFctr;
306  Vector<Complex> xPol;
307**/
308
309#ifdef HAS_ALMA
310  Double min = 0.0;
311  Double max = nInDataRow;
312  ProgressMeter fillpm(min, max, "Data importing progress");
313#endif
314  PKSrecord pksrec;
315  pksrec.srcType=-1;
316  int n = 0;
317  bool isGBTFITS = false ;
318  if ((header_->antennaname.find( "GBT" ) != String::npos) && File(filename_).isRegular()) {
319    FILE *fp = fopen( filename_.c_str(), "r" ) ;
320    fseek( fp, 640, SEEK_SET ) ;
321    char buf[81] ;
322    size_t tmpret = fread( buf, 80, 1, fp ) ;
323    (void) tmpret;
324    buf[80] = '\0' ;
325    if ( strstr( buf, "NRAO_GBT" ) != NULL ) {
326      isGBTFITS = true ;
327    }
328    fclose( fp ) ;
329  }
330  while ( status == 0 ) {
331    status = reader_->read(pksrec);
332    if ( status != 0 ) break;
333    n += 1;
334    Regex filterrx(".*[SL|PA]$");
335    Regex obsrx("^AT.+");
336    if ( header_->antennaname.matches(obsrx) &&
337         pksrec.obsType.matches(filterrx)) {
338        //cerr << "ignoring paddle scan" << endl;
339        continue;
340    }
341    TableRow row(table_->table());
342    TableRecord& rec = row.record();
343    // fields that don't get used and are just passed through asap
344    RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
345    // MRC changed type from double to float
346    Vector<Double> sratedbl(pksrec.scanRate.nelements());
347    convertArray(sratedbl, pksrec.scanRate);
348    *srateCol = sratedbl;
349    RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION");
350    *spmCol = pksrec.srcPM;
351    RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION");
352    *sdirCol = pksrec.srcDir;
353    RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY");
354    *svelCol = pksrec.srcVel;
355    // the real stuff
356    RecordFieldPtr<Int> fitCol(rec, "FIT_ID");
357    *fitCol = -1;
358    RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
359    *scanoCol = pksrec.scanNo-1;
360    RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
361    *cyclenoCol = pksrec.cycleNo-1;
362    RecordFieldPtr<Double> mjdCol(rec, "TIME");
363    *mjdCol = pksrec.mjd;
364    RecordFieldPtr<Double> intCol(rec, "INTERVAL");
365    *intCol = pksrec.interval;
366    RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
367    RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
368    RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
369    *fieldnCol = pksrec.fieldName;
370    // try to auto-identify if it is on or off.
371    Regex rx(refRx_);
372    Regex rx2("_S$");
373    Int match = pksrec.srcName.matches(rx);
374    if (match) {
375      *srcnCol = pksrec.srcName;
376    } else {
377      *srcnCol = pksrec.srcName.before(rx2);
378    }
379    //*srcnCol = pksrec.srcName;//.before(rx2);
380    *srctCol = match;
381    if ( pksrec.srcType != -1 ) {
382      *srctCol = pksrec.srcType ;
383    }
384    RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
385    *beamCol = pksrec.beamNo-beamOffset_-1;
386    RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
387    Int rb = -1;
388    if (nBeam_ > 1 ) rb = pksrec.refBeam-1;
389    *rbCol = rb;
390    RecordFieldPtr<uInt> ifCol(rec, "IFNO");
391    *ifCol = pksrec.IFno-ifOffset_- 1;
392    uInt id = table_->frequencies().addEntry(Double(pksrec.spectra.nrow()/2),
393                                             pksrec.refFreq, pksrec.freqInc);
394    RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
395    *mfreqidCol = id;
396    //*ifCol = id;
397
398    id = table_->molecules().addEntry(pksrec.restFreq);
399    RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
400    *molidCol = id;
401
402    id = table_->tcal().addEntry(pksrec.tcalTime, pksrec.tcal);
403    RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
404    *mcalidCol = id;
405    id = table_->weather().addEntry(pksrec.temperature, pksrec.pressure,
406                                    pksrec.humidity, pksrec.windSpeed,
407                                    pksrec.windAz);
408    RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
409    *mweatheridCol = id;
410
411    RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
412    id = table_->focus().addEntry(pksrec.parAngle, pksrec.focusAxi,
413                                  pksrec.focusTan, pksrec.focusRot);
414    *mfocusidCol = id;
415    RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
416    *dirCol = pksrec.direction;
417    RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
418    *azCol = pksrec.azimuth;
419    RecordFieldPtr<Float> elCol(rec, "ELEVATION");
420    *elCol = pksrec.elevation;
421
422    RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
423    RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
424    RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
425
426    RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
427    // Turn the (nchan,npol) matrix and possible complex xPol vector
428    // into 2-4 rows in the scantable
429    Vector<Float> tsysvec(1);
430    // Why is spectra.ncolumn() == 3 for haveXPol_ == True
431    uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2);
432    for ( uInt i=0; i< npol; ++i ) {
433      tsysvec = pksrec.tsys(i);
434      *tsysCol = tsysvec;
435      if (isGBTFITS)
436        *polnoCol = pksrec.polNo ;
437      else
438        *polnoCol = i;
439
440      *specCol = pksrec.spectra.column(i);
441      *flagCol = pksrec.flagged.column(i);
442      table_->table().addRow();
443      row.put(table_->table().nrow()-1, rec);
444    }
445
446    RecordFieldPtr< uInt > flagrowCol(rec, "FLAGROW");
447    *flagrowCol = pksrec.flagrow;
448
449    if ( haveXPol_[0] ) {
450      // no tsys given for xpol, so emulate it
451      tsysvec = sqrt(pksrec.tsys[0]*pksrec.tsys[1]);
452      *tsysCol = tsysvec;
453      // add real part of cross pol
454      *polnoCol = 2;
455      Vector<Float> r(real(pksrec.xPol));
456      *specCol = r;
457      // make up flags from linears
458      /// @fixme this has to be a bitwise or of both pols
459      *flagCol = pksrec.flagged.column(0);// | pksrec.flagged.column(1);
460      table_->table().addRow();
461      row.put(table_->table().nrow()-1, rec);
462      // ad imaginary part of cross pol
463      *polnoCol = 3;
464      Vector<Float> im(imag(pksrec.xPol));
465      *specCol = im;
466      table_->table().addRow();
467      row.put(table_->table().nrow()-1, rec);
468    }
469#ifdef HAS_ALMA
470    fillpm._update(n);
471#endif
472  }
473  if (status > 0) {
474    close();
475    throw(AipsError("Reading error occured, data possibly corrupted."));
476  }
477#ifdef HAS_ALMA
478  fillpm.done();
479#endif
480  return status;
481}
482
483/**
484 * For NRO data
485 *
486 * 2008/11/11 Takeshi Nakazato
487 **/
488void STFiller::openNRO( int whichIF, int whichBeam )
489{
490  // open file
491  // DEBUG
492  time_t t0 ;
493  time( &t0 ) ;
494  tm *ttm = localtime( &t0 ) ;
495  LogIO os( LogOrigin( "STFiller", "openNRO()", WHERE ) ) ;
496//   cout << "STFiller::openNRO()  Start time = " << t0
497//        << " ("
498//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
499//        << " "
500//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
501//        << ")" << endl ;
502  os << "Start time = " << t0
503     << " ("
504     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
505     << " "
506     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
507     << ")" << LogIO::POST ;
508
509  // fill STHeader
510  header_ = new STHeader() ;
511
512  if ( nreader_->getHeaderInfo( header_->nchan,
513                                header_->npol,
514                                nIF_,
515                                nBeam_,
516                                header_->observer,
517                                header_->project,
518                                header_->obstype,
519                                header_->antennaname,
520                                header_->antennaposition,
521                                header_->equinox,
522                                header_->freqref,
523                                header_->reffreq,
524                                header_->bandwidth,
525                                header_->utc,
526                                header_->fluxunit,
527                                header_->epoch,
528                                header_->poltype ) ) {
529//     cout << "STFiller::openNRO()  Failed to get header information." << endl ;
530//     return ;
531    throw( AipsError("Failed to get header information.") ) ;
532  }
533
534  // set FRAME and BASEFRAME keyword of FREQUENCIES table
535  if ( header_->freqref != "TOPO" ) {
536    table_->frequencies().setFrame( header_->freqref, false ) ;
537    table_->frequencies().setFrame( header_->freqref, true ) ;
538  }
539
540  ifOffset_ = 0;
541  vector<Bool> ifs = nreader_->getIFs() ;
542  if ( whichIF >= 0 ) {
543    if ( whichIF >= 0 && whichIF < nIF_ ) {
544      for ( int i = 0 ; i < nIF_ ; i++ )
545        ifs[i] = False ;
546      ifs[whichIF] = True ;
547      header_->nif = 1;
548      nIF_ = 1;
549      ifOffset_ = whichIF;
550    } else {
551      delete reader_;
552      reader_ = 0;
553      delete header_;
554      header_ = 0;
555      throw(AipsError("Illegal IF selection"));
556    }
557  }
558
559  beamOffset_ = 0;
560  vector<Bool> beams = nreader_->getBeams() ;
561  if (whichBeam>=0) {
562    if (whichBeam>=0 && whichBeam<nBeam_) {
563      for ( int i = 0 ; i < nBeam_ ; i++ )
564        beams[i] = False ;
565      beams[whichBeam] = True;
566      header_->nbeam = 1;
567      nBeam_ = 1;
568      beamOffset_ = whichBeam;
569    } else {
570      delete reader_;
571      reader_ = 0;
572      delete header_;
573      header_ = 0;
574      throw(AipsError("Illegal Beam selection"));
575    }
576  }
577
578  header_->nbeam = nBeam_ ;
579  header_->nif = nIF_ ;
580
581  // set header
582  table_->setHeader( *header_ ) ;
583
584  // DEBUG
585  time_t t1 ;
586  time( &t1 ) ;
587  ttm = localtime( &t1 ) ;
588//   cout << "STFiller::openNRO()  End time = " << t1
589//        << " ("
590//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
591//        << " "
592//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
593//        << ")" << endl ;
594//   cout << "STFiller::openNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
595  os << "End time = " << t1
596     << " ("
597     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
598     << " "
599     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
600     << ")" << endl ;
601  os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
602  os.post() ;
603  //
604
605  return ;
606}
607
608int STFiller::readNRO()
609{
610  // DEBUG
611  time_t t0 ;
612  time( &t0 ) ;
613  tm *ttm = localtime( &t0 ) ;
614  LogIO os( LogOrigin( "STFiller", "readNRO()", WHERE ) ) ;
615//   cout << "STFiller::readNRO()  Start time = " << t0
616//        << " ("
617//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
618//        << " "
619//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
620//        << ")" << endl ;
621  os << "Start time = " << t0
622     << " ("
623     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
624     << " "
625     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
626     << ")" << LogIO::POST ;
627  //
628
629  // fill row
630  uInt id ;
631  uInt imax = nreader_->getRowNum() ;
632  vector< vector<double > > freqs ;
633  uInt i = 0 ;
634  int count = 0 ;
635  uInt scanno ;
636  uInt cycleno ;
637  uInt iftmp ;
638  uInt beamno ;
639  uInt polno ;
640  vector<double> fqs ;
641  Vector<Double> restfreq ;
642  uInt refbeamno ;
643  Double scantime ;
644  Double interval ;
645  String srcname ;
646  String fieldname ;
647  Vector<Float> spectra ;
648  Vector<uChar> flagtra ;
649  Vector<Float> tsys ;
650  Vector<Double> direction ;
651  Float azimuth ;
652  Float elevation ;
653  Float parangle ;
654  Float opacity ;
655  uInt tcalid ;
656  Int fitid ;
657  uInt focusid ;
658  Float temperature ;
659  Float pressure ;
660  Float humidity ;
661  Float windvel ;
662  Float winddir ;
663  Double srcvel ;
664  Vector<Double> propermotion( 2, 0.0 ) ;
665  Vector<Double> srcdir ;
666  Vector<Double> scanrate( 2, 0.0 ) ;
667  for ( i = 0 ; i < imax ; i++ ) {
668    string scanType = nreader_->getScanType( i ) ;
669    Int srcType = -1 ;
670    if ( scanType.compare( 0, 2, "ON") == 0 ) {
671      // os << "ON srcType: " << i << LogIO::POST ;
672      srcType = 0 ;
673    }
674    else if ( scanType.compare( 0, 3, "OFF" ) == 0 ) {
675      //os << "OFF srcType: " << i << LogIO::POST ;
676      srcType = 1 ;
677    }
678    else if ( scanType.compare( 0, 4, "ZERO" ) == 0 ) {
679      //os << "ZERO srcType: " << i << LogIO::POST ;
680      srcType = 2 ;
681    }
682    else {
683      //os << "Undefined srcType: " << i << LogIO::POST ;
684      srcType = 3 ;
685    }
686
687    // if srcType is 2 (ZERO scan), ignore scan
688    if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
689      TableRow row( table_->table() ) ;
690      TableRecord& rec = row.record();
691
692      if ( nreader_->getScanInfo( i,
693                                  scanno,
694                                  cycleno,
695                                  iftmp,
696                                  beamno,
697                                  polno,
698                                  fqs,
699                                  restfreq,
700                                  refbeamno,
701                                  scantime,
702                                  interval,
703                                  srcname,
704                                  fieldname,
705                                  spectra,
706                                  flagtra,
707                                  tsys,
708                                  direction,
709                                  azimuth,
710                                  elevation,
711                                  parangle,
712                                  opacity,
713                                  tcalid,
714                                  fitid,
715                                  focusid,
716                                  temperature,
717                                  pressure,
718                                  humidity,
719                                  windvel,
720                                  winddir,
721                                  srcvel,
722                                  propermotion,
723                                  srcdir,
724                                  scanrate ) ) {
725//         cerr << "STFiller::readNRO()  Failed to get scan information." << endl ;
726//         return 1 ;
727        throw( AipsError("Failed to get scan information.") ) ;
728      }
729
730      RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
731      *scannoCol = scanno ;
732      RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO") ;
733      *cyclenoCol = cycleno ;
734      RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
735      *beamCol = beamno ;
736      RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
737      RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
738      *polnoCol = polno ;
739      RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
740      if ( freqs.size() == 0 ) {
741        id = table_->frequencies().addEntry( Double( fqs[0] ),
742                                             Double( fqs[1] ),
743                                             Double( fqs[2] ) ) ;
744        *mfreqidCol = id ;
745        *ifCol = id ;
746        freqs.push_back( fqs ) ;
747      }
748      else {
749        int iadd = -1 ;
750        for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) {
751          //os << "freqs[" << iif << "][1] = " << freqs[iif][1] << LogIO::POST ;
752          double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ;
753          //os << "fdiff = " << fdiff << LogIO::POST ;
754          if ( fdiff < 1.0e-8 ) {
755            iadd = iif ;
756            break ;
757          }
758        }
759        if ( iadd == -1 ) {
760          id = table_->frequencies().addEntry( Double( fqs[0] ),
761                                               Double( fqs[1] ),
762                                               Double( fqs[2] ) ) ;
763          *mfreqidCol = id ;
764          *ifCol = id ;
765          freqs.push_back( fqs ) ;
766        }
767        else {
768          *mfreqidCol = iadd ;
769          *ifCol = iadd ;
770        }
771      }
772      RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
773      id = table_->molecules().addEntry( restfreq ) ;
774      *molidCol = id ;
775      RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
776      *rbCol = refbeamno ;
777      RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
778      *mjdCol = scantime ;
779      RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
780      *intervalCol = interval ;
781      RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
782      *srcnCol = srcname ;
783      RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
784      *srctCol = srcType ;
785      RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
786      *fieldnCol = fieldname ;
787      RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
788      *specCol = spectra ;
789      RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
790      *flagCol = flagtra ;
791      RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
792      *tsysCol = tsys ;
793      RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
794      *dirCol = direction ;
795      RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
796      *azCol = azimuth ;
797      RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
798      *elCol = elevation ;
799      //RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
800      //*parCol = parangle ;
801      RecordFieldPtr<Float> tauCol(rec, "OPACITY") ;
802      *tauCol = opacity ;
803      id = table_->tcal().addEntry("",Vector<Float>(1,0.0));
804      RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
805      *mcalidCol = id;
806      //RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
807      //*mcalidCol = tcalid ;
808      RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
809      *fitCol = fitid ;
810      //RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
811      //*mfocusidCol = focusid ;
812      RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
813      id = table_->focus().addEntry(parangle, 0.0, 0.0, 0.0 ) ;
814      *mfocusidCol = id;
815      RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
816      id = table_->weather().addEntry( temperature,
817                                       pressure,
818                                       humidity,
819                                       windvel,
820                                       winddir ) ;
821      *mweatheridCol = id ;
822      RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
823      *svelCol = srcvel ;
824      RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
825      *spmCol = propermotion ;
826      RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
827      *sdirCol = srcdir ;
828      RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
829      *srateCol = scanrate ;
830
831      table_->table().addRow() ;
832      row.put(table_->table().nrow()-1, rec) ;
833    }
834    else {
835      count++ ;
836    }
837    // DEBUG
838    //int rownum = nreader_->getRowNum() ;
839    //os << "Finished row " << i << "/" << rownum << LogIO::POST ;
840    //
841  }
842
843  // DEBUG
844  time_t t1 ;
845  time( &t1 ) ;
846  ttm = localtime( &t1 ) ;
847//   cout << "STFiller::readNRO()  Processed " << i << " rows" << endl ;
848//   cout << "STFiller::readNRO()  Added " << i - count << " rows (ignored "
849//        << count << " \"ZERO\" scans)" << endl ;
850//   cout << "STFiller::readNRO()  End time = " << t1
851//        << " ("
852//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
853//        << " "
854//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
855//        << ")" << endl ;
856//   cout << "STFiller::readNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
857  os << "Processed " << i << " rows" << endl ;
858  os << "Added " << i - count << " rows (ignored "
859     << count << " \"ZERO\" scans)" << endl ;
860  os.post() ;
861  os << "End time = " << t1
862     << " ("
863     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
864     << " "
865     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
866     << ")" << endl ;
867  os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
868  os.post() ;
869  //
870
871  return 0 ;
872}
873
874Bool STFiller::fileCheck()
875{
876  bool bval = false ;
877
878  // if filename_ is directory, return false
879  File inFile( filename_ ) ;
880  if ( inFile.isDirectory() )
881    return bval ;
882
883  // if beginning of header data is "RW", return true
884  // otherwise, return false ;
885  FILE *fp = fopen( filename_.c_str(), "r" ) ;
886  char buf[9] ;
887  char buf2[80] ;
888  size_t tmpret;
889  tmpret = fread( buf, 4, 1, fp ) ;
890  buf[4] = '\0' ;
891  fseek( fp, 640, SEEK_SET ) ;
892  tmpret = fread( buf2, 80, 1, fp ) ;
893  (void) tmpret; //suppress unused warning
894  if ( ( strncmp( buf, "RW", 2 ) == 0 ) || ( strstr( buf2, "NRO45M" ) != NULL ) ) {
895    bval = true ;
896  }
897  fclose( fp ) ;
898  return bval ;
899}
900
901}//namespace asap
Note: See TracBrowser for help on using the repository browser.