source: trunk/external/atnf/PKSIO/SDFITSwriter.cc@ 1722

Last change on this file since 1722 was 1720, checked in by Malte Marquarding, 15 years ago

Update from livedata CVS repository
File size: 25.8 KB
Line 
1//#---------------------------------------------------------------------------
2//# SDFITSwriter.cc: ATNF CFITSIO interface class for SDFITS output.
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: SDFITSwriter.cc,v 19.18 2009-09-29 07:33:39 cal103 Exp $
32//#---------------------------------------------------------------------------
33//# Original: 2000/07/24, Mark Calabretta, ATNF
34//#---------------------------------------------------------------------------
35
36#include <atnf/PKSIO/MBrecord.h>
37#include <atnf/PKSIO/PKSmsg.h>
38#include <atnf/PKSIO/SDFITSwriter.h>
39
40#include <casa/iostream.h>
41
42#include <algorithm>
43#include <math.h>
44
45using namespace std;
46
47// Numerical constants.
48const double PI = 3.141592653589793238462643;
49
50// Factor to convert radians to degrees.
51const double R2D = 180.0 / PI;
52
53//------------------------------------------------- SDFITSwriter::SDFITSwriter
54
55SDFITSwriter::SDFITSwriter()
56{
57 // Default constructor.
58 cSDptr = 0x0;
59
60 // By default, messages are written to stderr.
61 initMsg();
62}
63
64//------------------------------------------------ SDFITSwriter::~SDFITSwriter
65
66SDFITSwriter::~SDFITSwriter()
67{
68 close();
69}
70
71//------------------------------------------------------- SDFITSwriter::create
72
73// Create the output SDFITS file.
74
75int SDFITSwriter::create(
76 char* sdName,
77 char* observer,
78 char* project,
79 char* telescope,
80 double antPos[3],
81 char* obsMode,
82 char* bunit,
83 float equinox,
84 char* dopplerFrame,
85 int nIF,
86 int* nChan,
87 int* nPol,
88 int* haveXPol,
89 int haveBase,
90 int extraSysCal)
91{
92 if (cSDptr) {
93 logMsg("ERROR: Output file already open, close it first.");
94 return 1;
95 }
96
97 // Clear the message stack.
98 clearMsg();
99
100 // Prepend an '!' to the output name to force it to be overwritten.
101 char sdname[80];
102 sdname[0] = '!';
103 strcpy(sdname+1, sdName);
104
105 // Create a new SDFITS file.
106 cStatus = 0;
107 if (fits_create_file(&cSDptr, sdname, &cStatus)) {
108 sprintf(cMsg, "ERROR: Failed to create SDFITS file\n %s", sdName);
109 logMsg(cMsg);
110 return cStatus;
111 }
112
113 cIsMX = strstr(obsMode, "MX") != 0;
114 cNIF = nIF;
115 cNChan = nChan;
116 cNPol = nPol;
117 cHaveXPol = haveXPol;
118 cHaveBase = haveBase;
119 cExtraSysCal = extraSysCal;
120
121 // Do all IFs have the same number of channels and polarizations?
122 cDoTDIM = 0;
123 int nprod = cNChan[0] * cNPol[0];
124 for (int iIF = 0; iIF < nIF; iIF++) {
125 if (cNChan[iIF]*cNPol[iIF] != nprod) {
126 // Need variable-length arrays as well as a TDIM column.
127 cDoTDIM = 2;
128 break;
129 }
130
131 if (cNChan[iIF] != cNChan[0] || cNPol[iIF] != cNPol[0]) {
132 // Varying channels and/or polarizations, need a TDIM column at least.
133 cDoTDIM = 1;
134 }
135 }
136
137 // Find the maximum number of polarizations in any IF.
138 int maxNPol = 0;
139 for (int iIF = 0; iIF < nIF; iIF++) {
140 if (cNPol[iIF] > maxNPol) maxNPol = cNPol[iIF];
141 }
142
143 // Do any IFs have cross-polarizations?
144 cDoXPol = 0;
145 for (int iIF = 0; iIF < nIF; iIF++) {
146 if (cHaveXPol[iIF]) {
147 cDoXPol = 1;
148 break;
149 }
150 }
151
152
153 cRow = 0;
154
155 // Write required primary header keywords.
156 if (fits_write_imghdr(cSDptr, 8, 0, 0, &cStatus)) {
157 logMsg("ERROR: Failed to write required primary header keywords.");
158 return cStatus;
159 }
160
161 // Identify the origin of the data.
162 fits_write_comment(cSDptr, " ", &cStatus);
163 fits_write_comment(cSDptr,
164 "This single dish FITS (SDFITS) file has a binary table extension which",
165 &cStatus);
166 fits_write_comment(cSDptr,
167 "contains data obtained from a telescope run by the Australia Telescope",
168 &cStatus);
169 fits_write_comment(cSDptr, "National Facility (ATNF).", &cStatus);
170 fits_write_comment(cSDptr, " ", &cStatus);
171
172 fits_write_date(cSDptr, &cStatus);
173
174 char text[72];
175 char version[7];
176 char date[11];
177 sscanf("$Revision: 19.18 $", "%*s%s", version);
178 sscanf("$Date: 2009-09-29 07:33:39 $", "%*s%s", date);
179 sprintf(text, "SDFITSwriter (v%s, %s)", version, date);
180 fits_write_key_str(cSDptr, "ORIGIN", text, "output class", &cStatus);
181
182 float cfvers;
183 fits_write_comment(cSDptr, "Written by Mark Calabretta "
184 "(mcalabre@atnf.csiro.au)", &cStatus);
185 sprintf(text, "using cfitsio v%.3f.", fits_get_version(&cfvers));
186 fits_write_comment(cSDptr, text, &cStatus);
187
188 if (cStatus) {
189 logMsg("ERROR: Failed in writing primary header.");
190 return cStatus;
191 }
192
193
194 // Create an SDFITS extension.
195 long nrow = 0;
196 int ncol = 0;
197 if (fits_create_tbl(cSDptr, BINARY_TBL, nrow, ncol, NULL, NULL, NULL,
198 "SINGLE DISH", &cStatus)) {
199 logMsg("ERROR: Failed to create a binary table extension.");
200 return 1;
201 }
202
203 char ttype[16];
204 char tform[9];
205 char tunit[9];
206
207 // NMATRIX (core, virtual).
208 fits_write_key_lng(cSDptr, "NMATRIX", 1l, "Number of DATA arrays",
209 &cStatus);
210
211 // OBSERVER (shared, virtual).
212 fits_write_key_str(cSDptr, "OBSERVER", observer, "Observer name(s)",
213 &cStatus);
214
215 // PROJID (shared, virtual).
216 fits_write_key_str(cSDptr, "PROJID", project, "Project name", &cStatus);
217
218 // TELESCOP (core, virtual).
219 fits_write_key_str(cSDptr, "TELESCOP", telescope, "Telescope name",
220 &cStatus);
221
222 // OBSGEO-X/Y/Z (additional, virtual).
223 fits_write_key_dbl(cSDptr, "OBSGEO-X", antPos[0], 9,
224 "[m] Antenna ITRF X-coordinate", &cStatus);
225 fits_write_key_dbl(cSDptr, "OBSGEO-Y", antPos[1], 9,
226 "[m] Antenna ITRF Y-coordinate", &cStatus);
227 fits_write_key_dbl(cSDptr, "OBSGEO-Z", antPos[2], 9,
228 "[m] Antenna ITRF Z-coordinate", &cStatus);
229
230 // SCAN (shared, real).
231 fits_insert_col(cSDptr, ++ncol, "SCAN", "1I", &cStatus);
232
233 // CYCLE (additional, real).
234 fits_insert_col(cSDptr, ++ncol, "CYCLE", "1J", &cStatus);
235
236 // DATE-OBS (core, real).
237 fits_insert_col(cSDptr, ++ncol, "DATE-OBS", "10A", &cStatus);
238
239 // TIME (core, real).
240 fits_insert_col(cSDptr, ++ncol, "TIME", "1D", &cStatus);
241 sprintf(tunit, "TUNIT%d", ncol);
242 fits_write_key_str(cSDptr, tunit, "s", "units of field", &cStatus);
243
244 // EXPOSURE (core, real).
245 fits_insert_col(cSDptr, ++ncol, "EXPOSURE", "1E", &cStatus);
246 sprintf(tunit, "TUNIT%d", ncol);
247 fits_write_key_str(cSDptr, tunit, "s", "units of field", &cStatus);
248
249 // OBJECT (core, real).
250 fits_insert_col(cSDptr, ++ncol, "OBJECT", "16A", &cStatus);
251
252 // OBJ-RA (additional, real).
253 fits_insert_col(cSDptr, ++ncol, "OBJ-RA", "1D", &cStatus);
254 sprintf(tunit, "TUNIT%d", ncol);
255 fits_write_key_str(cSDptr, tunit, "deg", "units of field", &cStatus);
256
257 // OBJ-DEC (additional, real).
258 fits_insert_col(cSDptr, ++ncol, "OBJ-DEC", "1D", &cStatus);
259 sprintf(tunit, "TUNIT%d", ncol);
260 fits_write_key_str(cSDptr, tunit, "deg", "units of field", &cStatus);
261
262 // RESTFRQ (additional, real).
263 fits_insert_col(cSDptr, ++ncol, "RESTFRQ", "1D", &cStatus);
264 sprintf(tunit, "TUNIT%d", ncol);
265 fits_write_key_str(cSDptr, tunit, "Hz", "units of field", &cStatus);
266
267 // OBSMODE (shared, real).
268 fits_insert_col(cSDptr, ++ncol, "OBSMODE", "16A", &cStatus);
269
270 // BEAM (additional, real).
271 fits_insert_col(cSDptr, ++ncol, "BEAM", "1I", &cStatus);
272
273 // IF (additional, real).
274 fits_insert_col(cSDptr, ++ncol, "IF", "1I", &cStatus);
275
276 // FREQRES (core, real).
277 fits_insert_col(cSDptr, ++ncol, "FREQRES", "1D", &cStatus);
278 sprintf(tunit, "TUNIT%d", ncol);
279 fits_write_key_str(cSDptr, tunit, "Hz", "units of field", &cStatus);
280
281 // BANDWID (core, real).
282 fits_insert_col(cSDptr, ++ncol, "BANDWID", "1D", &cStatus);
283 sprintf(tunit, "TUNIT%d", ncol);
284 fits_write_key_str(cSDptr, tunit, "Hz", "units of field", &cStatus);
285
286 // CTYPE1 (core, virtual).
287 fits_write_key_str(cSDptr, "CTYPE1", "FREQ",
288 "DATA array axis 1: frequency in Hz.", &cStatus);
289
290 // CRPIX1 (core, real).
291 fits_insert_col(cSDptr, ++ncol, "CRPIX1", "1E", &cStatus);
292
293 // CRVAL1 (core, real).
294 fits_insert_col(cSDptr, ++ncol, "CRVAL1", "1D", &cStatus);
295 sprintf(tunit, "TUNIT%d", ncol);
296 fits_write_key_str(cSDptr, tunit, "Hz", "units of field", &cStatus);
297
298 // CDELT1 (core, real).
299 fits_insert_col(cSDptr, ++ncol, "CDELT1", "1D", &cStatus);
300 sprintf(tunit, "TUNIT%d", ncol);
301 fits_write_key_str(cSDptr, tunit, "Hz", "units of field", &cStatus);
302
303
304 // CTYPE2 (core, virtual).
305 fits_write_key_str(cSDptr, "CTYPE2", "STOKES",
306 "DATA array axis 2: polarization code", &cStatus);
307
308 // CRPIX2 (core, virtual).
309 fits_write_key_flt(cSDptr, "CRPIX2", 1.0f, 1,
310 "Polarization code reference pixel", &cStatus);
311
312 // CRVAL2 (core, virtual).
313 fits_write_key_dbl(cSDptr, "CRVAL2", -5.0, 1,
314 "Polarization code at reference pixel (XX)", &cStatus);
315
316 // CDELT2 (core, virtual).
317 fits_write_key_dbl(cSDptr, "CDELT2", -1.0, 1,
318 "Polarization code axis increment", &cStatus);
319
320
321 // CTYPE3 (core, virtual).
322 fits_write_key_str(cSDptr, "CTYPE3", "RA",
323 "DATA array axis 3 (degenerate): RA (mid-int)",
324 &cStatus);
325
326 // CRPIX3 (core, virtual).
327 fits_write_key_flt(cSDptr, "CRPIX3", 1.0f, 1, "RA reference pixel",
328 &cStatus);
329
330 // CRVAL3 (core, real).
331 fits_insert_col(cSDptr, ++ncol, "CRVAL3", "1D", &cStatus);
332 sprintf(tunit, "TUNIT%d", ncol);
333 fits_write_key_str(cSDptr, tunit, "deg", "units of field", &cStatus);
334
335 // CDELT3 (core, virtual).
336 fits_write_key_dbl(cSDptr, "CDELT3", -1.0, 1, "RA axis increment", &cStatus);
337
338 // CTYPE4 (core, virtual).
339 fits_write_key_str(cSDptr, "CTYPE4", "DEC",
340 "DATA array axis 4 (degenerate): Dec (mid-int)",
341 &cStatus);
342
343 // CRPIX4 (core, virtual).
344 fits_write_key_flt(cSDptr, "CRPIX4", 1.0f, 1, "Dec reference pixel",
345 &cStatus);
346
347 // CRVAL4 (core, real).
348 fits_insert_col(cSDptr, ++ncol, "CRVAL4", "1D", &cStatus);
349 sprintf(tunit, "TUNIT%d", ncol);
350 fits_write_key_str(cSDptr, tunit, "deg", "units of field", &cStatus);
351
352 // CDELT4 (core, virtual).
353 fits_write_key_dbl(cSDptr, "CDELT4", 1.0, 1, "Dec axis increment", &cStatus);
354
355 // SCANRATE (additional, real).
356 fits_insert_col(cSDptr, ++ncol, "SCANRATE", "2E", &cStatus);
357 sprintf(tunit, "TUNIT%d", ncol);
358 fits_write_key_str(cSDptr, tunit, "deg/s", "units of field", &cStatus);
359
360 // SPECSYS (additional, virtual).
361 fits_write_key_str(cSDptr, "SPECSYS", dopplerFrame,
362 "Doppler reference frame (transformed)", &cStatus);
363
364 // SSYSOBS (additional, virtual).
365 fits_write_key_str(cSDptr, "SSYSOBS", "TOPOCENT",
366 "Doppler reference frame of observation", &cStatus);
367
368 // EQUINOX (shared, virtual).
369 fits_write_key_flt(cSDptr, "EQUINOX", equinox, 1,
370 "Equinox of equatorial coordinates", &cStatus);
371
372 // RADESYS (additional, virtual).
373 fits_write_key_str(cSDptr, "RADESYS", "FK5", "Equatorial coordinate frame",
374 &cStatus);
375
376 // TSYS (core, real).
377 sprintf(tform, "%dE", maxNPol);
378 fits_insert_col(cSDptr, ++ncol, "TSYS", tform, &cStatus);
379 sprintf(tunit, "TUNIT%d", ncol);
380 fits_write_key_str(cSDptr, tunit, bunit, "units of field", &cStatus);
381
382 // CALFCTR (additional, real).
383 sprintf(tform, "%dE", maxNPol);
384 fits_insert_col(cSDptr, ++ncol, "CALFCTR", tform, &cStatus);
385
386 if (cHaveBase) {
387 // BASELIN (additional, real).
388 sprintf(tform, "%dE", 2*maxNPol);
389 fits_insert_col(cSDptr, ++ncol, "BASELIN", tform, &cStatus);
390 long tdim[] = {2, maxNPol};
391 fits_write_tdim(cSDptr, ncol, 2, tdim, &cStatus);
392
393 // BASESUB (additional, real).
394 sprintf(tform, "%dE", 24*maxNPol);
395 fits_insert_col(cSDptr, ++ncol, "BASESUB", tform, &cStatus);
396 tdim[0] = 24;
397 fits_write_tdim(cSDptr, ncol, 2, tdim, &cStatus);
398 }
399
400 // DATA (core, real).
401 if (cDoTDIM < 2) {
402 // IFs all have the same number of products, use fixed-length arrays.
403 sprintf(tform, "%dE", cNChan[0]*cNPol[0]);
404 } else {
405 // IFs have a differing number of products, use variable-length arrays.
406 strcpy(tform, "1PE");
407 }
408 fits_insert_col(cSDptr, ++ncol, "DATA", tform, &cStatus);
409
410 if (cDoTDIM) {
411 // TDIMn varies with IF, write a TDIM column.
412 sprintf(ttype, "TDIM%d", ncol);
413 fits_insert_col(cSDptr, ++ncol, ttype, "16A", &cStatus);
414 } else {
415 // TDIMn fixed for each IF, write a TDIM keyword.
416 long tdim[] = {cNChan[0], cNPol[0], 1, 1};
417 fits_write_tdim(cSDptr, ncol, 4, tdim, &cStatus);
418 }
419
420 sprintf(tunit, "TUNIT%d", ncol);
421 fits_write_key_str(cSDptr, tunit, bunit, "units of field", &cStatus);
422
423 // FLAGGED (additional, logical).
424 if (cDoTDIM < 2) {
425 // IFs all have the same number of products, use fixed-length arrays.
426 sprintf(tform, "%dB", cNChan[0]*cNPol[0]);
427 } else {
428 // IFs have a differing number of products, use variable-length arrays.
429 strcpy(tform, "1PB");
430 }
431 fits_insert_col(cSDptr, ++ncol, "FLAGGED", tform, &cStatus);
432
433 if (cDoTDIM) {
434 // TDIMn varies with IF, write a TDIM column.
435 sprintf(ttype, "TDIM%d", ncol);
436 fits_insert_col(cSDptr, ++ncol, ttype, "16A", &cStatus);
437 } else {
438 // TDIMn fixed for each IF, write a TDIM keyword.
439 long tdim[] = {cNChan[0], cNPol[0], 1, 1};
440 fits_write_tdim(cSDptr, ncol, 4, tdim, &cStatus);
441 }
442
443 if (cDoXPol) {
444 // XCALFCTR (additional, real).
445 sprintf(tform, "%dE", 2);
446 fits_insert_col(cSDptr, ++ncol, "XCALFCTR", tform, &cStatus);
447
448 // XPOLDATA (additional, real).
449 if (cDoTDIM < 2) {
450 // IFs all have the same number of products, use fixed-length arrays.
451 sprintf(tform, "%dE", 2*cNChan[0]);
452 } else {
453 // IFs have a differing number of products, use variable-length arrays.
454 strcpy(tform, "1PE");
455 }
456 fits_insert_col(cSDptr, ++ncol, "XPOLDATA", tform, &cStatus);
457
458 if (cDoTDIM) {
459 // TDIMn varies with IF, write a TDIM column.
460 sprintf(ttype, "TDIM%d", ncol);
461 fits_insert_col(cSDptr, ++ncol, ttype, "16A", &cStatus);
462 } else {
463 // TDIMn fixed for each IF, write a TDIM keyword.
464 long tdim[] = {2, cNChan[0]};
465 fits_write_tdim(cSDptr, ncol, 2, tdim, &cStatus);
466 }
467
468 sprintf(tunit, "TUNIT%d", ncol);
469 fits_write_key_str(cSDptr, tunit, bunit, "units of field", &cStatus);
470 }
471
472 if (cExtraSysCal) {
473 if (cIsMX) {
474 // REFBEAM (additional, real).
475 fits_insert_col(cSDptr, ++ncol, "REFBEAM", "1I", &cStatus);
476 }
477
478 // TCAL (shared, real).
479 sprintf(tform, "%dE", min(maxNPol,2));
480 fits_insert_col(cSDptr, ++ncol, "TCAL", tform, &cStatus);
481 sprintf(tunit, "TUNIT%d", ncol);
482 fits_write_key_str(cSDptr, tunit, "Jy", "units of field", &cStatus);
483
484 // TCALTIME (additional, real).
485 fits_insert_col(cSDptr, ++ncol, "TCALTIME", "16A", &cStatus);
486
487 // AZIMUTH (shared, real).
488 fits_insert_col(cSDptr, ++ncol, "AZIMUTH", "1E", &cStatus);
489 sprintf(tunit, "TUNIT%d", ncol);
490 fits_write_key_str(cSDptr, tunit, "deg", "units of field", &cStatus);
491
492 // ELEVATIO (shared, real).
493 fits_insert_col(cSDptr, ++ncol, "ELEVATIO", "1E", &cStatus);
494 sprintf(tunit, "TUNIT%d", ncol);
495 fits_write_key_str(cSDptr, tunit, "deg", "units of field", &cStatus);
496
497 // PARANGLE (additional, real).
498 fits_insert_col(cSDptr, ++ncol, "PARANGLE", "1E", &cStatus);
499 sprintf(tunit, "TUNIT%d", ncol);
500 fits_write_key_str(cSDptr, tunit, "deg", "units of field", &cStatus);
501
502 // FOCUSAXI (additional, real).
503 fits_insert_col(cSDptr, ++ncol, "FOCUSAXI", "1E", &cStatus);
504 sprintf(tunit, "TUNIT%d", ncol);
505 fits_write_key_str(cSDptr, tunit, "m", "units of field", &cStatus);
506
507 // FOCUSTAN (additional, real).
508 fits_insert_col(cSDptr, ++ncol, "FOCUSTAN", "1E", &cStatus);
509 sprintf(tunit, "TUNIT%d", ncol);
510 fits_write_key_str(cSDptr, tunit, "m", "units of field", &cStatus);
511
512 // FOCUSROT (additional, real).
513 fits_insert_col(cSDptr, ++ncol, "FOCUSROT", "1E", &cStatus);
514 sprintf(tunit, "TUNIT%d", ncol);
515 fits_write_key_str(cSDptr, tunit, "deg", "units of field", &cStatus);
516
517 // TAMBIENT (shared, real).
518 fits_insert_col(cSDptr, ++ncol, "TAMBIENT", "1E", &cStatus);
519 sprintf(tunit, "TUNIT%d", ncol);
520 fits_write_key_str(cSDptr, tunit, "C", "units of field", &cStatus);
521
522 // PRESSURE (shared, real).
523 fits_insert_col(cSDptr, ++ncol, "PRESSURE", "1E", &cStatus);
524 sprintf(tunit, "TUNIT%d", ncol);
525 fits_write_key_str(cSDptr, tunit, "Pa", "units of field", &cStatus);
526
527 // HUMIDITY (shared, real).
528 fits_insert_col(cSDptr, ++ncol, "HUMIDITY", "1E", &cStatus);
529 sprintf(tunit, "TUNIT%d", ncol);
530 fits_write_key_str(cSDptr, tunit, "%", "units of field", &cStatus);
531
532 // WINDSPEE (shared, real).
533 fits_insert_col(cSDptr, ++ncol, "WINDSPEE", "1E", &cStatus);
534 sprintf(tunit, "TUNIT%d", ncol);
535 fits_write_key_str(cSDptr, tunit, "m/s", "units of field", &cStatus);
536
537 // WINDDIRE (shared, real).
538 fits_insert_col(cSDptr, ++ncol, "WINDDIRE", "1E", &cStatus);
539 sprintf(tunit, "TUNIT%d", ncol);
540 fits_write_key_str(cSDptr, tunit, "deg", "units of field", &cStatus);
541 }
542
543 // Set scaling parameters.
544 for (int j = 1; j <= ncol; j++) {
545 fits_set_tscale(cSDptr, j, 1.0, 0.0, &cStatus);
546 }
547
548 if (cStatus) {
549 logMsg("ERROR: Failed in writing binary table header.");
550 }
551
552 return cStatus;
553}
554
555//-------------------------------------------------------- SDFITSwriter::write
556
557// Write a record to the SDFITS file.
558
559int SDFITSwriter::write(MBrecord &mbrec)
560{
561 char *cptr;
562
563 // Check parameters.
564 int IFno = mbrec.IFno[0];
565 if (IFno < 1 || cNIF < IFno) {
566 cerr << "SDFITSwriter::write: "
567 << "Invalid IF number " << IFno
568 << " (maximum " << cNIF << ")." << endl;
569 return 1;
570 }
571
572 int iIF = IFno - 1;
573 int nChan = cNChan[iIF];
574 if (mbrec.nChan[0] != nChan) {
575 cerr << "SDFITSriter::write: "
576 << "Wrong number of channels for IF " << IFno << "," << endl
577 << " "
578 << "got " << nChan << " should be " << mbrec.nChan[0] << "." << endl;
579 return 1;
580 }
581
582 int nPol = cNPol[iIF];
583 if (mbrec.nPol[0] != nPol) {
584 cerr << "SDFITSriter::write: "
585 << "Wrong number of polarizations for IF " << IFno << "," << endl
586 << " "
587 << "got " << nPol << " should be " << mbrec.nPol[0] << "." << endl;
588 return 1;
589 }
590
591
592 // Next row.
593 cRow++;
594
595 int icol = 0;
596
597 // SCAN.
598 fits_write_col_int(cSDptr, ++icol, cRow, 1, 1, &mbrec.scanNo, &cStatus);
599
600 // CYCLE.
601 fits_write_col_int(cSDptr, ++icol, cRow, 1, 1, &mbrec.cycleNo, &cStatus);
602
603 // DATE-OBS.
604 cptr = mbrec.datobs;
605 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
606
607 // TIME.
608 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &mbrec.utc, &cStatus);
609
610 // EXPOSURE.
611 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.exposure, &cStatus);
612
613 // OBJECT.
614 cptr = mbrec.srcName;
615 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
616
617 // OBJ-RA.
618 double srcRA = mbrec.srcRA * R2D;
619 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &srcRA, &cStatus);
620
621 // OBJ-DEC.
622 double srcDec = mbrec.srcDec * R2D;
623 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &srcDec, &cStatus);
624
625 // RESTFRQ.
626 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &mbrec.restFreq, &cStatus);
627
628 // OBJECT.
629 cptr = mbrec.obsType;
630 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
631
632 // BEAM.
633 fits_write_col_sht(cSDptr, ++icol, cRow, 1, 1, &mbrec.beamNo, &cStatus);
634
635 // IF.
636 fits_write_col_sht(cSDptr, ++icol, cRow, 1, 1, &mbrec.IFno[0], &cStatus);
637
638 // FREQRES.
639 double freqRes = fabs(mbrec.fqDelt[0]);
640 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &freqRes, &cStatus);
641
642 // BANDWID.
643 double bandwidth = freqRes * nChan;
644 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &bandwidth, &cStatus);
645
646 // CRPIX1.
647 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.fqRefPix[0],
648 &cStatus);
649
650 // CRVAL1.
651 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &mbrec.fqRefVal[0],
652 &cStatus);
653
654 // CDELT1.
655 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &mbrec.fqDelt[0], &cStatus);
656
657 // CRVAL3.
658 double ra = mbrec.ra * R2D;
659 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &ra, &cStatus);
660
661 // CRVAL4.
662 double dec = mbrec.dec * R2D;
663 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &dec, &cStatus);
664
665 // SCANRATE.
666 float scanrate[2];
667 scanrate[0] = mbrec.raRate * R2D;
668 scanrate[1] = mbrec.decRate * R2D;
669 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 2, scanrate, &cStatus);
670
671 // TSYS.
672 fits_write_col_flt(cSDptr, ++icol, cRow, 1, nPol, mbrec.tsys[0], &cStatus);
673
674 // CALFCTR.
675 fits_write_col_flt(cSDptr, ++icol, cRow, 1, nPol, mbrec.calfctr[0],
676 &cStatus);
677
678 if (cHaveBase) {
679 // BASELIN.
680 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 2*nPol, mbrec.baseLin[0][0],
681 &cStatus);
682
683 // BASESUB.
684 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 24*nPol, mbrec.baseSub[0][0],
685 &cStatus);
686 }
687
688 // DATA.
689 fits_write_col_flt(cSDptr, ++icol, cRow, 1, nChan*nPol, mbrec.spectra[0],
690 &cStatus);
691
692 if (cDoTDIM) {
693 // TDIM(DATA).
694 char tdim[16];
695 sprintf(tdim, "(%d,%d,1,1)", nChan, nPol);
696 cptr = tdim;
697 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
698 }
699
700 // FLAGGED.
701 fits_write_col_byt(cSDptr, ++icol, cRow, 1, nChan*nPol, mbrec.flagged[0],
702 &cStatus);
703
704 if (cDoTDIM) {
705 // TDIM(FLAGGED).
706 char tdim[16];
707 sprintf(tdim, "(%d,%d,1,1)", nChan, nPol);
708 cptr = tdim;
709 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
710 }
711
712 if (cDoXPol) {
713 if (cHaveXPol[iIF] && mbrec.xpol[0]) {
714 // XCALFCTR.
715 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 2, mbrec.xcalfctr[0],
716 &cStatus);
717
718 // XPOLDATA.
719 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 2*nChan, mbrec.xpol[0],
720 &cStatus);
721
722 if (cDoTDIM) {
723 // TDIM(XPOLDATA).
724 char tdim[16];
725 sprintf(tdim, "(2,%d)", nChan);
726 cptr = tdim;
727 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
728 }
729
730 } else {
731 // Skip columns.
732 icol += cDoTDIM ? 3 : 2;
733 }
734 }
735
736
737 // Extra system calibration quantities from Parkes.
738 if (cExtraSysCal) {
739 if (cIsMX) {
740 fits_write_col_sht(cSDptr, ++icol, cRow, 1, 1, &mbrec.refBeam, &cStatus);
741 }
742
743 fits_write_col_flt(cSDptr, ++icol, cRow, 1, min(nPol,2), mbrec.tcal[0],
744 &cStatus);
745 cptr = mbrec.tcalTime;
746 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
747
748 float azimuth = mbrec.azimuth * R2D;
749 float elevation = mbrec.elevation * R2D;
750 float parAngle = mbrec.parAngle * R2D;
751 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &azimuth, &cStatus);
752 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &elevation, &cStatus);
753 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &parAngle, &cStatus);
754
755 float focusRot = mbrec.focusRot * R2D;
756 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.focusAxi, &cStatus);
757 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.focusTan, &cStatus);
758 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &focusRot, &cStatus);
759
760 float windAz = mbrec.windAz * R2D;
761 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.temp, &cStatus);
762 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.pressure, &cStatus);
763 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.humidity, &cStatus);
764 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.windSpeed, &cStatus);
765 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &windAz, &cStatus);
766 }
767
768 if (cStatus) {
769 logMsg("ERROR: Failed in writing binary table entry.");
770 }
771
772 return cStatus;
773}
774
775
776//------------------------------------------------------ SDFITSwriter::history
777
778// Write a history record.
779
780int SDFITSwriter::history(char *text)
781
782{
783 if (!cSDptr) {
784 return 1;
785 }
786
787 if (fits_write_history(cSDptr, text, &cStatus)) {
788 logMsg("ERROR: Failed in writing HISTORY records.");
789 }
790
791 return cStatus;
792}
793
794//-------------------------------------------------------- SDFITSwriter::close
795
796// Close the SDFITS file.
797
798void SDFITSwriter::close()
799{
800 if (cSDptr) {
801 cStatus = 0;
802 if (fits_close_file(cSDptr, &cStatus)) {
803 logMsg("ERROR: Failed to close file.");
804 }
805
806 cSDptr = 0;
807 }
808}
809
810//--------------------------------------------------- SDFITSwriter::deleteFile
811
812// Delete the SDFITS file.
813
814void SDFITSwriter::deleteFile()
815{
816 if (cSDptr) {
817 cStatus = 0;
818 if (fits_delete_file(cSDptr, &cStatus)) {
819 logMsg("ERROR: Failed to close and delete file.");
820 }
821
822 cSDptr = 0;
823 }
824}
825
826//------------------------------------------------------- SDFITSwriter::logMsg
827
828// Log a message. If the current CFITSIO status value is non-zero, also log
829// the corresponding error message and dump the CFITSIO message stack.
830
831void SDFITSwriter::logMsg(const char *msg)
832{
833 PKSmsg::logMsg(msg);
834
835 if (cStatus) {
836 fits_get_errstatus(cStatus, cMsg);
837 PKSmsg::logMsg(cMsg);
838
839 while (fits_read_errmsg(cMsg)) {
840 PKSmsg::logMsg(cMsg);
841 }
842 }
843}
Note: See TracBrowser for help on using the repository browser.