source: trunk/src/STFiller.cpp @ 2163

Last change on this file since 2163 was 2163, checked in by Malte Marquarding, 13 years ago

Remove various compiler warnings
  • 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 beamno ;
638  uInt polno ;
639  vector<double> fqs ;
640  Vector<Double> restfreq ;
641  uInt refbeamno ;
642  Double scantime ;
643  Double interval ;
644  String srcname ;
645  String fieldname ;
646  Array<Float> spectra ;
647  Array<uChar> flagtra ;
648  Array<Float> tsys ;
649  Array<Double> direction ;
650  Float azimuth ;
651  Float elevation ;
652  Float parangle ;
653  Float opacity ;
654  uInt tcalid ;
655  Int fitid ;
656  uInt focusid ;
657  Float temperature ;
658  Float pressure ;
659  Float humidity ;
660  Float windvel ;
661  Float winddir ;
662  Double srcvel ;
663  Array<Double> propermotion ;
664  Vector<Double> srcdir ;
665  Array<Double> scanrate ;
666  for ( i = 0 ; i < imax ; i++ ) {
667    string scanType = nreader_->getScanType( i ) ;
668    Int srcType = -1 ;
669    if ( scanType.compare( 0, 2, "ON") == 0 ) {
670      // os << "ON srcType: " << i << LogIO::POST ;
671      srcType = 0 ;
672    }
673    else if ( scanType.compare( 0, 3, "OFF" ) == 0 ) {
674      //os << "OFF srcType: " << i << LogIO::POST ;
675      srcType = 1 ;
676    }
677    else if ( scanType.compare( 0, 4, "ZERO" ) == 0 ) {
678      //os << "ZERO srcType: " << i << LogIO::POST ;
679      srcType = 2 ;
680    }
681    else {
682      //os << "Undefined srcType: " << i << LogIO::POST ;
683      srcType = 3 ;
684    }
685
686    // if srcType is 2 (ZERO scan), ignore scan
687    if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
688      TableRow row( table_->table() ) ;
689      TableRecord& rec = row.record();
690
691      if ( nreader_->getScanInfo( i,
692                                  scanno,
693                                  cycleno,
694                                  beamno,
695                                  polno,
696                                  fqs,
697                                  restfreq,
698                                  refbeamno,
699                                  scantime,
700                                  interval,
701                                  srcname,
702                                  fieldname,
703                                  spectra,
704                                  flagtra,
705                                  tsys,
706                                  direction,
707                                  azimuth,
708                                  elevation,
709                                  parangle,
710                                  opacity,
711                                  tcalid,
712                                  fitid,
713                                  focusid,
714                                  temperature,
715                                  pressure,
716                                  humidity,
717                                  windvel,
718                                  winddir,
719                                  srcvel,
720                                  propermotion,
721                                  srcdir,
722                                  scanrate ) ) {
723//         cerr << "STFiller::readNRO()  Failed to get scan information." << endl ;
724//         return 1 ;
725        throw( AipsError("Failed to get scan information.") ) ;
726      }
727
728      RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
729      *scannoCol = scanno ;
730      RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO") ;
731      *cyclenoCol = cycleno ;
732      RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
733      *beamCol = beamno ;
734      RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
735      RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
736      *polnoCol = polno ;
737      RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
738      if ( freqs.size() == 0 ) {
739        id = table_->frequencies().addEntry( Double( fqs[0] ),
740                                             Double( fqs[1] ),
741                                             Double( fqs[2] ) ) ;
742        *mfreqidCol = id ;
743        *ifCol = id ;
744        freqs.push_back( fqs ) ;
745      }
746      else {
747        int iadd = -1 ;
748        for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) {
749          //os << "freqs[" << iif << "][1] = " << freqs[iif][1] << LogIO::POST ;
750          double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ;
751          //os << "fdiff = " << fdiff << LogIO::POST ;
752          if ( fdiff < 1.0e-8 ) {
753            iadd = iif ;
754            break ;
755          }
756        }
757        if ( iadd == -1 ) {
758          id = table_->frequencies().addEntry( Double( fqs[0] ),
759                                               Double( fqs[1] ),
760                                               Double( fqs[2] ) ) ;
761          *mfreqidCol = id ;
762          *ifCol = id ;
763          freqs.push_back( fqs ) ;
764        }
765        else {
766          *mfreqidCol = iadd ;
767          *ifCol = iadd ;
768        }
769      }
770      RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
771      id = table_->molecules().addEntry( restfreq ) ;
772      *molidCol = id ;
773      RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
774      *rbCol = refbeamno ;
775      RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
776      *mjdCol = scantime ;
777      RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
778      *intervalCol = interval ;
779      RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
780      *srcnCol = srcname ;
781      RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
782      *srctCol = srcType ;
783      RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
784      *fieldnCol = fieldname ;
785      RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
786      *specCol = spectra ;
787      RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
788      *flagCol = flagtra ;
789      RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
790      *tsysCol = tsys ;
791      RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
792      *dirCol = direction ;
793      RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
794      *azCol = azimuth ;
795      RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
796      *elCol = elevation ;
797      //RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
798      //*parCol = parangle ;
799      RecordFieldPtr<Float> tauCol(rec, "OPACITY") ;
800      *tauCol = opacity ;
801      id = table_->tcal().addEntry("",Vector<Float>(1,0.0));
802      RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
803      *mcalidCol = id;
804      //RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
805      //*mcalidCol = tcalid ;
806      RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
807      *fitCol = fitid ;
808      //RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
809      //*mfocusidCol = focusid ;
810      RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
811      id = table_->focus().addEntry(parangle, 0.0, 0.0, 0.0 ) ;
812      *mfocusidCol = id;
813      RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
814      id = table_->weather().addEntry( temperature,
815                                       pressure,
816                                       humidity,
817                                       windvel,
818                                       winddir ) ;
819      *mweatheridCol = id ;
820      RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
821      *svelCol = srcvel ;
822      RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
823      *spmCol = propermotion ;
824      RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
825      *sdirCol = srcdir ;
826      RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
827      *srateCol = scanrate ;
828
829      table_->table().addRow() ;
830      row.put(table_->table().nrow()-1, rec) ;
831    }
832    else {
833      count++ ;
834    }
835    // DEBUG
836    //int rownum = nreader_->getRowNum() ;
837    //os << "Finished row " << i << "/" << rownum << LogIO::POST ;
838    //
839  }
840
841  // DEBUG
842  time_t t1 ;
843  time( &t1 ) ;
844  ttm = localtime( &t1 ) ;
845//   cout << "STFiller::readNRO()  Processed " << i << " rows" << endl ;
846//   cout << "STFiller::readNRO()  Added " << i - count << " rows (ignored "
847//        << count << " \"ZERO\" scans)" << endl ;
848//   cout << "STFiller::readNRO()  End time = " << t1
849//        << " ("
850//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
851//        << " "
852//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
853//        << ")" << endl ;
854//   cout << "STFiller::readNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
855  os << "Processed " << i << " rows" << endl ;
856  os << "Added " << i - count << " rows (ignored "
857     << count << " \"ZERO\" scans)" << endl ;
858  os.post() ;
859  os << "End time = " << t1
860     << " ("
861     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
862     << " "
863     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
864     << ")" << endl ;
865  os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
866  os.post() ;
867  //
868
869  return 0 ;
870}
871
872Bool STFiller::fileCheck()
873{
874  bool bval = false ;
875
876  // if filename_ is directory, return false
877  File inFile( filename_ ) ;
878  if ( inFile.isDirectory() )
879    return bval ;
880
881  // if beginning of header data is "RW", return true
882  // otherwise, return false ;
883  FILE *fp = fopen( filename_.c_str(), "r" ) ;
884  char buf[9] ;
885  char buf2[80] ;
886  size_t tmpret;
887  tmpret = fread( buf, 4, 1, fp ) ;
888  buf[4] = '\0' ;
889  fseek( fp, 640, SEEK_SET ) ;
890  tmpret = fread( buf2, 80, 1, fp ) ;
891  (void) tmpret; //suppress unused warning
892  if ( ( strncmp( buf, "RW", 2 ) == 0 ) || ( strstr( buf2, "NRO45M" ) != NULL ) ) {
893    bval = true ;
894  }
895  fclose( fp ) ;
896  return bval ;
897}
898
899}//namespace asap
Note: See TracBrowser for help on using the repository browser.