source: trunk/external/atnf/PKSIO/PKSFITSreader.cc @ 1399

Last change on this file since 1399 was 1399, checked in by Malte Marquarding, 17 years ago

Mark C added brightness unit to getHeader()

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