source: trunk/external-alma/atnf/PKSIO/SDFITSreader.cc@ 2140

Last change on this file since 2140 was 1757, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

Module(s): all the CASA sd tools and tasks are affected.

Description: Merged ATNF-ASAP 3.0.0 developments to CASA (alma) branch.

Note you also need to update casa/code/atnf.
File size: 66.3 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDFITSreader.cc: ATNF interface class for SDFITS input using CFITSIO.
3//#---------------------------------------------------------------------------
4//# livedata - processing pipeline for single-dish, multibeam spectral data.
5//# Copyright (C) 2000-2009, Australia Telescope National Facility, CSIRO
6//#
7//# This file is part of livedata.
8//#
9//# livedata is free software: you can redistribute it and/or modify it under
10//# the terms of the GNU General Public License as published by the Free
11//# Software Foundation, either version 3 of the License, or (at your option)
12//# any later version.
13//#
14//# livedata is distributed in the hope that it will be useful, but WITHOUT
15//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
17//# more details.
18//#
19//# You should have received a copy of the GNU General Public License along
20//# with livedata. If not, see <http://www.gnu.org/licenses/>.
21//#
22//# Correspondence concerning livedata may be directed to:
23//# Internet email: mcalabre@atnf.csiro.au
24//# Postal address: Dr. Mark Calabretta
25//# Australia Telescope National Facility, CSIRO
26//# PO Box 76
27//# Epping NSW 1710
28//# AUSTRALIA
29//#
30//# http://www.atnf.csiro.au/computing/software/livedata.html
31//# $Id: SDFITSreader.cc,v 19.45 2009-09-30 07:23:48 cal103 Exp $
32//#---------------------------------------------------------------------------
33//# The SDFITSreader class reads single dish FITS files such as those written
34//# by SDFITSwriter containing Parkes Multibeam data.
35//#
36//# Original: 2000/08/09, Mark Calabretta, ATNF
37//#---------------------------------------------------------------------------
38
39#include <atnf/pks/pks_maths.h>
40#include <atnf/PKSIO/MBrecord.h>
41#include <atnf/PKSIO/SDFITSreader.h>
42
43#include <casa/Logging/LogIO.h>
44#include <casa/Quanta/MVTime.h>
45#include <casa/math.h>
46#include <casa/stdio.h>
47
48#include <algorithm>
49#include <strings.h>
50#include <cstring>
51
52class FITSparm
53{
54 public:
55 char *name; // Keyword or column name.
56 int type; // Expected keyvalue or column data type.
57 int colnum; // Column number; 0 for keyword; -1 absent.
58 int coltype; // Column data type, as found.
59 long nelem; // Column data repeat count; < 0 for vardim.
60 int tdimcol; // TDIM column number; 0 for keyword; -1 absent.
61 char units[32]; // Units from TUNITn keyword.
62};
63
64// Numerical constants.
65const double PI = 3.141592653589793238462643;
66
67// Factor to convert radians to degrees.
68const double D2R = PI / 180.0;
69
70// Class name
71const string className = "SDFITSreader" ;
72
73//---------------------------------------------------- SDFITSreader::(statics)
74
75int SDFITSreader::sInit = 1;
76int SDFITSreader::sReset = 0;
77int (*SDFITSreader::sALFAcalNon)[2] = (int (*)[2])(new float[16]);
78int (*SDFITSreader::sALFAcalNoff)[2] = (int (*)[2])(new float[16]);
79float (*SDFITSreader::sALFAcalOn)[2] = (float (*)[2])(new float[16]);
80float (*SDFITSreader::sALFAcalOff)[2] = (float (*)[2])(new float[16]);
81float (*SDFITSreader::sALFAcal)[2] = (float (*)[2])(new float[16]);
82
83//------------------------------------------------- SDFITSreader::SDFITSreader
84
85SDFITSreader::SDFITSreader()
86{
87 // Default constructor.
88 cSDptr = 0x0;
89
90 // Allocate space for data descriptors.
91 cData = new FITSparm[NDATA];
92
93 for (int iData = 0; iData < NDATA; iData++) {
94 cData[iData].colnum = -1;
95 }
96
97 // Initialize pointers.
98 cBeams = 0x0;
99 cIFs = 0x0;
100 cStartChan = 0x0;
101 cEndChan = 0x0;
102 cRefChan = 0x0;
103 cPols = 0x0;
104}
105
106//------------------------------------------------ SDFITSreader::~SDFITSreader
107
108SDFITSreader::~SDFITSreader()
109{
110 close();
111
112 delete [] cData;
113}
114
115//--------------------------------------------------------- SDFITSreader::open
116
117// Open an SDFITS file for reading.
118
119int SDFITSreader::open(
120 char* sdName,
121 int &nBeam,
122 int* &beams,
123 int &nIF,
124 int* &IFs,
125 int* &nChan,
126 int* &nPol,
127 int* &haveXPol,
128 int &haveBase,
129 int &haveSpectra,
130 int &extraSysCal)
131{
132 const string methodName = "open()" ;
133
134 if (cSDptr) {
135 close();
136 }
137
138 // Open the SDFITS file.
139 cStatus = 0;
140 if (fits_open_file(&cSDptr, sdName, READONLY, &cStatus)) {
141 sprintf(cMsg, "ERROR: Failed to open SDFITS file\n %s", sdName);
142 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, cMsg);
143 return 1;
144 }
145
146 // Move to the SDFITS extension.
147 cALFA = cALFA_BD = cALFA_CIMA = 0;
148 if (fits_movnam_hdu(cSDptr, BINARY_TBL, "SINGLE DISH", 0, &cStatus)) {
149 // No SDFITS table, look for BDFITS or CIMAFITS.
150 cStatus = 0;
151 if (fits_movnam_hdu(cSDptr, BINARY_TBL, "BDFITS", 0, &cStatus) == 0) {
152 cALFA_BD = 1;
153
154 } else {
155 cStatus = 0;
156 if (fits_movnam_hdu(cSDptr, BINARY_TBL, "CIMAFITS", 0, &cStatus) == 0) {
157 cALFA_CIMA = 1;
158
159 // Check for later versions of CIMAFITS.
160 float version;
161 readParm("VERSION", TFLOAT, &version);
162 if (version >= 2.0f) cALFA_CIMA = int(version);
163
164 } else {
165 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed to locate SDFITS binary table.");
166 close();
167 return 1;
168 }
169 }
170
171 // Arecibo ALFA data of some kind.
172 cALFA = 1;
173 if (sInit) {
174 for (int iBeam = 0; iBeam < 8; iBeam++) {
175 for (int iPol = 0; iPol < 2; iPol++) {
176 sALFAcalOn[iBeam][iPol] = 0.0f;
177 sALFAcalOff[iBeam][iPol] = 0.0f;
178
179 // Nominal factor to calibrate spectra in Jy.
180 sALFAcal[iBeam][iPol] = 3.0f;
181 }
182 }
183
184 sInit = 0;
185 }
186 }
187
188 // GBT data.
189 char telescope[32];
190 readParm("TELESCOP", TSTRING, telescope); // Core.
191 cGBT = strncmp(telescope, "GBT", 3) == 0 ||
192 strncmp(telescope, "NRAO_GBT", 8) == 0;
193
194
195 // Check that the DATA array column is present.
196 findData(DATA, "DATA", TFLOAT);
197 haveSpectra = cHaveSpectra = cData[DATA].colnum > 0;
198
199 cNAxisTime = 0;
200 if (cHaveSpectra) {
201 // Find the number of data axes (must be the same for each IF).
202 cNAxes = 5;
203 if (readDim(DATA, 1, &cNAxes, cNAxis)) {
204 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
205 close();
206 return 1;
207 }
208
209 if (cALFA_BD) {
210 // ALFA BDFITS: variable length arrays don't actually vary and there is
211 // no TDIM (or MAXISn) card; use the LAGS_IN value.
212 cNAxes = 5;
213 readParm("LAGS_IN", TLONG, cNAxis);
214 cNAxis[1] = 1;
215 cNAxis[2] = 1;
216 cNAxis[3] = 1;
217 cNAxis[4] = 1;
218 cData[DATA].nelem = cNAxis[0];
219 }
220
221 if (cNAxes < 4) {
222 // Need at least four axes (for now).
223 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "DATA array contains fewer than four axes.");
224 close();
225 return 1;
226 } else if (cNAxes > 5) {
227 // We support up to five axes.
228 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "DATA array contains more than five axes.");
229 close();
230 return 1;
231 }
232
233 findData(FLAGGED, "FLAGGED", TBYTE);
234
235 } else {
236 // DATA column not present, check for a DATAXED keyword.
237 findData(DATAXED, "DATAXED", TSTRING);
238 if (cData[DATAXED].colnum < 0) {
239 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "DATA array column absent from binary table.");
240 close();
241 return 1;
242 }
243
244 // Determine the number of axes and their length.
245 char dataxed[32];
246 readParm("DATAXED", TSTRING, dataxed);
247
248 for (int iaxis = 0; iaxis < 5; iaxis++) cNAxis[iaxis] = 0;
249 sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAxis, cNAxis+1, cNAxis+2,
250 cNAxis+3, cNAxis+4);
251 for (int iaxis = 4; iaxis > -1; iaxis--) {
252 if (cNAxis[iaxis] == 0) cNAxes = iaxis;
253 }
254 }
255
256 char *CTYPE[5] = {"CTYPE1", "CTYPE2", "CTYPE3", "CTYPE4", "CTYPE5"};
257 char *CRPIX[5] = {"CRPIX1", "CRPIX2", "CRPIX3", "CRPIX4", "CRPIX5"};
258 char *CRVAL[5] = {"CRVAL1", "CRVAL2", "CRVAL3", "CRVAL4", "CRVAL5"};
259 char *CDELT[5] = {"CDELT1", "CDELT2", "CDELT3", "CDELT4", "CDELT5"};
260
261 // Find required DATA array axes.
262 char ctype[5][72];
263 for (int iaxis = 0; iaxis < cNAxes; iaxis++) {
264 strcpy(ctype[iaxis], "");
265 readParm(CTYPE[iaxis], TSTRING, ctype[iaxis]); // Core.
266 }
267
268 if (cStatus) {
269 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
270 close();
271 return 1;
272 }
273
274 char *fqCRVAL = 0;
275 char *fqCDELT = 0;
276 char *fqCRPIX = 0;
277 char *raCRVAL = 0;
278 char *decCRVAL = 0;
279 char *timeCRVAL = 0;
280 char *timeCDELT = 0;
281 char *timeCRPIX = 0;
282 char *beamCRVAL = 0;
283 char *polCRVAL = 0;
284
285 cFreqAxis = -1;
286 cStokesAxis = -1;
287 cRaAxis = -1;
288 cDecAxis = -1;
289 cTimeAxis = -1;
290 cBeamAxis = -1;
291
292 for (int iaxis = 0; iaxis < cNAxes; iaxis++) {
293 if (strncmp(ctype[iaxis], "FREQ", 4) == 0) {
294 cFreqAxis = iaxis;
295 fqCRVAL = CRVAL[iaxis];
296 fqCDELT = CDELT[iaxis];
297 fqCRPIX = CRPIX[iaxis];
298
299 } else if (strncmp(ctype[iaxis], "STOKES", 6) == 0) {
300 cStokesAxis = iaxis;
301 polCRVAL = CRVAL[iaxis];
302
303 } else if (strncmp(ctype[iaxis], "RA", 2) == 0) {
304 cRaAxis = iaxis;
305 raCRVAL = CRVAL[iaxis];
306
307 } else if (strncmp(ctype[iaxis], "DEC", 3) == 0) {
308 cDecAxis = iaxis;
309 decCRVAL = CRVAL[iaxis];
310
311 } else if (strcmp(ctype[iaxis], "TIME") == 0) {
312 // TIME (UTC seconds since midnight); axis type, if present, takes
313 // precedence over keyword.
314 cTimeAxis = iaxis;
315 timeCRVAL = CRVAL[iaxis];
316
317 // Check for non-degeneracy.
318 if ((cNAxisTime = cNAxis[iaxis]) > 1) {
319 timeCDELT = CDELT[iaxis];
320 timeCRPIX = CRPIX[iaxis];
321 sprintf(cMsg, "DATA array contains a TIME axis of length %ld.",
322 cNAxisTime);
323 //logMsg(cMsg);
324 log(LogOrigin( className, methodName, WHERE ), LogIO::NORMAL, cMsg);
325 }
326
327 } else if (strcmp(ctype[iaxis], "BEAM") == 0) {
328 // BEAM can be a keyword or axis type.
329 cBeamAxis = iaxis;
330 beamCRVAL = CRVAL[iaxis];
331 }
332 }
333
334 if (cALFA_BD) {
335 // Fixed in ALFA CIMAFITS.
336 cRaAxis = 2;
337 raCRVAL = "CRVAL2A";
338
339 cDecAxis = 3;
340 decCRVAL = "CRVAL3A";
341 }
342
343
344 // Check that required axes are present.
345 if (cFreqAxis < 0 || cStokesAxis < 0 || cRaAxis < 0 || cDecAxis < 0) {
346 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Could not find required DATA array axes.");
347 close();
348 return 1;
349 }
350
351 // Set up machinery for data retrieval.
352 findData(SCAN, "SCAN", TINT); // Shared.
353 findData(CYCLE, "CYCLE", TINT); // Additional.
354 findData(DATE_OBS, "DATE-OBS", TSTRING); // Core.
355
356 if (cTimeAxis >= 0) {
357 // The DATA array has a TIME axis.
358 if (cNAxisTime > 1) {
359 // Non-degenerate.
360 findData(TimeRefVal, timeCRVAL, TDOUBLE); // Time reference value.
361 findData(TimeDelt, timeCDELT, TDOUBLE); // Time increment.
362 findData(TimeRefPix, timeCRPIX, TFLOAT); // Time reference pixel.
363 } else {
364 // Degenerate, treat its like a simple TIME keyword.
365 findData(TIME, timeCRVAL, TDOUBLE);
366 }
367
368 } else {
369 findData(TIME, "TIME", TDOUBLE); // Core.
370 }
371
372 findData(EXPOSURE, "EXPOSURE", TFLOAT); // Core.
373 findData(OBJECT, "OBJECT", TSTRING); // Core.
374 findData(OBJ_RA, "OBJ-RA", TDOUBLE); // Additional.
375 findData(OBJ_DEC, "OBJ-DEC", TDOUBLE); // Additional.
376 findData(RESTFRQ, "RESTFRQ", TDOUBLE); // Additional.
377 findData(OBSMODE, "OBSMODE", TSTRING); // Shared.
378
379 findData(BEAM, "BEAM", TSHORT); // Additional.
380 findData(IF, "IF", TSHORT); // Additional.
381 findData(FqRefVal, fqCRVAL, TDOUBLE); // Frequency reference value.
382 findData(FqDelt, fqCDELT, TDOUBLE); // Frequency increment.
383 findData(FqRefPix, fqCRPIX, TFLOAT); // Frequency reference pixel.
384 findData(RA, raCRVAL, TDOUBLE); // Right ascension.
385 findData(DEC, decCRVAL, TDOUBLE); // Declination.
386 findData(SCANRATE, "SCANRATE", TFLOAT); // Additional.
387
388 findData(TSYS, "TSYS", TFLOAT); // Core.
389 findData(CALFCTR, "CALFCTR", TFLOAT); // Additional.
390 findData(XCALFCTR, "XCALFCTR", TFLOAT); // Additional.
391 findData(BASELIN, "BASELIN", TFLOAT); // Additional.
392 findData(BASESUB, "BASESUB", TFLOAT); // Additional.
393 findData(XPOLDATA, "XPOLDATA", TFLOAT); // Additional.
394
395 findData(REFBEAM, "REFBEAM", TSHORT); // Additional.
396 findData(TCAL, "TCAL", TFLOAT); // Shared.
397 findData(TCALTIME, "TCALTIME", TSTRING); // Additional.
398 findData(AZIMUTH, "AZIMUTH", TFLOAT); // Shared.
399 findData(ELEVATIO, "ELEVATIO", TFLOAT); // Shared.
400 findData(PARANGLE, "PARANGLE", TFLOAT); // Additional.
401 findData(FOCUSAXI, "FOCUSAXI", TFLOAT); // Additional.
402 findData(FOCUSTAN, "FOCUSTAN", TFLOAT); // Additional.
403 findData(FOCUSROT, "FOCUSROT", TFLOAT); // Additional.
404 findData(TAMBIENT, "TAMBIENT", TFLOAT); // Shared.
405 findData(PRESSURE, "PRESSURE", TFLOAT); // Shared.
406 findData(HUMIDITY, "HUMIDITY", TFLOAT); // Shared.
407 findData(WINDSPEE, "WINDSPEE", TFLOAT); // Shared.
408 findData(WINDDIRE, "WINDDIRE", TFLOAT); // Shared.
409
410 findData(STOKES, polCRVAL, TINT);
411 findData(SIG, "SIG", TSTRING);
412 findData(CAL, "CAL", TSTRING);
413
414 findData(RVSYS, "RVSYS", TDOUBLE);
415 findData(VFRAME, "VFRAME", TDOUBLE);
416 findData(VELDEF, "VELDEF", TSTRING);
417
418 if (cStatus) {
419 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
420 close();
421 return 1;
422 }
423
424
425 // Check for alternative column names.
426 if (cALFA) {
427 // ALFA data.
428 cALFAscan = 0;
429 cScanNo = 0;
430 if (cALFA_CIMA) {
431 findData(SCAN, "SCAN_ID", TINT);
432 if (cALFA_CIMA > 1) {
433 // Note that RECNUM increases by cNAxisTime per row.
434 findData(CYCLE, "RECNUM", TINT);
435 } else {
436 findData(CYCLE, "SUBSCAN", TINT);
437 }
438 } else if (cALFA_BD) {
439 findData(SCAN, "SCAN_NUMBER", TINT);
440 findData(CYCLE, "PATTERN_NUMBER", TINT);
441 }
442 } else {
443 readData(SCAN, 1, &cFirstScanNo);
444 }
445
446 cCycleNo = 0;
447 cLastUTC = 0.0;
448 for ( int i = 0 ; i < 4 ; i++ ) {
449 cGLastUTC[i] = 0.0 ;
450 cGLastScan[i] = -1 ;
451 cGCycleNo[i] = 0 ;
452 }
453
454 // Beam number, 1-relative by default.
455 cBeam_1rel = 1;
456 if (cALFA) {
457 // ALFA INPUT_ID, 0-relative (overrides BEAM column if present).
458 findData(BEAM, "INPUT_ID", TSHORT);
459 cBeam_1rel = 0;
460
461 } else if (cData[BEAM].colnum < 0) {
462 if (beamCRVAL) {
463 // There is a BEAM axis.
464 findData(BEAM, beamCRVAL, TDOUBLE);
465 } else {
466 // ms2sdfits output, 0-relative "feed" number.
467 findData(BEAM, "MAIN_FEED1", TSHORT);
468 cBeam_1rel = 0;
469 }
470 }
471
472 // IF number, 1-relative by default.
473 cIF_1rel = 1;
474 if (cALFA && cData[IF].colnum < 0) {
475 // ALFA data, 0-relative.
476 if (cALFA_CIMA > 1) {
477 findData(IF, "IFN", TSHORT);
478 } else {
479 findData(IF, "IFVAL", TSHORT);
480 }
481 cIF_1rel = 0;
482 }
483
484 // ms2sdfits writes a scalar "TSYS" column that averages the polarizations.
485 int colnum;
486 findCol("SYSCAL_TSYS", &colnum);
487 if (colnum > 0) {
488 // This contains the vector Tsys.
489 findData(TSYS, "SYSCAL_TSYS", TFLOAT);
490 }
491
492 // XPOLDATA?
493
494 if (cData[SCANRATE].colnum < 0) {
495 findData(SCANRATE, "FIELD_POINTING_DIR_RATE", TFLOAT);
496 }
497
498 if (cData[RESTFRQ].colnum < 0) {
499 findData(RESTFRQ, "RESTFREQ", TDOUBLE);
500 if (cData[RESTFRQ].colnum < 0) {
501 findData(RESTFRQ, "SPECTRAL_WINDOW_REST_FREQUENCY", TDOUBLE);
502 }
503 }
504
505 if (cData[OBJ_RA].colnum < 0) {
506 findData(OBJ_RA, "SOURCE_DIRECTION", TDOUBLE);
507 }
508 if (cData[OBJ_DEC].colnum < 0) {
509 findData(OBJ_DEC, "SOURCE_DIRECTION", TDOUBLE);
510 }
511
512 // REFBEAM?
513
514 if (cData[TCAL].colnum < 0) {
515 findData(TCAL, "SYSCAL_TCAL", TFLOAT);
516 } else if (cALFA_BD) {
517 // ALFA BDFITS has a different TCAL with 64 elements - kill it!
518 findData(TCAL, "NO NO NO", TFLOAT);
519 }
520
521 if (cALFA_BD) {
522 // ALFA BDFITS.
523 findData(AZIMUTH, "CRVAL2B", TFLOAT);
524 findData(ELEVATIO, "CRVAL3B", TFLOAT);
525 }
526
527 if (cALFA) {
528 // ALFA data.
529 findData(PARANGLE, "PARA_ANG", TFLOAT);
530 }
531
532 if (cData[TAMBIENT].colnum < 0) {
533 findData(TAMBIENT, "WEATHER_TEMPERATURE", TFLOAT);
534 }
535
536 if (cData[PRESSURE].colnum < 0) {
537 findData(PRESSURE, "WEATHER_PRESSURE", TFLOAT);
538 }
539
540 if (cData[HUMIDITY].colnum < 0) {
541 findData(HUMIDITY, "WEATHER_REL_HUMIDITY", TFLOAT);
542 }
543
544 if (cData[WINDSPEE].colnum < 0) {
545 findData(WINDSPEE, "WEATHER_WIND_SPEED", TFLOAT);
546 }
547
548 if (cData[WINDDIRE].colnum < 0) {
549 findData(WINDDIRE, "WEATHER_WIND_DIRECTION", TFLOAT);
550 }
551
552
553 // Find the number of rows.
554 fits_get_num_rows(cSDptr, &cNRow, &cStatus);
555 if (!cNRow) {
556 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Table contains no entries.");
557 close();
558 return 1;
559 }
560
561
562 // Determine which beams are present in the data.
563 if (cData[BEAM].colnum > 0) {
564 short *beamCol = new short[cNRow];
565 short beamNul = 1;
566 int anynul;
567 if (fits_read_col(cSDptr, TSHORT, cData[BEAM].colnum, 1, 1, cNRow,
568 &beamNul, beamCol, &anynul, &cStatus)) {
569 delete [] beamCol;
570 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
571 close();
572 return 1;
573 }
574
575 // Find the maximum beam number.
576 cNBeam = cBeam_1rel - 1;
577 for (int irow = 0; irow < cNRow; irow++) {
578 if (beamCol[irow] > cNBeam) {
579 cNBeam = beamCol[irow];
580 }
581
582 // Check validity.
583 if (beamCol[irow] < cBeam_1rel) {
584 delete [] beamCol;
585 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "SDFITS file contains invalid beam number.");
586 close();
587 return 1;
588 }
589 }
590
591 if (!cBeam_1rel) cNBeam++;
592
593 // Find all beams present in the data.
594 cBeams = new int[cNBeam];
595 for (int ibeam = 0; ibeam < cNBeam; ibeam++) {
596 cBeams[ibeam] = 0;
597 }
598
599 for (int irow = 0; irow < cNRow; irow++) {
600 cBeams[beamCol[irow] - cBeam_1rel] = 1;
601 }
602
603 delete [] beamCol;
604
605 } else {
606 // No BEAM column.
607 cNBeam = 1;
608 cBeams = new int[1];
609 cBeams[0] = 1;
610 }
611
612 // Passing back the address of the array allows PKSFITSreader::select() to
613 // modify its elements directly.
614 nBeam = cNBeam;
615 beams = cBeams;
616
617
618 // Determine which IFs are present in the data.
619 if (cData[IF].colnum > 0) {
620 short *IFCol = new short[cNRow];
621 short IFNul = 1;
622 int anynul;
623 if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
624 &IFNul, IFCol, &anynul, &cStatus)) {
625 delete [] IFCol;
626 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
627 close();
628 return 1;
629 }
630
631 // Find the maximum IF number.
632 cNIF = cIF_1rel - 1;
633 for (int irow = 0; irow < cNRow; irow++) {
634 if (IFCol[irow] > cNIF) {
635 cNIF = IFCol[irow];
636 }
637
638 // Check validity.
639 if (IFCol[irow] < cIF_1rel) {
640 delete [] IFCol;
641 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "SDFITS file contains invalid IF number.");
642 close();
643 return 1;
644 }
645 }
646
647 if (!cIF_1rel) cNIF++;
648
649 // Find all IFs present in the data.
650 cIFs = new int[cNIF];
651 cNChan = new int[cNIF];
652 cNPol = new int[cNIF];
653 cHaveXPol = new int[cNIF];
654 cGetXPol = 0;
655
656 for (int iIF = 0; iIF < cNIF; iIF++) {
657 cIFs[iIF] = 0;
658 cNChan[iIF] = 0;
659 cNPol[iIF] = 0;
660 cHaveXPol[iIF] = 0;
661 }
662
663 for (int irow = 0; irow < cNRow; irow++) {
664 int iIF = IFCol[irow] - cIF_1rel;
665 if (cIFs[iIF] == 0) {
666 cIFs[iIF] = 1;
667
668 // Find the axis lengths.
669 if (cHaveSpectra) {
670 if (cData[DATA].nelem < 0) {
671 // Variable dimension array.
672 if (readDim(DATA, irow+1, &cNAxes, cNAxis)) {
673 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
674 close();
675 return 1;
676 }
677 }
678
679 } else {
680 if (cData[DATAXED].colnum > 0) {
681 char dataxed[32];
682 readParm("DATAXED", TSTRING, dataxed);
683
684 sscanf(dataxed, "(%ld,%ld,%ld,%ld,%ld)", cNAxis, cNAxis+1,
685 cNAxis+2, cNAxis+3, cNAxis+4);
686 }
687 }
688
689 // Number of channels and polarizations.
690 cNChan[iIF] = cNAxis[cFreqAxis];
691 cNPol[iIF] = cNAxis[cStokesAxis];
692 cHaveXPol[iIF] = 0;
693
694 // Is cross-polarization data present?
695 if (cData[XPOLDATA].colnum > 0) {
696 // Check that it conforms.
697 int nAxis;
698 long nAxes[2];
699
700 if (readDim(XPOLDATA, irow+1, &nAxis, nAxes)) {
701 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE );
702 close();
703 return 1;
704 }
705
706 // Default is to get it if we have it.
707 if (nAxis == 2 &&
708 nAxes[0] == 2 &&
709 nAxes[1] == cNChan[iIF]) {
710 cGetXPol = cHaveXPol[iIF] = 1;
711 }
712 }
713 }
714 }
715
716 delete [] IFCol;
717
718 } else {
719 // No IF column.
720 cNIF = 1;
721 cIFs = new int[1];
722 cIFs[0] = 1;
723
724 cNChan = new int[1];
725 cNPol = new int[1];
726 cHaveXPol = new int[1];
727 cGetXPol = 0;
728
729 // Number of channels and polarizations.
730 cNChan[0] = cNAxis[cFreqAxis];
731 cNPol[0] = cNAxis[cStokesAxis];
732 cHaveXPol[0] = 0;
733 }
734
735 if (cALFA && cALFA_CIMA < 2) {
736 // Older ALFA data labels each polarization as a separate IF.
737 cNPol[0] = cNIF;
738 cNIF = 1;
739 }
740
741 // For GBT data that stores spectra for each polarization in separate rows
742 if ( cData[STOKES].colnum > 0 ) {
743 int *stokesCol = new int[cNRow];
744 int stokesNul = 1;
745 int anynul;
746 if (fits_read_col(cSDptr, TINT, cData[STOKES].colnum, 1, 1, cNRow,
747 &stokesNul, stokesCol, &anynul, &cStatus)) {
748 delete [] stokesCol;
749 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
750 close();
751 return 1;
752 }
753
754 vector<int> pols ;
755 pols.push_back( stokesCol[0] ) ;
756 for ( int i = 0 ; i < cNRow ; i++ ) {
757 bool pmatch = false ;
758 for ( uint j = 0 ; j < pols.size() ; j++ ) {
759 if ( stokesCol[i] == pols[j] ) {
760 pmatch = true ;
761 break ;
762 }
763 }
764 if ( !pmatch ) {
765 pols.push_back( stokesCol[i] ) ;
766 }
767 }
768
769 cPols = new int[pols.size()] ;
770 for ( uint i = 0 ; i < pols.size() ; i++ ) {
771 cPols[i] = pols[i] ;
772 }
773
774 for ( int i = 0 ; i < cNIF ; i++ ) {
775 cNPol[i] = pols.size() ;
776 }
777
778 delete [] stokesCol ;
779 }
780
781 // Passing back the address of the array allows PKSFITSreader::select() to
782 // modify its elements directly.
783 nIF = cNIF;
784 IFs = cIFs;
785
786 nChan = cNChan;
787 nPol = cNPol;
788 haveXPol = cHaveXPol;
789
790
791 // Default channel range selection.
792 cStartChan = new int[cNIF];
793 cEndChan = new int[cNIF];
794 cRefChan = new int[cNIF];
795
796 for (int iIF = 0; iIF < cNIF; iIF++) {
797 cStartChan[iIF] = 1;
798 cEndChan[iIF] = cNChan[iIF];
799 cRefChan[iIF] = cNChan[iIF]/2 + 1;
800 }
801
802 // Default is to get it if we have it.
803 cGetSpectra = cHaveSpectra;
804
805
806 // Are baseline parameters present?
807 cHaveBase = 0;
808 if (cData[BASELIN].colnum) {
809 // Check that it conforms.
810 int nAxis, status = 0;
811 long nAxes[2];
812
813 if (fits_read_tdim(cSDptr, cData[BASELIN].colnum, 2, &nAxis, nAxes,
814 &status) == 0) {
815 cHaveBase = (nAxis == 2);
816 }
817 }
818 haveBase = cHaveBase;
819
820
821 // Is extra system calibration data available?
822 cExtraSysCal = 0;
823 for (int iparm = REFBEAM; iparm < NDATA; iparm++) {
824 if (cData[iparm].colnum >= 0) {
825 cExtraSysCal = 1;
826 break;
827 }
828 }
829
830 extraSysCal = cExtraSysCal;
831
832
833 // Extras for ALFA data.
834 cALFAacc = 0.0f;
835 if (cALFA_CIMA > 1) {
836 // FFTs per second when the Mock correlator operates in RFI blanking mode.
837 readData("PHFFTACC", TFLOAT, 0, &cALFAacc);
838 }
839
840
841 cRow = 0;
842 cTimeIdx = cNAxisTime;
843
844 return 0;
845}
846
847//---------------------------------------------------- SDFITSreader::getHeader
848
849// Get parameters describing the data.
850
851int SDFITSreader::getHeader(
852 char observer[32],
853 char project[32],
854 char telescope[32],
855 double antPos[3],
856 char obsMode[32],
857 char bunit[32],
858 float &equinox,
859 char radecsys[32],
860 char dopplerFrame[32],
861 char datobs[32],
862 double &utc,
863 double &refFreq,
864 double &bandwidth)
865{
866 const string methodName = "getHeader()" ;
867
868 // Has the file been opened?
869 if (!cSDptr) {
870 return 1;
871 }
872
873 // Read parameter values.
874 readParm("OBSERVER", TSTRING, observer); // Shared.
875 readParm("PROJID", TSTRING, project); // Shared.
876 readParm("TELESCOP", TSTRING, telescope); // Core.
877
878 antPos[0] = 0.0;
879 antPos[1] = 0.0;
880 antPos[2] = 0.0;
881 if (readParm("ANTENNA_POSITION", TDOUBLE, antPos)) {
882 readParm("OBSGEO-X", TDOUBLE, antPos); // Additional.
883 readParm("OBSGEO-Y", TDOUBLE, antPos + 1); // Additional.
884 readParm("OBSGEO-Z", TDOUBLE, antPos + 2); // Additional.
885 }
886
887 if (antPos[0] == 0.0) {
888 if (strncmp(telescope, "ATPKS", 5) == 0) {
889 // Parkes coordinates.
890 antPos[0] = -4554232.087;
891 antPos[1] = 2816759.046;
892 antPos[2] = -3454035.950;
893 } else if (strncmp(telescope, "ATMOPRA", 7) == 0) {
894 // Mopra coordinates.
895 antPos[0] = -4682768.630;
896 antPos[1] = 2802619.060;
897 antPos[2] = -3291759.900;
898 } else if (strncmp(telescope, "ARECIBO", 7) == 0) {
899 // Arecibo coordinates.
900 antPos[0] = 2390486.900;
901 antPos[1] = -5564731.440;
902 antPos[2] = 1994720.450;
903 }
904 }
905
906 readData(OBSMODE, 1, obsMode); // Shared.
907
908 // Brightness unit.
909 if (cData[DATAXED].colnum >= 0) {
910 strcpy(bunit, "Jy");
911 } else {
912 strcpy(bunit, cData[DATA].units);
913 }
914
915 if (strcmp(bunit, "JY") == 0) {
916 bunit[1] = 'y';
917 } else if (strcmp(bunit, "JY/BEAM") == 0) {
918 strcpy(bunit, "Jy/beam");
919 }
920
921 readParm("EQUINOX", TFLOAT, &equinox); // Shared.
922 if (cStatus == 405) {
923 // EQUINOX was written as string value in early versions.
924 cStatus = 0;
925 char strtmp[32];
926 readParm("EQUINOX", TSTRING, strtmp);
927 sscanf(strtmp, "%f", &equinox);
928 }
929
930 if (readParm("RADESYS", TSTRING, radecsys)) { // Additional.
931 if (readParm("RADECSYS", TSTRING, radecsys)) { // Additional.
932 strcpy(radecsys, "");
933 }
934 }
935
936 if (readParm("SPECSYS", TSTRING, dopplerFrame)) { // Additional.
937 // Fallback value.
938 strcpy(dopplerFrame, "TOPOCENT");
939
940 // Look for VELFRAME, written by earlier versions of Livedata.
941 //
942 // Added few more codes currently (as of 2009 Oct) used in the GBT
943 // SDFITS (based io_sdfits_define.pro of GBTIDL). - TT
944 if (readParm("VELFRAME", TSTRING, dopplerFrame)) { // Additional.
945 // No, try digging it out of the CTYPE card (AIPS convention).
946 char keyw[9], ctype[9];
947 sprintf(keyw, "CTYPE%ld", cFreqAxis+1);
948 readParm(keyw, TSTRING, ctype);
949
950 if (strncmp(ctype, "FREQ-", 5) == 0) {
951 strcpy(dopplerFrame, ctype+5);
952 if (strcmp(dopplerFrame, "LSR") == 0) {
953 // LSR unqualified usually means LSR (kinematic).
954 strcpy(dopplerFrame, "LSRK");
955 } else if (strcmp(dopplerFrame, "LSD") == 0) {
956 // LSR as a dynamical defintion
957 strcpy(dopplerFrame, "LSRD");
958 } else if (strcmp(dopplerFrame, "HEL") == 0) {
959 // Almost certainly barycentric.
960 strcpy(dopplerFrame, "BARYCENT");
961 } else if (strcmp(dopplerFrame, "BAR") == 0) {
962 // barycentric.
963 strcpy(dopplerFrame, "BARYCENT");
964 } else if (strcmp(dopplerFrame, "OBS") == 0) {
965 // observed or topocentric.
966 strcpy(dopplerFrame, "TOPO");
967 } else if (strcmp(dopplerFrame, "GEO") == 0) {
968 // geocentric
969 strcpy(dopplerFrame, "GEO");
970 } else if (strcmp(dopplerFrame, "GAL") == 0) {
971 // galactic
972 strcpy(dopplerFrame, "GAL");
973 } else if (strcmp(dopplerFrame, "LGR") == 0) {
974 // Local group
975 strcpy(dopplerFrame, "LGROUP");
976 } else if (strcmp(dopplerFrame, "CMB") == 0) {
977 // Cosimic Microwave Backgroup
978 strcpy(dopplerFrame, "CMB");
979 }
980 } else {
981 strcpy(dopplerFrame, "");
982 }
983 }
984 // Translate to FITS standard names.
985 if (strncmp(dopplerFrame, "TOP", 3) == 0) {
986 strcpy(dopplerFrame, "TOPOCENT");
987 } else if (strncmp(dopplerFrame, "GEO", 3) == 0) {
988 strcpy(dopplerFrame, "GEOCENTR");
989 } else if (strncmp(dopplerFrame, "HEL", 3) == 0) {
990 strcpy(dopplerFrame, "HELIOCEN");
991 } else if (strncmp(dopplerFrame, "BARY", 4) == 0) {
992 strcpy(dopplerFrame, "BARYCENT");
993 } else if (strncmp(dopplerFrame, "GAL", 3) == 0) {
994 strcpy(dopplerFrame, "GALACTOC");
995 } else if (strncmp(dopplerFrame, "LGROUP", 6) == 0) {
996 strcpy(dopplerFrame, "LOCALGRP");
997 } else if (strncmp(dopplerFrame, "CMB", 3) == 0) {
998 strcpy(dopplerFrame, "CMBDIPOL");
999 }
1000 }
1001
1002 if (cStatus) {
1003 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1004 return 1;
1005 }
1006
1007 // Get parameters from first row of table.
1008 readTime(1, 1, datobs, utc);
1009 readData(FqRefVal, 1, &refFreq);
1010 readParm("BANDWID", TDOUBLE, &bandwidth); // Core.
1011
1012 if (cStatus) {
1013 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1014 return 1;
1015 }
1016
1017 return 0;
1018}
1019
1020//-------------------------------------------------- SDFITSreader::getFreqInfo
1021
1022// Get frequency parameters for each IF.
1023
1024int SDFITSreader::getFreqInfo(
1025 int &nIF,
1026 double* &startFreq,
1027 double* &endFreq)
1028{
1029 const string methodName = "getFreqInfo()" ;
1030
1031 float fqRefPix;
1032 double fqDelt, fqRefVal;
1033
1034 nIF = cNIF;
1035 startFreq = new double[nIF];
1036 endFreq = new double[nIF];
1037
1038 if (cData[IF].colnum > 0) {
1039 short *IFCol = new short[cNRow];
1040 short IFNul = 1;
1041 int anynul;
1042 if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
1043 &IFNul, IFCol, &anynul, &cStatus)) {
1044 delete [] IFCol;
1045 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1046 close();
1047 return 1;
1048 }
1049
1050 for (int iIF = 0; iIF < nIF; iIF++) {
1051 if (cIFs[iIF]) {
1052 // Find the first occurrence of this IF in the table.
1053 int IFno = iIF + cIF_1rel;
1054 for (int irow = 0; irow < cNRow;) {
1055 if (IFCol[irow++] == IFno) {
1056 readData(FqRefPix, irow, &fqRefPix);
1057 readData(FqRefVal, irow, &fqRefVal);
1058 readData(FqDelt, irow, &fqDelt);
1059
1060 if (cALFA_BD) {
1061 unsigned char invert;
1062 readData("UPPERSB", TBYTE, irow, &invert);
1063
1064 if (invert) {
1065 fqDelt = -fqDelt;
1066 }
1067 }
1068
1069 startFreq[iIF] = fqRefVal + ( 1 - fqRefPix) * fqDelt;
1070 endFreq[iIF] = fqRefVal + (cNChan[iIF] - fqRefPix) * fqDelt;
1071
1072 break;
1073 }
1074 }
1075
1076 } else {
1077 startFreq[iIF] = 0.0;
1078 endFreq[iIF] = 0.0;
1079 }
1080 }
1081
1082 delete [] IFCol;
1083
1084 } else {
1085 // No IF column, read the first table entry.
1086 readData(FqRefPix, 1, &fqRefPix);
1087 readData(FqRefVal, 1, &fqRefVal);
1088 readData(FqDelt, 1, &fqDelt);
1089
1090 startFreq[0] = fqRefVal + ( 1 - fqRefPix) * fqDelt;
1091 endFreq[0] = fqRefVal + (cNChan[0] - fqRefPix) * fqDelt;
1092 }
1093
1094 return cStatus;
1095}
1096
1097//---------------------------------------------------- SDFITSreader::findRange
1098
1099// Find the range of the data in time and position.
1100
1101int SDFITSreader::findRange(
1102 int &nRow,
1103 int &nSel,
1104 char dateSpan[2][32],
1105 double utcSpan[2],
1106 double* &positions)
1107{
1108 const string methodName = "findRange()" ;
1109
1110 // Has the file been opened?
1111 if (!cSDptr) {
1112 return 1;
1113 }
1114
1115 nRow = cNRow;
1116
1117 // Find the number of rows selected.
1118 short *sel = new short[cNRow];
1119 for (int irow = 0; irow < cNRow; irow++) {
1120 sel[irow] = 1;
1121 }
1122
1123 int anynul;
1124 if (cData[BEAM].colnum > 0) {
1125 short *beamCol = new short[cNRow];
1126 short beamNul = 1;
1127 if (fits_read_col(cSDptr, TSHORT, cData[BEAM].colnum, 1, 1, cNRow,
1128 &beamNul, beamCol, &anynul, &cStatus)) {
1129 delete [] beamCol;
1130 delete [] sel;
1131 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1132 return 1;
1133 }
1134
1135 for (int irow = 0; irow < cNRow; irow++) {
1136 if (!cBeams[beamCol[irow]-cBeam_1rel]) {
1137 sel[irow] = 0;
1138 }
1139 }
1140
1141 delete [] beamCol;
1142 }
1143
1144 if (cData[IF].colnum > 0) {
1145 short *IFCol = new short[cNRow];
1146 short IFNul = 1;
1147 if (fits_read_col(cSDptr, TSHORT, cData[IF].colnum, 1, 1, cNRow,
1148 &IFNul, IFCol, &anynul, &cStatus)) {
1149 delete [] IFCol;
1150 delete [] sel;
1151 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1152 return 1;
1153 }
1154
1155 for (int irow = 0; irow < cNRow; irow++) {
1156 if (!cIFs[IFCol[irow]-cIF_1rel]) {
1157 sel[irow] = 0;
1158 }
1159 }
1160
1161 delete [] IFCol;
1162 }
1163
1164 nSel = 0;
1165 for (int irow = 0; irow < cNRow; irow++) {
1166 nSel += sel[irow];
1167 }
1168
1169
1170 // Find the time range assuming the data is in chronological order.
1171 readTime(1, 1, dateSpan[0], utcSpan[0]);
1172 readTime(cNRow, cNAxisTime, dateSpan[1], utcSpan[1]);
1173
1174
1175 // Retrieve positions for selected data.
1176 int isel = 0;
1177 positions = new double[2*nSel];
1178
1179 if (cCoordSys == 1) {
1180 // Horizontal (Az,El).
1181 if (cData[AZIMUTH].colnum < 0 ||
1182 cData[ELEVATIO].colnum < 0) {
1183 log(LogOrigin( className, methodName, WHERE ), LogIO::WARN, "Azimuth/elevation information absent.");
1184 cStatus = -1;
1185
1186 } else {
1187 float *az = new float[cNRow];
1188 float *el = new float[cNRow];
1189 readCol(AZIMUTH, az);
1190 readCol(ELEVATIO, el);
1191
1192 if (!cStatus) {
1193 for (int irow = 0; irow < cNRow; irow++) {
1194 if (sel[irow]) {
1195 positions[isel++] = az[irow] * D2R;
1196 positions[isel++] = el[irow] * D2R;
1197 }
1198 }
1199 }
1200
1201 delete [] az;
1202 delete [] el;
1203 }
1204
1205 } else if (cCoordSys == 3) {
1206 // ZPA-EL.
1207 if (cData[BEAM].colnum < 0 ||
1208 cData[FOCUSROT].colnum < 0 ||
1209 cData[ELEVATIO].colnum < 0) {
1210 log(LogOrigin( className, methodName, WHERE ), LogIO::WARN, "ZPA/elevation information absent.");
1211 cStatus = -1;
1212
1213 } else {
1214 short *beam = new short[cNRow];
1215 float *rot = new float[cNRow];
1216 float *el = new float[cNRow];
1217 readCol(BEAM, beam);
1218 readCol(FOCUSROT, rot);
1219 readCol(ELEVATIO, el);
1220
1221 if (!cStatus) {
1222 for (int irow = 0; irow < cNRow; irow++) {
1223 if (sel[irow]) {
1224 Int beamNo = beam[irow];
1225 Double zpa = rot[irow];
1226 if (beamNo > 1) {
1227 // Beam geometry for the Parkes multibeam.
1228 if (beamNo < 8) {
1229 zpa += -60.0 + 60.0*(beamNo-2);
1230 } else {
1231 zpa += -90.0 + 60.0*(beamNo-8);
1232 }
1233
1234 if (zpa < -180.0) {
1235 zpa += 360.0;
1236 } else if (zpa > 180.0) {
1237 zpa -= 360.0;
1238 }
1239 }
1240
1241 positions[isel++] = zpa * D2R;
1242 positions[isel++] = el[irow] * D2R;
1243 }
1244 }
1245 }
1246
1247 delete [] beam;
1248 delete [] rot;
1249 delete [] el;
1250 }
1251
1252 } else {
1253 double *ra = new double[cNRow];
1254 double *dec = new double[cNRow];
1255 readCol(RA, ra);
1256 readCol(DEC, dec);
1257
1258 if (cStatus) {
1259 delete [] ra;
1260 delete [] dec;
1261 goto cleanup;
1262 }
1263
1264 if (cALFA_BD) {
1265 for (int irow = 0; irow < cNRow; irow++) {
1266 // Convert hours to degrees.
1267 ra[irow] *= 15.0;
1268 }
1269 }
1270
1271 if (cCoordSys == 0) {
1272 // Equatorial (RA,Dec).
1273 for (int irow = 0; irow < cNRow; irow++) {
1274 if (sel[irow]) {
1275 positions[isel++] = ra[irow] * D2R;
1276 positions[isel++] = dec[irow] * D2R;
1277 }
1278 }
1279
1280 } else if (cCoordSys == 2) {
1281 // Feed-plane.
1282 if (cData[OBJ_RA].colnum < 0 ||
1283 cData[OBJ_DEC].colnum < 0 ||
1284 cData[PARANGLE].colnum < 0 ||
1285 cData[FOCUSROT].colnum < 0) {
1286 log( LogOrigin( className, methodName, WHERE ), LogIO::WARN,
1287 "Insufficient information to compute feed-plane\n"
1288 " coordinates.");
1289 cStatus = -1;
1290
1291 } else {
1292 double *srcRA = new double[cNRow];
1293 double *srcDec = new double[cNRow];
1294 float *par = new float[cNRow];
1295 float *rot = new float[cNRow];
1296
1297 readCol(OBJ_RA, srcRA);
1298 readCol(OBJ_DEC, srcDec);
1299 readCol(PARANGLE, par);
1300 readCol(FOCUSROT, rot);
1301
1302 if (!cStatus) {
1303 for (int irow = 0; irow < cNRow; irow++) {
1304 if (sel[irow]) {
1305 // Convert to feed-plane coordinates.
1306 Double dist, pa;
1307 distPA(ra[irow]*D2R, dec[irow]*D2R, srcRA[irow]*D2R,
1308 srcDec[irow]*D2R, dist, pa);
1309
1310 Double spin = (par[irow] + rot[irow])*D2R - pa;
1311 if (spin > 2.0*PI) spin -= 2.0*PI;
1312 Double squint = PI/2.0 - dist;
1313
1314 positions[isel++] = spin;
1315 positions[isel++] = squint;
1316 }
1317 }
1318 }
1319
1320 delete [] srcRA;
1321 delete [] srcDec;
1322 delete [] par;
1323 delete [] rot;
1324 }
1325 }
1326
1327 delete [] ra;
1328 delete [] dec;
1329 }
1330
1331cleanup:
1332 delete [] sel;
1333
1334 if (cStatus) {
1335 nSel = 0;
1336 delete [] positions;
1337 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1338 cStatus = 0;
1339 return 1;
1340 }
1341
1342 return 0;
1343}
1344
1345
1346//--------------------------------------------------------- SDFITSreader::read
1347
1348// Read the next data record.
1349
1350int SDFITSreader::read(
1351 MBrecord &mbrec)
1352{
1353 const string methodName = "read()" ;
1354
1355 // Has the file been opened?
1356 if (!cSDptr) {
1357 return 1;
1358 }
1359 // Find the next selected beam and IF.
1360 short iBeam = 0, iIF = 0;
1361 int iPol = -1 ;
1362 while (1) {
1363 if (++cTimeIdx > cNAxisTime) {
1364 if (++cRow > cNRow) break;
1365 cTimeIdx = 1;
1366 }
1367
1368 if (cData[BEAM].colnum > 0) {
1369 readData(BEAM, cRow, &iBeam);
1370
1371 // Convert to 0-relative.
1372 if (cBeam_1rel) iBeam--;
1373 }
1374
1375
1376 if (cBeams[iBeam]) {
1377 if (cData[IF].colnum > 0) {
1378 readData(IF, cRow, &iIF);
1379
1380 // Convert to 0-relative.
1381 if (cIF_1rel) iIF--;
1382 }
1383
1384 if (cIFs[iIF]) {
1385 if (cALFA) {
1386 // ALFA data, check for calibration data.
1387 char chars[32];
1388 readData(OBSMODE, cRow, chars);
1389 if (strcmp(chars, "DROP") == 0) {
1390 // Completely flagged integration.
1391 continue;
1392
1393 } else if (strcmp(chars, "CAL") == 0) {
1394 sReset = 1;
1395 if (cALFA_CIMA > 1) {
1396 for (short iPol = 0; iPol < cNPol[iIF]; iPol++) {
1397 alfaCal(iBeam, iIF, iPol);
1398 }
1399 continue;
1400 } else {
1401 // iIF is really the polarization in older ALFA data.
1402 alfaCal(iBeam, 0, iIF);
1403 continue;
1404 }
1405
1406 } else {
1407 // Reset for the next CAL record.
1408 if (sReset) {
1409 for (short iPol = 0; iPol < cNPol[iIF]; iPol++) {
1410 sALFAcalNon[iBeam][iPol] = 0;
1411 sALFAcalNoff[iBeam][iPol] = 0;
1412 sALFAcalOn[iBeam][iPol] = 0.0f;
1413 sALFAcalOff[iBeam][iPol] = 0.0f;
1414 }
1415 sReset = 0;
1416
1417 sprintf(cMsg, "ALFA cal factors for beam %d: %.3e, %.3e",
1418 iBeam+1, sALFAcal[iBeam][0], sALFAcal[iBeam][1]);
1419 log(LogOrigin( className, methodName, WHERE ), LogIO::NORMAL, cMsg);
1420 //logMsg(cMsg);
1421 }
1422 }
1423 }
1424
1425 // for GBT SDFITS
1426 if (cData[STOKES].colnum > 0 ) {
1427 readData(STOKES, cRow, &iPol ) ;
1428 for ( int i = 0 ; i < cNPol[iIF] ; i++ ) {
1429 if ( cPols[i] == iPol ) {
1430 iPol = i ;
1431 break ;
1432 }
1433 }
1434 }
1435 break;
1436 }
1437 }
1438 }
1439
1440 // EOF?
1441 if (cRow > cNRow) {
1442 return -1;
1443 }
1444
1445
1446 if (cALFA) {
1447 int scanNo;
1448 readData(SCAN, cRow, &scanNo);
1449 if (scanNo != cALFAscan) {
1450 cScanNo++;
1451 cALFAscan = scanNo;
1452 }
1453 mbrec.scanNo = cScanNo;
1454
1455 } else {
1456 readData(SCAN, cRow, &mbrec.scanNo);
1457
1458 // Ensure that scan number is 1-relative.
1459 mbrec.scanNo -= (cFirstScanNo - 1);
1460 }
1461
1462 // Times.
1463 char datobs[32];
1464 readTime(cRow, cTimeIdx, datobs, mbrec.utc);
1465 strcpy(mbrec.datobs, datobs);
1466
1467 if (cData[CYCLE].colnum > 0) {
1468 readData(CYCLE, cRow, &mbrec.cycleNo);
1469 mbrec.cycleNo += cTimeIdx - 1;
1470 if (cALFA_BD) mbrec.cycleNo++;
1471 } else {
1472 // Cycle number not recorded, must do our own bookkeeping.
1473 if (mbrec.utc != cLastUTC) {
1474 mbrec.cycleNo = ++cCycleNo;
1475 cLastUTC = mbrec.utc;
1476 }
1477 }
1478
1479 if ( iPol != -1 ) {
1480 if ( mbrec.scanNo != cGLastScan[iPol] ) {
1481 cGLastScan[iPol] = mbrec.scanNo ;
1482 cGCycleNo[iPol] = 0 ;
1483 mbrec.cycleNo = ++cGCycleNo[iPol] ;
1484 }
1485 else {
1486 mbrec.cycleNo = ++cGCycleNo[iPol] ;
1487 }
1488 }
1489
1490 readData(EXPOSURE, cRow, &mbrec.exposure);
1491
1492 // Source identification.
1493 readData(OBJECT, cRow, mbrec.srcName);
1494
1495 if ( iPol != -1 ) {
1496 char obsmode[32] ;
1497 readData( OBSMODE, cRow, obsmode ) ;
1498 char sig[1] ;
1499 char cal[1] ;
1500 readData( SIG, cRow, sig ) ;
1501 readData( CAL, cRow, cal ) ;
1502 if ( strstr( obsmode, "PSWITCH" ) != NULL ) {
1503 // position switch
1504 strcat( mbrec.srcName, "_p" ) ;
1505 if ( strstr( obsmode, "PSWITCHON" ) != NULL ) {
1506 strcat( mbrec.srcName, "s" ) ;
1507 }
1508 else if ( strstr( obsmode, "PSWITCHOFF" ) != NULL ) {
1509 strcat( mbrec.srcName, "r" ) ;
1510 }
1511 }
1512 else if ( strstr( obsmode, "Nod" ) != NULL ) {
1513 // nod
1514 strcat( mbrec.srcName, "_n" ) ;
1515 if ( sig[0] == 'T' ) {
1516 strcat( mbrec.srcName, "s" ) ;
1517 }
1518 else {
1519 strcat( mbrec.srcName, "r" ) ;
1520 }
1521 }
1522 else if ( strstr( obsmode, "FSWITCH" ) != NULL ) {
1523 // frequency switch
1524 strcat( mbrec.srcName, "_f" ) ;
1525 if ( sig[0] == 'T' ) {
1526 strcat( mbrec.srcName, "s" ) ;
1527 }
1528 else {
1529 strcat( mbrec.srcName, "r" ) ;
1530 }
1531 }
1532 if ( cal[0] == 'T' ) {
1533 strcat( mbrec.srcName, "c" ) ;
1534 }
1535 else {
1536 strcat( mbrec.srcName, "o" ) ;
1537 }
1538 }
1539
1540 readData(OBJ_RA, cRow, &mbrec.srcRA);
1541 if (strcmp(cData[OBJ_RA].name, "OBJ-RA") == 0) {
1542 mbrec.srcRA *= D2R;
1543 }
1544
1545 if (strcmp(cData[OBJ_DEC].name, "OBJ-DEC") == 0) {
1546 readData(OBJ_DEC, cRow, &mbrec.srcDec);
1547 mbrec.srcDec *= D2R;
1548 }
1549
1550 // Line rest frequency (Hz).
1551 readData(RESTFRQ, cRow, &mbrec.restFreq);
1552 if (mbrec.restFreq == 0.0 && cALFA_BD) {
1553 mbrec.restFreq = 1420.40575e6;
1554 }
1555
1556 // Observation mode.
1557 readData(OBSMODE, cRow, mbrec.obsType);
1558
1559 // Beam-dependent parameters.
1560 mbrec.beamNo = iBeam + 1;
1561
1562 readData(RA, cRow, &mbrec.ra);
1563 readData(DEC, cRow, &mbrec.dec);
1564 mbrec.ra *= D2R;
1565 mbrec.dec *= D2R;
1566
1567 if (cALFA_BD) mbrec.ra *= 15.0;
1568
1569 float scanrate[2];
1570 readData(SCANRATE, cRow, &scanrate);
1571 if (strcmp(cData[SCANRATE].name, "SCANRATE") == 0) {
1572 mbrec.raRate = scanrate[0] * D2R;
1573 mbrec.decRate = scanrate[1] * D2R;
1574 }
1575 mbrec.paRate = 0.0f;
1576
1577 // IF-dependent parameters.
1578 int startChan = cStartChan[iIF];
1579 int endChan = cEndChan[iIF];
1580 int refChan = cRefChan[iIF];
1581
1582 // Allocate data storage.
1583 int nChan = abs(endChan - startChan) + 1;
1584 int nPol = cNPol[iIF];
1585
1586 if ( cData[STOKES].colnum > 0 )
1587 nPol = 1 ;
1588
1589 if (cGetSpectra || cGetXPol) {
1590 int nxpol = cGetXPol ? 2*nChan : 0;
1591 mbrec.allocate(0, nChan*nPol, nxpol);
1592 }
1593
1594 mbrec.nIF = 1;
1595 mbrec.IFno[0] = iIF + 1;
1596 mbrec.nChan[0] = nChan;
1597 mbrec.nPol[0] = nPol;
1598 mbrec.polNo = iPol ;
1599
1600 readData(FqRefPix, cRow, mbrec.fqRefPix);
1601 readData(FqRefVal, cRow, mbrec.fqRefVal);
1602 readData(FqDelt, cRow, mbrec.fqDelt);
1603
1604 if (cALFA_BD) {
1605 unsigned char invert;
1606 int anynul, colnum;
1607 findCol("UPPERSB", &colnum);
1608 fits_read_col(cSDptr, TBYTE, colnum, cRow, 1, 1, 0, &invert, &anynul,
1609 &cStatus);
1610
1611 if (invert) {
1612 mbrec.fqDelt[0] = -mbrec.fqDelt[0];
1613 }
1614 }
1615
1616 if (cStatus) {
1617 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1618 return 1;
1619 }
1620
1621 // Adjust for channel selection.
1622 if (mbrec.fqRefPix[0] != refChan) {
1623 mbrec.fqRefVal[0] += (refChan - mbrec.fqRefPix[0]) * mbrec.fqDelt[0];
1624 mbrec.fqRefPix[0] = refChan;
1625 }
1626
1627 if (endChan < startChan) {
1628 mbrec.fqDelt[0] = -mbrec.fqDelt[0];
1629 }
1630
1631 // The data may only have a scalar Tsys value.
1632 mbrec.tsys[0][0] = 0.0f;
1633 mbrec.tsys[0][1] = 0.0f;
1634 if (cData[TSYS].nelem >= nPol) {
1635 readData(TSYS, cRow, mbrec.tsys[0]);
1636 }
1637
1638 for (int j = 0; j < 2; j++) {
1639 mbrec.calfctr[0][j] = 0.0f;
1640 }
1641 if (cData[CALFCTR].colnum > 0) {
1642 readData(CALFCTR, cRow, mbrec.calfctr);
1643 }
1644
1645 if (cHaveBase) {
1646 mbrec.haveBase = 1;
1647 readData(BASELIN, cRow, mbrec.baseLin);
1648 readData(BASESUB, cRow, mbrec.baseSub);
1649 } else {
1650 mbrec.haveBase = 0;
1651 }
1652
1653 if (cStatus) {
1654 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1655 return 1;
1656 }
1657
1658 // Read data, sectioning and transposing it in the process.
1659 long *blc = new long[cNAxes+1];
1660 long *trc = new long[cNAxes+1];
1661 long *inc = new long[cNAxes+1];
1662 for (int iaxis = 0; iaxis <= cNAxes; iaxis++) {
1663 blc[iaxis] = 1;
1664 trc[iaxis] = 1;
1665 inc[iaxis] = 1;
1666 }
1667
1668 blc[cFreqAxis] = std::min(startChan, endChan);
1669 trc[cFreqAxis] = std::max(startChan, endChan);
1670 if (cTimeAxis >= 0) {
1671 blc[cTimeAxis] = cTimeIdx;
1672 trc[cTimeAxis] = cTimeIdx;
1673 }
1674 blc[cNAxes] = cRow;
1675 trc[cNAxes] = cRow;
1676
1677 mbrec.haveSpectra = cGetSpectra;
1678 if (cGetSpectra) {
1679 int anynul;
1680
1681 for (int iPol = 0; iPol < nPol; iPol++) {
1682 blc[cStokesAxis] = iPol+1;
1683 trc[cStokesAxis] = iPol+1;
1684
1685 if (cALFA && cALFA_CIMA < 2) {
1686 // ALFA data: polarizations are stored in successive rows.
1687 blc[cStokesAxis] = 1;
1688 trc[cStokesAxis] = 1;
1689
1690 if (iPol) {
1691 if (++cRow > cNRow) {
1692 return -1;
1693 }
1694
1695 blc[cNAxes] = cRow;
1696 trc[cNAxes] = cRow;
1697 }
1698
1699 } else if (cData[DATA].nelem < 0) {
1700 // Variable dimension array; get axis lengths.
1701 int naxes = 5, status;
1702
1703 if ((status = readDim(DATA, cRow, &naxes, cNAxis))) {
1704 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1705
1706 } else if ((status = (naxes != cNAxes))) {
1707 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "DATA array dimensions changed.");
1708 }
1709
1710 if (status) {
1711 delete [] blc;
1712 delete [] trc;
1713 delete [] inc;
1714 return 1;
1715 }
1716 }
1717
1718 if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxes, cNAxis,
1719 blc, trc, inc, 0, mbrec.spectra[0] + iPol*nChan, &anynul,
1720 &cStatus)) {
1721 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1722 delete [] blc;
1723 delete [] trc;
1724 delete [] inc;
1725 return 1;
1726 }
1727
1728 if (endChan < startChan) {
1729 // Reverse the spectrum.
1730 float *iptr = mbrec.spectra[0] + iPol*nChan;
1731 float *jptr = iptr + nChan - 1;
1732 float *mid = iptr + nChan/2;
1733 while (iptr < mid) {
1734 float tmp = *iptr;
1735 *(iptr++) = *jptr;
1736 *(jptr--) = tmp;
1737 }
1738 }
1739
1740 if (cALFA) {
1741 // ALFA data, rescale the spectrum.
1742 float el, zd;
1743 readData(ELEVATIO, cRow, &el);
1744 zd = 90.0f - el;
1745
1746 float factor = sALFAcal[iBeam][iPol] / alfaGain(zd);
1747
1748 if (cALFA_CIMA > 1) {
1749 // Rescale according to the number of unblanked accumulations.
1750 int colnum, naccum;
1751 findCol("STAT", &colnum);
1752 fits_read_col(cSDptr, TINT, colnum, cRow, 10*(cTimeIdx-1)+2, 1, 0,
1753 &naccum, &anynul, &cStatus);
1754 factor *= cALFAacc / naccum;
1755 }
1756
1757 float *chan = mbrec.spectra[0] + iPol*nChan;
1758 float *chanN = chan + nChan;
1759 while (chan < chanN) {
1760 // Approximate conversion to Jy.
1761 *(chan++) *= factor;
1762 }
1763 }
1764
1765 if (mbrec.tsys[0][iPol] == 0.0) {
1766 // Compute Tsys as the average across the spectrum.
1767 float *chan = mbrec.spectra[0] + iPol*nChan;
1768 float *chanN = chan + nChan;
1769 float *tsys = mbrec.tsys[0] + iPol;
1770 while (chan < chanN) {
1771 *tsys += *(chan++);
1772 }
1773
1774 *tsys /= nChan;
1775 }
1776
1777 // Read data flags.
1778 if (cData[FLAGGED].colnum > 0) {
1779 if (fits_read_subset_byt(cSDptr, cData[FLAGGED].colnum, cNAxes,
1780 cNAxis, blc, trc, inc, 0, mbrec.flagged[0] + iPol*nChan, &anynul,
1781 &cStatus)) {
1782 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1783 delete [] blc;
1784 delete [] trc;
1785 delete [] inc;
1786 return 1;
1787 }
1788
1789 if (endChan < startChan) {
1790 // Reverse the flag vector.
1791 unsigned char *iptr = mbrec.flagged[0] + iPol*nChan;
1792 unsigned char *jptr = iptr + nChan - 1;
1793 for (int ichan = 0; ichan < nChan/2; ichan++) {
1794 unsigned char tmp = *iptr;
1795 *(iptr++) = *jptr;
1796 *(jptr--) = tmp;
1797 }
1798 }
1799
1800 } else {
1801 // All channels are unflagged by default.
1802 unsigned char *iptr = mbrec.flagged[0] + iPol*nChan;
1803 for (int ichan = 0; ichan < nChan; ichan++) {
1804 *(iptr++) = 0;
1805 }
1806 }
1807 }
1808 }
1809
1810
1811 // Read cross-polarization data.
1812 if (cGetXPol) {
1813 int anynul;
1814 for (int j = 0; j < 2; j++) {
1815 mbrec.xcalfctr[0][j] = 0.0f;
1816 }
1817 if (cData[XCALFCTR].colnum > 0) {
1818 readData(XCALFCTR, cRow, mbrec.xcalfctr);
1819 }
1820
1821 blc[0] = 1;
1822 trc[0] = 2;
1823 blc[1] = std::min(startChan, endChan);
1824 trc[1] = std::max(startChan, endChan);
1825 blc[2] = cRow;
1826 trc[2] = cRow;
1827
1828 int nAxis = 2;
1829 long nAxes[] = {2, nChan};
1830
1831 if (fits_read_subset_flt(cSDptr, cData[XPOLDATA].colnum, nAxis, nAxes,
1832 blc, trc, inc, 0, mbrec.xpol[0], &anynul, &cStatus)) {
1833 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1834 delete [] blc;
1835 delete [] trc;
1836 delete [] inc;
1837 return 1;
1838 }
1839
1840 if (endChan < startChan) {
1841 // Invert the cross-polarization spectrum.
1842 float *iptr = mbrec.xpol[0];
1843 float *jptr = iptr + nChan - 2;
1844 for (int ichan = 0; ichan < nChan/2; ichan++) {
1845 float tmp = *iptr;
1846 *iptr = *jptr;
1847 *jptr = tmp;
1848
1849 tmp = *(iptr+1);
1850 *(iptr+1) = *(jptr+1);
1851 *(jptr+1) = tmp;
1852
1853 iptr += 2;
1854 jptr -= 2;
1855 }
1856 }
1857 }
1858
1859 delete [] blc;
1860 delete [] trc;
1861 delete [] inc;
1862
1863 if (cStatus) {
1864 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1865 return 1;
1866 }
1867
1868 mbrec.extraSysCal = cExtraSysCal;
1869 readData(REFBEAM, cRow, &mbrec.refBeam);
1870 readData(TCAL, cRow, &mbrec.tcal[0]);
1871 readData(TCALTIME, cRow, mbrec.tcalTime);
1872
1873 readData(AZIMUTH, cRow, &mbrec.azimuth);
1874 readData(ELEVATIO, cRow, &mbrec.elevation);
1875 readData(PARANGLE, cRow, &mbrec.parAngle);
1876
1877 readData(FOCUSAXI, cRow, &mbrec.focusAxi);
1878 readData(FOCUSTAN, cRow, &mbrec.focusTan);
1879 readData(FOCUSROT, cRow, &mbrec.focusRot);
1880
1881 readData(TAMBIENT, cRow, &mbrec.temp);
1882 readData(PRESSURE, cRow, &mbrec.pressure);
1883 readData(HUMIDITY, cRow, &mbrec.humidity);
1884 readData(WINDSPEE, cRow, &mbrec.windSpeed);
1885 readData(WINDDIRE, cRow, &mbrec.windAz);
1886
1887 if (cALFA_BD) {
1888 // ALFA BDFITS stores zenith angle rather than elevation.
1889 mbrec.elevation = 90.0 - mbrec.elevation;
1890 }
1891
1892 mbrec.azimuth *= D2R;
1893 mbrec.elevation *= D2R;
1894 mbrec.parAngle *= D2R;
1895 mbrec.focusRot *= D2R;
1896 mbrec.windAz *= D2R;
1897
1898 // For GBT data, source velocity can be evaluated
1899 if ( cData[RVSYS].colnum > 0 && cData[VFRAME].colnum > 0 ) {
1900 float vframe;
1901 readData(VFRAME, cRow, &vframe);
1902 float rvsys;
1903 readData(RVSYS, cRow, &rvsys);
1904 //mbrec.srcVelocity = rvsys - vframe ;
1905 mbrec.srcVelocity = rvsys ;
1906 }
1907
1908 if (cStatus) {
1909 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
1910 return 1;
1911 }
1912
1913 return 0;
1914}
1915
1916//-------------------------------------------------------- SDFITSreader::close
1917
1918// Close the SDFITS file.
1919
1920void SDFITSreader::close()
1921{
1922 if (cSDptr) {
1923 int status = 0;
1924 fits_close_file(cSDptr, &status);
1925 cSDptr = 0x0;
1926
1927 if (cBeams) delete [] cBeams;
1928 if (cIFs) delete [] cIFs;
1929 if (cStartChan) delete [] cStartChan;
1930 if (cEndChan) delete [] cEndChan;
1931 if (cRefChan) delete [] cRefChan;
1932 }
1933}
1934
1935//------------------------------------------------------- SDFITSreader::log
1936
1937// Log a message. If the current CFITSIO status value is non-zero, also log
1938// the corresponding error message and the CFITSIO message stack.
1939
1940void SDFITSreader::log(LogOrigin origin, LogIO::Command cmd, const char *msg)
1941{
1942 LogIO os( origin ) ;
1943
1944 os << msg << endl ;
1945
1946 if (cStatus > 0) {
1947 fits_get_errstatus(cStatus, cMsg);
1948 os << cMsg << endl ;
1949
1950 while (fits_read_errmsg(cMsg)) {
1951 os << cMsg << endl ;
1952 }
1953 }
1954 os << LogIO::POST ;
1955}
1956
1957//----------------------------------------------------- SDFITSreader::findData
1958
1959// Locate a data item in the SDFITS file.
1960
1961void SDFITSreader::findData(
1962 int iData,
1963 char *name,
1964 int type)
1965{
1966 cData[iData].name = name;
1967 cData[iData].type = type;
1968
1969 int colnum;
1970 findCol(name, &colnum);
1971 cData[iData].colnum = colnum;
1972
1973 // Determine the number of data elements.
1974 if (colnum > 0) {
1975 int coltype;
1976 long nelem, width;
1977 fits_get_coltype(cSDptr, colnum, &coltype, &nelem, &width, &cStatus);
1978 fits_get_bcolparms(cSDptr, colnum, 0x0, cData[iData].units, 0x0, 0x0, 0x0,
1979 0x0, 0x0, 0x0, &cStatus);
1980
1981 // Look for a TDIMnnn keyword or column.
1982 char tdim[8];
1983 sprintf(tdim, "TDIM%d", colnum);
1984 findCol(tdim, &cData[iData].tdimcol);
1985
1986 if (coltype < 0) {
1987 // CFITSIO returns coltype < 0 for variable length arrays.
1988 cData[iData].coltype = -coltype;
1989 cData[iData].nelem = -nelem;
1990
1991 } else {
1992 cData[iData].coltype = coltype;
1993
1994 // Is there a TDIMnnn column?
1995 if (cData[iData].tdimcol > 0) {
1996 // Yes, dimensions of the fixed-length array could still vary.
1997 cData[iData].nelem = -nelem;
1998 } else {
1999 cData[iData].nelem = nelem;
2000 }
2001 }
2002
2003 } else if (colnum == 0) {
2004 // Keyword.
2005 cData[iData].coltype = 0;
2006 cData[iData].nelem = 1;
2007 cData[iData].tdimcol = -1;
2008 }
2009}
2010
2011//------------------------------------------------------ SDFITSreader::findCol
2012
2013// Locate a parameter in the SDFITS file.
2014
2015void SDFITSreader::findCol(
2016 char *name,
2017 int *colnum)
2018{
2019 *colnum = 0;
2020 int status = 0;
2021 fits_get_colnum(cSDptr, CASESEN, name, colnum, &status);
2022
2023 if (status) {
2024 // Not a real column - maybe it's virtual.
2025 char card[81];
2026
2027 status = 0;
2028 fits_read_card(cSDptr, name, card, &status);
2029 if (status) {
2030 // Not virtual either.
2031 *colnum = -1;
2032 }
2033
2034 // Clear error messages.
2035 fits_clear_errmsg();
2036 }
2037}
2038
2039//------------------------------------------------------ SDFITSreader::readDim
2040
2041// Determine the dimensions of an array in the SDFITS file.
2042
2043int SDFITSreader::readDim(
2044 int iData,
2045 long iRow,
2046 int *naxes,
2047 long naxis[])
2048{
2049 int colnum = cData[iData].colnum;
2050 if (colnum <= 0) {
2051 return 1;
2052 }
2053
2054 int maxdim = *naxes;
2055 if (cData[iData].tdimcol < 0) {
2056 // No TDIMnnn column for this array.
2057 if (cData[iData].nelem < 0) {
2058 // Variable length array; read the array descriptor.
2059 *naxes = 1;
2060 long dummy;
2061 if (fits_read_descript(cSDptr, colnum, iRow, naxis, &dummy, &cStatus)) {
2062 return 1;
2063 }
2064
2065 } else {
2066 // Read the repeat count from TFORMnnn.
2067 if (fits_read_tdim(cSDptr, colnum, maxdim, naxes, naxis, &cStatus)) {
2068 return 1;
2069 }
2070 }
2071
2072 } else {
2073 // Read the TDIMnnn value from the header or table.
2074 char tdim[8], tdimval[64];
2075 sprintf(tdim, "TDIM%d", colnum);
2076 readData(tdim, TSTRING, iRow, tdimval);
2077
2078 // fits_decode_tdim() checks that the TDIMnnn value is within the length
2079 // of the array in the specified column number but unfortunately doesn't
2080 // recognize variable-length arrays. Hence we must decode it here.
2081 char *tp = tdimval;
2082 if (*tp != '(') return 1;
2083
2084 tp++;
2085 *naxes = 0;
2086 for (size_t j = 1; j < strlen(tdimval); j++) {
2087 if (tdimval[j] == ',' || tdimval[j] == ')') {
2088 sscanf(tp, "%ld", naxis + (*naxes)++);
2089 if (tdimval[j] == ')') break;
2090 tp = tdimval + j + 1;
2091 }
2092 }
2093 }
2094
2095 return 0;
2096}
2097
2098//----------------------------------------------------- SDFITSreader::readParm
2099
2100// Read a parameter value from the SDFITS file.
2101
2102int SDFITSreader::readParm(
2103 char *name,
2104 int type,
2105 void *value)
2106{
2107 return readData(name, type, 1, value);
2108}
2109
2110//----------------------------------------------------- SDFITSreader::readData
2111
2112// Read a data value from the SDFITS file.
2113
2114int SDFITSreader::readData(
2115 char *name,
2116 int type,
2117 long iRow,
2118 void *value)
2119{
2120 int colnum;
2121 findCol(name, &colnum);
2122
2123 if (colnum > 0 && iRow > 0) {
2124 // Read the first value from the specified row of the table.
2125 int coltype;
2126 long nelem, width;
2127 fits_get_coltype(cSDptr, colnum, &coltype, &nelem, &width, &cStatus);
2128
2129 int anynul;
2130 if (type == TSTRING) {
2131 if (nelem) {
2132 fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, &value, &anynul,
2133 &cStatus);
2134 } else {
2135 strcpy((char *)value, "");
2136 }
2137
2138 } else {
2139 if (nelem) {
2140 fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, value, &anynul,
2141 &cStatus);
2142 } else {
2143 if (type == TSHORT) {
2144 *((short *)value) = 0;
2145 } else if (type == TINT) {
2146 *((int *)value) = 0;
2147 } else if (type == TFLOAT) {
2148 *((float *)value) = 0.0f;
2149 } else if (type == TDOUBLE) {
2150 *((double *)value) = 0.0;
2151 }
2152 }
2153 }
2154
2155 } else if (colnum == 0) {
2156 // Read keyword value.
2157 fits_read_key(cSDptr, type, name, value, 0, &cStatus);
2158
2159 } else {
2160 // Not present.
2161 if (type == TSTRING) {
2162 strcpy((char *)value, "");
2163 } else if (type == TSHORT) {
2164 *((short *)value) = 0;
2165 } else if (type == TINT) {
2166 *((int *)value) = 0;
2167 } else if (type == TFLOAT) {
2168 *((float *)value) = 0.0f;
2169 } else if (type == TDOUBLE) {
2170 *((double *)value) = 0.0;
2171 }
2172 }
2173
2174 return colnum < 0;
2175}
2176
2177//----------------------------------------------------- SDFITSreader::readData
2178
2179// Read data from the SDFITS file.
2180
2181int SDFITSreader::readData(
2182 int iData,
2183 long iRow,
2184 void *value)
2185{
2186 int type = cData[iData].type;
2187 int colnum = cData[iData].colnum;
2188
2189 if (colnum > 0 && iRow > 0) {
2190 // Read the required number of values from the specified row of the table.
2191 long nelem = cData[iData].nelem;
2192 int anynul;
2193 if (type == TSTRING) {
2194 if (nelem) {
2195 fits_read_col(cSDptr, type, colnum, iRow, 1, 1, 0, &value, &anynul,
2196 &cStatus);
2197 } else {
2198 strcpy((char *)value, "");
2199 }
2200
2201 } else {
2202 if (nelem) {
2203 fits_read_col(cSDptr, type, colnum, iRow, 1, abs(nelem), 0, value,
2204 &anynul, &cStatus);
2205 } else {
2206 if (type == TSHORT) {
2207 *((short *)value) = 0;
2208 } else if (type == TINT) {
2209 *((int *)value) = 0;
2210 } else if (type == TFLOAT) {
2211 *((float *)value) = 0.0f;
2212 } else if (type == TDOUBLE) {
2213 *((double *)value) = 0.0;
2214 }
2215 }
2216 }
2217
2218 } else if (colnum == 0) {
2219 // Read keyword value.
2220 char *name = cData[iData].name;
2221 fits_read_key(cSDptr, type, name, value, 0, &cStatus);
2222
2223 } else {
2224 // Not present.
2225 if (type == TSTRING) {
2226 strcpy((char *)value, "");
2227 } else if (type == TSHORT) {
2228 *((short *)value) = 0;
2229 } else if (type == TINT) {
2230 *((int *)value) = 0;
2231 } else if (type == TFLOAT) {
2232 *((float *)value) = 0.0f;
2233 } else if (type == TDOUBLE) {
2234 *((double *)value) = 0.0;
2235 }
2236 }
2237
2238 return colnum < 0;
2239}
2240
2241//------------------------------------------------------ SDFITSreader::readCol
2242
2243// Read a scalar column from the SDFITS file.
2244
2245int SDFITSreader::readCol(
2246 int iData,
2247 void *value)
2248{
2249 int type = cData[iData].type;
2250
2251 if (cData[iData].colnum > 0) {
2252 // Table column.
2253 int anynul;
2254 fits_read_col(cSDptr, type, cData[iData].colnum, 1, 1, cNRow, 0,
2255 value, &anynul, &cStatus);
2256
2257 } else {
2258 // Header keyword.
2259 readData(iData, 0, value);
2260 for (int irow = 1; irow < cNRow; irow++) {
2261 if (type == TSHORT) {
2262 ((short *)value)[irow] = *((short *)value);
2263 } else if (type == TINT) {
2264 ((int *)value)[irow] = *((int *)value);
2265 } else if (type == TFLOAT) {
2266 ((float *)value)[irow] = *((float *)value);
2267 } else if (type == TDOUBLE) {
2268 ((double *)value)[irow] = *((double *)value);
2269 }
2270 }
2271 }
2272
2273 return cData[iData].colnum < 0;
2274}
2275
2276//----------------------------------------------------- SDFITSreader::readTime
2277
2278// Read the time from the SDFITS file.
2279
2280int SDFITSreader::readTime(
2281 long iRow,
2282 int iPix,
2283 char *datobs,
2284 double &utc)
2285{
2286 readData(DATE_OBS, iRow, datobs);
2287 if (cData[TIME].colnum >= 0) {
2288 readData(TIME, iRow, &utc);
2289 } else if (cGBT) {
2290 Int yy, mm ;
2291 Double dd, hour, min, sec ;
2292 sscanf( datobs, "%d-%d-%lfT%lf:%lf:%lf", &yy, &mm, &dd, &hour, &min, &sec ) ;
2293 dd = dd + ( hour * 3600.0 + min * 60.0 + sec ) / 86400.0 ;
2294 MVTime mvt( yy, mm, dd ) ;
2295 dd = mvt.day() ;
2296 utc = fmod( dd, 1.0 ) * 86400.0 ;
2297 } else if (cNAxisTime > 1) {
2298 double timeDelt, timeRefPix, timeRefVal;
2299 readData(TimeRefVal, iRow, &timeRefVal);
2300 readData(TimeDelt, iRow, &timeDelt);
2301 readData(TimeRefPix, iRow, &timeRefPix);
2302 utc = timeRefVal + (iPix - timeRefPix) * timeDelt;
2303 }
2304
2305 if (cALFA_BD) utc *= 3600.0;
2306
2307 // Check DATE-OBS format.
2308 if (datobs[2] == '/') {
2309 // Translate an old-format DATE-OBS.
2310 datobs[9] = datobs[1];
2311 datobs[8] = datobs[0];
2312 datobs[2] = datobs[6];
2313 datobs[5] = datobs[3];
2314 datobs[3] = datobs[7];
2315 datobs[6] = datobs[4];
2316 datobs[7] = '-';
2317 datobs[4] = '-';
2318 datobs[1] = '9';
2319 datobs[0] = '1';
2320
2321 } else if (datobs[10] == 'T' && cData[TIME].colnum < 0) {
2322 // Dig UTC out of a new-format DATE-OBS.
2323 int hh, mm;
2324 float ss;
2325 sscanf(datobs+11, "%d:%d:%f", &hh, &mm, &ss);
2326 utc = (hh*60 + mm)*60 + ss;
2327 }
2328
2329 datobs[10] = '\0';
2330
2331 return 0;
2332}
2333
2334//------------------------------------------------------ SDFITSreader::alfaCal
2335
2336// Process ALFA calibration data.
2337
2338int SDFITSreader::alfaCal(
2339 short iBeam,
2340 short iIF,
2341 short iPol)
2342{
2343 const string methodName = "alfaCal()" ;
2344
2345 int calOn;
2346 char chars[32];
2347 if (cALFA_BD) {
2348 readData("OBS_NAME", TSTRING, cRow, chars);
2349 } else {
2350 readData("SCANTYPE", TSTRING, cRow, chars);
2351 }
2352
2353 if (strcmp(chars, "ON") == 0) {
2354 calOn = 1;
2355 } else if (strcmp(chars, "OFF") == 0) {
2356 calOn = 0;
2357 } else {
2358 return 1;
2359 }
2360
2361 // Read cal data.
2362 long *blc = new long[cNAxes+1];
2363 long *trc = new long[cNAxes+1];
2364 long *inc = new long[cNAxes+1];
2365 for (int iaxis = 0; iaxis <= cNAxes; iaxis++) {
2366 blc[iaxis] = 1;
2367 trc[iaxis] = 1;
2368 inc[iaxis] = 1;
2369 }
2370
2371 // User channel selection.
2372 int startChan = cStartChan[iIF];
2373 int endChan = cEndChan[iIF];
2374
2375 blc[cFreqAxis] = std::min(startChan, endChan);
2376 trc[cFreqAxis] = std::max(startChan, endChan);
2377 if (cALFA_CIMA > 1) {
2378 // CIMAFITS 2.x has a legitimate STOKES axis...
2379 blc[cStokesAxis] = iPol+1;
2380 trc[cStokesAxis] = iPol+1;
2381 } else {
2382 // ...older ALFA data does not.
2383 blc[cStokesAxis] = 1;
2384 trc[cStokesAxis] = 1;
2385 }
2386 if (cTimeAxis >= 0) {
2387 blc[cTimeAxis] = cTimeIdx;
2388 trc[cTimeAxis] = cTimeIdx;
2389 }
2390 blc[cNAxes] = cRow;
2391 trc[cNAxes] = cRow;
2392
2393 float spectrum[endChan];
2394 int anynul;
2395 if (fits_read_subset_flt(cSDptr, cData[DATA].colnum, cNAxes, cNAxis,
2396 blc, trc, inc, 0, spectrum, &anynul, &cStatus)) {
2397 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE);
2398 delete [] blc;
2399 delete [] trc;
2400 delete [] inc;
2401 return 1;
2402 }
2403
2404 // Factor to rescale according to the number of unblanked accumulations.
2405 float factor = 1.0f;
2406 if (cALFA_CIMA > 1) {
2407 int colnum, naccum;
2408 findCol("STAT", &colnum);
2409 fits_read_col(cSDptr, TINT, colnum, cRow, 2, 1, 0, &naccum, &anynul,
2410 &cStatus);
2411 factor = cALFAacc / naccum;
2412 }
2413
2414 // Average the spectrum.
2415 float mean = 1e9f;
2416 for (int k = 0; k < 2; k++) {
2417 float discrim = 2.0f * mean;
2418
2419 int nChan = 0;
2420 float sum = 0.0f;
2421
2422 float *chanN = spectrum + abs(endChan - startChan) + 1;
2423 for (float *chan = spectrum; chan < chanN; chan++) {
2424 // Simple discriminant that eliminates strong radar interference.
2425 if (*chan < discrim) {
2426 nChan++;
2427 sum += *chan * factor;
2428 }
2429 }
2430
2431 mean = sum / nChan;
2432 }
2433
2434 if (calOn) {
2435 sALFAcalOn[iBeam][iPol] *= sALFAcalNon[iBeam][iPol];
2436 sALFAcalOn[iBeam][iPol] += mean;
2437 sALFAcalOn[iBeam][iPol] /= ++sALFAcalNon[iBeam][iPol];
2438 } else {
2439 sALFAcalOff[iBeam][iPol] *= sALFAcalNoff[iBeam][iPol];
2440 sALFAcalOff[iBeam][iPol] += mean;
2441 sALFAcalOff[iBeam][iPol] /= ++sALFAcalNoff[iBeam][iPol];
2442 }
2443
2444 if (sALFAcalNon[iBeam][iPol] && sALFAcalNoff[iBeam][iPol]) {
2445 // Tcal should come from the TCAL table, it varies weakly with beam,
2446 // polarization, and frequency. However, TCAL is not written properly.
2447 float Tcal = 12.0f;
2448 sALFAcal[iBeam][iPol] = Tcal / (sALFAcalOn[iBeam][iPol] -
2449 sALFAcalOff[iBeam][iPol]);
2450
2451 // Scale from K to Jy; the gain also varies weakly with beam,
2452 // polarization, frequency, and zenith angle.
2453 float fluxCal = 10.0f;
2454 sALFAcal[iBeam][iPol] /= fluxCal;
2455 }
2456
2457 return 0;
2458}
2459
2460//----------------------------------------------------- SDFITSreader::alfaGain
2461
2462// ALFA gain factor.
2463
2464float SDFITSreader::alfaGain(
2465 float zd)
2466{
2467 // Gain vs zenith distance table from Robert Minchin, 2008/12/08.
2468 const int nZD = 37;
2469 const float zdLim[] = {1.5f, 19.5f};
2470 const float zdInc = (nZD - 1) / (zdLim[1] - zdLim[0]);
2471 float zdGain[] = { 1.00723708,
2472 1.16644573, 1.15003645, 1.07117307, 1.02532673,
2473 1.01788402, 1.01369524, 1.00000000, 0.989855111,
2474 0.990888834, 0.993996620, 0.989964068, 0.982213855,
2475 0.978662670, 0.979349494, 0.978478372, 0.974631131,
2476 0.972126007, 0.972835243, 0.972742677, 0.968671739,
2477 0.963891327, 0.963452935, 0.966831207, 0.969585896,
2478 0.970700860, 0.972644389, 0.973754644, 0.967344403,
2479 0.952168941, 0.937160134, 0.927843094, 0.914048433,
2480 0.886700928, 0.864701211, 0.869126320, 0.854309499};
2481
2482 float gain;
2483 // Do table lookup by linear interpolation.
2484 float lambda = zdInc * (zd - zdLim[0]);
2485 int j = int(lambda);
2486 if (j < 0) {
2487 gain = zdGain[0];
2488 } else if (j >= nZD-1) {
2489 gain = zdGain[nZD-1];
2490 } else {
2491 gain = zdGain[j] + (lambda - j) * (zdGain[j+1] - zdGain[j]);
2492 }
2493
2494 return gain;
2495}
2496
Note: See TracBrowser for help on using the repository browser.