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

Last change on this file since 1635 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.