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

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

New Development: No

JIRA Issue: Yes CAS-1043

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Added lines to set FRAME attribute of FREQUENCIES

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