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

Last change on this file since 1399 was 1399, checked in by Malte Marquarding, 17 years ago

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