source: branches/mergetest/src/STFiller.cpp @ 1779

Last change on this file since 1779 was 1779, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: No (test merging alma branch)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

Description:


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