[1325] | 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 | |
---|
| 44 | using namespace std; |
---|
| 45 | |
---|
| 46 | // Numerical constants. |
---|
| 47 | const double PI = 3.141592653589793238462643; |
---|
| 48 | |
---|
| 49 | // Factor to convert radians to degrees. |
---|
| 50 | const double R2D = 180.0 / PI; |
---|
| 51 | |
---|
| 52 | //------------------------------------------------- SDFITSwriter::SDFITSwriter |
---|
| 53 | |
---|
| 54 | SDFITSwriter::SDFITSwriter() |
---|
| 55 | { |
---|
| 56 | // Default constructor. |
---|
| 57 | cSDptr = 0; |
---|
| 58 | } |
---|
| 59 | |
---|
| 60 | //------------------------------------------------ SDFITSwriter::~SDFITSwriter |
---|
| 61 | |
---|
| 62 | SDFITSwriter::~SDFITSwriter() |
---|
| 63 | { |
---|
| 64 | close(); |
---|
| 65 | } |
---|
| 66 | |
---|
| 67 | //------------------------------------------------------- SDFITSwriter::create |
---|
| 68 | |
---|
| 69 | // Create the output SDFITS file. |
---|
| 70 | |
---|
| 71 | int 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 | |
---|
| 532 | int 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 | |
---|
| 750 | void SDFITSwriter::reportError() |
---|
| 751 | { |
---|
| 752 | fits_report_error(stderr, cStatus); |
---|
| 753 | } |
---|
| 754 | |
---|
| 755 | //-------------------------------------------------------- SDFITSwriter::close |
---|
| 756 | |
---|
| 757 | // Close the SDFITS file. |
---|
| 758 | |
---|
| 759 | void 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 | |
---|
| 772 | void SDFITSwriter::deleteFile() |
---|
| 773 | { |
---|
| 774 | if (cSDptr) { |
---|
| 775 | cStatus = 0; |
---|
| 776 | fits_delete_file(cSDptr, &cStatus); |
---|
| 777 | cSDptr = 0; |
---|
| 778 | } |
---|
| 779 | } |
---|