source: branches/alma/src/STFiller.cpp @ 1654

Last change on this file since 1654 was 1654, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-1683

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: read/write GBT SDFITS data

Put in Release Notes: No

Module(s): atnf

Description: Describe your changes here...

Changes concurrent with modifications in atnf PKSIO modules, which
is several bug fixes for reading/writing GBT SDFITS data.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.3 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, 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, 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  table_->frequencies().setFrame(freqFrame);
267     
268}
269
270void STFiller::close( )
271{
272  delete reader_;reader_=0;
273  delete nreader_;nreader_=0;
274  delete header_;header_=0;
275  table_ = 0;
276}
277
278int asap::STFiller::read( )
279{
280  int status = 0;
281
282  //
283  // for NRO data
284  //
285  // 2008/11/12 Takeshi Nakazato
286  if ( isNRO_ ) {
287    status = readNRO() ;
288    return status ;
289  }
290  //
291
292/**
293  Int    beamNo, IFno, refBeam, scanNo, cycleNo;
294  Float  azimuth, elevation, focusAxi, focusRot, focusTan,
295    humidity, parAngle, pressure, temperature, windAz, windSpeed;
296  Double bandwidth, freqInc, interval, mjd, refFreq, srcVel;
297  String          fieldName, srcName, tcalTime, obsType;
298  Vector<Float>   calFctr, sigma, tcal, tsys;
299  Matrix<Float>   baseLin, baseSub;
300  Vector<Double>  direction(2), scanRate(2), srcDir(2), srcPM(2), restFreq(1);
301  Matrix<Float>   spectra;
302  Matrix<uChar>   flagtra;
303  Complex         xCalFctr;
304  Vector<Complex> xPol;
305**/
306
307  Double min = 0.0;
308  Double max = nInDataRow;
309#ifdef HAS_ALMA
310  ProgressMeter fillpm(min, max, "Data importing progress");
311#endif
312  PKSrecord pksrec;
313  int n = 0;
314  bool isGBTFITS = false ;
315  if ((header_->antennaname.find( "GBT" ) != String::npos) && File(filename_).isRegular()) {
316    FILE *fp = fopen( filename_.c_str(), "r" ) ;
317    fseek( fp, 640, SEEK_SET ) ;
318    char buf[81] ;
319    fread( buf, 80, 1, fp ) ;
320    buf[80] = '\0' ;
321    if ( strstr( buf, "NRAO_GBT" ) != NULL ) {
322      isGBTFITS = true ;
323    }
324    fclose( fp ) ;
325  }
326  while ( status == 0 ) {
327    status = reader_->read(pksrec);
328    if ( status != 0 ) break;
329    n += 1;
330
331    Regex filterrx(".*[SL|PA]$");
332    Regex obsrx("^AT.+");
333    if ( header_->antennaname.matches(obsrx) &&
334         pksrec.obsType.matches(filterrx)) {
335        //cerr << "ignoring paddle scan" << endl;
336        continue;
337    }
338    TableRow row(table_->table());
339    TableRecord& rec = row.record();
340    // fields that don't get used and are just passed through asap
341    RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
342    // MRC changed type from double to float
343    Vector<Double> sratedbl(pksrec.scanRate.nelements());
344    convertArray(sratedbl, pksrec.scanRate);
345    *srateCol = sratedbl;
346    RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION");
347    *spmCol = pksrec.srcPM;
348    RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION");
349    *sdirCol = pksrec.srcDir;
350    RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY");
351    *svelCol = pksrec.srcVel;
352    // the real stuff
353    RecordFieldPtr<Int> fitCol(rec, "FIT_ID");
354    *fitCol = -1;
355    RecordFieldPtr<uInt> scanoCol(rec, "SCANNO");
356    *scanoCol = pksrec.scanNo-1;
357    RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
358    *cyclenoCol = pksrec.cycleNo-1;
359    RecordFieldPtr<Double> mjdCol(rec, "TIME");
360    *mjdCol = pksrec.mjd;
361    RecordFieldPtr<Double> intCol(rec, "INTERVAL");
362    *intCol = pksrec.interval;
363    RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
364    RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
365    RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
366    *fieldnCol = pksrec.fieldName;
367    // try to auto-identify if it is on or off.
368    Regex rx(refRx_);
369    Regex rx2("_S$");
370    Int match = pksrec.srcName.matches(rx);
371    if (match) {
372      *srcnCol = pksrec.srcName;
373    } else {
374      *srcnCol = pksrec.srcName.before(rx2);
375    }
376    //*srcnCol = pksrec.srcName;//.before(rx2);
377    *srctCol = match;
378    RecordFieldPtr<uInt> beamCol(rec, "BEAMNO");
379    *beamCol = pksrec.beamNo-beamOffset_-1;
380    RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO");
381    Int rb = -1;
382    if (nBeam_ > 1 ) rb = pksrec.refBeam-1;
383    *rbCol = rb;
384    RecordFieldPtr<uInt> ifCol(rec, "IFNO");
385    *ifCol = pksrec.IFno-ifOffset_- 1;
386    uInt id;
387    /// @todo this has to change when nchan isn't global anymore
388    id = table_->frequencies().addEntry(Double(header_->nchan/2),
389                                            pksrec.refFreq, pksrec.freqInc);
390    RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID");
391    *mfreqidCol = id;
392
393    id = table_->molecules().addEntry(pksrec.restFreq);
394    RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID");
395    *molidCol = id;
396
397    id = table_->tcal().addEntry(pksrec.tcalTime, pksrec.tcal);
398    RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID");
399    *mcalidCol = id;
400    id = table_->weather().addEntry(pksrec.temperature, pksrec.pressure,
401                                    pksrec.humidity, pksrec.windSpeed,
402                                    pksrec.windAz);
403    RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID");
404    *mweatheridCol = id;
405    RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID");
406    id = table_->focus().addEntry(pksrec.focusAxi, pksrec.focusTan,
407                                  pksrec.focusRot);
408    *mfocusidCol = id;
409    RecordFieldPtr<Array<Double> > dirCol(rec, "DIRECTION");
410    *dirCol = pksrec.direction;
411    RecordFieldPtr<Float> azCol(rec, "AZIMUTH");
412    *azCol = pksrec.azimuth;
413    RecordFieldPtr<Float> elCol(rec, "ELEVATION");
414    *elCol = pksrec.elevation;
415
416    RecordFieldPtr<Float> parCol(rec, "PARANGLE");
417    *parCol = pksrec.parAngle;
418
419    RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA");
420    RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA");
421    RecordFieldPtr< uInt > polnoCol(rec, "POLNO");
422
423    RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS");
424    // Turn the (nchan,npol) matrix and possible complex xPol vector
425    // into 2-4 rows in the scantable
426    Vector<Float> tsysvec(1);
427    // Why is spectra.ncolumn() == 3 for haveXPol_ == True
428    uInt npol = (pksrec.spectra.ncolumn()==1 ? 1: 2);
429    for ( uInt i=0; i< npol; ++i ) {
430      tsysvec = pksrec.tsys(i);
431      *tsysCol = tsysvec;
432      if (isGBTFITS)
433        *polnoCol = pksrec.polNo ;
434      else
435        *polnoCol = i;
436
437      *specCol = pksrec.spectra.column(i);
438      *flagCol = pksrec.flagged.column(i);
439      table_->table().addRow();
440      row.put(table_->table().nrow()-1, rec);
441    }
442    if ( haveXPol_[0] ) {
443      // no tsys given for xpol, so emulate it
444      tsysvec = sqrt(pksrec.tsys[0]*pksrec.tsys[1]);
445      *tsysCol = tsysvec;
446      // add real part of cross pol
447      *polnoCol = 2;
448      Vector<Float> r(real(pksrec.xPol));
449      *specCol = r;
450      // make up flags from linears
451      /// @fixme this has to be a bitwise or of both pols
452      *flagCol = pksrec.flagged.column(0);// | pksrec.flagged.column(1);
453      table_->table().addRow();
454      row.put(table_->table().nrow()-1, rec);
455      // ad imaginary part of cross pol
456      *polnoCol = 3;
457      Vector<Float> im(imag(pksrec.xPol));
458      *specCol = im;
459      table_->table().addRow();
460      row.put(table_->table().nrow()-1, rec);
461    }
462#ifdef HAS_ALMA
463    fillpm._update(n);
464#endif
465  }
466  if (status > 0) {
467    close();
468    throw(AipsError("Reading error occured, data possibly corrupted."));
469  }
470#ifdef HAS_ALMA
471  fillpm.done();
472#endif
473  return status;
474}
475
476/**
477 * For NRO data
478 *
479 * 2008/11/11 Takeshi Nakazato
480 **/
481void STFiller::openNRO( int whichIF, int whichBeam )
482{
483  // open file
484  // DEBUG
485  time_t t0 ;
486  time( &t0 ) ;
487  tm *ttm = localtime( &t0 ) ;
488  LogIO os( LogOrigin( "STFiller", "openNRO()", WHERE ) ) ;
489//   cout << "STFiller::openNRO()  Start time = " << t0
490//        << " ("
491//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
492//        << " "
493//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
494//        << ")" << endl ;
495  os << "Start time = " << t0
496     << " ("
497     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
498     << " "
499     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
500     << ")" << LogIO::POST ;
501
502  // fill STHeader
503  header_ = new STHeader() ;
504
505  if ( nreader_->getHeaderInfo( header_->nchan,
506                                header_->npol,
507                                nIF_,
508                                nBeam_,
509                                header_->observer,
510                                header_->project,
511                                header_->obstype,
512                                header_->antennaname,
513                                header_->antennaposition,
514                                header_->equinox,
515                                header_->freqref,
516                                header_->reffreq,
517                                header_->bandwidth,
518                                header_->utc,
519                                header_->fluxunit,
520                                header_->epoch,
521                                header_->poltype ) ) {
522//     cout << "STFiller::openNRO()  Failed to get header information." << endl ;
523//     return ;
524    throw( AipsError("Failed to get header information.") ) ;
525  }
526
527  // set frame keyword of FREQUENCIES table
528  if ( header_->freqref != "TOPO" ) {
529    table_->frequencies().setFrame( header_->freqref, false ) ;
530  }
531
532  ifOffset_ = 0;
533  vector<Bool> ifs = nreader_->getIFs() ;
534  if ( whichIF >= 0 ) {
535    if ( whichIF >= 0 && whichIF < nIF_ ) {
536      for ( int i = 0 ; i < nIF_ ; i++ )
537        ifs[i] = False ;
538      ifs[whichIF] = True ;
539      header_->nif = 1;
540      nIF_ = 1;
541      ifOffset_ = whichIF;
542    } else {
543      delete reader_;
544      reader_ = 0;
545      delete header_;
546      header_ = 0;
547      throw(AipsError("Illegal IF selection"));
548    }
549  }
550
551  beamOffset_ = 0;
552  vector<Bool> beams = nreader_->getBeams() ;
553  if (whichBeam>=0) {
554    if (whichBeam>=0 && whichBeam<nBeam_) {
555      for ( int i = 0 ; i < nBeam_ ; i++ )
556        beams[i] = False ;
557      beams[whichBeam] = True;
558      header_->nbeam = 1;
559      nBeam_ = 1;
560      beamOffset_ = whichBeam;
561    } else {
562      delete reader_;
563      reader_ = 0;
564      delete header_;
565      header_ = 0;
566      throw(AipsError("Illegal Beam selection"));
567    }
568  }
569
570  header_->nbeam = nBeam_ ;
571  header_->nif = nIF_ ;
572
573  // set header
574  table_->setHeader( *header_ ) ;
575
576  // DEBUG
577  time_t t1 ;
578  time( &t1 ) ;
579  ttm = localtime( &t1 ) ;
580//   cout << "STFiller::openNRO()  End time = " << t1
581//        << " ("
582//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
583//        << " "
584//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
585//        << ")" << endl ;
586//   cout << "STFiller::openNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
587  os << "End time = " << t1
588     << " ("
589     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
590     << " "
591     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
592     << ")" << endl ;
593  os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
594  os.post() ;
595  //
596
597  return ;
598}
599
600int STFiller::readNRO()
601{
602  // DEBUG
603  time_t t0 ;
604  time( &t0 ) ;
605  tm *ttm = localtime( &t0 ) ;
606  LogIO os( LogOrigin( "STFiller", "readNRO()", WHERE ) ) ;
607//   cout << "STFiller::readNRO()  Start time = " << t0
608//        << " ("
609//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
610//        << " "
611//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
612//        << ")" << endl ;
613  os << "Start time = " << t0
614     << " ("
615     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
616     << " "
617     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
618     << ")" << LogIO::POST ;
619  //
620
621  // fill row
622  uInt id ;
623  uInt imax = nreader_->getRowNum() ;
624  vector< vector<double > > freqs ;
625  uInt i = 0 ;
626  int count = 0 ;
627  uInt scanno ;
628  uInt cycleno ;
629  uInt beamno ;
630  uInt polno ;
631  vector<double> fqs ;
632  Vector<Double> restfreq ;
633  uInt refbeamno ;
634  Double scantime ;
635  Double interval ;
636  String srcname ;
637  String fieldname ;
638  Array<Float> spectra ;
639  Array<uChar> flagtra ;
640  Array<Float> tsys ;
641  Array<Double> direction ;
642  Float azimuth ;
643  Float elevation ;
644  Float parangle ;
645  Float opacity ;
646  uInt tcalid ;
647  Int fitid ;
648  uInt focusid ;
649  Float temperature ;
650  Float pressure ;
651  Float humidity ;
652  Float windvel ;
653  Float winddir ;
654  Double srcvel ;
655  Array<Double> propermotion ;
656  Vector<Double> srcdir ;
657  Array<Double> scanrate ;
658  for ( i = 0 ; i < imax ; i++ ) {
659    string scanType = nreader_->getScanType( i ) ;
660    Int srcType = -1 ;
661    if ( scanType.compare( 0, 2, "ON") == 0 ) {
662      // os << "ON srcType: " << i << LogIO::POST ;
663      srcType = 0 ;
664    }
665    else if ( scanType.compare( 0, 3, "OFF" ) == 0 ) {
666      //os << "OFF srcType: " << i << LogIO::POST ;
667      srcType = 1 ;
668    }
669    else if ( scanType.compare( 0, 4, "ZERO" ) == 0 ) {
670      //os << "ZERO srcType: " << i << LogIO::POST ;
671      srcType = 2 ;
672    }
673    else {
674      //os << "Undefined srcType: " << i << LogIO::POST ;
675      srcType = 3 ;
676    }
677 
678    // if srcType is 2 (ZERO scan), ignore scan
679    if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
680      TableRow row( table_->table() ) ;
681      TableRecord& rec = row.record();
682
683      if ( nreader_->getScanInfo( i,
684                                  scanno,
685                                  cycleno,
686                                  beamno,
687                                  polno,
688                                  fqs,
689                                  restfreq,
690                                  refbeamno,
691                                  scantime,
692                                  interval,
693                                  srcname,
694                                  fieldname,
695                                  spectra,
696                                  flagtra,
697                                  tsys,
698                                  direction,
699                                  azimuth,
700                                  elevation,
701                                  parangle,
702                                  opacity,
703                                  tcalid,
704                                  fitid,
705                                  focusid,
706                                  temperature,
707                                  pressure,
708                                  humidity,
709                                  windvel,
710                                  winddir,
711                                  srcvel,
712                                  propermotion,
713                                  srcdir,
714                                  scanrate ) ) {
715//         cerr << "STFiller::readNRO()  Failed to get scan information." << endl ;
716//         return 1 ;
717        throw( AipsError("Failed to get scan information.") ) ;
718      }
719
720      RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
721      *scannoCol = scanno ;
722      RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO") ;
723      *cyclenoCol = cycleno ;
724      RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
725      *beamCol = beamno ;
726      RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
727      RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
728      *polnoCol = polno ;
729      RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
730      if ( freqs.size() == 0 ) {
731        id = table_->frequencies().addEntry( Double( fqs[0] ),
732                                             Double( fqs[1] ),
733                                             Double( fqs[2] ) ) ;
734        *mfreqidCol = id ;
735        *ifCol = id ;
736        freqs.push_back( fqs ) ;
737      }
738      else {
739        int iadd = -1 ;
740        for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) {
741          //os << "freqs[" << iif << "][1] = " << freqs[iif][1] << LogIO::POST ;
742          double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ;
743          //os << "fdiff = " << fdiff << LogIO::POST ;
744          if ( fdiff < 1.0e-8 ) {
745            iadd = iif ;
746            break ;
747          }
748        }
749        if ( iadd == -1 ) {
750          id = table_->frequencies().addEntry( Double( fqs[0] ),
751                                               Double( fqs[1] ),
752                                               Double( fqs[2] ) ) ;
753          *mfreqidCol = id ;
754          *ifCol = id ;
755          freqs.push_back( fqs ) ;
756        }
757        else {
758          *mfreqidCol = iadd ;
759          *ifCol = iadd ;
760        }
761      }
762      RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
763      id = table_->molecules().addEntry( restfreq ) ;
764      *molidCol = id ;
765      RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
766      *rbCol = refbeamno ;
767      RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
768      *mjdCol = scantime ;
769      RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
770      *intervalCol = interval ;
771      RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
772      *srcnCol = srcname ;
773      RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
774      *srctCol = srcType ;     
775      RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
776      *fieldnCol = fieldname ;
777      RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
778      *specCol = spectra ;
779      RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
780      *flagCol = flagtra ;
781      RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
782      *tsysCol = tsys ;
783      RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
784      *dirCol = direction ;
785      RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
786      *azCol = azimuth ;
787      RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
788      *elCol = elevation ;
789      RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
790      *parCol = parangle ;
791      RecordFieldPtr<Float> tauCol(rec, "OPACITY") ;
792      *tauCol = opacity ;
793      RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
794      *mcalidCol = tcalid ;
795      RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
796      *fitCol = fitid ;
797      RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
798      *mfocusidCol = focusid ;
799      RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
800      id = table_->weather().addEntry( temperature,
801                                       pressure,
802                                       humidity,
803                                       windvel,
804                                       winddir ) ;
805      *mweatheridCol = id ;         
806      RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
807      *svelCol = srcvel ;
808      RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
809      *spmCol = propermotion ;
810      RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
811      *sdirCol = srcdir ;
812      RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
813      *srateCol = scanrate ;
814
815      table_->table().addRow() ;
816      row.put(table_->table().nrow()-1, rec) ;
817    }
818    else {
819      count++ ;
820    }
821    // DEBUG
822    //int rownum = nreader_->getRowNum() ;
823    //os << "Finished row " << i << "/" << rownum << LogIO::POST ;
824    //
825  }
826
827  // DEBUG
828  time_t t1 ;
829  time( &t1 ) ;
830  ttm = localtime( &t1 ) ;
831//   cout << "STFiller::readNRO()  Processed " << i << " rows" << endl ;
832//   cout << "STFiller::readNRO()  Added " << i - count << " rows (ignored "
833//        << count << " \"ZERO\" scans)" << endl ;
834//   cout << "STFiller::readNRO()  End time = " << t1
835//        << " ("
836//        << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
837//        << " "
838//        << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
839//        << ")" << endl ;
840//   cout << "STFiller::readNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
841  os << "Processed " << i << " rows" << endl ;
842  os << "Added " << i - count << " rows (ignored "
843     << count << " \"ZERO\" scans)" << endl ;
844  os.post() ;
845  os << "End time = " << t1
846     << " ("
847     << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
848     << " "
849     << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
850     << ")" << endl ;
851  os << "Elapsed time = " << t1 - t0 << " sec" << endl ;
852  os.post() ;
853  //
854
855  return 0 ;
856}
857
858Bool STFiller::fileCheck()
859{
860  bool bval = false ;
861
862  // if filename_ is directory, return false
863  File inFile( filename_ ) ;
864  if ( inFile.isDirectory() )
865    return bval ;
866 
867  // if beginning of header data is "RW", return true
868  // otherwise, return false ;
869  FILE *fp = fopen( filename_.c_str(), "r" ) ;
870  char buf[9] ;
871  char buf2[80] ;
872  fread( buf, 4, 1, fp ) ;
873  buf[4] = '\0' ;
874  fseek( fp, 640, SEEK_SET ) ;
875  fread( buf2, 80, 1, fp ) ;
876  if ( ( strncmp( buf, "RW", 2 ) == 0 ) || ( strstr( buf2, "NRO45M" ) != NULL ) ) {
877    bval = true ;
878  }
879  fclose( fp ) ;
880  return bval ;
881}
882
883}//namespace asap
Note: See TracBrowser for help on using the repository browser.