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

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