source: tags/asap2.3.1/external/atnf/PKSIO/SDFITSwriter.cc @ 2192

Last change on this file since 2192 was 1466, checked in by Malte Marquarding, 15 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.