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

Last change on this file since 1495 was 1495, checked in by TakTsutsumi, 15 years ago

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: a new parameter in open method

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Added a boolean, getPt in open()

to accommodate change in PKSIO
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.7 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//   if ( nreader_->open() != 0 ) {
459//     throw( AipsError( "Error while opening file "+filename_ ) ) ;
460//   }
461
462  isNRO_ = true ;
463
464  // store data into  NROHeader and NRODataset
465  //if ( nreader_->readHeader() != 0 ) {
466  //throw( AipsError( "Error while reading file "+filename_ ) ) ;
467  //}
468
469  // get header data
470  NROHeader *nheader =  nreader_->getHeader() ;
471
472  // fill STHeader
473  header_ = new STHeader() ;
474
475  header_->nchan = nheader->getNUMCH() ;
476  header_->npol = nreader_->getPolarizationNum() ;
477  header_->observer = nheader->getOBSVR() ;
478  header_->project = nheader->getPROJ() ;
479  header_->obstype = nheader->getSWMOD() ;
480  header_->antennaname = nheader->getSITE() ;
481  // TODO: should be investigated antenna position since there are
482  //       no corresponding information in the header
483  // 2008/11/13 Takeshi Nakazato
484  //
485  // INFO: tentative antenna posiiton is obtained for NRO 45m from ITRF website
486  // 2008/11/26 Takeshi Nakazato
487  vector<double> pos = nreader_->getAntennaPosition() ;
488  header_->antennaposition = pos ;
489  char *eq = nheader->getEPOCH() ;
490  if ( strncmp( eq, "B1950", 5 ) == 0 )
491    header_->equinox = 1950.0 ;
492  else if ( strncmp( eq, "J2000", 5 ) == 0 )
493    header_->equinox = 2000.0 ;
494//   char *vref = nheader->getVREF() ;
495//   if ( strncmp( vref, "LSR", 3 ) == 0 ) {
496//     strcat( vref, "K" ) ;
497//   }
498//   header_->freqref = vref ;
499  nreader_->getData( 0 ) ;
500  header_->reffreq = nreader_->getData()->FREQ0 ;
501  header_->bandwidth = nheader->getBEBW()[0] ;
502  header_->utc = nreader_->getStartTime() ;
503  header_->fluxunit = "K" ;
504  header_->epoch = "UTC" ; 
505  char *poltp = nheader->getPOLTP()[0] ;
506  if ( strcmp( poltp, "" ) == 0 )
507    poltp = "None" ;
508  header_->poltype = poltp ;
509  // DEBUG
510  cout << "STFiller::openNRO()  poltype = " << header_->poltype << endl ;
511  //
512
513  vector<Bool> ifs = nreader_->getIFs() ;
514  ifOffset_ = 0;
515  nIF_ = ifs.size() ;
516  if ( whichIF >= 0 ) {
517    if ( whichIF >= 0 && whichIF < nIF_ ) {
518      for ( int i = 0 ; i < nIF_ ; i++ )
519        ifs[i] = False ;
520      ifs[whichIF] = True ;
521      header_->nif = 1;
522      nIF_ = 1;
523      ifOffset_ = whichIF;
524    } else {
525      delete reader_;
526      reader_ = 0;
527      delete header_;
528      header_ = 0;
529      throw(AipsError("Illegal IF selection"));
530    }
531  }
532
533  beamOffset_ = 0;
534  vector<Bool> beams = nreader_->getBeams() ;
535  nBeam_ = beams.size() ;
536  if (whichBeam>=0) {
537    if (whichBeam>=0 && whichBeam<nBeam_) {
538      for ( int i = 0 ; i < nBeam_ ; i++ )
539        beams[i] = False ;
540      beams[whichBeam] = True;
541      header_->nbeam = 1;
542      nBeam_ = 1;
543      beamOffset_ = whichBeam;
544    } else {
545      delete reader_;
546      reader_ = 0;
547      delete header_;
548      header_ = 0;
549      throw(AipsError("Illegal Beam selection"));
550    }
551  }
552  header_->nbeam = nBeam_ ;
553  header_->nif = nIF_ ;
554
555  // set header
556  table_->setHeader( *header_ ) ;
557
558  // DEBUG
559  //cout << "STFiller::openNRO() Velocity Definition = " << nheader->getVDEF() << endl ;
560
561  // DEBUG
562  time_t t1 ;
563  time( &t1 ) ;
564  ttm = localtime( &t1 ) ;
565  cout << "STFiller::openNRO()  End time = " << t1
566       << " ("
567       << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
568       << " "
569       << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
570       << ")" << endl ;
571  cout << "STFiller::openNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
572  //
573
574  return ;
575}
576
577int STFiller::readNRO()
578{
579  // DEBUG
580  time_t t0 ;
581  time( &t0 ) ;
582  tm *ttm = localtime( &t0 ) ;
583  cout << "STFiller::readNRO()  Start time = " << t0
584       << " ("
585       << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
586       << " "
587       << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
588       << ")" << endl ;
589  //
590
591  // get header
592  NROHeader *h = nreader_->getHeader() ;
593
594  // fill row
595  uInt id ;
596  uInt imax = nreader_->getRowNum() ;
597  vector< vector<double > > freqs ;
598  uInt i = 0 ;
599  int count = 0 ;
600  for ( i = 0 ; i < imax ; i++ ) {
601    if( nreader_->getData( i ) != 0 ) {
602      cerr << "STFiller::readNRO()  error while reading row " << i << endl ;
603      return -1 ;
604    }
605    NRODataset *d = nreader_->getData() ;
606
607    char *scanType = d->SCANTP ;
608    Int srcType = -1 ;
609    if ( strncmp( scanType, "ON", 2 ) == 0 ) {
610      srcType = 0 ;
611    }
612    else if ( strncmp( scanType, "OFF", 3 ) == 0 ) {
613      //cout << "OFF srcType: " << i << endl ;
614      srcType = 1 ;
615    }
616    else if ( strncmp( scanType, "ZERO", 4 ) == 0 ) {
617      //cout << "ZERO srcType: " << i << endl ;
618      srcType = 2 ;
619    }
620    else {
621      //cout << "Undefined srcType: " << i << endl ;
622      srcType = 3 ;
623    }
624 
625    // if srcType is 2 (ZERO scan), ignore scan
626    if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
627      TableRow row( table_->table() ) ;
628      TableRecord& rec = row.record();
629
630      RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
631      *srctCol = srcType ;     
632      RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
633      Array<Double> srcarr( IPosition( 1, 2 ) ) ;
634      srcarr = 0.0 ;
635      *srateCol = srcarr ;
636      RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
637      *spmCol = srcarr ;
638      RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
639      *sdirCol = nreader_->getSourceDirection() ;
640      RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
641      *svelCol = h->getURVEL() ;   // [m/s]
642      RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
643      *fitCol = -1 ;
644      RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
645      //*scannoCol = (uInt)(d->getISCAN()) ;
646      *scannoCol = (uInt)(d->ISCAN) ;
647      RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
648      RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
649      *mjdCol = Double( nreader_->getStartIntTime( i ) ) ;
650      RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
651      *intervalCol = Double( h->getIPTIM() ) ;
652      RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
653      *srcnCol = String( h->getOBJ() ) ;
654      RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
655      *fieldnCol = String( h->getOBJ() ) ;
656      // BEAMNO is 0-base
657      RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
658      string arryt = string( d->ARRYT ) ;
659      string sbeamno = arryt.substr( 1, arryt.size()-1 ) ;
660      uInt ibeamno = atoi( sbeamno.c_str() ) ;
661      *beamCol = ibeamno - 1 ;
662      RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
663      RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
664      RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
665      vector<double> fqs = nreader_->getFrequencies( i ) ;
666      //cout << "STFiller::readNRO()  fqs[1] = " << fqs[1] << endl ;
667      if ( freqs.size() == 0 ) {
668        id = table_->frequencies().addEntry( Double( fqs[0] ),
669                                             Double( fqs[1] ),
670                                             Double( fqs[2] ) ) ;
671        *mfreqidCol = id ;
672        *ifCol = id ;
673        freqs.push_back( fqs ) ;
674      }
675      else {
676        int iadd = -1 ;
677        for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) {
678          //cout << "STFiller::readNRO()  freqs[" << iif << "][1] = " << freqs[iif][1] << endl ;
679          double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ;
680          //cout << "STFiller::readNRO()  fdiff = " << fdiff << endl ;
681          if ( fdiff < 1.0e-8 ) {
682            iadd = iif ;
683            break ;
684          }
685        }
686        if ( iadd == -1 ) {
687          id = table_->frequencies().addEntry( Double( fqs[0] ),
688                                               Double( fqs[1] ),
689                                               Double( fqs[2] ) ) ;
690          *mfreqidCol = id ;
691          *ifCol = id ;
692          freqs.push_back( fqs ) ;
693        }
694        else {
695          *mfreqidCol = iadd ;
696          *ifCol = iadd ;
697        }
698      }
699      RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
700      Vector<Double> restfreq( IPosition( 1, 1 ) ) ;
701      restfreq( 0 ) = d->FREQ0 ;
702      id = table_->molecules().addEntry( restfreq ) ;
703      *molidCol = id ;
704      RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
705      //
706      // No Tcal information in the data
707      //
708      // 2008/11/20 Takeshi Nakazato
709      //
710      *mcalidCol = 0 ;
711      RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
712
713      id = table_->weather().addEntry( Float( d->TEMP ),
714                                       Float( d->PATM ),
715                                       Float( d->PH2O ),
716                                       Float( d->VWIND ),
717                                       Float( d->DWIND ) ) ;
718
719      *mweatheridCol = id ;         
720      RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
721      RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
722      *dirCol = nreader_->getDirection( i ) ;
723      RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
724      *azCol = d->RAZ ;
725      RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
726      *elCol = d->REL ;
727      RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
728      RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
729      vector<double> spec = nreader_->getSpectrum( i ) ;
730      Array<Float> sp( IPosition( 1, spec.size() ) ) ;
731      int index = 0 ;
732      for ( Array<Float>::iterator itr = sp.begin() ; itr != sp.end() ; itr++ ) {
733        *itr = spec[index++] ;
734      }
735      *specCol = sp ;
736      RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
737      Array<uChar> flag( sp.shape() ) ;
738      flag.set( 0 ) ;
739      *flagCol = flag ;
740      RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
741      *polnoCol = 0 ;
742      RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
743      Array<Float> tsys( IPosition( 1, 1 ), d->TSYS ) ;
744      *tsysCol = tsys ;
745
746      table_->table().addRow() ;
747      row.put(table_->table().nrow()-1, rec) ;
748    }
749    else {
750      count++ ;
751    }
752    // DEBUG
753    //int rownum = nreader_->getRowNum() ;
754    //cout << "STFiller::readNRO() Finished row " << i << "/" << rownum << endl ;
755    //
756  }
757
758  // DEBUG
759  time_t t1 ;
760  time( &t1 ) ;
761  ttm = localtime( &t1 ) ;
762  cout << "STFiller::readNRO()  Processed " << i << " rows" << endl ;
763  cout << "STFiller::readNRO()  Added " << i - count << " rows (ignored "
764       << count << " \"ZERO\" scans)" << endl ;
765  cout << "STFiller::readNRO()  End time = " << t1
766       << " ("
767       << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
768       << " "
769       << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
770       << ")" << endl ;
771  cout << "STFiller::readNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
772  //
773
774  return 0 ;
775}
776
777Bool STFiller::fileCheck()
778{
779  // if filename_ is directory, return false
780  File inFile( filename_ ) ;
781  if ( inFile.isDirectory() )
782    return false ;
783 
784  // if beginning of header data is "RW", return true
785  // otherwise, return false ;
786  FILE *fp = fopen( filename_.c_str(), "r" ) ;
787  char buf[9] ;
788  fread( buf, 4, 1, fp ) ;
789  buf[4] = '\0' ;
790  if ( ( strncmp( buf, "RW", 2 ) == 0 ) )
791    return true ;
792
793  return false ;
794}
795
796}//namespace asap
Note: See TracBrowser for help on using the repository browser.