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

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

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