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

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

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