source: branches/parallel/external-alma/atnf/PKSIO/SDFITSwriter.cc@ 2287

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

New Development: Yes

JIRA Issue: Yes (CAS-2211)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: ASAP 3.0.0 interface changes

Test Programs:

Put in Release Notes: Yes

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

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

Note you also need to update casa/code/atnf.
File size: 26.7 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/SDFITSwriter.h>
38
39#include <casa/Logging/LogIO.h>
40
41#include <casa/iostream.h>
42
43#include <algorithm>
44#include <math.h>
45#include <cstring>
46
47using namespace std;
48
49// Numerical constants.
50const double PI = 3.141592653589793238462643;
51
52// Factor to convert radians to degrees.
53const double R2D = 180.0 / PI;
54
55// Class name
56const string className = "SDFITSwriter" ;
57
58//------------------------------------------------- SDFITSwriter::SDFITSwriter
59
60SDFITSwriter::SDFITSwriter()
61{
62 // Default constructor.
63 cSDptr = 0x0;
64}
65
66//------------------------------------------------ SDFITSwriter::~SDFITSwriter
67
68SDFITSwriter::~SDFITSwriter()
69{
70 close();
71}
72
73//------------------------------------------------------- SDFITSwriter::create
74
75// Create the output SDFITS file.
76
77int SDFITSwriter::create(
78 char* sdName,
79 char* observer,
80 char* project,
81 char* telescope,
82 double antPos[3],
83 char* obsMode,
84 char* bunit,
85 float equinox,
86 char* dopplerFrame,
87 int nIF,
88 int* nChan,
89 int* nPol,
90 int* haveXPol,
91 int haveBase,
92 int extraSysCal)
93{
94 const string methodName = "create()" ;
95
96 if (cSDptr) {
97 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Output file already open, close it first.");
98 return 1;
99 }
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, "Failed to create SDFITS file\n %s", sdName);
110 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, 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 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "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 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "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 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "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 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "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 const string methodName = "write()" ;
563 LogIO os( LogOrigin( className, methodName, WHERE ) ) ;
564
565 char *cptr;
566
567 // Check parameters.
568 int IFno = mbrec.IFno[0];
569 if (IFno < 1 || cNIF < IFno) {
570 os << LogIO::WARN
571 << "SDFITSwriter::write: "
572 << "Invalid IF number " << IFno
573 << " (maximum " << cNIF << ")." << LogIO::POST ;
574 return 1;
575 }
576
577 int iIF = IFno - 1;
578 int nChan = cNChan[iIF];
579 if (mbrec.nChan[0] != nChan) {
580 os << LogIO::WARN
581 << "SDFITSriter::write: "
582 << "Wrong number of channels for IF " << IFno << "," << endl
583 << " "
584 << "got " << nChan << " should be " << mbrec.nChan[0] << "." << endl;
585 os << LogIO::POST ;
586 return 1;
587 }
588
589 int nPol = cNPol[iIF];
590 if (mbrec.nPol[0] != nPol) {
591 os << LogIO::WARN
592 << "SDFITSriter::write: "
593 << "Wrong number of polarizations for IF " << IFno << "," << endl
594 << " "
595 << "got " << nPol << " should be " << mbrec.nPol[0] << "." << endl;
596 os << LogIO::POST ;
597 return 1;
598 }
599
600
601 // Next row.
602 cRow++;
603
604 int icol = 0;
605
606 // SCAN.
607 fits_write_col_int(cSDptr, ++icol, cRow, 1, 1, &mbrec.scanNo, &cStatus);
608
609 // CYCLE.
610 fits_write_col_int(cSDptr, ++icol, cRow, 1, 1, &mbrec.cycleNo, &cStatus);
611
612 // DATE-OBS.
613 cptr = mbrec.datobs;
614 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
615
616 // TIME.
617 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &mbrec.utc, &cStatus);
618
619 // EXPOSURE.
620 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.exposure, &cStatus);
621
622 // OBJECT.
623 cptr = mbrec.srcName;
624 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
625
626 // OBJ-RA.
627 double srcRA = mbrec.srcRA * R2D;
628 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &srcRA, &cStatus);
629
630 // OBJ-DEC.
631 double srcDec = mbrec.srcDec * R2D;
632 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &srcDec, &cStatus);
633
634 // RESTFRQ.
635 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &mbrec.restFreq, &cStatus);
636
637 // OBJECT.
638 cptr = mbrec.obsType;
639 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
640
641 // BEAM.
642 fits_write_col_sht(cSDptr, ++icol, cRow, 1, 1, &mbrec.beamNo, &cStatus);
643
644 // IF.
645 fits_write_col_sht(cSDptr, ++icol, cRow, 1, 1, &mbrec.IFno[0], &cStatus);
646
647 // FREQRES.
648 double freqRes = fabs(mbrec.fqDelt[0]);
649 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &freqRes, &cStatus);
650
651 // BANDWID.
652 double bandwidth = freqRes * nChan;
653 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &bandwidth, &cStatus);
654
655 // CRPIX1.
656 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.fqRefPix[0],
657 &cStatus);
658
659 // CRVAL1.
660 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &mbrec.fqRefVal[0],
661 &cStatus);
662
663 // CDELT1.
664 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &mbrec.fqDelt[0], &cStatus);
665
666 // CRVAL3.
667 double ra = mbrec.ra * R2D;
668 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &ra, &cStatus);
669
670 // CRVAL4.
671 double dec = mbrec.dec * R2D;
672 fits_write_col_dbl(cSDptr, ++icol, cRow, 1, 1, &dec, &cStatus);
673
674 // SCANRATE.
675 float scanrate[2];
676 scanrate[0] = mbrec.raRate * R2D;
677 scanrate[1] = mbrec.decRate * R2D;
678 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 2, scanrate, &cStatus);
679
680 // TSYS.
681 fits_write_col_flt(cSDptr, ++icol, cRow, 1, nPol, mbrec.tsys[0], &cStatus);
682
683 // CALFCTR.
684 fits_write_col_flt(cSDptr, ++icol, cRow, 1, nPol, mbrec.calfctr[0],
685 &cStatus);
686
687 if (cHaveBase) {
688 // BASELIN.
689 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 2*nPol, mbrec.baseLin[0][0],
690 &cStatus);
691
692 // BASESUB.
693 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 24*nPol, mbrec.baseSub[0][0],
694 &cStatus);
695 }
696
697 // DATA.
698 fits_write_col_flt(cSDptr, ++icol, cRow, 1, nChan*nPol, mbrec.spectra[0],
699 &cStatus);
700
701 if (cDoTDIM) {
702 // TDIM(DATA).
703 char tdim[16];
704 sprintf(tdim, "(%d,%d,1,1)", nChan, nPol);
705 cptr = tdim;
706 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
707 }
708
709 // FLAGGED.
710 fits_write_col_byt(cSDptr, ++icol, cRow, 1, nChan*nPol, mbrec.flagged[0],
711 &cStatus);
712
713 if (cDoTDIM) {
714 // TDIM(FLAGGED).
715 char tdim[16];
716 sprintf(tdim, "(%d,%d,1,1)", nChan, nPol);
717 cptr = tdim;
718 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
719 }
720
721 if (cDoXPol) {
722 if (cHaveXPol[iIF] && mbrec.xpol[0]) {
723 // XCALFCTR.
724 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 2, mbrec.xcalfctr[0],
725 &cStatus);
726
727 // XPOLDATA.
728 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 2*nChan, mbrec.xpol[0],
729 &cStatus);
730
731 if (cDoTDIM) {
732 // TDIM(XPOLDATA).
733 char tdim[16];
734 sprintf(tdim, "(2,%d)", nChan);
735 cptr = tdim;
736 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
737 }
738
739 } else {
740 // Skip columns.
741 icol += cDoTDIM ? 3 : 2;
742 }
743 }
744
745
746 // Extra system calibration quantities from Parkes.
747 if (cExtraSysCal) {
748 if (cIsMX) {
749 fits_write_col_sht(cSDptr, ++icol, cRow, 1, 1, &mbrec.refBeam, &cStatus);
750 }
751
752 fits_write_col_flt(cSDptr, ++icol, cRow, 1, min(nPol,2), mbrec.tcal[0],
753 &cStatus);
754 cptr = mbrec.tcalTime;
755 fits_write_col_str(cSDptr, ++icol, cRow, 1, 1, &cptr, &cStatus);
756
757 float azimuth = mbrec.azimuth * R2D;
758 float elevation = mbrec.elevation * R2D;
759 float parAngle = mbrec.parAngle * R2D;
760 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &azimuth, &cStatus);
761 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &elevation, &cStatus);
762 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &parAngle, &cStatus);
763
764 float focusRot = mbrec.focusRot * R2D;
765 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.focusAxi, &cStatus);
766 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.focusTan, &cStatus);
767 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &focusRot, &cStatus);
768
769 float windAz = mbrec.windAz * R2D;
770 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.temp, &cStatus);
771 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.pressure, &cStatus);
772 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.humidity, &cStatus);
773 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &mbrec.windSpeed, &cStatus);
774 fits_write_col_flt(cSDptr, ++icol, cRow, 1, 1, &windAz, &cStatus);
775 }
776
777 if (cStatus) {
778 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed in writing binary table entry.");
779 }
780
781 return cStatus;
782}
783
784
785//------------------------------------------------------ SDFITSwriter::history
786
787// Write a history record.
788
789int SDFITSwriter::history(char *text)
790
791{
792 const string methodName = "history()" ;
793
794 if (!cSDptr) {
795 return 1;
796 }
797
798 if (fits_write_history(cSDptr, text, &cStatus)) {
799 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed in writing HISTORY records.");
800 }
801
802 return cStatus;
803}
804
805//-------------------------------------------------------- SDFITSwriter::close
806
807// Close the SDFITS file.
808
809void SDFITSwriter::close()
810{
811 const string methodName = "close()" ;
812
813 if (cSDptr) {
814 cStatus = 0;
815 if (fits_close_file(cSDptr, &cStatus)) {
816 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed to close file.");
817 }
818
819 cSDptr = 0;
820 }
821}
822
823//--------------------------------------------------- SDFITSwriter::deleteFile
824
825// Delete the SDFITS file.
826
827void SDFITSwriter::deleteFile()
828{
829 const string methodName = "deleteFile()" ;
830
831 if (cSDptr) {
832 cStatus = 0;
833 if (fits_delete_file(cSDptr, &cStatus)) {
834 log(LogOrigin( className, methodName, WHERE ), LogIO::SEVERE, "Failed to close and delete file.");
835 }
836
837 cSDptr = 0;
838 }
839}
840
841//------------------------------------------------------- SDFITSwriter::log
842
843// Log a message. If the current CFITSIO status value is non-zero, also log
844// the corresponding error message and dump the CFITSIO message stack.
845
846void SDFITSwriter::log(LogOrigin origin, LogIO::Command cmd, const char *msg)
847{
848 LogIO os( origin ) ;
849
850 os << cmd << msg << endl ;
851
852 if (cStatus) {
853 fits_get_errstatus(cStatus, cMsg);
854 os << cMsg << endl ;
855
856 while (fits_read_errmsg(cMsg)) {
857 os << cMsg << endl ;
858 }
859 }
860
861 os << LogIO::POST ;
862}
Note: See TracBrowser for help on using the repository browser.