source: branches/alma/external/atnf/PKSIO/PKSFITSreader.cc @ 1453

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

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: many

Test Programs: sd.scantable(), sd.scantable.save()

Put in Release Notes: N/A

Description: copied from current casapy code tree


File size: 14.8 KB
Line 
1//#---------------------------------------------------------------------------
2//# PKSFITSreader.cc: Class to read Parkes multibeam data from a FITS file.
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2000-2006
5//# Associated Universities, Inc. Washington DC, USA.
6//#
7//# This library is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU Library General Public License as published by
9//# the Free Software Foundation; either version 2 of the License, or (at your
10//# option) any later version.
11//#
12//# This library is distributed in the hope that it will be useful, but
13//# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14//# or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
15//# License for more details.
16//#
17//# You should have received a copy of the GNU Library General Public License
18//# along with this library; if not, write to the Free Software Foundation,
19//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning AIPS++ should be addressed as follows:
22//#        Internet email: aips2-request@nrao.edu.
23//#        Postal address: AIPS++ Project Office
24//#                        National Radio Astronomy Observatory
25//#                        520 Edgemont Road
26//#                        Charlottesville, VA 22903-2475 USA
27//#
28//# $Id$
29//#---------------------------------------------------------------------------
30//# Original: 2000/08/02, Mark Calabretta, ATNF
31//#---------------------------------------------------------------------------
32
33#include <atnf/PKSIO/MBFITSreader.h>
34#include <atnf/PKSIO/SDFITSreader.h>
35#include <atnf/PKSIO/PKSFITSreader.h>
36#include <atnf/PKSIO/PKSMBrecord.h>
37
38#include <casa/Arrays/Array.h>
39#include <casa/BasicMath/Math.h>
40#include <casa/Quanta/MVTime.h>
41
42#include <casa/stdio.h>
43
44//----------------------------------------------- PKSFITSreader::PKSFITSreader
45
46// Constructor sets the method of position interpolation.
47
48PKSFITSreader::PKSFITSreader(
49        const String fitsType,
50        const Int    retry,
51        const Bool   interpolate)
52{
53  cMBrec.setNIFs(1);
54
55  if (fitsType == "SDFITS") {
56    cReader = new SDFITSreader();
57  } else {
58    cReader = new MBFITSreader(retry, interpolate ? 1 : 0);
59  }
60}
61
62//---------------------------------------------- PKSFITSreader::~PKSFITSreader
63
64// Destructor.
65
66PKSFITSreader::~PKSFITSreader()
67{
68  close();
69  delete cReader;
70}
71
72//-------------------------------------------------------- PKSFITSreader::open
73
74// Open the FITS file for reading.
75
76Int PKSFITSreader::open(
77        const String fitsName,
78        Vector<Bool> &beams,
79        Vector<Bool> &IFs,
80        Vector<uInt> &nChan,
81        Vector<uInt> &nPol,
82        Vector<Bool> &haveXPol,
83        Bool   &haveBase,
84        Bool   &haveSpectra)
85{
86  int    extraSysCal, haveBase_, *haveXPol_, haveSpectra_, nBeam, *nChan_,
87         nIF, *nPol_, status;
88  if ((status = cReader->open((char *)fitsName.chars(), nBeam, cBeams, nIF,
89                              cIFs, nChan_, nPol_, haveXPol_, haveBase_,
90                              haveSpectra_, extraSysCal))) {
91    return status;
92  }
93
94  // Beams present in data.
95  beams.resize(nBeam);
96  for (Int ibeam = 0; ibeam < nBeam; ibeam++) {
97    beams(ibeam) = cBeams[ibeam];
98  }
99
100  // IFs, channels, and polarizations present in data.
101  IFs.resize(nIF);
102  nChan.resize(nIF);
103  nPol.resize(nIF);
104  haveXPol.resize(nIF);
105
106  for (Int iIF = 0; iIF < nIF; iIF++) {
107    IFs(iIF)   = cIFs[iIF];
108    nChan(iIF) = nChan_[iIF];
109    nPol(iIF)  = nPol_[iIF];
110
111    // Cross-polarization data present?
112    haveXPol(iIF) = haveXPol_[iIF];
113  }
114
115  cNBeam = nBeam;
116  cNIF   = nIF;
117  cNChan.assign(nChan);
118  cNPol.assign(nPol);
119  cHaveXPol.assign(haveXPol);
120
121  // Baseline parameters present?
122  haveBase = haveBase_;
123
124  // Spectral data present?
125  haveSpectra = haveSpectra_;
126
127  return 0;
128}
129
130//--------------------------------------------------- PKSFITSreader::getHeader
131
132// Get parameters describing the data.
133
134Int PKSFITSreader::getHeader(
135        String &observer,
136        String &project,
137        String &antName,
138        Vector<Double> &antPosition,
139        String &obsType,
140        Float  &equinox,
141        String &dopplerFrame,
142        Double &mjd,
143        Double &refFreq,
144        Double &bandwidth,
145        String &fluxunit)
146{
147  char   datobs[32], dopplerFrame_[32], observer_[32], obsType_[32],
148         project_[32], radecsys[32], telescope[32];
149  float  equinox_;
150  double antPos[3], utc;
151
152  if (cReader->getHeader(observer_, project_, telescope, antPos, obsType_,
153                         equinox_, radecsys, dopplerFrame_, datobs, utc,
154                         refFreq, bandwidth)) {
155    return 1;
156  }
157
158  fluxunit = "";
159  observer = trim(observer_);
160  project  = trim(project_);
161  antName  = trim(telescope);
162  antPosition.resize(3);
163  antPosition(0) = antPos[0];
164  antPosition(1) = antPos[1];
165  antPosition(2) = antPos[2];
166  obsType = trim(obsType_);
167  equinox = equinox_;
168  dopplerFrame = trim(dopplerFrame_);
169
170  // Get UTC in MJD form.
171  Int day, month, year;
172  sscanf(datobs, "%4d-%2d-%2d", &year, &month, &day);
173  MVTime date(year, month, Double(day));
174  mjd = date.day() + utc/86400.0;
175
176  return 0;
177}
178
179//------------------------------------------------- PKSFITSreader::getFreqInfo
180
181// Get frequency parameters for each IF.
182
183Int PKSFITSreader::getFreqInfo(
184        Vector<Double> &startFreq,
185        Vector<Double> &endFreq)
186{
187  int     nIF;
188  double *startfreq, *endfreq;
189
190  Int status;
191  if (!(status = cReader->getFreqInfo(nIF, startfreq, endfreq))) {
192    startFreq.takeStorage(IPosition(1,nIF), startfreq, TAKE_OVER);
193    endFreq.takeStorage(IPosition(1,nIF), endfreq, TAKE_OVER);
194  }
195
196  return status;
197}
198
199//------------------------------------------------------ PKSFITSreader::select
200
201// Set data selection by beam number and channel.
202
203uInt PKSFITSreader::select(
204        const Vector<Bool> beamSel,
205        const Vector<Bool> IFsel,
206        const Vector<Int>  startChan,
207        const Vector<Int>  endChan,
208        const Vector<Int>  refChan,
209        const Bool getSpectra,
210        const Bool getXPol,
211        const Bool getFeedPos)
212{
213  // Apply beam selection.
214  uInt nBeamSel = beamSel.nelements();
215  for (uInt ibeam = 0; ibeam < cNBeam; ibeam++) {
216    // This modifies FITSreader::cBeams[].
217    if (ibeam < nBeamSel) {
218      cBeams[ibeam] = cBeams[ibeam] && beamSel(ibeam);
219    } else {
220      cBeams[ibeam] = 0;
221    }
222  }
223
224  // Apply IF selection.
225  int *end   = new int[cNIF];
226  int *ref   = new int[cNIF];
227  int *start = new int[cNIF];
228
229  for (uInt iIF = 0; iIF < cNIF; iIF++) {
230    // This modifies FITSreader::cIFs[].
231    if (iIF < IFsel.nelements()) {
232      cIFs[iIF] = cIFs[iIF] && IFsel(iIF);
233    } else {
234      cIFs[iIF] = 0;
235    }
236
237    if (cIFs[iIF]) {
238      if (iIF < startChan.nelements()) {
239        start[iIF] = startChan(iIF);
240
241        if (start[iIF] <= 0) {
242          start[iIF] += cNChan(iIF);
243        } else if (start[iIF] > Int(cNChan(iIF))) {
244          start[iIF]  = cNChan(iIF);
245        }
246
247      } else {
248        start[iIF] = 1;
249      }
250
251      if (iIF < endChan.nelements()) {
252        end[iIF] = endChan(iIF);
253
254        if (end[iIF] <= 0) {
255          end[iIF] += cNChan(iIF);
256        } else if (end[iIF] > Int(cNChan(iIF))) {
257          end[iIF]  = cNChan(iIF);
258        }
259
260      } else {
261        end[iIF] = cNChan(iIF);
262      }
263
264      if (iIF < refChan.nelements()) {
265        ref[iIF] = refChan(iIF);
266      } else {
267        if (start[iIF] <= end[iIF]) {
268          ref[iIF] = start[iIF] + (end[iIF] - start[iIF] + 1)/2;
269        } else {
270          ref[iIF] = start[iIF] - (start[iIF] - end[iIF] + 1)/2;
271        }
272      }
273    }
274  }
275
276  cGetSpectra = getSpectra;
277  cGetXPol    = getXPol;
278  cGetFeedPos = getFeedPos;
279
280  uInt maxNChan = cReader->select(start, end, ref, cGetSpectra, cGetXPol,
281                                  cGetFeedPos);
282
283  delete [] end;
284  delete [] ref;
285  delete [] start;
286
287  return maxNChan;
288}
289
290//--------------------------------------------------- PKSFITSreader::findRange
291
292// Read the TIME column.
293
294Int PKSFITSreader::findRange(
295        Int    &nRow,
296        Int    &nSel,
297        Vector<Double> &timeSpan,
298        Matrix<Double> &positions)
299{
300  char    dateSpan[2][32];
301  double  utcSpan[2];
302  double* posns;
303
304  Int status;
305  if (!(status = cReader->findRange(nRow, nSel, dateSpan, utcSpan, posns))) {
306    timeSpan.resize(2);
307
308    Int day, month, year;
309    sscanf(dateSpan[0], "%4d-%2d-%2d", &year, &month, &day);
310    timeSpan(0) = MVTime(year, month, Double(day)).second() + utcSpan[0];
311    sscanf(dateSpan[1], "%4d-%2d-%2d", &year, &month, &day);
312    timeSpan(1) = MVTime(year, month, Double(day)).second() + utcSpan[1];
313
314    positions.takeStorage(IPosition(2,2,nSel), posns, TAKE_OVER);
315  }
316
317  return status;
318}
319
320//-------------------------------------------------------- PKSFITSreader::read
321
322// Read the next data record.
323
324Int PKSFITSreader::read(
325        Int             &scanNo,
326        Int             &cycleNo,
327        Double          &mjd,
328        Double          &interval,
329        String          &fieldName,
330        String          &srcName,
331        Vector<Double>  &srcDir,
332        Vector<Double>  &srcPM,
333        Double          &srcVel,
334        String          &obsType,
335        Int             &IFno,
336        Double          &refFreq,
337        Double          &bandwidth,
338        Double          &freqInc,
339        Vector<Double>  &restFreq,
340        Vector<Float>   &tcal,
341        String          &tcalTime,
342        Float           &azimuth,
343        Float           &elevation,
344        Float           &parAngle,
345        Float           &focusAxi,
346        Float           &focusTan,
347        Float           &focusRot,
348        Float           &temperature,
349        Float           &pressure,
350        Float           &humidity,
351        Float           &windSpeed,
352        Float           &windAz,
353        Int             &refBeam,
354        Int             &beamNo,
355        Vector<Double>  &direction,
356        Vector<Double>  &scanRate,
357        Vector<Float>   &tsys,
358        Vector<Float>   &sigma,
359        Vector<Float>   &calFctr,
360        Matrix<Float>   &baseLin,
361        Matrix<Float>   &baseSub,
362        Matrix<Float>   &spectra,
363        Matrix<uChar>   &flagged,
364        Complex         &xCalFctr,
365        Vector<Complex> &xPol)
366{
367  Int status;
368
369  if ((status = cReader->read(cMBrec))) {
370    if (status != -1) {
371      status = 1;
372    }
373
374    return status;
375  }
376
377
378  uInt nChan = cMBrec.nChan[0];
379  uInt nPol  = cMBrec.nPol[0];
380
381  scanNo  = cMBrec.scanNo;
382  cycleNo = cMBrec.cycleNo;
383
384  // Extract MJD.
385  Int day, month, year;
386  sscanf(cMBrec.datobs, "%4d-%2d-%2d", &year, &month, &day);
387  mjd = MVTime(year, month, Double(day)).day() + cMBrec.utc/86400.0;
388
389  interval  = cMBrec.exposure;
390
391  fieldName = trim(cMBrec.srcName);
392  srcName   = fieldName;
393  srcDir(0) = cMBrec.srcRA;
394  srcDir(1) = cMBrec.srcDec;
395  srcPM(0)  = 0.0;
396  srcPM(1)  = 0.0;
397  srcVel    = 0.0;
398  obsType   = trim(cMBrec.obsType);
399
400  IFno = cMBrec.IFno[0];
401  Double chanWidth = fabs(cMBrec.fqDelt[0]);
402  refFreq   = cMBrec.fqRefVal[0];
403  bandwidth = chanWidth * nChan;
404  freqInc   = cMBrec.fqDelt[0];
405  //restFreq  = cMBrec.restFreq;
406  restFreq.resize(1);
407  restFreq(0)  = cMBrec.restFreq;
408
409  tcal.resize(nPol);
410  for (uInt ipol = 0; ipol < nPol; ipol++) {
411    tcal(ipol) = cMBrec.tcal[0][ipol];
412  }
413  tcalTime  = trim(cMBrec.tcalTime);
414  azimuth   = cMBrec.azimuth;
415  elevation = cMBrec.elevation;
416  parAngle  = cMBrec.parAngle;
417  focusAxi  = cMBrec.focusAxi;
418  focusTan  = cMBrec.focusTan;
419  focusRot  = cMBrec.focusRot;
420
421  temperature = cMBrec.temp;
422  pressure    = cMBrec.pressure;
423  humidity    = cMBrec.humidity;
424  windSpeed   = cMBrec.windSpeed;
425  windAz      = cMBrec.windAz;
426
427  refBeam = cMBrec.refBeam;
428  beamNo  = cMBrec.beamNo;
429
430  direction(0) = cMBrec.ra;
431  direction(1) = cMBrec.dec;
432  scanRate(0)  = cMBrec.raRate;
433  scanRate(1)  = cMBrec.decRate;
434
435  tsys.resize(nPol);
436  sigma.resize(nPol);
437  calFctr.resize(nPol);
438  for (uInt ipol = 0; ipol < nPol; ipol++) {
439    tsys(ipol)  = cMBrec.tsys[0][ipol];
440    sigma(ipol) = tsys(ipol) / 0.81 / sqrt(interval * chanWidth);
441    calFctr(ipol) = cMBrec.calfctr[0][ipol];
442  }
443
444  if (cMBrec.haveBase) {
445    baseLin.resize(2,nPol);
446    baseSub.resize(9,nPol);
447
448    for (uInt ipol = 0; ipol < nPol; ipol++) {
449      baseLin(0,ipol) = cMBrec.baseLin[0][ipol][0];
450      baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
451
452      for (uInt j = 0; j < 9; j++) {
453        baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
454      }
455    }
456
457  } else {
458    baseLin.resize(0,0);
459    baseSub.resize(0,0);
460  }
461
462  if (cGetSpectra && cMBrec.haveSpectra) {
463    spectra.resize(nChan,nPol);
464    spectra.takeStorage(IPosition(2,nChan,nPol), cMBrec.spectra[0], SHARE);
465
466    flagged.resize(nChan,nPol);
467    flagged.takeStorage(IPosition(2,nChan,nPol), cMBrec.flagged[0], SHARE);
468
469  } else {
470    spectra.resize(0,0);
471    flagged.resize(0,0);
472  }
473
474  if (cGetXPol) {
475    xCalFctr = Complex(cMBrec.xcalfctr[0][0], cMBrec.xcalfctr[0][1]);
476    xPol.resize(nChan);
477    xPol.takeStorage(IPosition(1,nChan), (Complex *)cMBrec.xpol[0], SHARE);
478  }
479
480  return 0;
481}
482
483//-------------------------------------------------------- PKSFITSreader::read
484
485// Read the next data record, just the basics.
486
487Int PKSFITSreader::read(
488        Int           &IFno,
489        Vector<Float> &tsys,
490        Vector<Float> &calFctr,
491        Matrix<Float> &baseLin,
492        Matrix<Float> &baseSub,
493        Matrix<Float> &spectra,
494        Matrix<uChar> &flagged)
495{
496  Int status;
497
498  if ((status = cReader->read(cMBrec))) {
499    if (status != -1) {
500      status = 1;
501    }
502
503    return status;
504  }
505
506  IFno = cMBrec.IFno[0];
507
508  uInt nChan = cMBrec.nChan[0];
509  uInt nPol  = cMBrec.nPol[0];
510
511  tsys.resize(nPol);
512  calFctr.resize(nPol);
513  for (uInt ipol = 0; ipol < nPol; ipol++) {
514    tsys(ipol) = cMBrec.tsys[0][ipol];
515    calFctr(ipol) = cMBrec.calfctr[0][ipol];
516  }
517
518  if (cMBrec.haveBase) {
519    baseLin.resize(2,nPol);
520    baseSub.resize(9,nPol);
521
522    for (uInt ipol = 0; ipol < nPol; ipol++) {
523      baseLin(0,ipol) = cMBrec.baseLin[0][ipol][0];
524      baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
525
526      for (uInt j = 0; j < 9; j++) {
527        baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
528      }
529    }
530
531  } else {
532    baseLin.resize(0,0);
533    baseSub.resize(0,0);
534  }
535
536  if (cGetSpectra && cMBrec.haveSpectra) {
537    spectra.resize(nChan,nPol);
538    spectra.takeStorage(IPosition(2,nChan,nPol), cMBrec.spectra[0], SHARE);
539
540    flagged.resize(nChan,nPol);
541    flagged.takeStorage(IPosition(2,nChan,nPol), cMBrec.flagged[0], SHARE);
542
543  } else {
544    spectra.resize(0,0);
545    flagged.resize(0,0);
546  }
547
548  return 0;
549}
550
551//------------------------------------------------------- PKSFITSreader::close
552
553// Close the FITS file.
554
555void PKSFITSreader::close()
556{
557  cReader->close();
558}
559
560//-------------------------------------------------------- PKSFITSreader::trim
561
562// Trim trailing blanks from a null-terminated character string.
563
564char* PKSFITSreader::trim(char *string)
565{
566  int j = 0, k = 0;
567  while (string[j] != '\0') {
568    if (string[j++] != ' ') {
569      k = j;
570    }
571  }
572
573  string[k] = '\0';
574
575  return string;
576}
Note: See TracBrowser for help on using the repository browser.