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

Last change on this file since 1325 was 1325, checked in by mar637, 17 years ago

Changes to use casacore instead of casa_asap/aips++\nAdded atnf PKSIO library snapshot to external and linking against this local copy
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: PKSFITSreader.cc,v 19.12 2006/07/05 04:49:57 mcalabre Exp $
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{
146  char   datobs[32], dopplerFrame_[32], observer_[32], obsType_[32],
147         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                         equinox_, radecsys, dopplerFrame_, datobs, utc,
153                         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  equinox = equinox_;
166  dopplerFrame = trim(dopplerFrame_);
167
168  // Get UTC in MJD form.
169  Int day, month, year;
170  sscanf(datobs, "%4d-%2d-%2d", &year, &month, &day);
171  MVTime date(year, month, Double(day));
172  mjd = date.day() + utc/86400.0;
173
174  return 0;
175}
176
177//------------------------------------------------- PKSFITSreader::getFreqInfo
178
179// Get frequency parameters for each IF.
180
181Int PKSFITSreader::getFreqInfo(
182        Vector<Double> &startFreq,
183        Vector<Double> &endFreq)
184{
185  int     nIF;
186  double *startfreq, *endfreq;
187
188  Int status;
189  if (!(status = cReader->getFreqInfo(nIF, startfreq, endfreq))) {
190    startFreq.takeStorage(IPosition(1,nIF), startfreq, TAKE_OVER);
191    endFreq.takeStorage(IPosition(1,nIF), endfreq, TAKE_OVER);
192  }
193
194  return status;
195}
196
197//------------------------------------------------------ PKSFITSreader::select
198
199// Set data selection by beam number and channel.
200
201uInt PKSFITSreader::select(
202        const Vector<Bool> beamSel,
203        const Vector<Bool> IFsel,
204        const Vector<Int>  startChan,
205        const Vector<Int>  endChan,
206        const Vector<Int>  refChan,
207        const Bool getSpectra,
208        const Bool getXPol,
209        const Bool getFeedPos)
210{
211  // Apply beam selection.
212  uInt nBeamSel = beamSel.nelements();
213  for (uInt ibeam = 0; ibeam < cNBeam; ibeam++) {
214    // This modifies FITSreader::cBeams[].
215    if (ibeam < nBeamSel) {
216      cBeams[ibeam] = cBeams[ibeam] && beamSel(ibeam);
217    } else {
218      cBeams[ibeam] = 0;
219    }
220  }
221
222  // Apply IF selection.
223  int *end   = new int[cNIF];
224  int *ref   = new int[cNIF];
225  int *start = new int[cNIF];
226
227  for (uInt iIF = 0; iIF < cNIF; iIF++) {
228    // This modifies FITSreader::cIFs[].
229    if (iIF < IFsel.nelements()) {
230      cIFs[iIF] = cIFs[iIF] && IFsel(iIF);
231    } else {
232      cIFs[iIF] = 0;
233    }
234
235    if (cIFs[iIF]) {
236      if (iIF < startChan.nelements()) {
237        start[iIF] = startChan(iIF);
238
239        if (start[iIF] <= 0) {
240          start[iIF] += cNChan(iIF);
241        } else if (start[iIF] > Int(cNChan(iIF))) {
242          start[iIF]  = cNChan(iIF);
243        }
244
245      } else {
246        start[iIF] = 1;
247      }
248
249      if (iIF < endChan.nelements()) {
250        end[iIF] = endChan(iIF);
251
252        if (end[iIF] <= 0) {
253          end[iIF] += cNChan(iIF);
254        } else if (end[iIF] > Int(cNChan(iIF))) {
255          end[iIF]  = cNChan(iIF);
256        }
257
258      } else {
259        end[iIF] = cNChan(iIF);
260      }
261
262      if (iIF < refChan.nelements()) {
263        ref[iIF] = refChan(iIF);
264      } else {
265        if (start[iIF] <= end[iIF]) {
266          ref[iIF] = start[iIF] + (end[iIF] - start[iIF] + 1)/2;
267        } else {
268          ref[iIF] = start[iIF] - (start[iIF] - end[iIF] + 1)/2;
269        }
270      }
271    }
272  }
273
274  cGetSpectra = getSpectra;
275  cGetXPol    = getXPol;
276  cGetFeedPos = getFeedPos;
277
278  uInt maxNChan = cReader->select(start, end, ref, cGetSpectra, cGetXPol,
279                                  cGetFeedPos);
280
281  delete [] end;
282  delete [] ref;
283  delete [] start;
284
285  return maxNChan;
286}
287
288//--------------------------------------------------- PKSFITSreader::findRange
289
290// Read the TIME column.
291
292Int PKSFITSreader::findRange(
293        Int    &nRow,
294        Int    &nSel,
295        Vector<Double> &timeSpan,
296        Matrix<Double> &positions)
297{
298  char    dateSpan[2][32];
299  double  utcSpan[2];
300  double* posns;
301
302  Int status;
303  if (!(status = cReader->findRange(nRow, nSel, dateSpan, utcSpan, posns))) {
304    timeSpan.resize(2);
305
306    Int day, month, year;
307    sscanf(dateSpan[0], "%4d-%2d-%2d", &year, &month, &day);
308    timeSpan(0) = MVTime(year, month, Double(day)).second() + utcSpan[0];
309    sscanf(dateSpan[1], "%4d-%2d-%2d", &year, &month, &day);
310    timeSpan(1) = MVTime(year, month, Double(day)).second() + utcSpan[1];
311
312    positions.takeStorage(IPosition(2,2,nSel), posns, TAKE_OVER);
313  }
314
315  return status;
316}
317
318//-------------------------------------------------------- PKSFITSreader::read
319
320// Read the next data record.
321
322Int PKSFITSreader::read(
323        Int             &scanNo,
324        Int             &cycleNo,
325        Double          &mjd,
326        Double          &interval,
327        String          &fieldName,
328        String          &srcName,
329        Vector<Double>  &srcDir,
330        Vector<Double>  &srcPM,
331        Double          &srcVel,
332        String          &obsType,
333        Int             &IFno,
334        Double          &refFreq,
335        Double          &bandwidth,
336        Double          &freqInc,
337        Double          &restFreq,
338        Vector<Float>   &tcal,
339        String          &tcalTime,
340        Float           &azimuth,
341        Float           &elevation,
342        Float           &parAngle,
343        Float           &focusAxi,
344        Float           &focusTan,
345        Float           &focusRot,
346        Float           &temperature,
347        Float           &pressure,
348        Float           &humidity,
349        Float           &windSpeed,
350        Float           &windAz,
351        Int             &refBeam,
352        Int             &beamNo,
353        Vector<Double>  &direction,
354        Vector<Double>  &scanRate,
355        Vector<Float>   &tsys,
356        Vector<Float>   &sigma,
357        Vector<Float>   &calFctr,
358        Matrix<Float>   &baseLin,
359        Matrix<Float>   &baseSub,
360        Matrix<Float>   &spectra,
361        Matrix<uChar>   &flagged,
362        Complex         &xCalFctr,
363        Vector<Complex> &xPol)
364{
365  Int status;
366
367  if ((status = cReader->read(cMBrec))) {
368    if (status != -1) {
369      status = 1;
370    }
371
372    return status;
373  }
374
375
376  uInt nChan = cMBrec.nChan[0];
377  uInt nPol  = cMBrec.nPol[0];
378
379  scanNo  = cMBrec.scanNo;
380  cycleNo = cMBrec.cycleNo;
381
382  // Extract MJD.
383  Int day, month, year;
384  sscanf(cMBrec.datobs, "%4d-%2d-%2d", &year, &month, &day);
385  mjd = MVTime(year, month, Double(day)).day() + cMBrec.utc/86400.0;
386
387  interval  = cMBrec.exposure;
388
389  fieldName = trim(cMBrec.srcName);
390  srcName   = fieldName;
391  srcDir(0) = cMBrec.srcRA;
392  srcDir(1) = cMBrec.srcDec;
393  srcPM(0)  = 0.0;
394  srcPM(1)  = 0.0;
395  srcVel    = 0.0;
396  obsType   = trim(cMBrec.obsType);
397
398  IFno = cMBrec.IFno[0];
399  Double chanWidth = fabs(cMBrec.fqDelt[0]);
400  refFreq   = cMBrec.fqRefVal[0];
401  bandwidth = chanWidth * nChan;
402  freqInc   = cMBrec.fqDelt[0];
403  restFreq  = cMBrec.restFreq;
404
405  tcal.resize(nPol);
406  for (uInt ipol = 0; ipol < nPol; ipol++) {
407    tcal(ipol) = cMBrec.tcal[0][ipol];
408  }
409  tcalTime  = trim(cMBrec.tcalTime);
410  azimuth   = cMBrec.azimuth;
411  elevation = cMBrec.elevation;
412  parAngle  = cMBrec.parAngle;
413  focusAxi  = cMBrec.focusAxi;
414  focusTan  = cMBrec.focusTan;
415  focusRot  = cMBrec.focusRot;
416
417  temperature = cMBrec.temp;
418  pressure    = cMBrec.pressure;
419  humidity    = cMBrec.humidity;
420  windSpeed   = cMBrec.windSpeed;
421  windAz      = cMBrec.windAz;
422
423  refBeam = cMBrec.refBeam;
424  beamNo  = cMBrec.beamNo;
425
426  direction(0) = cMBrec.ra;
427  direction(1) = cMBrec.dec;
428  scanRate(0)  = cMBrec.raRate;
429  scanRate(1)  = cMBrec.decRate;
430
431  tsys.resize(nPol);
432  sigma.resize(nPol);
433  calFctr.resize(nPol);
434  for (uInt ipol = 0; ipol < nPol; ipol++) {
435    tsys(ipol)  = cMBrec.tsys[0][ipol];
436    sigma(ipol) = tsys(ipol) / 0.81 / sqrt(interval * chanWidth);
437    calFctr(ipol) = cMBrec.calfctr[0][ipol];
438  }
439
440  if (cMBrec.haveBase) {
441    baseLin.resize(2,nPol);
442    baseSub.resize(9,nPol);
443
444    for (uInt ipol = 0; ipol < nPol; ipol++) {
445      baseLin(0,ipol) = cMBrec.baseLin[0][ipol][0];
446      baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
447
448      for (uInt j = 0; j < 9; j++) {
449        baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
450      }
451    }
452
453  } else {
454    baseLin.resize(0,0);
455    baseSub.resize(0,0);
456  }
457
458  if (cGetSpectra && cMBrec.haveSpectra) {
459    spectra.resize(nChan,nPol);
460    spectra.takeStorage(IPosition(2,nChan,nPol), cMBrec.spectra[0], SHARE);
461
462    flagged.resize(nChan,nPol);
463    flagged.takeStorage(IPosition(2,nChan,nPol), cMBrec.flagged[0], SHARE);
464
465  } else {
466    spectra.resize(0,0);
467    flagged.resize(0,0);
468  }
469
470  if (cGetXPol) {
471    xCalFctr = Complex(cMBrec.xcalfctr[0][0], cMBrec.xcalfctr[0][1]);
472    xPol.resize(nChan);
473    xPol.takeStorage(IPosition(1,nChan), (Complex *)cMBrec.xpol[0], SHARE);
474  }
475
476  return 0;
477}
478
479//-------------------------------------------------------- PKSFITSreader::read
480
481// Read the next data record, just the basics.
482
483Int PKSFITSreader::read(
484        Int           &IFno,
485        Vector<Float> &tsys,
486        Vector<Float> &calFctr,
487        Matrix<Float> &baseLin,
488        Matrix<Float> &baseSub,
489        Matrix<Float> &spectra,
490        Matrix<uChar> &flagged)
491{
492  Int status;
493
494  if ((status = cReader->read(cMBrec))) {
495    if (status != -1) {
496      status = 1;
497    }
498
499    return status;
500  }
501
502  IFno = cMBrec.IFno[0];
503
504  uInt nChan = cMBrec.nChan[0];
505  uInt nPol  = cMBrec.nPol[0];
506
507  tsys.resize(nPol);
508  calFctr.resize(nPol);
509  for (uInt ipol = 0; ipol < nPol; ipol++) {
510    tsys(ipol) = cMBrec.tsys[0][ipol];
511    calFctr(ipol) = cMBrec.calfctr[0][ipol];
512  }
513
514  if (cMBrec.haveBase) {
515    baseLin.resize(2,nPol);
516    baseSub.resize(9,nPol);
517
518    for (uInt ipol = 0; ipol < nPol; ipol++) {
519      baseLin(0,ipol) = cMBrec.baseLin[0][ipol][0];
520      baseLin(1,ipol) = cMBrec.baseLin[0][ipol][1];
521
522      for (uInt j = 0; j < 9; j++) {
523        baseSub(j,ipol) = cMBrec.baseSub[0][ipol][j];
524      }
525    }
526
527  } else {
528    baseLin.resize(0,0);
529    baseSub.resize(0,0);
530  }
531
532  if (cGetSpectra && cMBrec.haveSpectra) {
533    spectra.resize(nChan,nPol);
534    spectra.takeStorage(IPosition(2,nChan,nPol), cMBrec.spectra[0], SHARE);
535
536    flagged.resize(nChan,nPol);
537    flagged.takeStorage(IPosition(2,nChan,nPol), cMBrec.flagged[0], SHARE);
538
539  } else {
540    spectra.resize(0,0);
541    flagged.resize(0,0);
542  }
543
544  return 0;
545}
546
547//------------------------------------------------------- PKSFITSreader::close
548
549// Close the FITS file.
550
551void PKSFITSreader::close()
552{
553  cReader->close();
554}
555
556//-------------------------------------------------------- PKSFITSreader::trim
557
558// Trim trailing blanks from a null-terminated character string.
559
560char* PKSFITSreader::trim(char *string)
561{
562  int j = 0, k = 0;
563  while (string[j] != '\0') {
564    if (string[j++] != ' ') {
565      k = j;
566    }
567  }
568
569  string[k] = '\0';
570
571  return string;
572}
Note: See TracBrowser for help on using the repository browser.