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

Last change on this file since 1492 was 1492, 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: Changed reference of the value for

FreqRefVal? field in the MAIN table.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 23.9 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 )
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);
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  // set velocity reference
559//   table_->frequencies().setFrame( vref ) ;
560
561  // DEBUG
562  //cout << "STFiller::openNRO() Velocity Definition = " << nheader->getVDEF() << endl ;
563
564  // DEBUG
565  time_t t1 ;
566  time( &t1 ) ;
567  ttm = localtime( &t1 ) ;
568  cout << "STFiller::openNRO()  End time = " << t1
569       << " ("
570       << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
571       << " "
572       << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
573       << ")" << endl ;
574  cout << "STFiller::openNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
575  //
576
577  return ;
578}
579
580int STFiller::readNRO()
581{
582  // DEBUG
583  time_t t0 ;
584  time( &t0 ) ;
585  tm *ttm = localtime( &t0 ) ;
586  cout << "STFiller::readNRO()  Start time = " << t0
587       << " ("
588       << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
589       << " "
590       << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
591       << ")" << endl ;
592  //
593
594  // get header
595  NROHeader *h = nreader_->getHeader() ;
596
597  // fill row
598  uInt id ;
599  uInt imax = nreader_->getRowNum() ;
600  vector< vector<double > > freqs ;
601  uInt i = 0 ;
602  int count = 0 ;
603  for ( i = 0 ; i < imax ; i++ ) {
604    if( nreader_->getData( i ) != 0 ) {
605      cerr << "STFiller::readNRO()  error while reading row " << i << endl ;
606      return -1 ;
607    }
608    NRODataset *d = nreader_->getData() ;
609
610    char *scanType = d->SCANTP ;
611    Int srcType = -1 ;
612    if ( strncmp( scanType, "ON", 2 ) == 0 ) {
613      srcType = 0 ;
614    }
615    else if ( strncmp( scanType, "OFF", 3 ) == 0 ) {
616      //cout << "OFF srcType: " << i << endl ;
617      srcType = 1 ;
618    }
619    else if ( strncmp( scanType, "ZERO", 4 ) == 0 ) {
620      //cout << "ZERO srcType: " << i << endl ;
621      srcType = 2 ;
622    }
623    else {
624      //cout << "Undefined srcType: " << i << endl ;
625      srcType = 3 ;
626    }
627 
628    // if srcType is 2 (ZERO scan), ignore scan
629    if ( srcType != 2 && srcType != -1 && srcType != 3 ) {
630      TableRow row( table_->table() ) ;
631      TableRecord& rec = row.record();
632
633      RecordFieldPtr<Int> srctCol(rec, "SRCTYPE") ;
634      *srctCol = srcType ;     
635      RecordFieldPtr<Array<Double> > srateCol(rec, "SCANRATE");
636      Array<Double> srcarr( IPosition( 1, 2 ) ) ;
637      srcarr = 0.0 ;
638      *srateCol = srcarr ;
639      RecordFieldPtr<Array<Double> > spmCol(rec, "SRCPROPERMOTION") ;
640      *spmCol = srcarr ;
641      RecordFieldPtr<Array<Double> > sdirCol(rec, "SRCDIRECTION") ;
642      *sdirCol = nreader_->getSourceDirection() ;
643      RecordFieldPtr<Double> svelCol(rec, "SRCVELOCITY") ;
644      *svelCol = h->getURVEL() ;   // [m/s]
645      //*svelCol = h->getURVEL() + d->VRAD ;  // [m/s]
646      //*svelCol = d->VRAD ;  // [m/s]
647      //*svelCol = 0.0 ;
648      RecordFieldPtr<Int> fitCol(rec, "FIT_ID") ;
649      *fitCol = -1 ;
650      RecordFieldPtr<uInt> scannoCol( rec, "SCANNO" ) ;
651      //*scannoCol = (uInt)(d->getISCAN()) ;
652      *scannoCol = (uInt)(d->ISCAN) ;
653      RecordFieldPtr<uInt> cyclenoCol(rec, "CYCLENO");
654      RecordFieldPtr<Double> mjdCol( rec, "TIME" ) ;
655      *mjdCol = Double( nreader_->getStartIntTime( i ) ) ;
656      RecordFieldPtr<Double> intervalCol( rec, "INTERVAL" ) ;
657      *intervalCol = Double( h->getIPTIM() ) ;
658      RecordFieldPtr<String> srcnCol(rec, "SRCNAME") ;
659      *srcnCol = String( h->getOBJ() ) ;
660      RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
661      *fieldnCol = String( h->getOBJ() ) ;
662      // BEAMNO is 0-base
663      RecordFieldPtr<uInt> beamCol(rec, "BEAMNO") ;
664      string arryt = string( d->ARRYT ) ;
665      string sbeamno = arryt.substr( 1, arryt.size()-1 ) ;
666      uInt ibeamno = atoi( sbeamno.c_str() ) ;
667      *beamCol = ibeamno - 1 ;
668      RecordFieldPtr<Int> rbCol(rec, "REFBEAMNO") ;
669      RecordFieldPtr<uInt> ifCol(rec, "IFNO") ;
670      RecordFieldPtr<uInt> mfreqidCol(rec, "FREQ_ID") ;
671      vector<double> fqs = nreader_->getFrequencies( i ) ;
672      //cout << "STFiller::readNRO()  fqs[1] = " << fqs[1] << endl ;
673      if ( freqs.size() == 0 ) {
674        id = table_->frequencies().addEntry( Double( fqs[0] ),
675                                             Double( fqs[1] ),
676                                             Double( fqs[2] ) ) ;
677        *mfreqidCol = id ;
678        *ifCol = id ;
679        freqs.push_back( fqs ) ;
680      }
681      else {
682        int iadd = -1 ;
683        for ( uInt iif = 0 ; iif < freqs.size() ; iif++ ) {
684          //cout << "STFiller::readNRO()  freqs[" << iif << "][1] = " << freqs[iif][1] << endl ;
685          double fdiff = abs( freqs[iif][1] - fqs[1] ) / freqs[iif][1] ;
686          //cout << "STFiller::readNRO()  fdiff = " << fdiff << endl ;
687          if ( fdiff < 1.0e-8 ) {
688            iadd = iif ;
689            break ;
690          }
691        }
692        if ( iadd == -1 ) {
693          id = table_->frequencies().addEntry( Double( fqs[0] ),
694                                               Double( fqs[1] ),
695                                               Double( fqs[2] ) ) ;
696          *mfreqidCol = id ;
697          *ifCol = id ;
698          freqs.push_back( fqs ) ;
699        }
700        else {
701          *mfreqidCol = iadd ;
702          *ifCol = iadd ;
703        }
704      }
705      RecordFieldPtr<uInt> molidCol(rec, "MOLECULE_ID") ;
706      Vector<Double> restfreq( IPosition( 1, 1 ) ) ;
707      restfreq( 0 ) = d->FREQ0 ;
708      id = table_->molecules().addEntry( restfreq ) ;
709      *molidCol = id ;
710      RecordFieldPtr<uInt> mcalidCol(rec, "TCAL_ID") ;
711      //
712      // No Tcal information in the data
713      //
714      // 2008/11/20 Takeshi Nakazato
715      //
716      *mcalidCol = 0 ;
717      RecordFieldPtr<uInt> mweatheridCol(rec, "WEATHER_ID") ;
718
719      id = table_->weather().addEntry( Float( d->TEMP ),
720                                       Float( d->PATM ),
721                                       Float( d->PH2O ),
722                                       Float( d->VWIND ),
723                                       Float( d->DWIND ) ) ;
724
725      *mweatheridCol = id ;         
726      RecordFieldPtr<uInt> mfocusidCol(rec, "FOCUS_ID") ;
727      RecordFieldPtr< Array<Double> > dirCol(rec, "DIRECTION") ;
728      *dirCol = nreader_->getDirection( i ) ;
729      RecordFieldPtr<Float> azCol(rec, "AZIMUTH") ;
730      *azCol = d->RAZ ;
731      RecordFieldPtr<Float> elCol(rec, "ELEVATION") ;
732      *elCol = d->REL ;
733      RecordFieldPtr<Float> parCol(rec, "PARANGLE") ;
734      RecordFieldPtr< Array<Float> > specCol(rec, "SPECTRA") ;
735      vector<double> spec = nreader_->getSpectrum( i ) ;
736      Array<Float> sp( IPosition( 1, spec.size() ) ) ;
737      int index = 0 ;
738      for ( Array<Float>::iterator itr = sp.begin() ; itr != sp.end() ; itr++ ) {
739        *itr = spec[index++] ;
740      }
741      *specCol = sp ;
742      RecordFieldPtr< Array<uChar> > flagCol(rec, "FLAGTRA") ;
743      Array<uChar> flag( sp.shape() ) ;
744      flag.set( 0 ) ;
745      *flagCol = flag ;
746      RecordFieldPtr< uInt > polnoCol(rec, "POLNO") ;
747      *polnoCol = 0 ;
748      RecordFieldPtr< Array<Float> > tsysCol(rec, "TSYS") ;
749      Array<Float> tsys( IPosition( 1, 1 ), d->TSYS ) ;
750      *tsysCol = tsys ;
751
752      table_->table().addRow() ;
753      row.put(table_->table().nrow()-1, rec) ;
754    }
755    else {
756      count++ ;
757    }
758    // DEBUG
759    //int rownum = nreader_->getRowNum() ;
760    //cout << "STFiller::readNRO() Finished row " << i << "/" << rownum << endl ;
761    //
762  }
763
764  // DEBUG
765  time_t t1 ;
766  time( &t1 ) ;
767  ttm = localtime( &t1 ) ;
768  cout << "STFiller::readNRO()  Processed " << i << " rows" << endl ;
769  cout << "STFiller::readNRO()  Added " << i - count << " rows (ignored "
770       << count << " \"ZERO\" scans)" << endl ;
771  cout << "STFiller::readNRO()  End time = " << t1
772       << " ("
773       << ttm->tm_year + 1900 << "/" << ttm->tm_mon + 1 << "/" << ttm->tm_mday
774       << " "
775       << ttm->tm_hour << ":" << ttm->tm_min << ":" << ttm->tm_sec
776       << ")" << endl ;
777  cout << "STFiller::readNRO()  Elapsed time = " << t1 - t0 << " sec" << endl ;
778  //
779
780  return 0 ;
781}
782
783Bool STFiller::fileCheck()
784{
785  // if filename_ is directory, return false
786  File inFile( filename_ ) ;
787  if ( inFile.isDirectory() )
788    return false ;
789 
790  // if beginning of header data is "RW", return true
791  // otherwise, return false ;
792  FILE *fp = fopen( filename_.c_str(), "r" ) ;
793  char buf[9] ;
794  fread( buf, 4, 1, fp ) ;
795  buf[4] = '\0' ;
796  if ( ( strncmp( buf, "RW", 2 ) == 0 ) )
797    return true ;
798
799  return false ;
800}
801
802}//namespace asap
Note: See TracBrowser for help on using the repository browser.