source: branches/mergetest/external/atnf/PKSIO/SDFITSwriter.cc @ 1779

Last change on this file since 1779 was 1779, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: No (test merging alma branch)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

Description:


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.