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

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

Mark C added brightness unit to getHeader()
File size: 50.9 KB
RevLine 
[1325]1//#---------------------------------------------------------------------------
2//# SDFITSreader.cc: ATNF CFITSIO interface class for SDFITS input.
3//#---------------------------------------------------------------------------
[1399]4//# Copyright (C) 2000-2007
[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//#
[1399]28//# $Id: SDFITSreader.cc,v 19.23 2007/11/12 03:37:56 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 }
1291
1292 // IF-dependent parameters.
1293 int startChan = cStartChan[iIF];
1294 int endChan = cEndChan[iIF];
1295 int refChan = cRefChan[iIF];
1296
1297 // Allocate data storage.
1298 int nChan = abs(endChan - startChan) + 1;
1299 int nPol = cNPol[iIF];
1300
1301 if (cGetSpectra || cGetXPol) {
1302 int nxpol = cGetXPol ? 2*nChan : 0;
1303 mbrec.allocate(0, nChan*nPol, nxpol);
1304 }
1305
1306 mbrec.nIF = 1;
1307 mbrec.IFno[0] = iIF + 1;
1308 mbrec.nChan[0] = nChan;
1309 mbrec.nPol[0] = nPol;
1310
1311 readData(FqRefPix, cRow, mbrec.fqRefPix);
1312 readData(FqRefVal, cRow, mbrec.fqRefVal);
1313 readData(FqDelt, cRow, mbrec.fqDelt);
1314
1315 if (cALFA_BD) {
1316 unsigned char invert;
1317 int anynul, colnum;
1318 findCol("UPPERSB", &colnum);
1319 fits_read_col(cSDptr, TBYTE, colnum, cRow, 1, 1, 0, &invert, &anynul,
1320 &cStatus);
1321
1322 if (invert) {
1323 mbrec.fqDelt[0] = -mbrec.fqDelt[0];
1324 }
1325 }
1326
1327 if (cStatus) {
1328 reportError();
1329 return 1;
1330 }
1331
1332 // Adjust for channel selection.
1333 if (mbrec.fqRefPix[0] != refChan) {
1334 mbrec.fqRefVal[0] += (refChan - mbrec.fqRefPix[0]) * mbrec.fqDelt[0];
1335 mbrec.fqRefPix[0] = refChan;
1336 }
1337
1338 if (endChan < startChan) {
1339 mbrec.fqDelt[0] = -mbrec.fqDelt[0];
1340 }
1341
1342 // The data may only have a scalar Tsys value.
1343 mbrec.tsys[0][0] = 0.0f;
1344 mbrec.tsys[0][1] = 0.0f;
1345 if (cData[TSYS].nelem >= nPol) {
1346 readData(TSYS, cRow, mbrec.tsys[0]);
1347 }
1348
1349 for (int j = 0; j < 2; j++) {
1350 mbrec.calfctr[0][j] = 0.0f;
1351 }
1352 if (cData[CALFCTR].colnum > 0) {
1353 readData(CALFCTR, cRow, mbrec.calfctr);
1354 }
1355
1356 if (cHaveBase) {
1357 mbrec.haveBase = 1;
1358 readData(BASELIN, cRow, mbrec.baseLin);
1359 readData(BASESUB, cRow, mbrec.baseSub);
1360 } else {
1361 mbrec.haveBase = 0;
1362 }
1363
1364 if (cStatus) {
1365 reportError();
1366 return 1;
1367 }
1368
1369 // Read data, sectioning and transposing it in the process.
1370 long *blc = new long[cNAxis+1];
1371 long *trc = new long[cNAxis+1];
1372 long *inc = new long[cNAxis+1];
1373 for (int iaxis = 0; iaxis <= cNAxis; iaxis++) {
1374 blc[iaxis] = 1;
1375 trc[iaxis] = 1;
1376 inc[iaxis] = 1;
1377 }
1378
1379 blc[cReqax[0]] = std::min(startChan, endChan);
1380 trc[cReqax[0]] = std::max(startChan, endChan);
1381 blc[cNAxis] = cRow;
1382 trc[cNAxis] = cRow;
1383
1384 mbrec.haveSpectra = cGetSpectra;
1385 if (cGetSpectra) {
1386 int anynul;
1387
1388 for (int ipol = 0; ipol < nPol; ipol++) {
1389 blc[cReqax[1]] = ipol+1;
1390 trc[cReqax[1]] = ipol+1;
1391
1392 if (cALFA) {
1393 // ALFA data: polarizations are stored in successive rows.
1394 blc[cReqax[1]] = 1;
1395 trc[cReqax[1]] = 1;
1396
1397 if (ipol) {
1398 if (++cRow > cNRow) {
1399 return -1;
1400 }
1401
1402 blc[cNAxis] = cRow;
1403 trc[cNAxis] = cRow;
1404 }
1405
1406 } else if (cData[DATA].nelem < 0) {
1407 // Variable dimension array; get axis lengths.
1408 int naxis = 5, status;
1409
1410 if ((status = readDim(DATA, cRow, &naxis, cNAxes))) {
1411 reportError();
1412
1413 } else if ((status = (naxis != cNAxis))) {
1414 cerr << "DATA array dimensions changed." << endl;
1415 }
1416
1417 if (status) {
1418 delete [] blc;
1419 delete [] trc;
1420 delete [] inc;
1421 return 1;
1422 }
1423 }
1424
1425 if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxis, cNAxes,
1426 blc, trc, inc, 0, mbrec.spectra[0] + ipol*nChan, &anynul,
1427 &cStatus)) {
1428 reportError();
1429 delete [] blc;
1430 delete [] trc;
1431 delete [] inc;
1432 return 1;
1433 }
1434
1435 if (endChan < startChan) {
1436 // Reverse the spectrum.
1437 float *iptr = mbrec.spectra[0] + ipol*nChan;
1438 float *jptr = iptr + nChan - 1;
1439 float *mid = iptr + nChan/2;
1440 while (iptr < mid) {
1441 float tmp = *iptr;
1442 *(iptr++) = *jptr;
1443 *(jptr--) = tmp;
1444 }
1445 }
1446
1447 if (cALFA) {
1448 // ALFA data, rescale the spectrum.
1449 float *chan = mbrec.spectra[0] + ipol*nChan;
1450 float *chanN = chan + nChan;
1451 while (chan < chanN) {
1452 // Approximate conversion to Jy.
1453 *(chan++) *= cALFAcal[iBeam][iIF];
1454 }
1455 }
1456
1457 if (mbrec.tsys[0][ipol] == 0.0) {
1458 // Compute Tsys as the average across the spectrum.
1459 float *chan = mbrec.spectra[0] + ipol*nChan;
1460 float *chanN = chan + nChan;
1461 float *tsys = mbrec.tsys[0] + ipol;
1462 while (chan < chanN) {
1463 *tsys += *(chan++);
1464 }
1465
1466 *tsys /= nChan;
1467 }
1468
1469 // Read data flags.
1470 if (cData[FLAGGED].colnum > 0) {
1471 if (fits_read_subset_byt(cSDptr, cData[FLAGGED].colnum, cNAxis,
1472 cNAxes, blc, trc, inc, 0, mbrec.flagged[0] + ipol*nChan, &anynul,
1473 &cStatus)) {
1474 reportError();
1475 delete [] blc;
1476 delete [] trc;
1477 delete [] inc;
1478 return 1;
1479 }
1480
1481 if (endChan < startChan) {
1482 // Reverse the flag vector.
1483 unsigned char *iptr = mbrec.flagged[0] + ipol*nChan;
1484 unsigned char *jptr = iptr + nChan - 1;
1485 for (int ichan = 0; ichan < nChan/2; ichan++) {
1486 unsigned char tmp = *iptr;
1487 *(iptr++) = *jptr;
1488 *(jptr--) = tmp;
1489 }
1490 }
1491
1492 } else {
1493 // All channels are unflagged by default.
1494 unsigned char *iptr = mbrec.flagged[0] + ipol*nChan;
1495 for (int ichan = 0; ichan < nChan; ichan++) {
1496 *(iptr++) = 0;
1497 }
1498 }
1499 }
1500 }
1501
1502
1503 // Read cross-polarization data.
1504 if (cGetXPol) {
1505 int anynul;
1506 for (int j = 0; j < 2; j++) {
1507 mbrec.xcalfctr[0][j] = 0.0f;
1508 }
1509 if (cData[XCALFCTR].colnum > 0) {
1510 readData(XCALFCTR, cRow, mbrec.xcalfctr);
1511 }
1512
1513 blc[0] = 1;
1514 trc[0] = 2;
1515 blc[1] = std::min(startChan, endChan);
1516 trc[1] = std::max(startChan, endChan);
1517 blc[2] = cRow;
1518 trc[2] = cRow;
1519
1520 int nAxis = 2;
1521 long nAxes[] = {2, nChan};
1522
1523 if (fits_read_subset_flt(cSDptr, cData[XPOLDATA].colnum, nAxis, nAxes,
1524 blc, trc, inc, 0, mbrec.xpol[0], &anynul, &cStatus)) {
1525 reportError();
1526 delete [] blc;
1527 delete [] trc;
1528 delete [] inc;
1529 return 1;
1530 }
1531
1532 if (endChan < startChan) {
1533 // Invert the cross-polarization spectrum.
1534 float *iptr = mbrec.xpol[0];
1535 float *jptr = iptr + nChan - 2;
1536 for (int ichan = 0; ichan < nChan/2; ichan++) {
1537 float tmp = *iptr;
1538 *iptr = *jptr;
1539 *jptr = tmp;
1540
1541 tmp = *(iptr+1);
1542 *(iptr+1) = *(jptr+1);
1543 *(jptr+1) = tmp;
1544
1545 iptr += 2;
1546 jptr -= 2;
1547 }
1548 }
1549 }
1550
1551 delete [] blc;
1552 delete [] trc;
1553 delete [] inc;
1554
1555 if (cStatus) {
1556 reportError();
1557 return 1;
1558 }
1559
1560 mbrec.extraSysCal = cExtraSysCal;
1561 readData(REFBEAM, cRow, &mbrec.refBeam);
1562 readData(TCAL, cRow, &mbrec.tcal[0]);
1563 readData(TCALTIME, cRow, mbrec.tcalTime);
1564 readData(AZIMUTH, cRow, &mbrec.azimuth);
1565 readData(ELEVATIO, cRow, &mbrec.elevation);
1566 readData(PARANGLE, cRow, &mbrec.parAngle);
1567 readData(FOCUSAXI, cRow, &mbrec.focusAxi);
1568 readData(FOCUSTAN, cRow, &mbrec.focusTan);
1569 readData(FOCUSROT, cRow, &mbrec.focusRot);
1570 readData(TAMBIENT, cRow, &mbrec.temp);
1571 readData(PRESSURE, cRow, &mbrec.pressure);
1572 readData(HUMIDITY, cRow, &mbrec.humidity);
1573 readData(WINDSPEE, cRow, &mbrec.windSpeed);
1574 readData(WINDDIRE, cRow, &mbrec.windAz);
1575
1576 if (cALFA_BD) {
1577 // ALFA BDFITS stores zenith angle rather than elevation.
1578 mbrec.elevation = 90.0 - mbrec.elevation;
1579 }
1580
1581 mbrec.azimuth *= D2R;
1582 mbrec.elevation *= D2R;
1583 mbrec.parAngle *= D2R;
1584 mbrec.focusRot *= D2R;
1585 mbrec.windAz *= D2R;
1586
1587 if (cStatus) {
1588 reportError();
1589 return 1;
1590 }
1591
1592 return 0;
1593}
1594
1595
1596//------------------------------------------------------ SDFITSreader::alfaCal
1597
1598// Process ALFA calibration data.
1599
1600int SDFITSreader::alfaCal(
1601 short iBeam,
1602 short iPol)
1603{
1604 int calOn;
1605 char chars[32];
1606 if (cALFA_BD) {
1607 readData("OBS_NAME", TSTRING, cRow, chars);
1608 } else {
1609 readData("SCANTYPE", TSTRING, cRow, chars);
1610 }
1611
1612 if (strcmp(chars, "ON") == 0) {
1613 calOn = 1;
1614 } else if (strcmp(chars, "OFF") == 0) {
1615 calOn = 0;
1616 } else {
1617 return 1;
1618 }
1619
1620 // Read cal data.
1621 long *blc = new long[cNAxis+1];
1622 long *trc = new long[cNAxis+1];
1623 long *inc = new long[cNAxis+1];
1624 for (int iaxis = 0; iaxis <= cNAxis; iaxis++) {
1625 blc[iaxis] = 1;
1626 trc[iaxis] = 1;
1627 inc[iaxis] = 1;
1628 }
1629
1630 // User channel selection.
1631 int startChan = cStartChan[0];
1632 int endChan = cEndChan[0];
1633
1634 blc[cNAxis] = cRow;
1635 trc[cNAxis] = cRow;
1636 blc[cReqax[0]] = std::min(startChan, endChan);
1637 trc[cReqax[0]] = std::max(startChan, endChan);
1638 blc[cReqax[1]] = 1;
1639 trc[cReqax[1]] = 1;
1640
1641 float spectrum[endChan];
1642 int anynul;
1643 if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxis, cNAxes,
1644 blc, trc, inc, 0, spectrum, &anynul, &cStatus)) {
1645 reportError();
1646 delete [] blc;
1647 delete [] trc;
1648 delete [] inc;
1649 return 1;
1650 }
1651
1652 // Average the spectrum.
1653 float mean = 1e9f;
1654 for (int k = 0; k < 2; k++) {
1655 float discrim = 2.0f * mean;
1656
1657 int nChan = 0;
1658 float sum = 0.0f;
1659
1660 float *chanN = spectrum + abs(endChan - startChan) + 1;
1661 for (float *chan = spectrum; chan < chanN; chan++) {
1662 // Simple discriminant that eliminates strong radar interference.
1663 if (*chan < discrim) {
1664 nChan++;
1665 sum += *chan;
1666 }
1667 }
1668
1669 mean = sum / nChan;
1670 }
1671
1672 if (calOn) {
1673 cALFAcalOn[iBeam][iPol] += mean;
1674 } else {
1675 cALFAcalOff[iBeam][iPol] += mean;
1676 }
1677
1678 if (cALFAcalOn[iBeam][iPol] != 0.0f &&
1679 cALFAcalOff[iBeam][iPol] != 0.0f) {
1680 // Tcal should come from the TCAL table, it varies weakly with beam,
1681 // polarization, and frequency. However, TCAL is not written properly.
1682 float Tcal = 12.0f;
1683 cALFAcal[iBeam][iPol] = Tcal / (cALFAcalOn[iBeam][iPol] -
1684 cALFAcalOff[iBeam][iPol]);
1685
1686 // Scale from K to Jy; the gain also varies weakly with beam,
1687 // polarization, frequency, and zenith angle.
1688 float fluxCal = 10.0f;
1689 cALFAcal[iBeam][iPol] /= fluxCal;
1690
1691 cALFAcalOn[iBeam][iPol] = 0.0f;
1692 cALFAcalOff[iBeam][iPol] = 0.0f;
1693 }
1694
1695 return 0;
1696}
1697
1698
1699//-------------------------------------------------- SDFITSreader::reportError
1700
1701// Print the error message corresponding to the input status value and all the
1702// messages on the CFITSIO error stack to stderr.
1703
1704void SDFITSreader::reportError()
1705{
1706 fits_report_error(stderr, cStatus);
1707}
1708
1709//-------------------------------------------------------- SDFITSreader::close
1710
1711// Close the SDFITS file.
1712
1713void SDFITSreader::close()
1714{
1715 if (cSDptr) {
1716 int status = 0;
1717 fits_close_file(cSDptr, &status);
1718 cSDptr = 0;
1719
1720 if (cBeams) delete [] cBeams;
1721 if (cIFs) delete [] cIFs;
1722 if (cStartChan) delete [] cStartChan;
1723 if (cEndChan) delete [] cEndChan;
1724 if (cRefChan) delete [] cRefChan;
1725 }
1726}
1727
1728//----------------------------------------------------- SDFITSreader::findData
1729
1730// Locate a data item in the SDFITS file.
1731
1732void SDFITSreader::findData(
1733 int iData,
1734 char *name,
1735 int type)
1736{
1737 cData[iData].name = name;
1738 cData[iData].type = type;
1739
1740 int colnum;
1741 findCol(name, &colnum);
1742 cData[iData].colnum = colnum;
1743
1744 // Determine the number of data elements.
1745 if (colnum > 0) {
1746 int coltype;
1747 long nelem, width;
1748 fits_get_coltype(cSDptr, colnum, &coltype, &nelem, &width, &cStatus);
[1399]1749 fits_get_bcolparms(cSDptr, colnum, 0x0, cData[iData].units, 0x0, 0x0, 0x0,
1750 0x0, 0x0, 0x0, &cStatus);
[1325]1751
1752 // Look for a TDIMnnn keyword or column.
1753 char tdim[8];
1754 sprintf(tdim, "TDIM%d", colnum);
1755 findCol(tdim, &cData[iData].tdimcol);
1756
1757 if (coltype < 0) {
1758 // CFITSIO returns coltype < 0 for variable length arrays.
1759 cData[iData].coltype = -coltype;
1760 cData[iData].nelem = -nelem;
1761
1762 } else {
1763 cData[iData].coltype = coltype;
1764
1765 // Is there a TDIMnnn column?
1766 if (cData[iData].tdimcol > 0) {
1767 // Yes, dimensions of the fixed-length array could still vary.
1768 cData[iData].nelem = -nelem;
1769 } else {
1770 cData[iData].nelem = nelem;
1771 }
1772 }
1773
1774 } else if (colnum == 0) {
1775 // Keyword.
1776 cData[iData].coltype = 0;
1777 cData[iData].nelem = 1;
1778 cData[iData].tdimcol = -1;
1779 }
1780}
1781
1782//------------------------------------------------------ SDFITSreader::readDim
1783
1784// Determine the dimensions of an array in the SDFITS file.
1785
1786int SDFITSreader::readDim(
1787 int iData,
1788 long iRow,
1789 int *naxis,
1790 long naxes[])
1791{
1792 int colnum = cData[iData].colnum;
1793 if (colnum <= 0) {
1794 return 1;
1795 }
1796
1797 int maxdim = *naxis;
1798 if (cData[iData].tdimcol < 0) {
1799 // No TDIMnnn column for this array.
1800 if (cData[iData].nelem < 0) {
1801 // Variable length array; read the array descriptor.
1802 *naxis = 1;
1803 long dummy;
1804 if (fits_read_descript(cSDptr, colnum, iRow, naxes, &dummy, &cStatus)) {
1805 return 1;
1806 }
1807
1808 } else {
1809 // Read the repeat count from TFORMnnn.
1810 if (fits_read_tdim(cSDptr, colnum, maxdim, naxis, naxes, &cStatus)) {
1811 return 1;
1812 }
1813 }
1814
1815 } else {
1816 // Read the TDIMnnn value from the header or table.
1817 char tdim[8], tdimval[64];
1818 sprintf(tdim, "TDIM%d", colnum);
1819 readData(tdim, TSTRING, iRow, tdimval);
1820
1821 // fits_decode_tdim() checks that the TDIMnnn value is within the length
1822 // of the array in the specified column number but unfortunately doesn't
1823 // recognize variable-length arrays. Hence we must decode it here.
1824 char *tp = tdimval;
1825 if (*tp != '(') return 1;
1826
1827 tp++;
1828 *naxis = 0;
1829 for (size_t j = 1; j < strlen(tdimval); j++) {
1830 if (tdimval[j] == ',' || tdimval[j] == ')') {
1831 sscanf(tp, "%ld", naxes + (*naxis)++);
1832 if (tdimval[j] == ')') break;
1833 tp = tdimval + j + 1;
1834 }
1835 }
1836 }
1837
1838 return 0;
1839}
1840
1841//----------------------------------------------------- SDFITSreader::readParm
1842
1843// Read a parameter value from the SDFITS file.
1844
1845int SDFITSreader::readParm(
1846 char *name,
1847 int type,
1848 void *value)
1849{
1850 return readData(name, type, 1, value);
1851}
1852
1853//----------------------------------------------------- SDFITSreader::readData
1854
1855// Read a data value from the SDFITS file.
1856
1857int SDFITSreader::readData(
1858 char *name,
1859 int type,
1860 long iRow,
1861 void *value)
1862{
1863 int colnum;
1864 findCol(name, &colnum);
1865
1866 if (colnum > 0) {
1867 // Read the first value from the specified row of the table.
1868 int coltype;
1869 long nelem, width;
1870 fits_get_coltype(cSDptr, colnum, &coltype, &nelem, &width, &cStatus);
1871
1872 int anynul;
1873 if (type == TSTRING) {
1874 if (nelem) {
1875 fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, &value, &anynul,
1876 &cStatus);
1877 } else {
1878 strcpy((char *)value, "");
1879 }
1880
1881 } else {
1882 if (nelem) {
1883 fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, value, &anynul,
1884 &cStatus);
1885 } else {
1886 if (type == TSHORT) {
1887 *((short *)value) = 0;
1888 } else if (type == TINT) {
1889 *((int *)value) = 0;
1890 } else if (type == TFLOAT) {
1891 *((float *)value) = 0.0f;
1892 } else if (type == TDOUBLE) {
1893 *((double *)value) = 0.0;
1894 }
1895 }
1896 }
1897
1898 } else if (colnum == 0) {
1899 // Read keyword value.
1900 fits_read_key(cSDptr, type, name, value, 0, &cStatus);
1901
1902 } else {
1903 // Not present.
1904 if (type == TSTRING) {
1905 strcpy((char *)value, "");
1906 } else if (type == TSHORT) {
1907 *((short *)value) = 0;
1908 } else if (type == TINT) {
1909 *((int *)value) = 0;
1910 } else if (type == TFLOAT) {
1911 *((float *)value) = 0.0f;
1912 } else if (type == TDOUBLE) {
1913 *((double *)value) = 0.0;
1914 }
1915 }
1916
1917 return colnum < 0;
1918}
1919
1920//----------------------------------------------------- SDFITSreader::readData
1921
1922// Read data from the SDFITS file.
1923
1924int SDFITSreader::readData(
1925 int iData,
1926 long iRow,
1927 void *value)
1928{
1929 char *name = cData[iData].name;
1930 int type = cData[iData].type;
1931 int colnum = cData[iData].colnum;
1932 long nelem = cData[iData].nelem;
1933
1934 if (colnum > 0) {
1935 // Read the required number of values from the specified row of the table.
1936 int anynul;
1937 if (type == TSTRING) {
1938 if (nelem) {
1939 fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, &value, &anynul,
1940 &cStatus);
1941 } else {
1942 strcpy((char *)value, "");
1943 }
1944
1945 } else {
1946 if (nelem) {
1947 fits_read_col(cSDptr, type, colnum, iRow, 1, abs(nelem), 0, value,
1948 &anynul, &cStatus);
1949 } else {
1950 if (type == TSHORT) {
1951 *((short *)value) = 0;
1952 } else if (type == TINT) {
1953 *((int *)value) = 0;
1954 } else if (type == TFLOAT) {
1955 *((float *)value) = 0.0f;
1956 } else if (type == TDOUBLE) {
1957 *((double *)value) = 0.0;
1958 }
1959 }
1960 }
1961
1962 } else if (colnum == 0) {
1963 // Read keyword value.
1964 fits_read_key(cSDptr, type, name, value, 0, &cStatus);
1965
1966 } else {
1967 // Not present.
1968 if (type == TSTRING) {
1969 strcpy((char *)value, "");
1970 } else if (type == TSHORT) {
1971 *((short *)value) = 0;
1972 } else if (type == TINT) {
1973 *((int *)value) = 0;
1974 } else if (type == TFLOAT) {
1975 *((float *)value) = 0.0f;
1976 } else if (type == TDOUBLE) {
1977 *((double *)value) = 0.0;
1978 }
1979 }
1980
1981 return colnum < 0;
1982}
1983
1984//------------------------------------------------------ SDFITSreader::findCol
1985
1986// Locate a parameter in the SDFITS file.
1987
1988void SDFITSreader::findCol(
1989 char *name,
1990 int *colnum)
1991{
1992 *colnum = 0;
1993 int status = 0;
1994 fits_get_colnum(cSDptr, CASESEN, name, colnum, &status);
1995
1996 if (status) {
1997 // Not a real column - maybe it's virtual.
1998 char card[81];
1999
2000 status = 0;
2001 fits_read_card(cSDptr, name, card, &status);
2002 if (status) {
2003 // Not virtual either.
2004 *colnum = -1;
2005 }
2006
2007 // Clear error messages.
2008 fits_clear_errmsg();
2009 }
2010}
Note: See TracBrowser for help on using the repository browser.