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

Last change on this file since 1430 was 1427, checked in by Malte Marquarding, 16 years ago

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