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

Last change on this file since 1325 was 1325, checked in by mar637, 17 years ago

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