source: trunk/external/atnf/PKSIO/SDFITSreader.cc@ 1361

Last change on this file since 1361 was 1325, checked in by mar637, 18 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: 50.5 KB
RevLine 
[1325]1//#---------------------------------------------------------------------------
2//# SDFITSreader.cc: ATNF CFITSIO interface class for SDFITS input.
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 WITHOUT
13//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14//# 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 this software 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: SDFITSreader.cc,v 19.22 2006/07/12 00:14:26 mcalabre Exp $
29//#---------------------------------------------------------------------------
30//# The SDFITSreader class reads single dish FITS files such as those written
31//# by SDFITSwriter containing Parkes Multibeam data.
32//#
33//# Original: 2000/08/09, Mark Calabretta, ATNF
34//#---------------------------------------------------------------------------
35
36#include <algorithm>
37#include <strings.h>
38
39// AIPS++ includes.
40#include <casa/iostream.h>
41#include <casa/math.h>
42#include <casa/stdio.h>
43
44// ATNF includes.
45#include <atnf/pks/pks_maths.h>
46#include <atnf/PKSIO/PKSMBrecord.h>
47#include <atnf/PKSIO/SDFITSreader.h>
48
49
50class FITSparm
51{
52 public:
53 char *name; // Keyword or column name.
54 int type; // Expected keyvalue or column data type.
55 int colnum; // Column number; 0 for keyword; -1 absent.
56 int coltype; // Column data type, as found.
57 long nelem; // Column data repeat count; < 0 for vardim.
58 int tdimcol; // TDIM column number; 0 for keyword; -1 absent.
59};
60
61// Numerical constants.
62const double PI = 3.141592653589793238462643;
63
64// Factor to convert radians to degrees.
65const double D2R = PI / 180.0;
66
67//------------------------------------------------- SDFITSreader::SDFITSreader
68
69SDFITSreader::SDFITSreader()
70{
71 // Default constructor.
72 cSDptr = 0;
73
74 // Allocate space for data descriptors.
75 cData = new FITSparm[NDATA];
76
77 for (int iData = 0; iData < NDATA; iData++) {
78 cData[iData].colnum = -1;
79 }
80
81 // Initialize pointers.
82 cBeams = 0x0;
83 cIFs = 0x0;
84 cStartChan = 0x0;
85 cEndChan = 0x0;
86 cRefChan = 0x0;
87}
88
89//------------------------------------------------ SDFITSreader::~SDFITSreader
90
91SDFITSreader::~SDFITSreader()
92{
93 close();
94
95 delete [] cData;
96}
97
98//--------------------------------------------------------- SDFITSreader::open
99
100// Open an SDFITS file for reading.
101
102int SDFITSreader::open(
103 char* sdName,
104 int &nBeam,
105 int* &beams,
106 int &nIF,
107 int* &IFs,
108 int* &nChan,
109 int* &nPol,
110 int* &haveXPol,
111 int &haveBase,
112 int &haveSpectra,
113 int &extraSysCal)
114{
115 if (cSDptr) {
116 close();
117 }
118
119 // Open the SDFITS file.
120 cStatus = 0;
121 if (fits_open_file(&cSDptr, sdName, READONLY, &cStatus)) {
122 cerr << "Failed to open SDFITS file: " << sdName << endl;
123 reportError();
124 return 1;
125 }
126
127 // Move to the SDFITS extension.
128 cALFA = cALFA_BD = cALFA_CIMA = 0;
129 if (fits_movnam_hdu(cSDptr, BINARY_TBL, "SINGLE DISH", 0, &cStatus)) {
130 // No SDFITS table, look for BDFITS or CIMAFITS.
131 cStatus = 0;
132 if (fits_movnam_hdu(cSDptr, BINARY_TBL, "BDFITS", 0, &cStatus) == 0) {
133 cALFA_BD = 1;
134
135 } else {
136 cStatus = 0;
137 if (fits_movnam_hdu(cSDptr, BINARY_TBL, "CIMAFITS", 0, &cStatus) == 0) {
138 cALFA_CIMA = 1;
139
140 } else {
141 cerr << "Failed to locate SDFITS binary table." << endl;
142 reportError();
143 close();
144 return 1;
145 }
146 }
147
148 // Arecibo ALFA data of some kind.
149 cALFA = 1;
150 for (int iBeam = 0; iBeam < 8; iBeam++) {
151 for (int iPol = 0; iPol < 2; iPol++) {
152 cALFAcalOn[iBeam][iPol] = 0.0f;
153 cALFAcalOff[iBeam][iPol] = 0.0f;
154
155 // Nominal factor to calibrate spectra in Jy.
156 cALFAcal[iBeam][iPol] = 3.0f;
157 }
158 }
159 }
160
161 // GBT data.
162 char telescope[32];
163 readParm("TELESCOP", TSTRING, telescope); // Core.
164 cGBT = strncmp(telescope, "GBT", 3) == 0 ||
165 strncmp(telescope, "NRAO_GBT", 8) == 0;
166
167 cRow = 0;
168
169
170 // Check that the DATA array column is present.
171 findData(DATA, "DATA", TFLOAT);
172 haveSpectra = cHaveSpectra = cData[DATA].colnum > 0;
173
174 if (cHaveSpectra) {
175 // Find the number of data axes (must be the same for each IF).
176 cNAxis = 5;
177 if (readDim(DATA, 1, &cNAxis, cNAxes)) {
178 reportError();
179 close();
180 return 1;
181 }
182
183 if (cALFA_BD) {
184 // ALFA BDFITS: variable length arrays don't actually vary and there is
185 // no TDIM (or MAXISn) card; use the LAGS_IN value.
186 cNAxis = 5;
187 readParm("LAGS_IN", TLONG, cNAxes);
188 cNAxes[1] = 1;
189 cNAxes[2] = 1;
190 cNAxes[3] = 1;
191 cNAxes[4] = 1;
192 cData[DATA].nelem = cNAxes[0];
193 }
194
195 if (cNAxis < 4) {
196 // Need at least four axes (for now).
197 cerr << "DATA array contains fewer than four axes." << endl;
198 close();
199 return 1;
200 } else if (cNAxis > 5) {
201 // We support up to five axes.
202 cerr << "DATA array contains more than five axes." << endl;
203 close();
204 return 1;
205 }
206
207 findData(FLAGGED, "FLAGGED", TBYTE);
208
209 } else {
210 // DATA column not present, check for a DATAXED keyword.
211 findData(DATAXED, "DATAXED", TSTRING);
212 if (cData[DATAXED].colnum < 0) {
213 cerr << "DATA array column absent from binary table." << endl;
214 close();
215 return 1;
216 }
217
218 // Determine the number of axes and their length.
219 char dataxed[32];
220 readParm("DATAXED", TSTRING, dataxed);
221
222 for (int iaxis = 0; iaxis < 5; iaxis++) cNAxes[iaxis] = 0;
223 sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAxes, cNAxes+1, cNAxes+2,
224 cNAxes+3, cNAxes+4);
225 for (int iaxis = 4; iaxis > -1; iaxis--) {
226 if (cNAxes[iaxis] == 0) cNAxis = iaxis;
227 }
228 }
229
230 char *CTYPE[5] = {"CTYPE1", "CTYPE2", "CTYPE3", "CTYPE4", "CTYPE5"};
231 char *CRPIX[5] = {"CRPIX1", "CRPIX2", "CRPIX3", "CRPIX4", "CRPIX5"};
232 char *CRVAL[5] = {"CRVAL1", "CRVAL2", "CRVAL3", "CRVAL4", "CRVAL5"};
233 char *CDELT[5] = {"CDELT1", "CDELT2", "CDELT3", "CDELT4", "CDELT5"};
234
235 // Find required DATA array axes.
236 char ctype[5][72];
237 for (int iaxis = 0; iaxis < cNAxis; iaxis++) {
238 strcpy(ctype[iaxis], "");
239 readParm(CTYPE[iaxis], TSTRING, ctype[iaxis]); // Core.
240 }
241
242 if (cStatus) {
243 reportError();
244 close();
245 return 1;
246 }
247
248 char *fqCRPIX = 0;
249 char *fqCRVAL = 0;
250 char *fqCDELT = 0;
251 char *raCRVAL = 0;
252 char *decCRVAL = 0;
253 char *timeCRVAL = 0;
254 char *beamCRVAL = 0;
255
256 for (int iaxis = 0; iaxis < cNAxis; iaxis++) {
257 if (strncmp(ctype[iaxis], "FREQ", 4) == 0) {
258 cReqax[0] = iaxis;
259 fqCRPIX = CRPIX[iaxis];
260 fqCRVAL = CRVAL[iaxis];
261 fqCDELT = CDELT[iaxis];
262
263 } else if (strncmp(ctype[iaxis], "STOKES", 6) == 0) {
264 cReqax[1] = iaxis;
265
266 } else if (strncmp(ctype[iaxis], "RA", 2) == 0) {
267 cReqax[2] = iaxis;
268 raCRVAL = CRVAL[iaxis];
269
270 } else if (strncmp(ctype[iaxis], "DEC", 3) == 0) {
271 cReqax[3] = iaxis;
272 decCRVAL = CRVAL[iaxis];
273
274 } else if (strcmp(ctype[iaxis], "TIME") == 0) {
275 // TIME (UTC seconds since midnight) can be a keyword or axis type.
276 timeCRVAL = CRVAL[iaxis];
277
278 } else if (strcmp(ctype[iaxis], "BEAM") == 0) {
279 // BEAM can be a keyword or axis type.
280 beamCRVAL = CRVAL[iaxis];
281 }
282 }
283
284 if (cALFA_BD) {
285 // Fixed in ALFA CIMAFITS.
286 cReqax[2] = 2;
287 raCRVAL = "CRVAL2A";
288
289 cReqax[3] = 3;
290 decCRVAL = "CRVAL3A";
291 }
292
293 // Check that all are present.
294 for (int iaxis = 0; iaxis < 4; iaxis++) {
295 if (cReqax[iaxis] < 0) {
296 cerr << "Could not find required DATA array axes." << endl;
297 close();
298 return 1;
299 }
300 }
301
302 // Set up machinery for data retrieval.
303 findData(SCAN, "SCAN", TINT); // Shared.
304 findData(CYCLE, "CYCLE", TINT); // Additional.
305 findData(DATE_OBS, "DATE-OBS", TSTRING); // Core.
306 findData(TIME, "TIME", TDOUBLE); // Core.
307 findData(EXPOSURE, "EXPOSURE", TFLOAT); // Core.
308 findData(OBJECT, "OBJECT", TSTRING); // Core.
309 findData(OBJ_RA, "OBJ-RA", TDOUBLE); // Additional.
310 findData(OBJ_DEC, "OBJ-DEC", TDOUBLE); // Additional.
311 findData(RESTFRQ, "RESTFRQ", TDOUBLE); // Additional.
312 findData(OBSMODE, "OBSMODE", TSTRING); // Shared.
313
314 findData(BEAM, "BEAM", TSHORT); // Additional.
315 findData(IF, "IF", TSHORT); // Additional.
316 findData(FqRefPix, fqCRPIX, TFLOAT); // Frequency reference pixel.
317 findData(FqRefVal, fqCRVAL, TDOUBLE); // Frequency reference value.
318 findData(FqDelt, fqCDELT, TDOUBLE); // Frequency increment.
319 findData(RA, raCRVAL, TDOUBLE); // Right ascension.
320 findData(DEC, decCRVAL, TDOUBLE); // Declination.
321 findData(SCANRATE, "SCANRATE", TFLOAT); // Additional.
322
323 findData(TSYS, "TSYS", TFLOAT); // Core.
324 findData(CALFCTR, "CALFCTR", TFLOAT); // Additional.
325 findData(XCALFCTR, "XCALFCTR", TFLOAT); // Additional.
326 findData(BASELIN, "BASELIN", TFLOAT); // Additional.
327 findData(BASESUB, "BASESUB", TFLOAT); // Additional.
328 findData(XPOLDATA, "XPOLDATA", TFLOAT); // Additional.
329
330 findData(REFBEAM, "REFBEAM", TSHORT); // Additional.
331 findData(TCAL, "TCAL", TFLOAT); // Shared.
332 findData(TCALTIME, "TCALTIME", TSTRING); // Additional.
333 findData(AZIMUTH, "AZIMUTH", TFLOAT); // Shared.
334 findData(ELEVATIO, "ELEVATIO", TFLOAT); // Shared.
335 findData(PARANGLE, "PARANGLE", TFLOAT); // Additional.
336 findData(FOCUSAXI, "FOCUSAXI", TFLOAT); // Additional.
337 findData(FOCUSTAN, "FOCUSTAN", TFLOAT); // Additional.
338 findData(FOCUSROT, "FOCUSROT", TFLOAT); // Additional.
339 findData(TAMBIENT, "TAMBIENT", TFLOAT); // Shared.
340 findData(PRESSURE, "PRESSURE", TFLOAT); // Shared.
341 findData(HUMIDITY, "HUMIDITY", TFLOAT); // Shared.
342 findData(WINDSPEE, "WINDSPEE", TFLOAT); // Shared.
343 findData(WINDDIRE, "WINDDIRE", TFLOAT); // Shared.
344
345 if (cStatus) {
346 reportError();
347 close();
348 return 1;
349 }
350
351
352 // Check for alternative column names.
353 if (cALFA) {
354 // ALFA data.
355 cALFAscan = 0;
356 cScanNo = 0;
357 if (cALFA_BD) {
358 findData(SCAN, "SCAN_NUMBER", TINT);
359 findData(CYCLE, "PATTERN_NUMBER", TINT);
360 } else if (cALFA_CIMA) {
361 findData(SCAN, "SCAN_ID", TINT);
362 findData(CYCLE, "SUBSCAN", TINT);
363 }
364 } else {
365 readData(SCAN, 1, &cFirstScanNo);
366 }
367
368 cCycleNo = 0;
369 cLastUTC = 0.0;
370
371 // Beam number, 1-relative by default.
372 cBeam_1rel = 1;
373 if (cData[BEAM].colnum < 0) {
374 if (beamCRVAL) {
375 // There is a BEAM axis.
376 findData(BEAM, beamCRVAL, TDOUBLE);
377
378 } else {
379 if (cALFA) {
380 // ALFA data, 0-relative.
381 findData(BEAM, "INPUT_ID", TSHORT);
382 } else {
383 // ms2sdfits output, 0-relative "feed" number.
384 findData(BEAM, "MAIN_FEED1", TSHORT);
385 }
386
387 cBeam_1rel = 0;
388 }
389 }
390
391 // IF number, 1-relative by default.
392 cIF_1rel = 1;
393 if (cALFA && cData[IF].colnum < 0) {
394 // ALFA data, 0-relative.
395 findData(IF, "IFVAL", TSHORT);
396 cIF_1rel = 0;
397 }
398
399 if (cData[TIME].colnum < 0) {
400 if (timeCRVAL) {
401 // There is a TIME axis.
402 findData(TIME, timeCRVAL, TDOUBLE);
403 }
404 }
405
406 // ms2sdfits writes a scalar "TSYS" column that averages the polarizations.
407 int colnum;
408 findCol("SYSCAL_TSYS", &colnum);
409 if (colnum > 0) {
410 // This contains the vector Tsys.
411 findData(TSYS, "SYSCAL_TSYS", TFLOAT);
412 }
413
414 // XPOLDATA?
415
416 if (cData[SCANRATE].colnum < 0) {
417 findData(SCANRATE, "FIELD_POINTING_DIR_RATE", TFLOAT);
418 }
419
420 if (cData[RESTFRQ].colnum < 0) {
421 findData(RESTFRQ, "RESTFREQ", TDOUBLE);
422 if (cData[RESTFRQ].colnum < 0) {
423 findData(RESTFRQ, "SPECTRAL_WINDOW_REST_FREQUENCY", TDOUBLE);
424 }
425 }
426
427 if (cData[OBJ_RA].colnum < 0) {
428 findData(OBJ_RA, "SOURCE_DIRECTION", TDOUBLE);
429 }
430 if (cData[OBJ_DEC].colnum < 0) {
431 findData(OBJ_DEC, "SOURCE_DIRECTION", TDOUBLE);
432 }
433
434 // REFBEAM?
435
436 if (cData[TCAL].colnum < 0) {
437 findData(TCAL, "SYSCAL_TCAL", TFLOAT);
438 } else if (cALFA_BD) {
439 // ALFA BDFITS has a different TCAL with 64 elements - kill it!
440 findData(TCAL, "NO NO NO", TFLOAT);
441 }
442
443 if (cALFA_BD) {
444 // ALFA BDFITS.
445 findData(AZIMUTH, "CRVAL2B", TFLOAT);
446 findData(ELEVATIO, "CRVAL3B", TFLOAT);
447 }
448
449 if (cALFA) {
450 // ALFA data.
451 findData(PARANGLE, "PARA_ANG", TFLOAT);
452 }
453
454 if (cData[TAMBIENT].colnum < 0) {
455 findData(TAMBIENT, "WEATHER_TEMPERATURE", TFLOAT);
456 }
457
458 if (cData[PRESSURE].colnum < 0) {
459 findData(PRESSURE, "WEATHER_PRESSURE", TFLOAT);
460 }
461
462 if (cData[HUMIDITY].colnum < 0) {
463 findData(HUMIDITY, "WEATHER_REL_HUMIDITY", TFLOAT);
464 }
465
466 if (cData[WINDSPEE].colnum < 0) {
467 findData(WINDSPEE, "WEATHER_WIND_SPEED", TFLOAT);
468 }
469
470 if (cData[WINDDIRE].colnum < 0) {
471 findData(WINDDIRE, "WEATHER_WIND_DIRECTION", TFLOAT);
472 }
473
474
475 // Find the number of rows.
476 fits_get_num_rows(cSDptr, &cNRow, &cStatus);
477 if (!cNRow) {
478 cerr << "Table contains no entries." << endl;
479 close();
480 return 1;
481 }
482
483
484 // Determine which beams are present in the data.
485 if (cData[BEAM].colnum > 0) {
486 short *beamCol = new short[cNRow];
487 short beamNul = 1;
488 int anynul;
489 if (fits_read_col(cSDptr, TSHORT, cData[BEAM].colnum, 1, 1, cNRow,
490 &beamNul, beamCol, &anynul, &cStatus)) {
491 reportError();
492 delete [] beamCol;
493 close();
494 return 1;
495 }
496
497 // Find the maximum beam number.
498 cNBeam = cBeam_1rel - 1;
499 for (int irow = 0; irow < cNRow; irow++) {
500 if (beamCol[irow] > cNBeam) {
501 cNBeam = beamCol[irow];
502 }
503
504 // Check validity.
505 if (beamCol[irow] < cBeam_1rel) {
506 cerr << "SDFITS file contains invalid beam number." << endl;
507 delete [] beamCol;
508 close();
509 return 1;
510 }
511 }
512
513 if (!cBeam_1rel) cNBeam++;
514
515 // Find all beams present in the data.
516 cBeams = new int[cNBeam];
517 for (int ibeam = 0; ibeam < cNBeam; ibeam++) {
518 cBeams[ibeam] = 0;
519 }
520
521 for (int irow = 0; irow < cNRow; irow++) {
522 cBeams[beamCol[irow] - cBeam_1rel] = 1;
523 }
524
525 delete [] beamCol;
526
527 } else {
528 // No BEAM column.
529 cNBeam = 1;
530 cBeams = new int[1];
531 cBeams[0] = 1;
532 }
533
534 // Passing back the address of the array allows PKSFITSreader::select() to
535 // modify its elements directly.
536 nBeam = cNBeam;
537 beams = cBeams;
538
539
540 // Determine which IFs are present in the data.
541 if (cData[IF].colnum > 0) {
542 short *IFCol = new short[cNRow];
543 short IFNul = 1;
544 int anynul;
545 if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
546 &IFNul, IFCol, &anynul, &cStatus)) {
547 reportError();
548 delete [] IFCol;
549 close();
550 return 1;
551 }
552
553 // Find the maximum IF number.
554 cNIF = cIF_1rel - 1;
555 for (int irow = 0; irow < cNRow; irow++) {
556 if (IFCol[irow] > cNIF) {
557 cNIF = IFCol[irow];
558 }
559
560 // Check validity.
561 if (IFCol[irow] < cIF_1rel) {
562 cerr << "SDFITS file contains invalid IF number." << endl;
563 delete [] IFCol;
564 close();
565 return 1;
566 }
567 }
568
569 if (!cIF_1rel) cNIF++;
570
571 // Find all IFs present in the data.
572 cIFs = new int[cNIF];
573 cNChan = new int[cNIF];
574 cNPol = new int[cNIF];
575 cHaveXPol = new int[cNIF];
576 cGetXPol = 0;
577
578 for (int iIF = 0; iIF < cNIF; iIF++) {
579 cIFs[iIF] = 0;
580 cNChan[iIF] = 0;
581 cNPol[iIF] = 0;
582 cHaveXPol[iIF] = 0;
583 }
584
585 for (int irow = 0; irow < cNRow; irow++) {
586 int iIF = IFCol[irow] - cIF_1rel;
587 if (cIFs[iIF] == 0) {
588 cIFs[iIF] = 1;
589
590 // Find the axis lengths.
591 if (cHaveSpectra) {
592 if (cData[DATA].nelem < 0) {
593 // Variable dimension array.
594 if (readDim(DATA, irow+1, &cNAxis, cNAxes)) {
595 reportError();
596 close();
597 return 1;
598 }
599 }
600
601 } else {
602 if (cData[DATAXED].colnum > 0) {
603 char dataxed[32];
604 readParm("DATAXED", TSTRING, dataxed);
605
606 sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAxes, cNAxes+1,
607 cNAxes+2, cNAxes+3, cNAxes+4);
608 }
609 }
610
611 // Number of channels and polarizations.
612 cNChan[iIF] = cNAxes[cReqax[0]];
613 cNPol[iIF] = cNAxes[cReqax[1]];
614 cHaveXPol[iIF] = 0;
615
616 // Is cross-polarization data present?
617 if (cData[XPOLDATA].colnum > 0) {
618 // Check that it conforms.
619 int nAxis;
620 long nAxes[2];
621
622 if (readDim(XPOLDATA, irow+1, &nAxis, nAxes)) {
623 reportError();
624 close();
625 return 1;
626 }
627
628 // Default is to get it if we have it.
629 if (nAxis == 2 &&
630 nAxes[0] == 2 &&
631 nAxes[1] == cNChan[iIF]) {
632 cGetXPol = cHaveXPol[iIF] = 1;
633 }
634 }
635 }
636 }
637
638 delete [] IFCol;
639
640 } else {
641 // No IF column.
642 cNIF = 1;
643 cIFs = new int[1];
644 cIFs[0] = 1;
645
646 cNChan = new int[1];
647 cNPol = new int[1];
648 cHaveXPol = new int[1];
649 cGetXPol = 0;
650
651 // Number of channels and polarizations.
652 cNChan[0] = cNAxes[cReqax[0]];
653 cNPol[0] = cNAxes[cReqax[1]];
654 cHaveXPol[0] = 0;
655 }
656
657 if (cALFA) {
658 // ALFA labels each polarization as a separate IF.
659 cNPol[0] = cNIF;
660 cNIF = 1;
661 }
662
663 // Passing back the address of the array allows PKSFITSreader::select() to
664 // modify its elements directly.
665 nIF = cNIF;
666 IFs = cIFs;
667
668 nChan = cNChan;
669 nPol = cNPol;
670 haveXPol = cHaveXPol;
671
672
673 // Default channel range selection.
674 cStartChan = new int[cNIF];
675 cEndChan = new int[cNIF];
676 cRefChan = new int[cNIF];
677
678 for (int iIF = 0; iIF < cNIF; iIF++) {
679 cStartChan[iIF] = 1;
680 cEndChan[iIF] = cNChan[iIF];
681 cRefChan[iIF] = cNChan[iIF]/2 + 1;
682 }
683
684 // Default is to get it if we have it.
685 cGetSpectra = cHaveSpectra;
686
687
688 // Are baseline parameters present?
689 cHaveBase = 0;
690 if (cData[BASELIN].colnum) {
691 // Check that it conforms.
692 int nAxis, status = 0;
693 long nAxes[2];
694
695 if (fits_read_tdim(cSDptr, cData[BASELIN].colnum, 2, &nAxis, nAxes,
696 &status) == 0) {
697 cHaveBase = (nAxis == 2);
698 }
699 }
700 haveBase = cHaveBase;
701
702
703 // Is extra system calibration data available?
704 cExtraSysCal = 0;
705 for (int iparm = REFBEAM; iparm < NDATA; iparm++) {
706 if (cData[iparm].colnum >= 0) {
707 cExtraSysCal = 1;
708 break;
709 }
710 }
711
712 extraSysCal = cExtraSysCal;
713
714 return 0;
715}
716
717//---------------------------------------------------- SDFITSreader::getHeader
718
719// Get parameters describing the data.
720
721int SDFITSreader::getHeader(
722 char observer[32],
723 char project[32],
724 char telescope[32],
725 double antPos[3],
726 char obsMode[32],
727 float &equinox,
728 char radecsys[32],
729 char dopplerFrame[32],
730 char datobs[32],
731 double &utc,
732 double &refFreq,
733 double &bandwidth)
734{
735 // Has the file been opened?
736 if (!cSDptr) {
737 return 1;
738 }
739
740 // Read parameter values.
741 readParm("OBSERVER", TSTRING, observer); // Shared.
742 readParm("PROJID", TSTRING, project); // Shared.
743 readParm("TELESCOP", TSTRING, telescope); // Core.
744
745 antPos[0] = 0.0;
746 antPos[1] = 0.0;
747 antPos[2] = 0.0;
748 if (readParm("ANTENNA_POSITION", TDOUBLE, antPos)) {
749 readParm("OBSGEO-X", TDOUBLE, antPos); // Additional.
750 readParm("OBSGEO-Y", TDOUBLE, antPos + 1); // Additional.
751 readParm("OBSGEO-Z", TDOUBLE, antPos + 2); // Additional.
752 }
753
754 if (antPos[0] == 0.0) {
755 if (strncmp(telescope, "ATPKS", 5) == 0) {
756 // Parkes coordinates.
757 antPos[0] = -4554232.087;
758 antPos[1] = 2816759.046;
759 antPos[2] = -3454035.950;
760 } else if (strncmp(telescope, "ATMOPRA", 7) == 0) {
761 // Mopra coordinates.
762 antPos[0] = -4682768.630;
763 antPos[1] = 2802619.060;
764 antPos[2] = -3291759.900;
765 } else if (strncmp(telescope, "ARECIBO", 7) == 0) {
766 // Arecibo coordinates.
767 antPos[0] = 2390486.900;
768 antPos[1] = -5564731.440;
769 antPos[2] = 1994720.450;
770 }
771 }
772
773 readData(OBSMODE, 1, obsMode); // Shared.
774
775 readParm("EQUINOX", TFLOAT, &equinox); // Shared.
776 if (cStatus == 405) {
777 // EQUINOX was written as string value in early versions.
778 cStatus = 0;
779 char strtmp[32];
780 readParm("EQUINOX", TSTRING, strtmp);
781 sscanf(strtmp, "%f", &equinox);
782 }
783
784 if (readParm("RADESYS", TSTRING, radecsys)) { // Additional.
785 if (readParm("RADECSYS", TSTRING, radecsys)) { // Additional.
786 strcpy(radecsys, "");
787 }
788 }
789
790 if (readParm("SPECSYS", TSTRING, dopplerFrame)) { // Additional.
791 // Fallback value.
792 strcpy(dopplerFrame, "TOPOCENT");
793
794 // Look for VELFRAME, written by earlier versions of Livedata.
795 if (readParm("VELFRAME", TSTRING, dopplerFrame)) { // Additional.
796 // No, try digging it out of the CTYPE card (AIPS convention).
797 char keyw[9], ctype[9];
798 sprintf(keyw, "CTYPE%ld", cReqax[0]+1);
799 readParm(keyw, TSTRING, ctype);
800
801 if (strncmp(ctype, "FREQ-", 5) == 0) {
802 strcpy(dopplerFrame, ctype+5);
803 if (strcmp(dopplerFrame, "LSR") == 0) {
804 // LSR unqualified usually means LSR (kinematic).
805 strcpy(dopplerFrame, "LSRK");
806 } else if (strcmp(dopplerFrame, "HEL") == 0) {
807 // Almost certainly barycentric.
808 strcpy(dopplerFrame, "BARYCENT");
809 }
810 } else {
811 strcpy(dopplerFrame, "");
812 }
813 }
814
815 // Translate to FITS standard names.
816 if (strncmp(dopplerFrame, "TOP", 3) == 0) {
817 strcpy(dopplerFrame, "TOPOCENT");
818 } else if (strncmp(dopplerFrame, "GEO", 3) == 0) {
819 strcpy(dopplerFrame, "GEOCENTR");
820 } else if (strncmp(dopplerFrame, "HEL", 3) == 0) {
821 strcpy(dopplerFrame, "HELIOCEN");
822 } else if (strncmp(dopplerFrame, "BARY", 4) == 0) {
823 strcpy(dopplerFrame, "BARYCENT");
824 }
825 }
826
827 if (cStatus) {
828 reportError();
829 return 1;
830 }
831
832 // Get parameters from first row of table.
833 readData(DATE_OBS, 1, datobs);
834 readData(TIME, 1, &utc);
835 readData(FqRefVal, 1, &refFreq);
836 readParm("BANDWID", TDOUBLE, &bandwidth); // Core.
837
838 if (cALFA_BD) utc *= 3600.0;
839
840 if (cStatus) {
841 reportError();
842 return 1;
843 }
844
845 // Check DATE-OBS format.
846 if (datobs[2] == '/') {
847 // Translate an old-format DATE-OBS.
848 datobs[9] = datobs[1];
849 datobs[8] = datobs[0];
850 datobs[2] = datobs[6];
851 datobs[5] = datobs[3];
852 datobs[3] = datobs[7];
853 datobs[6] = datobs[4];
854 datobs[7] = '-';
855 datobs[4] = '-';
856 datobs[1] = '9';
857 datobs[0] = '1';
858 datobs[10] = '\0';
859
860 } else if (datobs[10] == 'T' && cData[TIME].colnum < 0) {
861 // Dig UTC out of a new-format DATE-OBS.
862 int hh, mm;
863 float ss;
864 sscanf(datobs+11, "%d:%d:%f", &hh, &mm, &ss);
865 utc = (hh*60 + mm)*60 + ss;
866 datobs[10] = '\0';
867 }
868
869 return 0;
870}
871
872//-------------------------------------------------- SDFITSreader::getFreqInfo
873
874// Get frequency parameters for each IF.
875
876int SDFITSreader::getFreqInfo(
877 int &nIF,
878 double* &startFreq,
879 double* &endFreq)
880{
881 float fqRefPix;
882 double fqDelt, fqRefVal;
883
884 nIF = cNIF;
885 startFreq = new double[nIF];
886 endFreq = new double[nIF];
887
888 if (cData[IF].colnum > 0) {
889 short *IFCol = new short[cNRow];
890 short IFNul = 1;
891 int anynul;
892 if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
893 &IFNul, IFCol, &anynul, &cStatus)) {
894 reportError();
895 delete [] IFCol;
896 close();
897 return 1;
898 }
899
900 for (int iIF = 0; iIF < nIF; iIF++) {
901 if (cIFs[iIF]) {
902 // Find the first occurrence of this IF in the table.
903 int IFno = iIF + cIF_1rel;
904 for (int irow = 0; irow < cNRow;) {
905 if (IFCol[irow++] == IFno) {
906 readData(FqRefPix, irow, &fqRefPix);
907 readData(FqRefVal, irow, &fqRefVal);
908 readData(FqDelt, irow, &fqDelt);
909
910 if (cALFA_BD) {
911 unsigned char invert;
912 readData("UPPERSB", TBYTE, irow, &invert);
913
914 if (invert) {
915 fqDelt = -fqDelt;
916 }
917 }
918
919 startFreq[iIF] = fqRefVal + ( 1 - fqRefPix) * fqDelt;
920 endFreq[iIF] = fqRefVal + (cNChan[iIF] - fqRefPix) * fqDelt;
921
922 break;
923 }
924 }
925
926 } else {
927 startFreq[iIF] = 0.0;
928 endFreq[iIF] = 0.0;
929 }
930 }
931
932 delete [] IFCol;
933
934 } else {
935 // No IF column, read the first table entry.
936 readData(FqRefPix, 1, &fqRefPix);
937 readData(FqRefVal, 1, &fqRefVal);
938 readData(FqDelt, 1, &fqDelt);
939
940 startFreq[0] = fqRefVal + ( 1 - fqRefPix) * fqDelt;
941 endFreq[0] = fqRefVal + (cNChan[0] - fqRefPix) * fqDelt;
942 }
943
944 return cStatus;
945}
946
947//---------------------------------------------------- SDFITSreader::findRange
948
949// Find the range of the data in time and position.
950
951int SDFITSreader::findRange(
952 int &nRow,
953 int &nSel,
954 char dateSpan[2][32],
955 double utcSpan[2],
956 double* &positions)
957{
958 int anynul;
959
960 // Has the file been opened?
961 if (!cSDptr) {
962 return 1;
963 }
964
965 nRow = cNRow;
966
967 // Find the number of rows selected.
968 short *sel = new short[nRow];
969 for (int irow = 0; irow < nRow; irow++) {
970 sel[irow] = 1;
971 }
972
973 if (cData[BEAM].colnum > 0) {
974 short *beamCol = new short[cNRow];
975 short beamNul = 1;
976 if (fits_read_col(cSDptr, TSHORT, cData[BEAM].colnum, 1, 1, cNRow,
977 &beamNul, beamCol, &anynul, &cStatus)) {
978 reportError();
979 delete [] beamCol;
980 delete [] sel;
981 return 1;
982 }
983
984 for (int irow = 0; irow < nRow; irow++) {
985 if (!cBeams[beamCol[irow]-cBeam_1rel]) {
986 sel[irow] = 0;
987 }
988 }
989
990 delete [] beamCol;
991 }
992
993 if (cData[IF].colnum > 0) {
994 short *IFCol = new short[cNRow];
995 short IFNul = 1;
996 if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
997 &IFNul, IFCol, &anynul, &cStatus)) {
998 reportError();
999 delete [] IFCol;
1000 delete [] sel;
1001 return 1;
1002 }
1003
1004 for (int irow = 0; irow < nRow; irow++) {
1005 if (!cIFs[IFCol[irow]-cIF_1rel]) {
1006 sel[irow] = 0;
1007 }
1008 }
1009
1010 delete [] IFCol;
1011 }
1012
1013 nSel = 0;
1014 for (int irow = 0; irow < nRow; irow++) {
1015 nSel += sel[irow];
1016 }
1017
1018
1019 // Find the time range assuming the data is in chronological order.
1020 readData(DATE_OBS, 1, dateSpan[0]);
1021 readData(DATE_OBS, nRow, dateSpan[1]);
1022 readData(TIME, 1, utcSpan);
1023 readData(TIME, nRow, utcSpan+1);
1024
1025 if (cALFA_BD) {
1026 utcSpan[0] *= 3600.0;
1027 utcSpan[1] *= 3600.0;
1028 }
1029
1030 // Check DATE-OBS format.
1031 for (int i = 0; i < 2; i++) {
1032 if (dateSpan[0][2] == '/') {
1033 // Translate an old-format DATE-OBS.
1034 dateSpan[i][9] = dateSpan[i][1];
1035 dateSpan[i][8] = dateSpan[i][0];
1036 dateSpan[i][2] = dateSpan[i][6];
1037 dateSpan[i][5] = dateSpan[i][3];
1038 dateSpan[i][3] = dateSpan[i][7];
1039 dateSpan[i][6] = dateSpan[i][4];
1040 dateSpan[i][7] = '-';
1041 dateSpan[i][4] = '-';
1042 dateSpan[i][1] = '9';
1043 dateSpan[i][0] = '1';
1044 dateSpan[i][10] = '\0';
1045 }
1046
1047 if (dateSpan[i][10] == 'T' && cData[TIME].colnum < 0) {
1048 // Dig UTC out of a new-format DATE-OBS.
1049 int hh, mm;
1050 float ss;
1051 sscanf(dateSpan[i]+11, "%d:%d:%f", &hh, &mm, &ss);
1052 utcSpan[i] = (hh*60 + mm)*60 + ss;
1053 }
1054 }
1055
1056
1057 // Retrieve positions for selected data.
1058 double *ra = new double[cNRow];
1059 double *dec = new double[cNRow];
1060 fits_read_col(cSDptr, TDOUBLE, cData[RA].colnum, 1, 1, nRow, 0, ra,
1061 &anynul, &cStatus);
1062 fits_read_col(cSDptr, TDOUBLE, cData[DEC].colnum, 1, 1, nRow, 0, dec,
1063 &anynul, &cStatus);
1064
1065 if (cALFA_BD) {
1066 for (int irow = 0; irow < nRow; irow++) {
1067 // Convert hours to degrees.
1068 ra[irow] *= 15.0;
1069 }
1070 }
1071
1072 int isel = 0;
1073 positions = new double[2*nSel];
1074
1075 // Parameters needed to compute feed-plane coordinates.
1076 double *srcRA, *srcDec;
1077 float *par, *rot;
1078 if (cGetFeedPos) {
1079 srcRA = new double[cNRow];
1080 srcDec = new double[cNRow];
1081 par = new float[cNRow];
1082 rot = new float[cNRow];
1083 fits_read_col(cSDptr, TDOUBLE, cData[OBJ_RA].colnum, 1, 1, nRow, 0,
1084 srcRA, &anynul, &cStatus);
1085 fits_read_col(cSDptr, TDOUBLE, cData[OBJ_DEC].colnum, 1, 1, nRow, 0,
1086 srcDec, &anynul, &cStatus);
1087 fits_read_col(cSDptr, TFLOAT, cData[PARANGLE].colnum, 1, 1, nRow, 0,
1088 par, &anynul, &cStatus);
1089 fits_read_col(cSDptr, TFLOAT, cData[FOCUSROT].colnum, 1, 1, nRow, 0,
1090 rot, &anynul, &cStatus);
1091
1092 for (int irow = 0; irow < nRow; irow++) {
1093 if (sel[irow]) {
1094 // Convert to feed-plane coordinates.
1095 Double dist, pa;
1096 distPA(ra[irow]*D2R, dec[irow]*D2R, srcRA[irow]*D2R, srcDec[irow]*D2R,
1097 dist, pa);
1098
1099 Double spin = (par[irow] + rot[irow])*D2R - pa + PI;
1100 if (spin > 2.0*PI) spin -= 2.0*PI;
1101 Double squint = PI/2.0 - dist;
1102
1103 positions[isel++] = spin;
1104 positions[isel++] = squint;
1105 }
1106 }
1107
1108 delete [] srcRA;
1109 delete [] srcDec;
1110 delete [] par;
1111 delete [] rot;
1112
1113 } else {
1114 for (int irow = 0; irow < nRow; irow++) {
1115 if (sel[irow]) {
1116 positions[isel++] = ra[irow] * D2R;
1117 positions[isel++] = dec[irow] * D2R;
1118 }
1119 }
1120 }
1121
1122 delete [] sel;
1123 delete [] ra;
1124 delete [] dec;
1125
1126 return cStatus;
1127}
1128
1129
1130//--------------------------------------------------------- SDFITSreader::read
1131
1132// Read the next data record.
1133
1134int SDFITSreader::read(
1135 PKSMBrecord &mbrec)
1136{
1137 // Has the file been opened?
1138 if (!cSDptr) {
1139 return 1;
1140 }
1141
1142 // Find the next selected beam and IF.
1143 short iBeam = 0, iIF = 0;
1144 while (++cRow <= cNRow) {
1145 if (cData[BEAM].colnum > 0) {
1146 readData(BEAM, cRow, &iBeam);
1147
1148 // Convert to 0-relative.
1149 if (cBeam_1rel) iBeam--;
1150 }
1151
1152
1153 if (cBeams[iBeam]) {
1154 if (cData[IF].colnum > 0) {
1155 readData(IF, cRow, &iIF);
1156
1157 // Convert to 0-relative.
1158 if (cIF_1rel) iIF--;
1159 }
1160
1161 if (cIFs[iIF]) {
1162 if (cALFA) {
1163 // ALFA data, check for calibration data.
1164 char chars[32];
1165 readData(OBSMODE, cRow, chars);
1166 if (strcmp(chars, "CAL") == 0) {
1167 // iIF is really the polarization in ALFA data.
1168 alfaCal(iBeam, iIF);
1169 continue;
1170 }
1171 }
1172
1173 break;
1174 }
1175 }
1176 }
1177
1178 // EOF?
1179 if (cRow > cNRow) {
1180 return -1;
1181 }
1182
1183
1184 if (cALFA) {
1185 int scanNo;
1186 readData(SCAN, cRow, &scanNo);
1187 if (scanNo != cALFAscan) {
1188 cScanNo++;
1189 cALFAscan = scanNo;
1190 }
1191 mbrec.scanNo = cScanNo;
1192
1193 } else {
1194 readData(SCAN, cRow, &mbrec.scanNo);
1195
1196 // Ensure that scan number is 1-relative.
1197 mbrec.scanNo -= (cFirstScanNo - 1);
1198 }
1199
1200 // Times.
1201 char datobs[32];
1202 readData(DATE_OBS, cRow, datobs);
1203 readData(TIME, cRow, &mbrec.utc);
1204 if (cALFA_BD) mbrec.utc *= 3600.0;
1205
1206 if (datobs[2] == '/') {
1207 // Translate an old-format DATE-OBS.
1208 datobs[9] = datobs[1];
1209 datobs[8] = datobs[0];
1210 datobs[2] = datobs[6];
1211 datobs[5] = datobs[3];
1212 datobs[3] = datobs[7];
1213 datobs[6] = datobs[4];
1214 datobs[7] = '-';
1215 datobs[4] = '-';
1216 datobs[1] = '9';
1217 datobs[0] = '1';
1218
1219 } else if (datobs[10] == 'T' && cData[TIME].colnum < 0) {
1220 // Dig UTC out of a new-format DATE-OBS.
1221 int hh, mm;
1222 float ss;
1223 sscanf(datobs+11, "%d:%d:%f", &hh, &mm, &ss);
1224 mbrec.utc = (hh*60 + mm)*60 + ss;
1225 }
1226
1227 datobs[10] = '\0';
1228 strcpy(mbrec.datobs, datobs);
1229
1230 if (cData[CYCLE].colnum > 0) {
1231 readData(CYCLE, cRow, &mbrec.cycleNo);
1232 if (cALFA_BD) mbrec.cycleNo++;
1233 } else {
1234 // Cycle number not recorded, must do our own bookkeeping.
1235 if (mbrec.utc != cLastUTC) {
1236 mbrec.cycleNo = ++cCycleNo;
1237 cLastUTC = mbrec.utc;
1238 }
1239 }
1240
1241 readData(EXPOSURE, cRow, &mbrec.exposure);
1242
1243 // Source identification.
1244 readData(OBJECT, cRow, mbrec.srcName);
1245
1246 readData(OBJ_RA, cRow, &mbrec.srcRA);
1247 if (strcmp(cData[OBJ_RA].name, "OBJ-RA") == 0) {
1248 mbrec.srcRA *= D2R;
1249 }
1250
1251 if (strcmp(cData[OBJ_DEC].name, "OBJ-DEC") == 0) {
1252 readData(OBJ_DEC, cRow, &mbrec.srcDec);
1253 mbrec.srcDec *= D2R;
1254 }
1255
1256 // Line rest frequency (Hz).
1257 readData(RESTFRQ, cRow, &mbrec.restFreq);
1258 if (mbrec.restFreq == 0.0 && cALFA_BD) {
1259 mbrec.restFreq = 1420.40575e6;
1260 }
1261
1262 // Observation mode.
1263 readData(OBSMODE, cRow, mbrec.obsType);
1264
1265 // Beam-dependent parameters.
1266 mbrec.beamNo = iBeam + 1;
1267
1268 readData(RA, cRow, &mbrec.ra);
1269 readData(DEC, cRow, &mbrec.dec);
1270 mbrec.ra *= D2R;
1271 mbrec.dec *= D2R;
1272
1273 if (cALFA_BD) mbrec.ra *= 15.0;
1274
1275 float scanrate[2];
1276 readData(SCANRATE, cRow, &scanrate);
1277 if (strcmp(cData[SCANRATE].name, "SCANRATE") == 0) {
1278 mbrec.raRate = scanrate[0] * D2R;
1279 mbrec.decRate = scanrate[1] * D2R;
1280 }
1281
1282 // IF-dependent parameters.
1283 int startChan = cStartChan[iIF];
1284 int endChan = cEndChan[iIF];
1285 int refChan = cRefChan[iIF];
1286
1287 // Allocate data storage.
1288 int nChan = abs(endChan - startChan) + 1;
1289 int nPol = cNPol[iIF];
1290
1291 if (cGetSpectra || cGetXPol) {
1292 int nxpol = cGetXPol ? 2*nChan : 0;
1293 mbrec.allocate(0, nChan*nPol, nxpol);
1294 }
1295
1296 mbrec.nIF = 1;
1297 mbrec.IFno[0] = iIF + 1;
1298 mbrec.nChan[0] = nChan;
1299 mbrec.nPol[0] = nPol;
1300
1301 readData(FqRefPix, cRow, mbrec.fqRefPix);
1302 readData(FqRefVal, cRow, mbrec.fqRefVal);
1303 readData(FqDelt, cRow, mbrec.fqDelt);
1304
1305 if (cALFA_BD) {
1306 unsigned char invert;
1307 int anynul, colnum;
1308 findCol("UPPERSB", &colnum);
1309 fits_read_col(cSDptr, TBYTE, colnum, cRow, 1, 1, 0, &invert, &anynul,
1310 &cStatus);
1311
1312 if (invert) {
1313 mbrec.fqDelt[0] = -mbrec.fqDelt[0];
1314 }
1315 }
1316
1317 if (cStatus) {
1318 reportError();
1319 return 1;
1320 }
1321
1322 // Adjust for channel selection.
1323 if (mbrec.fqRefPix[0] != refChan) {
1324 mbrec.fqRefVal[0] += (refChan - mbrec.fqRefPix[0]) * mbrec.fqDelt[0];
1325 mbrec.fqRefPix[0] = refChan;
1326 }
1327
1328 if (endChan < startChan) {
1329 mbrec.fqDelt[0] = -mbrec.fqDelt[0];
1330 }
1331
1332 // The data may only have a scalar Tsys value.
1333 mbrec.tsys[0][0] = 0.0f;
1334 mbrec.tsys[0][1] = 0.0f;
1335 if (cData[TSYS].nelem >= nPol) {
1336 readData(TSYS, cRow, mbrec.tsys[0]);
1337 }
1338
1339 for (int j = 0; j < 2; j++) {
1340 mbrec.calfctr[0][j] = 0.0f;
1341 }
1342 if (cData[CALFCTR].colnum > 0) {
1343 readData(CALFCTR, cRow, mbrec.calfctr);
1344 }
1345
1346 if (cHaveBase) {
1347 mbrec.haveBase = 1;
1348 readData(BASELIN, cRow, mbrec.baseLin);
1349 readData(BASESUB, cRow, mbrec.baseSub);
1350 } else {
1351 mbrec.haveBase = 0;
1352 }
1353
1354 if (cStatus) {
1355 reportError();
1356 return 1;
1357 }
1358
1359 // Read data, sectioning and transposing it in the process.
1360 long *blc = new long[cNAxis+1];
1361 long *trc = new long[cNAxis+1];
1362 long *inc = new long[cNAxis+1];
1363 for (int iaxis = 0; iaxis <= cNAxis; iaxis++) {
1364 blc[iaxis] = 1;
1365 trc[iaxis] = 1;
1366 inc[iaxis] = 1;
1367 }
1368
1369 blc[cReqax[0]] = std::min(startChan, endChan);
1370 trc[cReqax[0]] = std::max(startChan, endChan);
1371 blc[cNAxis] = cRow;
1372 trc[cNAxis] = cRow;
1373
1374 mbrec.haveSpectra = cGetSpectra;
1375 if (cGetSpectra) {
1376 int anynul;
1377
1378 for (int ipol = 0; ipol < nPol; ipol++) {
1379 blc[cReqax[1]] = ipol+1;
1380 trc[cReqax[1]] = ipol+1;
1381
1382 if (cALFA) {
1383 // ALFA data: polarizations are stored in successive rows.
1384 blc[cReqax[1]] = 1;
1385 trc[cReqax[1]] = 1;
1386
1387 if (ipol) {
1388 if (++cRow > cNRow) {
1389 return -1;
1390 }
1391
1392 blc[cNAxis] = cRow;
1393 trc[cNAxis] = cRow;
1394 }
1395
1396 } else if (cData[DATA].nelem < 0) {
1397 // Variable dimension array; get axis lengths.
1398 int naxis = 5, status;
1399
1400 if ((status = readDim(DATA, cRow, &naxis, cNAxes))) {
1401 reportError();
1402
1403 } else if ((status = (naxis != cNAxis))) {
1404 cerr << "DATA array dimensions changed." << endl;
1405 }
1406
1407 if (status) {
1408 delete [] blc;
1409 delete [] trc;
1410 delete [] inc;
1411 return 1;
1412 }
1413 }
1414
1415 if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxis, cNAxes,
1416 blc, trc, inc, 0, mbrec.spectra[0] + ipol*nChan, &anynul,
1417 &cStatus)) {
1418 reportError();
1419 delete [] blc;
1420 delete [] trc;
1421 delete [] inc;
1422 return 1;
1423 }
1424
1425 if (endChan < startChan) {
1426 // Reverse the spectrum.
1427 float *iptr = mbrec.spectra[0] + ipol*nChan;
1428 float *jptr = iptr + nChan - 1;
1429 float *mid = iptr + nChan/2;
1430 while (iptr < mid) {
1431 float tmp = *iptr;
1432 *(iptr++) = *jptr;
1433 *(jptr--) = tmp;
1434 }
1435 }
1436
1437 if (cALFA) {
1438 // ALFA data, rescale the spectrum.
1439 float *chan = mbrec.spectra[0] + ipol*nChan;
1440 float *chanN = chan + nChan;
1441 while (chan < chanN) {
1442 // Approximate conversion to Jy.
1443 *(chan++) *= cALFAcal[iBeam][iIF];
1444 }
1445 }
1446
1447 if (mbrec.tsys[0][ipol] == 0.0) {
1448 // Compute Tsys as the average across the spectrum.
1449 float *chan = mbrec.spectra[0] + ipol*nChan;
1450 float *chanN = chan + nChan;
1451 float *tsys = mbrec.tsys[0] + ipol;
1452 while (chan < chanN) {
1453 *tsys += *(chan++);
1454 }
1455
1456 *tsys /= nChan;
1457 }
1458
1459 // Read data flags.
1460 if (cData[FLAGGED].colnum > 0) {
1461 if (fits_read_subset_byt(cSDptr, cData[FLAGGED].colnum, cNAxis,
1462 cNAxes, blc, trc, inc, 0, mbrec.flagged[0] + ipol*nChan, &anynul,
1463 &cStatus)) {
1464 reportError();
1465 delete [] blc;
1466 delete [] trc;
1467 delete [] inc;
1468 return 1;
1469 }
1470
1471 if (endChan < startChan) {
1472 // Reverse the flag vector.
1473 unsigned char *iptr = mbrec.flagged[0] + ipol*nChan;
1474 unsigned char *jptr = iptr + nChan - 1;
1475 for (int ichan = 0; ichan < nChan/2; ichan++) {
1476 unsigned char tmp = *iptr;
1477 *(iptr++) = *jptr;
1478 *(jptr--) = tmp;
1479 }
1480 }
1481
1482 } else {
1483 // All channels are unflagged by default.
1484 unsigned char *iptr = mbrec.flagged[0] + ipol*nChan;
1485 for (int ichan = 0; ichan < nChan; ichan++) {
1486 *(iptr++) = 0;
1487 }
1488 }
1489 }
1490 }
1491
1492
1493 // Read cross-polarization data.
1494 if (cGetXPol) {
1495 int anynul;
1496 for (int j = 0; j < 2; j++) {
1497 mbrec.xcalfctr[0][j] = 0.0f;
1498 }
1499 if (cData[XCALFCTR].colnum > 0) {
1500 readData(XCALFCTR, cRow, mbrec.xcalfctr);
1501 }
1502
1503 blc[0] = 1;
1504 trc[0] = 2;
1505 blc[1] = std::min(startChan, endChan);
1506 trc[1] = std::max(startChan, endChan);
1507 blc[2] = cRow;
1508 trc[2] = cRow;
1509
1510 int nAxis = 2;
1511 long nAxes[] = {2, nChan};
1512
1513 if (fits_read_subset_flt(cSDptr, cData[XPOLDATA].colnum, nAxis, nAxes,
1514 blc, trc, inc, 0, mbrec.xpol[0], &anynul, &cStatus)) {
1515 reportError();
1516 delete [] blc;
1517 delete [] trc;
1518 delete [] inc;
1519 return 1;
1520 }
1521
1522 if (endChan < startChan) {
1523 // Invert the cross-polarization spectrum.
1524 float *iptr = mbrec.xpol[0];
1525 float *jptr = iptr + nChan - 2;
1526 for (int ichan = 0; ichan < nChan/2; ichan++) {
1527 float tmp = *iptr;
1528 *iptr = *jptr;
1529 *jptr = tmp;
1530
1531 tmp = *(iptr+1);
1532 *(iptr+1) = *(jptr+1);
1533 *(jptr+1) = tmp;
1534
1535 iptr += 2;
1536 jptr -= 2;
1537 }
1538 }
1539 }
1540
1541 delete [] blc;
1542 delete [] trc;
1543 delete [] inc;
1544
1545 if (cStatus) {
1546 reportError();
1547 return 1;
1548 }
1549
1550 mbrec.extraSysCal = cExtraSysCal;
1551 readData(REFBEAM, cRow, &mbrec.refBeam);
1552 readData(TCAL, cRow, &mbrec.tcal[0]);
1553 readData(TCALTIME, cRow, mbrec.tcalTime);
1554 readData(AZIMUTH, cRow, &mbrec.azimuth);
1555 readData(ELEVATIO, cRow, &mbrec.elevation);
1556 readData(PARANGLE, cRow, &mbrec.parAngle);
1557 readData(FOCUSAXI, cRow, &mbrec.focusAxi);
1558 readData(FOCUSTAN, cRow, &mbrec.focusTan);
1559 readData(FOCUSROT, cRow, &mbrec.focusRot);
1560 readData(TAMBIENT, cRow, &mbrec.temp);
1561 readData(PRESSURE, cRow, &mbrec.pressure);
1562 readData(HUMIDITY, cRow, &mbrec.humidity);
1563 readData(WINDSPEE, cRow, &mbrec.windSpeed);
1564 readData(WINDDIRE, cRow, &mbrec.windAz);
1565
1566 if (cALFA_BD) {
1567 // ALFA BDFITS stores zenith angle rather than elevation.
1568 mbrec.elevation = 90.0 - mbrec.elevation;
1569 }
1570
1571 mbrec.azimuth *= D2R;
1572 mbrec.elevation *= D2R;
1573 mbrec.parAngle *= D2R;
1574 mbrec.focusRot *= D2R;
1575 mbrec.windAz *= D2R;
1576
1577 if (cStatus) {
1578 reportError();
1579 return 1;
1580 }
1581
1582 return 0;
1583}
1584
1585
1586//------------------------------------------------------ SDFITSreader::alfaCal
1587
1588// Process ALFA calibration data.
1589
1590int SDFITSreader::alfaCal(
1591 short iBeam,
1592 short iPol)
1593{
1594 int calOn;
1595 char chars[32];
1596 if (cALFA_BD) {
1597 readData("OBS_NAME", TSTRING, cRow, chars);
1598 } else {
1599 readData("SCANTYPE", TSTRING, cRow, chars);
1600 }
1601
1602 if (strcmp(chars, "ON") == 0) {
1603 calOn = 1;
1604 } else if (strcmp(chars, "OFF") == 0) {
1605 calOn = 0;
1606 } else {
1607 return 1;
1608 }
1609
1610 // Read cal data.
1611 long *blc = new long[cNAxis+1];
1612 long *trc = new long[cNAxis+1];
1613 long *inc = new long[cNAxis+1];
1614 for (int iaxis = 0; iaxis <= cNAxis; iaxis++) {
1615 blc[iaxis] = 1;
1616 trc[iaxis] = 1;
1617 inc[iaxis] = 1;
1618 }
1619
1620 // User channel selection.
1621 int startChan = cStartChan[0];
1622 int endChan = cEndChan[0];
1623
1624 blc[cNAxis] = cRow;
1625 trc[cNAxis] = cRow;
1626 blc[cReqax[0]] = std::min(startChan, endChan);
1627 trc[cReqax[0]] = std::max(startChan, endChan);
1628 blc[cReqax[1]] = 1;
1629 trc[cReqax[1]] = 1;
1630
1631 float spectrum[endChan];
1632 int anynul;
1633 if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxis, cNAxes,
1634 blc, trc, inc, 0, spectrum, &anynul, &cStatus)) {
1635 reportError();
1636 delete [] blc;
1637 delete [] trc;
1638 delete [] inc;
1639 return 1;
1640 }
1641
1642 // Average the spectrum.
1643 float mean = 1e9f;
1644 for (int k = 0; k < 2; k++) {
1645 float discrim = 2.0f * mean;
1646
1647 int nChan = 0;
1648 float sum = 0.0f;
1649
1650 float *chanN = spectrum + abs(endChan - startChan) + 1;
1651 for (float *chan = spectrum; chan < chanN; chan++) {
1652 // Simple discriminant that eliminates strong radar interference.
1653 if (*chan < discrim) {
1654 nChan++;
1655 sum += *chan;
1656 }
1657 }
1658
1659 mean = sum / nChan;
1660 }
1661
1662 if (calOn) {
1663 cALFAcalOn[iBeam][iPol] += mean;
1664 } else {
1665 cALFAcalOff[iBeam][iPol] += mean;
1666 }
1667
1668 if (cALFAcalOn[iBeam][iPol] != 0.0f &&
1669 cALFAcalOff[iBeam][iPol] != 0.0f) {
1670 // Tcal should come from the TCAL table, it varies weakly with beam,
1671 // polarization, and frequency. However, TCAL is not written properly.
1672 float Tcal = 12.0f;
1673 cALFAcal[iBeam][iPol] = Tcal / (cALFAcalOn[iBeam][iPol] -
1674 cALFAcalOff[iBeam][iPol]);
1675
1676 // Scale from K to Jy; the gain also varies weakly with beam,
1677 // polarization, frequency, and zenith angle.
1678 float fluxCal = 10.0f;
1679 cALFAcal[iBeam][iPol] /= fluxCal;
1680
1681 cALFAcalOn[iBeam][iPol] = 0.0f;
1682 cALFAcalOff[iBeam][iPol] = 0.0f;
1683 }
1684
1685 return 0;
1686}
1687
1688
1689//-------------------------------------------------- SDFITSreader::reportError
1690
1691// Print the error message corresponding to the input status value and all the
1692// messages on the CFITSIO error stack to stderr.
1693
1694void SDFITSreader::reportError()
1695{
1696 fits_report_error(stderr, cStatus);
1697}
1698
1699//-------------------------------------------------------- SDFITSreader::close
1700
1701// Close the SDFITS file.
1702
1703void SDFITSreader::close()
1704{
1705 if (cSDptr) {
1706 int status = 0;
1707 fits_close_file(cSDptr, &status);
1708 cSDptr = 0;
1709
1710 if (cBeams) delete [] cBeams;
1711 if (cIFs) delete [] cIFs;
1712 if (cStartChan) delete [] cStartChan;
1713 if (cEndChan) delete [] cEndChan;
1714 if (cRefChan) delete [] cRefChan;
1715 }
1716}
1717
1718//----------------------------------------------------- SDFITSreader::findData
1719
1720// Locate a data item in the SDFITS file.
1721
1722void SDFITSreader::findData(
1723 int iData,
1724 char *name,
1725 int type)
1726{
1727 cData[iData].name = name;
1728 cData[iData].type = type;
1729
1730 int colnum;
1731 findCol(name, &colnum);
1732 cData[iData].colnum = colnum;
1733
1734 // Determine the number of data elements.
1735 if (colnum > 0) {
1736 int coltype;
1737 long nelem, width;
1738 fits_get_coltype(cSDptr, colnum, &coltype, &nelem, &width, &cStatus);
1739
1740 // Look for a TDIMnnn keyword or column.
1741 char tdim[8];
1742 sprintf(tdim, "TDIM%d", colnum);
1743 findCol(tdim, &cData[iData].tdimcol);
1744
1745 if (coltype < 0) {
1746 // CFITSIO returns coltype < 0 for variable length arrays.
1747 cData[iData].coltype = -coltype;
1748 cData[iData].nelem = -nelem;
1749
1750 } else {
1751 cData[iData].coltype = coltype;
1752
1753 // Is there a TDIMnnn column?
1754 if (cData[iData].tdimcol > 0) {
1755 // Yes, dimensions of the fixed-length array could still vary.
1756 cData[iData].nelem = -nelem;
1757 } else {
1758 cData[iData].nelem = nelem;
1759 }
1760 }
1761
1762 } else if (colnum == 0) {
1763 // Keyword.
1764 cData[iData].coltype = 0;
1765 cData[iData].nelem = 1;
1766 cData[iData].tdimcol = -1;
1767 }
1768}
1769
1770//------------------------------------------------------ SDFITSreader::readDim
1771
1772// Determine the dimensions of an array in the SDFITS file.
1773
1774int SDFITSreader::readDim(
1775 int iData,
1776 long iRow,
1777 int *naxis,
1778 long naxes[])
1779{
1780 int colnum = cData[iData].colnum;
1781 if (colnum <= 0) {
1782 return 1;
1783 }
1784
1785 int maxdim = *naxis;
1786 if (cData[iData].tdimcol < 0) {
1787 // No TDIMnnn column for this array.
1788 if (cData[iData].nelem < 0) {
1789 // Variable length array; read the array descriptor.
1790 *naxis = 1;
1791 long dummy;
1792 if (fits_read_descript(cSDptr, colnum, iRow, naxes, &dummy, &cStatus)) {
1793 return 1;
1794 }
1795
1796 } else {
1797 // Read the repeat count from TFORMnnn.
1798 if (fits_read_tdim(cSDptr, colnum, maxdim, naxis, naxes, &cStatus)) {
1799 return 1;
1800 }
1801 }
1802
1803 } else {
1804 // Read the TDIMnnn value from the header or table.
1805 char tdim[8], tdimval[64];
1806 sprintf(tdim, "TDIM%d", colnum);
1807 readData(tdim, TSTRING, iRow, tdimval);
1808
1809 // fits_decode_tdim() checks that the TDIMnnn value is within the length
1810 // of the array in the specified column number but unfortunately doesn't
1811 // recognize variable-length arrays. Hence we must decode it here.
1812 char *tp = tdimval;
1813 if (*tp != '(') return 1;
1814
1815 tp++;
1816 *naxis = 0;
1817 for (size_t j = 1; j < strlen(tdimval); j++) {
1818 if (tdimval[j] == ',' || tdimval[j] == ')') {
1819 sscanf(tp, "%ld", naxes + (*naxis)++);
1820 if (tdimval[j] == ')') break;
1821 tp = tdimval + j + 1;
1822 }
1823 }
1824 }
1825
1826 return 0;
1827}
1828
1829//----------------------------------------------------- SDFITSreader::readParm
1830
1831// Read a parameter value from the SDFITS file.
1832
1833int SDFITSreader::readParm(
1834 char *name,
1835 int type,
1836 void *value)
1837{
1838 return readData(name, type, 1, value);
1839}
1840
1841//----------------------------------------------------- SDFITSreader::readData
1842
1843// Read a data value from the SDFITS file.
1844
1845int SDFITSreader::readData(
1846 char *name,
1847 int type,
1848 long iRow,
1849 void *value)
1850{
1851 int colnum;
1852 findCol(name, &colnum);
1853
1854 if (colnum > 0) {
1855 // Read the first value from the specified row of the table.
1856 int coltype;
1857 long nelem, width;
1858 fits_get_coltype(cSDptr, colnum, &coltype, &nelem, &width, &cStatus);
1859
1860 int anynul;
1861 if (type == TSTRING) {
1862 if (nelem) {
1863 fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, &value, &anynul,
1864 &cStatus);
1865 } else {
1866 strcpy((char *)value, "");
1867 }
1868
1869 } else {
1870 if (nelem) {
1871 fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, value, &anynul,
1872 &cStatus);
1873 } else {
1874 if (type == TSHORT) {
1875 *((short *)value) = 0;
1876 } else if (type == TINT) {
1877 *((int *)value) = 0;
1878 } else if (type == TFLOAT) {
1879 *((float *)value) = 0.0f;
1880 } else if (type == TDOUBLE) {
1881 *((double *)value) = 0.0;
1882 }
1883 }
1884 }
1885
1886 } else if (colnum == 0) {
1887 // Read keyword value.
1888 fits_read_key(cSDptr, type, name, value, 0, &cStatus);
1889
1890 } else {
1891 // Not present.
1892 if (type == TSTRING) {
1893 strcpy((char *)value, "");
1894 } else if (type == TSHORT) {
1895 *((short *)value) = 0;
1896 } else if (type == TINT) {
1897 *((int *)value) = 0;
1898 } else if (type == TFLOAT) {
1899 *((float *)value) = 0.0f;
1900 } else if (type == TDOUBLE) {
1901 *((double *)value) = 0.0;
1902 }
1903 }
1904
1905 return colnum < 0;
1906}
1907
1908//----------------------------------------------------- SDFITSreader::readData
1909
1910// Read data from the SDFITS file.
1911
1912int SDFITSreader::readData(
1913 int iData,
1914 long iRow,
1915 void *value)
1916{
1917 char *name = cData[iData].name;
1918 int type = cData[iData].type;
1919 int colnum = cData[iData].colnum;
1920 long nelem = cData[iData].nelem;
1921
1922 if (colnum > 0) {
1923 // Read the required number of values from the specified row of the table.
1924 int anynul;
1925 if (type == TSTRING) {
1926 if (nelem) {
1927 fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, &value, &anynul,
1928 &cStatus);
1929 } else {
1930 strcpy((char *)value, "");
1931 }
1932
1933 } else {
1934 if (nelem) {
1935 fits_read_col(cSDptr, type, colnum, iRow, 1, abs(nelem), 0, value,
1936 &anynul, &cStatus);
1937 } else {
1938 if (type == TSHORT) {
1939 *((short *)value) = 0;
1940 } else if (type == TINT) {
1941 *((int *)value) = 0;
1942 } else if (type == TFLOAT) {
1943 *((float *)value) = 0.0f;
1944 } else if (type == TDOUBLE) {
1945 *((double *)value) = 0.0;
1946 }
1947 }
1948 }
1949
1950 } else if (colnum == 0) {
1951 // Read keyword value.
1952 fits_read_key(cSDptr, type, name, value, 0, &cStatus);
1953
1954 } else {
1955 // Not present.
1956 if (type == TSTRING) {
1957 strcpy((char *)value, "");
1958 } else if (type == TSHORT) {
1959 *((short *)value) = 0;
1960 } else if (type == TINT) {
1961 *((int *)value) = 0;
1962 } else if (type == TFLOAT) {
1963 *((float *)value) = 0.0f;
1964 } else if (type == TDOUBLE) {
1965 *((double *)value) = 0.0;
1966 }
1967 }
1968
1969 return colnum < 0;
1970}
1971
1972//------------------------------------------------------ SDFITSreader::findCol
1973
1974// Locate a parameter in the SDFITS file.
1975
1976void SDFITSreader::findCol(
1977 char *name,
1978 int *colnum)
1979{
1980 *colnum = 0;
1981 int status = 0;
1982 fits_get_colnum(cSDptr, CASESEN, name, colnum, &status);
1983
1984 if (status) {
1985 // Not a real column - maybe it's virtual.
1986 char card[81];
1987
1988 status = 0;
1989 fits_read_card(cSDptr, name, card, &status);
1990 if (status) {
1991 // Not virtual either.
1992 *colnum = -1;
1993 }
1994
1995 // Clear error messages.
1996 fits_clear_errmsg();
1997 }
1998}
Note: See TracBrowser for help on using the repository browser.