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

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

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