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

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

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