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

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

New Development: Yes

JIRA Issue: No (merge)

Ready for Test: Yes

Interface Changes: No

Description:

Merged changes -r1774:1817 in newfiller branch to alma branch


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