source: branches/alma/external/atnf/PKSIO/SDFITSwriter.cc@ 1393

Last change on this file since 1393 was 1325, checked in by mar637, 18 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.