[1325] | 1 | //#--------------------------------------------------------------------------- |
---|
| 2 | //# MBFITSreader.cc: ATNF single-dish RPFITS reader. |
---|
| 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 | //# |
---|
[1372] | 29 | //# $Id: MBFITSreader.cc,v 19.34 2007/07/02 06:12:18 cal103 Exp $ |
---|
[1325] | 30 | //#--------------------------------------------------------------------------- |
---|
| 31 | //# The MBFITSreader class reads single dish RPFITS files (such as Parkes |
---|
| 32 | //# Multibeam MBFITS files). |
---|
| 33 | //# |
---|
| 34 | //# Original: 2000/07/28 Mark Calabretta |
---|
| 35 | //#--------------------------------------------------------------------------- |
---|
| 36 | |
---|
| 37 | #include <atnf/PKSIO/MBFITSreader.h> |
---|
| 38 | #include <atnf/PKSIO/PKSMBrecord.h> |
---|
| 39 | |
---|
| 40 | #include <RPFITS.h> |
---|
| 41 | |
---|
| 42 | #include <casa/math.h> |
---|
| 43 | #include <casa/iostream.h> |
---|
| 44 | #include <casa/stdio.h> |
---|
| 45 | #include <casa/stdlib.h> |
---|
| 46 | #include <casa/string.h> |
---|
| 47 | #include <unistd.h> |
---|
| 48 | |
---|
| 49 | using namespace std; |
---|
| 50 | |
---|
| 51 | // Numerical constants. |
---|
| 52 | const double PI = 3.141592653589793238462643; |
---|
| 53 | const double TWOPI = 2.0 * PI; |
---|
| 54 | |
---|
| 55 | //------------------------------------------------- MBFITSreader::MBFITSreader |
---|
| 56 | |
---|
| 57 | // Default constructor. |
---|
| 58 | |
---|
| 59 | MBFITSreader::MBFITSreader( |
---|
| 60 | const int retry, |
---|
| 61 | const int interpolate) |
---|
| 62 | { |
---|
| 63 | cRetry = retry; |
---|
| 64 | if (cRetry > 10) { |
---|
| 65 | cRetry = 10; |
---|
| 66 | } |
---|
| 67 | |
---|
| 68 | cInterp = interpolate; |
---|
| 69 | if (cInterp < 0 || cInterp > 2) { |
---|
| 70 | cInterp = 1; |
---|
| 71 | } |
---|
| 72 | |
---|
| 73 | // Initialize pointers. |
---|
| 74 | cBeams = 0x0; |
---|
| 75 | cIFs = 0x0; |
---|
| 76 | cNChan = 0x0; |
---|
| 77 | cNPol = 0x0; |
---|
| 78 | cHaveXPol = 0x0; |
---|
| 79 | cStartChan = 0x0; |
---|
| 80 | cEndChan = 0x0; |
---|
| 81 | cRefChan = 0x0; |
---|
| 82 | |
---|
| 83 | cVis = 0x0; |
---|
| 84 | cWgt = 0x0; |
---|
| 85 | |
---|
| 86 | cBeamSel = 0x0; |
---|
| 87 | cIFSel = 0x0; |
---|
| 88 | cChanOff = 0x0; |
---|
| 89 | cXpolOff = 0x0; |
---|
| 90 | cBuffer = 0x0; |
---|
| 91 | cPosUTC = 0x0; |
---|
| 92 | |
---|
| 93 | cMBopen = 0; |
---|
| 94 | } |
---|
| 95 | |
---|
| 96 | //------------------------------------------------ MBFITSreader::~MBFITSreader |
---|
| 97 | |
---|
| 98 | // Destructor. |
---|
| 99 | |
---|
| 100 | MBFITSreader::~MBFITSreader() |
---|
| 101 | { |
---|
| 102 | close(); |
---|
| 103 | } |
---|
| 104 | |
---|
| 105 | //--------------------------------------------------------- MBFITSreader::open |
---|
| 106 | |
---|
| 107 | // Open the RPFITS file for reading. |
---|
| 108 | |
---|
| 109 | int MBFITSreader::open( |
---|
| 110 | char *rpname, |
---|
| 111 | int &nBeam, |
---|
| 112 | int* &beams, |
---|
| 113 | int &nIF, |
---|
| 114 | int* &IFs, |
---|
| 115 | int* &nChan, |
---|
| 116 | int* &nPol, |
---|
| 117 | int* &haveXPol, |
---|
| 118 | int &haveBase, |
---|
| 119 | int &haveSpectra, |
---|
| 120 | int &extraSysCal) |
---|
| 121 | { |
---|
| 122 | if (cMBopen) { |
---|
| 123 | close(); |
---|
| 124 | } |
---|
| 125 | |
---|
| 126 | strcpy(names_.file, rpname); |
---|
| 127 | |
---|
| 128 | // Open the RPFITS file. |
---|
| 129 | int jstat = -3; |
---|
| 130 | rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag, |
---|
| 131 | &cBin, &cIFno, &cSrcNo); |
---|
| 132 | |
---|
| 133 | if (jstat) { |
---|
[1372] | 134 | fprintf(stderr, "ERROR, failed to open MBFITS file: %s\n", rpname); |
---|
[1325] | 135 | return 1; |
---|
| 136 | } |
---|
| 137 | |
---|
| 138 | cMBopen = 1; |
---|
| 139 | |
---|
| 140 | // Tell RPFITSIN that we want the OBSTYPE card. |
---|
| 141 | int j; |
---|
| 142 | param_.ncard = 1; |
---|
| 143 | for (j = 0; j < 80; j++) { |
---|
| 144 | names_.card[j] = ' '; |
---|
| 145 | } |
---|
| 146 | strncpy(names_.card, "OBSTYPE", 7); |
---|
| 147 | |
---|
| 148 | // Read the first header. |
---|
| 149 | jstat = -1; |
---|
| 150 | rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag, |
---|
| 151 | &cBin, &cIFno, &cSrcNo); |
---|
| 152 | |
---|
| 153 | if (jstat) { |
---|
[1372] | 154 | fprintf(stderr, "ERROR, failed to read MBFITS header: %s\n", rpname); |
---|
[1325] | 155 | close(); |
---|
| 156 | return 1; |
---|
| 157 | } |
---|
| 158 | |
---|
| 159 | // Mopra data has some peculiarities. |
---|
| 160 | cMopra = strncmp(names_.instrument, "ATMOPRA", 7) == 0; |
---|
| 161 | |
---|
[1372] | 162 | // Non-ATNF data may not store the position in (u,v,w). |
---|
| 163 | if (strncmp(names_.sta, "tid", 3) == 0) { |
---|
| 164 | fprintf(stderr, "WARNING, found Tidbinbilla data"); |
---|
| 165 | cSUpos = 1; |
---|
| 166 | } else if (strncmp(names_.sta, "HOB", 3) == 0) { |
---|
| 167 | fprintf(stderr, "WARNING, found Hobart data"); |
---|
| 168 | cSUpos = 1; |
---|
| 169 | } else if (strncmp(names_.sta, "CED", 3) == 0) { |
---|
| 170 | fprintf(stderr, "WARNING, found Ceduna data"); |
---|
| 171 | cSUpos = 1; |
---|
| 172 | } else { |
---|
| 173 | cSUpos = 0; |
---|
| 174 | } |
---|
| 175 | |
---|
| 176 | if (cSUpos) { |
---|
| 177 | fprintf(stderr, ", using telescope position from SU table.\n"); |
---|
[1325] | 178 | cInterp = 0; |
---|
| 179 | } |
---|
| 180 | |
---|
| 181 | |
---|
| 182 | // Find the maximum beam number. |
---|
| 183 | cNBeam = 0; |
---|
| 184 | for (int iBeam = 0; iBeam < anten_.nant; iBeam++) { |
---|
| 185 | if (anten_.ant_num[iBeam] > cNBeam) { |
---|
| 186 | cNBeam = anten_.ant_num[iBeam]; |
---|
| 187 | } |
---|
| 188 | } |
---|
| 189 | |
---|
| 190 | if (cNBeam <= 0) { |
---|
[1372] | 191 | fprintf(stderr, "ERROR, couldn't determine number of beams.\n"); |
---|
[1325] | 192 | close(); |
---|
| 193 | return 1; |
---|
| 194 | } |
---|
| 195 | |
---|
| 196 | // Construct the beam mask. |
---|
| 197 | cBeams = new int[cNBeam]; |
---|
| 198 | for (int iBeam = 0; iBeam < cNBeam; iBeam++) { |
---|
| 199 | cBeams[iBeam] = 0; |
---|
| 200 | } |
---|
| 201 | |
---|
| 202 | // ...beams present in the data. |
---|
| 203 | for (int iBeam = 0; iBeam < anten_.nant; iBeam++) { |
---|
| 204 | cBeams[anten_.ant_num[iBeam] - 1] = 1; |
---|
| 205 | } |
---|
| 206 | |
---|
| 207 | // Passing back the address of the array allows PKSFITSreader::select() to |
---|
| 208 | // modify its elements directly. |
---|
| 209 | nBeam = cNBeam; |
---|
| 210 | beams = cBeams; |
---|
| 211 | |
---|
| 212 | |
---|
| 213 | // Number of IFs. |
---|
| 214 | cNIF = if_.n_if; |
---|
| 215 | cIFs = new int[cNIF]; |
---|
| 216 | for (int iIF = 0; iIF < cNIF; iIF++) { |
---|
| 217 | cIFs[iIF] = 1; |
---|
| 218 | } |
---|
| 219 | |
---|
| 220 | // Passing back the address of the array allows PKSFITSreader::select() to |
---|
| 221 | // modify its elements directly. |
---|
| 222 | nIF = cNIF; |
---|
| 223 | IFs = cIFs; |
---|
| 224 | |
---|
| 225 | |
---|
| 226 | // Number of channels and polarizations. |
---|
| 227 | cNChan = new int[cNIF]; |
---|
| 228 | cNPol = new int[cNIF]; |
---|
| 229 | cHaveXPol = new int[cNIF]; |
---|
| 230 | cGetXPol = 0; |
---|
| 231 | |
---|
| 232 | int maxProd = 0; |
---|
| 233 | for (int iIF = 0; iIF < cNIF; iIF++) { |
---|
| 234 | cNChan[iIF] = if_.if_nfreq[iIF]; |
---|
| 235 | cNPol[iIF] = if_.if_nstok[iIF]; |
---|
| 236 | cNChan[iIF] -= cNChan[iIF]%2; |
---|
| 237 | |
---|
| 238 | // Do we have cross-polarization data? |
---|
| 239 | if ((cHaveXPol[iIF] = cNPol[iIF] > 2)) { |
---|
| 240 | // Cross-polarization data is handled separately. |
---|
| 241 | cNPol[iIF] = 2; |
---|
| 242 | |
---|
| 243 | // Default is to get it if we have it. |
---|
| 244 | cGetXPol = 1; |
---|
| 245 | } |
---|
| 246 | |
---|
| 247 | // Maximum number of spectral products in any IF. |
---|
| 248 | int nProd = if_.if_nfreq[iIF] * if_.if_nstok[iIF]; |
---|
| 249 | if (maxProd < nProd) maxProd = nProd; |
---|
| 250 | } |
---|
| 251 | |
---|
| 252 | // Allocate memory for RPFITSIN subroutine arguments. |
---|
| 253 | if (cVis) delete [] cVis; |
---|
| 254 | if (cWgt) delete [] cWgt; |
---|
| 255 | cVis = new float[2*maxProd]; |
---|
| 256 | cWgt = new float[maxProd]; |
---|
| 257 | |
---|
| 258 | nChan = cNChan; |
---|
| 259 | nPol = cNPol; |
---|
| 260 | haveXPol = cHaveXPol; |
---|
| 261 | |
---|
| 262 | |
---|
| 263 | // Default channel range selection. |
---|
| 264 | cStartChan = new int[cNIF]; |
---|
| 265 | cEndChan = new int[cNIF]; |
---|
| 266 | cRefChan = new int[cNIF]; |
---|
| 267 | |
---|
| 268 | for (int iIF = 0; iIF < cNIF; iIF++) { |
---|
| 269 | cStartChan[iIF] = 1; |
---|
| 270 | cEndChan[iIF] = cNChan[iIF]; |
---|
| 271 | cRefChan[iIF] = cNChan[iIF]/2 + 1; |
---|
| 272 | } |
---|
| 273 | |
---|
| 274 | cGetSpectra = 1; |
---|
| 275 | |
---|
| 276 | |
---|
| 277 | // No baseline parameters in MBFITS. |
---|
| 278 | haveBase = 0; |
---|
| 279 | |
---|
| 280 | // Always have spectra in MBFITS. |
---|
| 281 | haveSpectra = cHaveSpectra = 1; |
---|
| 282 | |
---|
| 283 | |
---|
| 284 | // Integration cycle time (s). |
---|
| 285 | cIntTime = param_.intime; |
---|
| 286 | |
---|
| 287 | // Can't deduce binning mode till later. |
---|
| 288 | cNBin = 0; |
---|
| 289 | |
---|
| 290 | |
---|
| 291 | // Read the first syscal record. |
---|
| 292 | if (rpget(1, cEOS)) { |
---|
[1372] | 293 | fprintf(stderr, "ERROR, failed to read first syscal record.\n"); |
---|
[1325] | 294 | close(); |
---|
| 295 | return 1; |
---|
| 296 | } |
---|
| 297 | |
---|
| 298 | // Additional information for Parkes Multibeam data? |
---|
| 299 | extraSysCal = (sc_.sc_ant > anten_.nant); |
---|
| 300 | |
---|
| 301 | |
---|
| 302 | cFirst = 1; |
---|
| 303 | cEOF = 0; |
---|
| 304 | cFlushing = 0; |
---|
| 305 | |
---|
| 306 | return 0; |
---|
| 307 | } |
---|
| 308 | |
---|
| 309 | //---------------------------------------------------- MBFITSreader::getHeader |
---|
| 310 | |
---|
| 311 | // Get parameters describing the data. |
---|
| 312 | |
---|
| 313 | int MBFITSreader::getHeader( |
---|
| 314 | char observer[32], |
---|
| 315 | char project[32], |
---|
| 316 | char telescope[32], |
---|
| 317 | double antPos[3], |
---|
| 318 | char obsType[32], |
---|
| 319 | float &equinox, |
---|
| 320 | char radecsys[32], |
---|
| 321 | char dopplerFrame[32], |
---|
| 322 | char datobs[32], |
---|
| 323 | double &utc, |
---|
| 324 | double &refFreq, |
---|
| 325 | double &bandwidth) |
---|
| 326 | { |
---|
| 327 | if (!cMBopen) { |
---|
[1372] | 328 | fprintf(stderr, "ERROR, an MBFITS file has not been opened.\n"); |
---|
[1325] | 329 | return 1; |
---|
| 330 | } |
---|
| 331 | |
---|
| 332 | sprintf(observer, "%-16.16s", names_.rp_observer); |
---|
| 333 | sprintf(project, "%-16.16s", names_.object); |
---|
| 334 | sprintf(telescope, "%-16.16s", names_.instrument); |
---|
| 335 | |
---|
| 336 | // Observatory coordinates (ITRF), in m. |
---|
| 337 | antPos[0] = doubles_.x[0]; |
---|
| 338 | antPos[1] = doubles_.y[0]; |
---|
| 339 | antPos[2] = doubles_.z[0]; |
---|
| 340 | |
---|
| 341 | // This is the only sure way to identify the telescope, maybe. |
---|
| 342 | if (strncmp(names_.sta, "MB0", 3) == 0) { |
---|
| 343 | // Parkes Multibeam. |
---|
| 344 | sprintf(telescope, "%-16.16s", "ATPKSMB"); |
---|
| 345 | antPos[0] = -4554232.087; |
---|
| 346 | antPos[1] = 2816759.046; |
---|
| 347 | antPos[2] = -3454035.950; |
---|
| 348 | } else if (strncmp(names_.sta, "HOH", 3) == 0) { |
---|
| 349 | // Parkes HOH receiver. |
---|
| 350 | sprintf(telescope, "%-16.16s", "ATPKSHOH"); |
---|
| 351 | antPos[0] = -4554232.087; |
---|
| 352 | antPos[1] = 2816759.046; |
---|
| 353 | antPos[2] = -3454035.950; |
---|
| 354 | } else if (strncmp(names_.sta, "CA0", 3) == 0) { |
---|
| 355 | // An ATCA antenna, use the array centre position. |
---|
| 356 | sprintf(telescope, "%-16.16s", "ATCA"); |
---|
| 357 | antPos[0] = -4750915.837; |
---|
| 358 | antPos[1] = 2792906.182; |
---|
| 359 | antPos[2] = -3200483.747; |
---|
| 360 | } else if (strncmp(names_.sta, "MOP", 3) == 0) { |
---|
| 361 | // Mopra. |
---|
| 362 | sprintf(telescope, "%-16.16s", "ATMOPRA"); |
---|
| 363 | antPos[0] = -4682768.630; |
---|
| 364 | antPos[1] = 2802619.060; |
---|
| 365 | antPos[2] = -3291759.900; |
---|
| 366 | } else if (strncmp(names_.sta, "HOB", 3) == 0) { |
---|
| 367 | // Hobart. |
---|
| 368 | sprintf(telescope, "%-16.16s", "HOBART"); |
---|
| 369 | antPos[0] = -3950236.735; |
---|
| 370 | antPos[1] = 2522347.567; |
---|
| 371 | antPos[2] = -4311562.569; |
---|
| 372 | } else if (strncmp(names_.sta, "CED", 3) == 0) { |
---|
| 373 | // Ceduna. |
---|
| 374 | sprintf(telescope, "%-16.16s", "CEDUNA"); |
---|
| 375 | antPos[0] = -3749943.657; |
---|
| 376 | antPos[1] = 3909017.709; |
---|
| 377 | antPos[2] = -3367518.309; |
---|
| 378 | } else if (strncmp(names_.sta, "tid", 3) == 0) { |
---|
| 379 | // DSS. |
---|
| 380 | sprintf(telescope, "%-16.16s", "DSS-43"); |
---|
| 381 | antPos[0] = -4460894.727; |
---|
| 382 | antPos[1] = 2682361.530; |
---|
| 383 | antPos[2] = -3674748.424; |
---|
| 384 | } |
---|
| 385 | |
---|
| 386 | // Observation type. |
---|
| 387 | int j; |
---|
| 388 | for (j = 0; j < 31; j++) { |
---|
| 389 | obsType[j] = names_.card[11+j]; |
---|
| 390 | if (obsType[j] == '\'') break; |
---|
| 391 | } |
---|
| 392 | obsType[j] = '\0'; |
---|
| 393 | |
---|
| 394 | // Coordinate frames. |
---|
| 395 | equinox = 2000.0f; |
---|
| 396 | strcpy(radecsys, "FK5"); |
---|
| 397 | strcpy(dopplerFrame, "TOPOCENT"); |
---|
| 398 | |
---|
| 399 | // Time at start of observation. |
---|
| 400 | sprintf(datobs, "%-10.10s", names_.datobs); |
---|
| 401 | utc = cUTC; |
---|
| 402 | |
---|
| 403 | // Spectral parameters. |
---|
| 404 | refFreq = doubles_.if_freq[0]; |
---|
| 405 | bandwidth = doubles_.if_bw[0]; |
---|
| 406 | |
---|
| 407 | return 0; |
---|
| 408 | } |
---|
| 409 | |
---|
| 410 | //-------------------------------------------------- MBFITSreader::getFreqInfo |
---|
| 411 | |
---|
| 412 | // Get frequency parameters for each IF. |
---|
| 413 | |
---|
| 414 | int MBFITSreader::getFreqInfo( |
---|
| 415 | int &nIF, |
---|
| 416 | double* &startFreq, |
---|
| 417 | double* &endFreq) |
---|
| 418 | { |
---|
| 419 | // This is RPFITS - can't do it! |
---|
| 420 | return 1; |
---|
| 421 | } |
---|
| 422 | |
---|
| 423 | //---------------------------------------------------- MBFITSreader::findRange |
---|
| 424 | |
---|
| 425 | // Find the range of the data selected in time and position. |
---|
| 426 | |
---|
| 427 | int MBFITSreader::findRange( |
---|
| 428 | int &nRow, |
---|
| 429 | int &nSel, |
---|
| 430 | char dateSpan[2][32], |
---|
| 431 | double utcSpan[2], |
---|
| 432 | double* &positions) |
---|
| 433 | { |
---|
| 434 | // This is RPFITS - can't do it! |
---|
| 435 | return 1; |
---|
| 436 | } |
---|
| 437 | |
---|
| 438 | //--------------------------------------------------------- MBFITSreader::read |
---|
| 439 | |
---|
| 440 | // Read the next data record. |
---|
| 441 | |
---|
| 442 | int MBFITSreader::read( |
---|
| 443 | PKSMBrecord &MBrec) |
---|
| 444 | { |
---|
| 445 | int beamNo = -1; |
---|
| 446 | int haveData, status; |
---|
| 447 | PKSMBrecord *iMBuff = 0x0; |
---|
| 448 | |
---|
| 449 | if (!cMBopen) { |
---|
[1372] | 450 | fprintf(stderr, "ERROR, an MBFITS file has not been opened.\n"); |
---|
[1325] | 451 | return 1; |
---|
| 452 | } |
---|
| 453 | |
---|
| 454 | // Positions recorded in the input records do not coincide with the midpoint |
---|
| 455 | // of the integration and hence the input must be buffered so that true |
---|
| 456 | // positions may be interpolated. |
---|
| 457 | // |
---|
| 458 | // On the first call nBeamSel buffers of length nBin, are allocated and |
---|
| 459 | // filled, where nBin is the number of time bins. |
---|
| 460 | // |
---|
| 461 | // The input records for binned, single beam data with multiple simultaneous |
---|
| 462 | // IFs are ordered by IF within each integration rather than by bin number |
---|
| 463 | // and hence are not in time order. No multibeam data exists with |
---|
| 464 | // nBin > 1 but the likelihood that the input records would be in beam/IF |
---|
| 465 | // order and the requirement that output records be in time order would |
---|
| 466 | // force an elaborate double-buffering system and we do not support it. |
---|
| 467 | // |
---|
| 468 | // Once all buffers are filled, the next record for each beam pertains to |
---|
| 469 | // the next integration and should contain new position information allowing |
---|
| 470 | // the proper position for each spectrum in the buffer to be interpolated. |
---|
| 471 | // The buffers are then flushed in time order. For single beam data there |
---|
| 472 | // is only one buffer and reads from the MBFITS file are suspended while the |
---|
| 473 | // flush is in progress. For multibeam data each buffer is of unit length |
---|
| 474 | // so the flush completes immediately and the new record takes its place. |
---|
| 475 | |
---|
| 476 | haveData = 0; |
---|
| 477 | while (!haveData) { |
---|
| 478 | int iBeamSel = -1, iIFSel = -1; |
---|
| 479 | |
---|
| 480 | if (!cFlushing) { |
---|
| 481 | if (cEOF) { |
---|
| 482 | return -1; |
---|
| 483 | } |
---|
| 484 | |
---|
| 485 | // Read the next record. |
---|
| 486 | if ((status = rpget(0, cEOS)) == -1) { |
---|
| 487 | // EOF. |
---|
| 488 | cEOF = 1; |
---|
| 489 | cFlushing = 1; |
---|
| 490 | cFlushBin = 0; |
---|
| 491 | cFlushIF = 0; |
---|
| 492 | |
---|
| 493 | #ifdef PKSIO_DEBUG |
---|
| 494 | printf("End-of-file detected, flushing last scan.\n"); |
---|
| 495 | #endif |
---|
| 496 | |
---|
| 497 | } else if (status) { |
---|
| 498 | // IO error. |
---|
| 499 | return 1; |
---|
| 500 | |
---|
| 501 | } else { |
---|
| 502 | if (cFirst) { |
---|
| 503 | // First data; cBeamSel[] stores the buffer index for each beam. |
---|
| 504 | cNBeamSel = 0; |
---|
| 505 | cBeamSel = new int[cNBeam]; |
---|
| 506 | |
---|
| 507 | for (int iBeam = 0; iBeam < cNBeam; iBeam++) { |
---|
| 508 | if (cBeams[iBeam]) { |
---|
| 509 | // Buffer offset for this beam. |
---|
| 510 | cBeamSel[iBeam] = cNBeamSel++; |
---|
| 511 | } else { |
---|
| 512 | // Signal that the beam is not selected. |
---|
| 513 | cBeamSel[iBeam] = -1; |
---|
| 514 | } |
---|
| 515 | } |
---|
| 516 | |
---|
| 517 | // Set up bookkeeping arrays for IFs. |
---|
| 518 | cIFSel = new int[cNIF]; |
---|
| 519 | cChanOff = new int[cNIF]; |
---|
| 520 | cXpolOff = new int[cNIF]; |
---|
| 521 | |
---|
| 522 | int simulIF = 0; |
---|
| 523 | int maxChan = 0; |
---|
| 524 | int maxXpol = 0; |
---|
| 525 | |
---|
| 526 | for (int iIF = 0; iIF < cNIF; iIF++) { |
---|
| 527 | if (cIFs[iIF]) { |
---|
| 528 | // Buffer index for each IF within each simultaneous set. |
---|
| 529 | cIFSel[iIF] = 0; |
---|
| 530 | |
---|
| 531 | // Array offsets for each IF within each simultaneous set. |
---|
| 532 | cChanOff[iIF] = 0; |
---|
| 533 | cXpolOff[iIF] = 0; |
---|
| 534 | |
---|
| 535 | // Look for earlier IFs in the same simultaneous set. |
---|
| 536 | for (int jIF = 0; jIF < iIF; jIF++) { |
---|
| 537 | if (!cIFs[jIF]) continue; |
---|
| 538 | |
---|
| 539 | if (if_.if_simul[jIF] == if_.if_simul[iIF]) { |
---|
| 540 | // Got one, increment indices. |
---|
| 541 | cIFSel[iIF]++; |
---|
| 542 | |
---|
| 543 | cChanOff[iIF] += cNChan[jIF] * cNPol[jIF]; |
---|
| 544 | if (cHaveXPol[jIF]) { |
---|
| 545 | cXpolOff[iIF] += 2 * cNChan[jIF]; |
---|
| 546 | } |
---|
| 547 | } |
---|
| 548 | } |
---|
| 549 | |
---|
| 550 | // Maximum number of selected IFs in any simultaneous set. |
---|
| 551 | simulIF = max(simulIF, cIFSel[iIF]+1); |
---|
| 552 | |
---|
| 553 | // Maximum memory required for any simultaneous set. |
---|
| 554 | maxChan = max(maxChan, cChanOff[iIF] + cNChan[iIF]*cNPol[iIF]); |
---|
| 555 | if (cHaveXPol[iIF]) { |
---|
| 556 | maxXpol = max(maxXpol, cXpolOff[iIF] + 2*cNChan[iIF]); |
---|
| 557 | } |
---|
| 558 | |
---|
| 559 | } else { |
---|
| 560 | // Signal that the IF is not selected. |
---|
| 561 | cIFSel[iIF] = -1; |
---|
| 562 | } |
---|
| 563 | } |
---|
| 564 | |
---|
| 565 | // Check for binning mode observations. |
---|
| 566 | if (param_.intbase > 0.0f) { |
---|
| 567 | cNBin = int((cIntTime / param_.intbase) + 0.5); |
---|
| 568 | |
---|
| 569 | // intbase sometimes contains rubbish. |
---|
| 570 | if (cNBin == 0) { |
---|
| 571 | cNBin = 1; |
---|
| 572 | } |
---|
| 573 | } else { |
---|
| 574 | cNBin = 1; |
---|
| 575 | } |
---|
| 576 | |
---|
| 577 | if (cNBin > 1 && cNBeamSel > 1) { |
---|
[1372] | 578 | fprintf(stderr, "ERROR, cannot handle binning mode for multiple " |
---|
[1325] | 579 | "beams.\n"); |
---|
| 580 | close(); |
---|
| 581 | return 1; |
---|
| 582 | } |
---|
| 583 | |
---|
| 584 | // Allocate buffer data storage; the PKSMBrecord constructor zeroes |
---|
| 585 | // class members such as cycleNo that are tested in the first pass |
---|
| 586 | // below. |
---|
| 587 | int nBuff = cNBeamSel * cNBin; |
---|
| 588 | cBuffer = new PKSMBrecord[nBuff]; |
---|
| 589 | |
---|
| 590 | // Allocate memory for spectral arrays. |
---|
| 591 | for (int ibuff = 0; ibuff < nBuff; ibuff++) { |
---|
| 592 | cBuffer[ibuff].setNIFs(simulIF); |
---|
| 593 | cBuffer[ibuff].allocate(0, maxChan, maxXpol); |
---|
| 594 | } |
---|
| 595 | |
---|
| 596 | cPosUTC = new double[cNBeamSel]; |
---|
| 597 | |
---|
| 598 | cFirst = 0; |
---|
| 599 | cScanNo = 1; |
---|
| 600 | cCycleNo = 0; |
---|
| 601 | cPrevUTC = 0.0; |
---|
| 602 | cStaleness = new int[cNBeamSel]; |
---|
| 603 | for (int iBeamSel = 0; iBeamSel < cNBeamSel; iBeamSel++) { |
---|
| 604 | cStaleness[iBeamSel] = 0; |
---|
| 605 | } |
---|
| 606 | } |
---|
| 607 | |
---|
| 608 | // Check for end-of-scan. |
---|
| 609 | if (cEOS) { |
---|
| 610 | cScanNo++; |
---|
| 611 | cCycleNo = 0; |
---|
| 612 | cPrevUTC = 0.0; |
---|
| 613 | } |
---|
| 614 | |
---|
| 615 | // Apply beam selection. |
---|
| 616 | beamNo = int(cBaseline / 256.0); |
---|
| 617 | iBeamSel = cBeamSel[beamNo-1]; |
---|
| 618 | if (iBeamSel < 0) continue; |
---|
| 619 | |
---|
| 620 | // Sanity check (mainly for MOPS). |
---|
| 621 | if (cIFno > cNIF) continue; |
---|
| 622 | |
---|
| 623 | // Apply IF selection. |
---|
| 624 | iIFSel = cIFSel[cIFno - 1]; |
---|
| 625 | if (iIFSel < 0) continue; |
---|
| 626 | |
---|
| 627 | sprintf(cDateObs, "%-10.10s", names_.datobs); |
---|
| 628 | |
---|
| 629 | // Check for change-of-day. |
---|
| 630 | if (cUTC < cPrevUTC - 85800.0) { |
---|
| 631 | cUTC += 86400.0; |
---|
| 632 | } |
---|
| 633 | |
---|
| 634 | if (cNBin > 1) { |
---|
| 635 | // Binning mode: correct the time. |
---|
| 636 | cUTC += param_.intbase * (cBin - (cNBin + 1)/2.0); |
---|
| 637 | } |
---|
| 638 | |
---|
| 639 | // New integration cycle? |
---|
| 640 | if (cUTC > cPrevUTC) { |
---|
| 641 | cCycleNo++; |
---|
| 642 | cPrevUTC = cUTC + 0.0001; |
---|
| 643 | } |
---|
| 644 | |
---|
| 645 | // Compute buffer number. |
---|
| 646 | iMBuff = cBuffer + iBeamSel; |
---|
| 647 | if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1); |
---|
| 648 | |
---|
| 649 | if (cCycleNo < iMBuff->cycleNo) { |
---|
| 650 | // Note that if the first beam and IF are not both selected cEOS |
---|
| 651 | // will be cleared by rpget() when the next beam/IF is read. |
---|
| 652 | cEOS = 1; |
---|
| 653 | } |
---|
| 654 | |
---|
| 655 | // Begin flush cycle? |
---|
| 656 | if (cEOS || (iMBuff->nIF && cUTC > iMBuff->utc + 0.0001)) { |
---|
| 657 | cFlushing = 1; |
---|
| 658 | cFlushBin = 0; |
---|
| 659 | cFlushIF = 0; |
---|
| 660 | } |
---|
| 661 | |
---|
| 662 | #ifdef PKSIO_DEBUG |
---|
| 663 | printf(" In:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, cIFno); |
---|
| 664 | if (cEOS) printf("Start of new scan, flushing previous scan.\n"); |
---|
| 665 | #endif |
---|
| 666 | } |
---|
| 667 | } |
---|
| 668 | |
---|
| 669 | |
---|
| 670 | if (cFlushing) { |
---|
| 671 | // Find the oldest integration to flush, noting that the last |
---|
| 672 | // integration cycle may be incomplete. |
---|
| 673 | beamNo = 0; |
---|
| 674 | int cycleNo = 0; |
---|
| 675 | for (; cFlushBin < cNBin; cFlushBin++) { |
---|
| 676 | for (iBeamSel = 0; iBeamSel < cNBeamSel; iBeamSel++) { |
---|
| 677 | iMBuff = cBuffer + iBeamSel + cNBeamSel*cFlushBin; |
---|
| 678 | |
---|
| 679 | // iMBuff->nIF is set to zero (below) to signal that all IFs in |
---|
| 680 | // an integration have been flushed. |
---|
| 681 | if (iMBuff->nIF) { |
---|
| 682 | if (cycleNo == 0 || iMBuff->cycleNo < cycleNo) { |
---|
| 683 | beamNo = iMBuff->beamNo; |
---|
| 684 | cycleNo = iMBuff->cycleNo; |
---|
| 685 | } |
---|
| 686 | } |
---|
| 687 | } |
---|
| 688 | |
---|
| 689 | if (beamNo) { |
---|
| 690 | // Found an integration to flush. |
---|
| 691 | break; |
---|
| 692 | } |
---|
| 693 | } |
---|
| 694 | |
---|
| 695 | if (beamNo) { |
---|
| 696 | iBeamSel = cBeamSel[beamNo-1]; |
---|
| 697 | iMBuff = cBuffer + iBeamSel + cNBeamSel*cFlushBin; |
---|
| 698 | |
---|
| 699 | // Find the IF to flush. |
---|
| 700 | for (; cFlushIF < iMBuff->nIF; cFlushIF++) { |
---|
| 701 | if (iMBuff->IFno[cFlushIF]) break; |
---|
| 702 | } |
---|
| 703 | |
---|
| 704 | } else { |
---|
| 705 | // Flush complete. |
---|
| 706 | cFlushing = 0; |
---|
| 707 | if (cEOF) { |
---|
| 708 | return -1; |
---|
| 709 | } |
---|
| 710 | |
---|
| 711 | // The last record read must have been the first of a new cycle. |
---|
| 712 | beamNo = int(cBaseline / 256.0); |
---|
| 713 | iBeamSel = cBeamSel[beamNo-1]; |
---|
| 714 | |
---|
| 715 | // Compute buffer number. |
---|
| 716 | iMBuff = cBuffer + iBeamSel; |
---|
| 717 | if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1); |
---|
| 718 | } |
---|
| 719 | } |
---|
| 720 | |
---|
| 721 | |
---|
| 722 | if (cFlushing && cFlushBin == 0 && cFlushIF == 0 && cInterp) { |
---|
| 723 | // Interpolate the beam position at the start of the flush cycle. |
---|
| 724 | #ifdef PKSIO_DEBUG |
---|
| 725 | printf("Doing position interpolation for beam %d.\n", iMBuff->beamNo); |
---|
| 726 | #endif |
---|
| 727 | |
---|
| 728 | double prevRA = iMBuff->ra; |
---|
| 729 | double prevDec = iMBuff->dec; |
---|
| 730 | double prevUTC = cPosUTC[iBeamSel]; |
---|
| 731 | |
---|
| 732 | if (!cEOF && !cEOS) { |
---|
| 733 | // The position is measured by the control system at a time returned |
---|
| 734 | // by RPFITSIN as the 'w' visibility coordinate. The ra and dec, |
---|
| 735 | // returned as the 'u' and 'v' visibility coordinates, must be |
---|
| 736 | // interpolated to the integration time which RPFITSIN returns as |
---|
| 737 | // 'cUTC', this usually being a second or two later. |
---|
| 738 | // |
---|
| 739 | // Note that the time recorded as the 'w' visibility coordinate |
---|
| 740 | // cycles through 86400 back to 0 at midnight, whereas that in 'cUTC' |
---|
| 741 | // continues to increase past 86400. |
---|
| 742 | |
---|
| 743 | double thisRA = cU; |
---|
| 744 | double thisDec = cV; |
---|
| 745 | double thisUTC = cW; |
---|
| 746 | |
---|
| 747 | if (thisUTC < prevUTC) { |
---|
| 748 | // Must have cycled through midnight. |
---|
| 749 | thisUTC += 86400.0; |
---|
| 750 | } |
---|
| 751 | |
---|
| 752 | // Guard against RA cycling through 24h in either direction. |
---|
| 753 | if (fabs(thisRA - prevRA) > PI) { |
---|
| 754 | if (thisRA < prevRA) { |
---|
| 755 | thisRA += TWOPI; |
---|
| 756 | } else { |
---|
| 757 | thisRA -= TWOPI; |
---|
| 758 | } |
---|
| 759 | } |
---|
| 760 | |
---|
| 761 | // The control system at Mopra typically does not update the |
---|
| 762 | // positions between successive integration cycles at the end of a |
---|
| 763 | // scan (nor are they flagged). In this case we use the previously |
---|
| 764 | // computed rates, even if from the previous scan since these are |
---|
| 765 | // likely to be a better guess than anything else. |
---|
| 766 | |
---|
| 767 | double dUTC = thisUTC - prevUTC; |
---|
| 768 | |
---|
| 769 | // Scan rate for this beam. |
---|
| 770 | if (dUTC > 0.0) { |
---|
| 771 | iMBuff->raRate = (thisRA - prevRA) / dUTC; |
---|
| 772 | iMBuff->decRate = (thisDec - prevDec) / dUTC; |
---|
| 773 | |
---|
| 774 | if (cInterp == 2) { |
---|
| 775 | // Use the same interpolation scheme as the original pksmbfits |
---|
| 776 | // client. This incorrectly assumed that (thisUTC - prevUTC) is |
---|
| 777 | // equal to the integration time and interpolated by computing a |
---|
| 778 | // weighted sum of the positions before and after the required |
---|
| 779 | // time. |
---|
| 780 | |
---|
| 781 | double utc = iMBuff->utc; |
---|
| 782 | if (utc - prevUTC > 100.0) { |
---|
| 783 | // Must have cycled through midnight. |
---|
| 784 | utc -= 86400.0; |
---|
| 785 | } |
---|
| 786 | |
---|
| 787 | double tw1 = 1.0 - (utc - prevUTC) / iMBuff->exposure; |
---|
| 788 | double tw2 = 1.0 - (thisUTC - utc) / iMBuff->exposure; |
---|
| 789 | double gamma = (tw2 / (tw1 + tw2)) * dUTC / (utc - prevUTC); |
---|
| 790 | |
---|
| 791 | iMBuff->raRate *= gamma; |
---|
| 792 | iMBuff->decRate *= gamma; |
---|
| 793 | } |
---|
| 794 | |
---|
| 795 | cStaleness[iBeamSel] = 0; |
---|
| 796 | |
---|
| 797 | } else { |
---|
| 798 | // Issue warnings. |
---|
| 799 | int nch = 0; |
---|
| 800 | fprintf(stderr, "WARNING, scan %d,%n cycle %d: Position ", |
---|
| 801 | iMBuff->scanNo, &nch, iMBuff->cycleNo); |
---|
| 802 | |
---|
| 803 | if (dUTC < 0.0) { |
---|
| 804 | fprintf(stderr, "timestamp went backwards!\n"); |
---|
| 805 | } else { |
---|
| 806 | if (thisRA != prevRA || thisDec != prevDec) { |
---|
| 807 | fprintf(stderr, "changed but timestamp unchanged!\n"); |
---|
| 808 | } else { |
---|
| 809 | fprintf(stderr, "and timestamp unchanged!\n"); |
---|
| 810 | } |
---|
| 811 | } |
---|
| 812 | |
---|
| 813 | cStaleness[iBeamSel]++; |
---|
| 814 | fprintf(stderr, "%-*s Using stale scan rate, staleness = %d " |
---|
| 815 | "cycle%s.\n", nch, "WARNING,", cStaleness[iBeamSel], |
---|
| 816 | (cStaleness[iBeamSel] == 1) ? "" : "s"); |
---|
| 817 | |
---|
| 818 | if (thisRA != prevRA || thisDec != prevDec) { |
---|
| 819 | if (iMBuff->raRate == 0.0 && iMBuff->decRate == 0.0) { |
---|
| 820 | fprintf(stderr, "%-*s But the previous rate was zero! " |
---|
| 821 | "Position will be inaccurate.\n", nch, "WARNING,"); |
---|
| 822 | } |
---|
| 823 | } |
---|
| 824 | } |
---|
| 825 | } |
---|
| 826 | |
---|
| 827 | // Compute the position of this beam for all bins. |
---|
| 828 | for (int idx = 0; idx < cNBin; idx++) { |
---|
| 829 | int jbuff = iBeamSel + cNBeamSel*idx; |
---|
| 830 | |
---|
| 831 | cBuffer[jbuff].raRate = iMBuff->raRate; |
---|
| 832 | cBuffer[jbuff].decRate = iMBuff->decRate; |
---|
| 833 | |
---|
| 834 | double dutc = cBuffer[jbuff].utc - prevUTC; |
---|
| 835 | if (dutc > 100.0) { |
---|
| 836 | // Must have cycled through midnight. |
---|
| 837 | dutc -= 86400.0; |
---|
| 838 | } |
---|
| 839 | |
---|
| 840 | cBuffer[jbuff].ra = prevRA + cBuffer[jbuff].raRate * dutc; |
---|
| 841 | cBuffer[jbuff].dec = prevDec + cBuffer[jbuff].decRate * dutc; |
---|
| 842 | if (cBuffer[jbuff].ra < 0.0) { |
---|
| 843 | cBuffer[jbuff].ra += TWOPI; |
---|
| 844 | } else if (cBuffer[jbuff].ra > TWOPI) { |
---|
| 845 | cBuffer[jbuff].ra -= TWOPI; |
---|
| 846 | } |
---|
| 847 | } |
---|
| 848 | } |
---|
| 849 | |
---|
| 850 | |
---|
| 851 | if (cFlushing) { |
---|
| 852 | // Copy buffer location out one IF at a time. |
---|
| 853 | MBrec.extract(*iMBuff, cFlushIF); |
---|
| 854 | haveData = 1; |
---|
| 855 | |
---|
| 856 | #ifdef PKSIO_DEBUG |
---|
| 857 | printf("Out:%4d%4d%3d%3d\n", MBrec.scanNo, MBrec.cycleNo, MBrec.beamNo, |
---|
| 858 | MBrec.IFno[0]); |
---|
| 859 | #endif |
---|
| 860 | |
---|
| 861 | // Signal that this IF in this buffer location has been flushed. |
---|
| 862 | iMBuff->IFno[cFlushIF] = 0; |
---|
| 863 | |
---|
| 864 | if (cFlushIF == iMBuff->nIF - 1) { |
---|
| 865 | // Signal that all IFs in this buffer location have been flushed. |
---|
| 866 | iMBuff->nIF = 0; |
---|
| 867 | |
---|
| 868 | // Stop cEOS being set when the next integration is read. |
---|
| 869 | iMBuff->cycleNo = 0; |
---|
| 870 | |
---|
| 871 | } else { |
---|
| 872 | // Carry on flushing the other IFs. |
---|
| 873 | continue; |
---|
| 874 | } |
---|
| 875 | |
---|
| 876 | // Has the whole buffer been flushed? |
---|
| 877 | if (cFlushBin == cNBin - 1) { |
---|
| 878 | if (cEOS || cEOF) { |
---|
| 879 | // Carry on flushing other buffers. |
---|
| 880 | cFlushIF = 0; |
---|
| 881 | continue; |
---|
| 882 | } |
---|
| 883 | |
---|
| 884 | cFlushing = 0; |
---|
| 885 | |
---|
| 886 | beamNo = int(cBaseline / 256.0); |
---|
| 887 | iBeamSel = cBeamSel[beamNo-1]; |
---|
| 888 | |
---|
| 889 | // Compute buffer number. |
---|
| 890 | iMBuff = cBuffer + iBeamSel; |
---|
| 891 | if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1); |
---|
| 892 | } |
---|
| 893 | } |
---|
| 894 | |
---|
| 895 | if (!cFlushing) { |
---|
| 896 | // Buffer this MBrec. |
---|
| 897 | if (cCycleNo == 1 && iMBuff->IFno[0]) { |
---|
| 898 | // Sanity check on the number of IFs in the new scan. |
---|
| 899 | if (if_.n_if != cNIF) { |
---|
| 900 | fprintf(stderr, "WARNING, scan %d has %d IFs instead of %d, " |
---|
| 901 | "continuing.\n", cScanNo, if_.n_if, cNIF); |
---|
| 902 | } |
---|
| 903 | } |
---|
| 904 | |
---|
[1372] | 905 | // Sanity check on incomplete integrations within a scan. |
---|
| 906 | if (iMBuff->nIF && (iMBuff->cycleNo != cCycleNo)) { |
---|
| 907 | // Force the incomplete integration to be flushed before proceeding. |
---|
| 908 | cFlushing = 1; |
---|
| 909 | continue; |
---|
| 910 | } |
---|
| 911 | |
---|
[1325] | 912 | iMBuff->scanNo = cScanNo; |
---|
| 913 | iMBuff->cycleNo = cCycleNo; |
---|
| 914 | |
---|
| 915 | // Times. |
---|
| 916 | strncpy(iMBuff->datobs, cDateObs, 10); |
---|
| 917 | iMBuff->utc = cUTC; |
---|
| 918 | iMBuff->exposure = param_.intbase; |
---|
| 919 | |
---|
| 920 | // Source identification. |
---|
| 921 | sprintf(iMBuff->srcName, "%-16.16s", |
---|
| 922 | names_.su_name + (cSrcNo-1)*16); |
---|
| 923 | iMBuff->srcRA = doubles_.su_ra[cSrcNo-1]; |
---|
| 924 | iMBuff->srcDec = doubles_.su_dec[cSrcNo-1]; |
---|
| 925 | |
---|
| 926 | // Rest frequency of the line of interest. |
---|
| 927 | iMBuff->restFreq = doubles_.rfreq; |
---|
| 928 | if (strncmp(names_.instrument, "ATPKSMB", 7) == 0) { |
---|
| 929 | // Fix the HI rest frequency recorded for Parkes multibeam data. |
---|
| 930 | double reffreq = doubles_.freq; |
---|
| 931 | double restfreq = doubles_.rfreq; |
---|
| 932 | if ((restfreq == 0.0 || fabs(restfreq - reffreq) == 0.0) && |
---|
| 933 | fabs(reffreq - 1420.40575e6) < 100.0) { |
---|
| 934 | iMBuff->restFreq = 1420.40575e6; |
---|
| 935 | } |
---|
| 936 | } |
---|
| 937 | |
---|
| 938 | // Observation type. |
---|
| 939 | int j; |
---|
| 940 | for (j = 0; j < 15; j++) { |
---|
| 941 | iMBuff->obsType[j] = names_.card[11+j]; |
---|
| 942 | if (iMBuff->obsType[j] == '\'') break; |
---|
| 943 | } |
---|
| 944 | iMBuff->obsType[j] = '\0'; |
---|
| 945 | |
---|
| 946 | // Beam-dependent parameters. |
---|
| 947 | iMBuff->beamNo = beamNo; |
---|
| 948 | |
---|
| 949 | // Beam position at the specified time. |
---|
[1372] | 950 | if (cSUpos) { |
---|
| 951 | // Non-ATNF data that does not store the position in (u,v,w). |
---|
[1325] | 952 | iMBuff->ra = doubles_.su_ra[cSrcNo-1]; |
---|
| 953 | iMBuff->dec = doubles_.su_dec[cSrcNo-1]; |
---|
| 954 | } else { |
---|
| 955 | iMBuff->ra = cU; |
---|
| 956 | iMBuff->dec = cV; |
---|
| 957 | } |
---|
| 958 | cPosUTC[iBeamSel] = cW; |
---|
| 959 | |
---|
| 960 | // IF-dependent parameters. |
---|
| 961 | int iIF = cIFno - 1; |
---|
| 962 | int startChan = cStartChan[iIF]; |
---|
| 963 | int endChan = cEndChan[iIF]; |
---|
| 964 | int refChan = cRefChan[iIF]; |
---|
| 965 | |
---|
| 966 | int nChan = abs(endChan - startChan) + 1; |
---|
| 967 | |
---|
| 968 | iIFSel = cIFSel[iIF]; |
---|
| 969 | iMBuff->nIF++; |
---|
| 970 | iMBuff->IFno[iIFSel] = cIFno; |
---|
| 971 | iMBuff->nChan[iIFSel] = nChan; |
---|
| 972 | iMBuff->nPol[iIFSel] = cNPol[iIF]; |
---|
| 973 | |
---|
| 974 | iMBuff->fqRefPix[iIFSel] = doubles_.if_ref[iIF]; |
---|
| 975 | iMBuff->fqRefVal[iIFSel] = doubles_.if_freq[iIF]; |
---|
| 976 | iMBuff->fqDelt[iIFSel] = |
---|
| 977 | if_.if_invert[iIF] * fabs(doubles_.if_bw[iIF] / |
---|
| 978 | (if_.if_nfreq[iIF] - 1)); |
---|
| 979 | |
---|
| 980 | // Adjust for channel selection. |
---|
| 981 | if (iMBuff->fqRefPix[iIFSel] != refChan) { |
---|
| 982 | iMBuff->fqRefVal[iIFSel] += |
---|
| 983 | (refChan - iMBuff->fqRefPix[iIFSel]) * |
---|
| 984 | iMBuff->fqDelt[iIFSel]; |
---|
| 985 | iMBuff->fqRefPix[iIFSel] = refChan; |
---|
| 986 | } |
---|
| 987 | |
---|
| 988 | if (endChan < startChan) { |
---|
| 989 | iMBuff->fqDelt[iIFSel] = -iMBuff->fqDelt[iIFSel]; |
---|
| 990 | } |
---|
| 991 | |
---|
| 992 | |
---|
| 993 | // System temperature. |
---|
| 994 | int iBeam = beamNo - 1; |
---|
| 995 | int scq = sc_.sc_q; |
---|
| 996 | float TsysPol1 = sc_.sc_cal[scq*iBeam + 3]; |
---|
| 997 | float TsysPol2 = sc_.sc_cal[scq*iBeam + 4]; |
---|
| 998 | iMBuff->tsys[iIFSel][0] = TsysPol1*TsysPol1; |
---|
| 999 | iMBuff->tsys[iIFSel][1] = TsysPol2*TsysPol2; |
---|
| 1000 | |
---|
| 1001 | // Calibration factor; may be changed later if the data is recalibrated. |
---|
| 1002 | if (scq > 14) { |
---|
| 1003 | // Will only be present for Parkes Multibeam or LBA data. |
---|
| 1004 | iMBuff->calfctr[iIFSel][0] = sc_.sc_cal[scq*iBeam + 14]; |
---|
| 1005 | iMBuff->calfctr[iIFSel][1] = sc_.sc_cal[scq*iBeam + 15]; |
---|
| 1006 | } else { |
---|
| 1007 | iMBuff->calfctr[iIFSel][0] = 0.0f; |
---|
| 1008 | iMBuff->calfctr[iIFSel][1] = 0.0f; |
---|
| 1009 | } |
---|
| 1010 | |
---|
| 1011 | // Cross-polarization calibration factor (unknown to MBFITS). |
---|
| 1012 | for (int j = 0; j < 2; j++) { |
---|
| 1013 | iMBuff->xcalfctr[iIFSel][j] = 0.0f; |
---|
| 1014 | } |
---|
| 1015 | |
---|
| 1016 | // Baseline parameters (unknown to MBFITS). |
---|
| 1017 | iMBuff->haveBase = 0; |
---|
| 1018 | |
---|
| 1019 | // Data (always present in MBFITS). |
---|
| 1020 | iMBuff->haveSpectra = 1; |
---|
| 1021 | |
---|
| 1022 | // Flag: bit 0 set if off source. |
---|
| 1023 | // bit 1 set if loss of sync in A polarization. |
---|
| 1024 | // bit 2 set if loss of sync in B polarization. |
---|
| 1025 | unsigned char rpflag = |
---|
| 1026 | (unsigned char)(sc_.sc_cal[scq*iBeam + 12] + 0.5f); |
---|
| 1027 | |
---|
| 1028 | // The baseline flag may be set independently. |
---|
| 1029 | if (rpflag == 0) rpflag = cFlag; |
---|
| 1030 | |
---|
| 1031 | // Copy and scale data. |
---|
| 1032 | int inc = 2 * if_.if_nstok[iIF]; |
---|
| 1033 | if (endChan < startChan) inc = -inc; |
---|
| 1034 | |
---|
| 1035 | float TsysF; |
---|
| 1036 | iMBuff->spectra[iIFSel] = iMBuff->spectra[0] + cChanOff[iIF]; |
---|
| 1037 | iMBuff->flagged[iIFSel] = iMBuff->flagged[0] + cChanOff[iIF]; |
---|
| 1038 | |
---|
| 1039 | float *spectra = iMBuff->spectra[iIFSel]; |
---|
| 1040 | unsigned char *flagged = iMBuff->flagged[iIFSel]; |
---|
| 1041 | for (int ipol = 0; ipol < cNPol[iIF]; ipol++) { |
---|
| 1042 | if (sc_.sc_cal[scq*iBeam + 3 + ipol] > 0.0f) { |
---|
| 1043 | // The correlator has already applied the calibration. |
---|
| 1044 | TsysF = 1.0f; |
---|
| 1045 | } else { |
---|
| 1046 | // The correlator has normalized cVis[k] to a Tsys of 500K. |
---|
| 1047 | TsysF = iMBuff->tsys[iIFSel][ipol] / 500.0f; |
---|
| 1048 | } |
---|
| 1049 | |
---|
| 1050 | int k = 2 * (if_.if_nstok[iIF]*(startChan - 1) + ipol); |
---|
| 1051 | for (int ichan = 0; ichan < nChan; ichan++) { |
---|
| 1052 | *(spectra++) = TsysF * cVis[k]; |
---|
| 1053 | *(flagged++) = rpflag; |
---|
| 1054 | k += inc; |
---|
| 1055 | } |
---|
| 1056 | } |
---|
| 1057 | |
---|
| 1058 | if (cHaveXPol[iIF]) { |
---|
| 1059 | int k = 2 * (3*(startChan - 1) + 2); |
---|
| 1060 | iMBuff->xpol[iIFSel] = iMBuff->xpol[0] + cXpolOff[iIF]; |
---|
| 1061 | float *xpol = iMBuff->xpol[iIFSel]; |
---|
| 1062 | for (int ichan = 0; ichan < nChan; ichan++) { |
---|
| 1063 | *(xpol++) = cVis[k]; |
---|
| 1064 | *(xpol++) = cVis[k+1]; |
---|
| 1065 | k += inc; |
---|
| 1066 | } |
---|
| 1067 | } |
---|
| 1068 | |
---|
| 1069 | |
---|
| 1070 | // Parallactic angle. |
---|
| 1071 | iMBuff->parAngle = sc_.sc_cal[scq*iBeam + 11]; |
---|
| 1072 | |
---|
| 1073 | // Calibration factor applied to the data by the correlator. |
---|
| 1074 | if (scq > 14) { |
---|
| 1075 | // Will only be present for Parkes Multibeam or LBA data. |
---|
| 1076 | iMBuff->tcal[iIFSel][0] = sc_.sc_cal[scq*iBeam + 14]; |
---|
| 1077 | iMBuff->tcal[iIFSel][1] = sc_.sc_cal[scq*iBeam + 15]; |
---|
| 1078 | } else { |
---|
| 1079 | iMBuff->tcal[iIFSel][0] = 0.0f; |
---|
| 1080 | iMBuff->tcal[iIFSel][1] = 0.0f; |
---|
| 1081 | } |
---|
| 1082 | |
---|
| 1083 | if (sc_.sc_ant <= anten_.nant) { |
---|
| 1084 | // No extra syscal information present. |
---|
| 1085 | iMBuff->extraSysCal = 0; |
---|
| 1086 | iMBuff->azimuth = 0.0f; |
---|
| 1087 | iMBuff->elevation = 0.0f; |
---|
| 1088 | iMBuff->parAngle = 0.0f; |
---|
| 1089 | iMBuff->focusAxi = 0.0f; |
---|
| 1090 | iMBuff->focusTan = 0.0f; |
---|
| 1091 | iMBuff->focusRot = 0.0f; |
---|
| 1092 | iMBuff->temp = 0.0f; |
---|
| 1093 | iMBuff->pressure = 0.0f; |
---|
| 1094 | iMBuff->humidity = 0.0f; |
---|
| 1095 | iMBuff->windSpeed = 0.0f; |
---|
| 1096 | iMBuff->windAz = 0.0f; |
---|
| 1097 | strcpy(iMBuff->tcalTime, " "); |
---|
| 1098 | iMBuff->refBeam = 0; |
---|
| 1099 | |
---|
| 1100 | } else { |
---|
| 1101 | // Additional information for Parkes Multibeam data. |
---|
| 1102 | int iOff = scq*(sc_.sc_ant - 1) - 1; |
---|
| 1103 | iMBuff->extraSysCal = 1; |
---|
| 1104 | iMBuff->azimuth = sc_.sc_cal[iOff + 2]; |
---|
| 1105 | iMBuff->elevation = sc_.sc_cal[iOff + 3]; |
---|
| 1106 | iMBuff->parAngle = sc_.sc_cal[iOff + 4]; |
---|
| 1107 | iMBuff->focusAxi = sc_.sc_cal[iOff + 5] * 1e-3; |
---|
| 1108 | iMBuff->focusTan = sc_.sc_cal[iOff + 6] * 1e-3; |
---|
| 1109 | iMBuff->focusRot = sc_.sc_cal[iOff + 7]; |
---|
| 1110 | iMBuff->temp = sc_.sc_cal[iOff + 8]; |
---|
| 1111 | iMBuff->pressure = sc_.sc_cal[iOff + 9]; |
---|
| 1112 | iMBuff->humidity = sc_.sc_cal[iOff + 10]; |
---|
| 1113 | iMBuff->windSpeed = sc_.sc_cal[iOff + 11]; |
---|
| 1114 | iMBuff->windAz = sc_.sc_cal[iOff + 12]; |
---|
| 1115 | |
---|
| 1116 | char *tcalTime = iMBuff->tcalTime; |
---|
| 1117 | sprintf(tcalTime, "%-16.16s", (char *)(&sc_.sc_cal[iOff+13])); |
---|
| 1118 | |
---|
| 1119 | #ifndef AIPS_LITTLE_ENDIAN |
---|
| 1120 | // Do byte swapping on the ASCII date string. |
---|
| 1121 | for (int j = 0; j < 16; j += 4) { |
---|
| 1122 | char ctmp; |
---|
| 1123 | ctmp = tcalTime[j]; |
---|
| 1124 | tcalTime[j] = tcalTime[j+3]; |
---|
| 1125 | tcalTime[j+3] = ctmp; |
---|
| 1126 | ctmp = tcalTime[j+1]; |
---|
| 1127 | tcalTime[j+1] = tcalTime[j+2]; |
---|
| 1128 | tcalTime[j+2] = ctmp; |
---|
| 1129 | } |
---|
| 1130 | #endif |
---|
| 1131 | |
---|
| 1132 | // Reference beam number. |
---|
| 1133 | float refbeam = sc_.sc_cal[iOff + 17]; |
---|
| 1134 | if (refbeam > 0.0f || refbeam < 100.0f) { |
---|
| 1135 | iMBuff->refBeam = int(refbeam); |
---|
| 1136 | } else { |
---|
| 1137 | iMBuff->refBeam = 0; |
---|
| 1138 | } |
---|
| 1139 | } |
---|
| 1140 | } |
---|
| 1141 | } |
---|
| 1142 | |
---|
| 1143 | return 0; |
---|
| 1144 | } |
---|
| 1145 | |
---|
| 1146 | //-------------------------------------------------------- MBFITSreader::rpget |
---|
| 1147 | |
---|
| 1148 | // Read the next data record from the RPFITS file. |
---|
| 1149 | |
---|
| 1150 | int MBFITSreader::rpget(int syscalonly, int &EOS) |
---|
| 1151 | { |
---|
| 1152 | EOS = 0; |
---|
| 1153 | |
---|
| 1154 | int retries = 0; |
---|
| 1155 | |
---|
| 1156 | // Allow 10 read errors. |
---|
| 1157 | int numErr = 0; |
---|
| 1158 | |
---|
| 1159 | int jstat = 0; |
---|
| 1160 | while (numErr < 10) { |
---|
| 1161 | int lastjstat = jstat; |
---|
| 1162 | rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag, |
---|
| 1163 | &cBin, &cIFno, &cSrcNo); |
---|
| 1164 | |
---|
| 1165 | switch(jstat) { |
---|
| 1166 | case -1: |
---|
| 1167 | // Read failed; retry. |
---|
| 1168 | numErr++; |
---|
| 1169 | fprintf(stderr, "RPFITS read failed - retrying.\n"); |
---|
| 1170 | jstat = 0; |
---|
| 1171 | break; |
---|
| 1172 | |
---|
| 1173 | case 0: |
---|
| 1174 | // Successful read. |
---|
| 1175 | if (lastjstat == 0) { |
---|
| 1176 | if (cBaseline == -1) { |
---|
| 1177 | // Syscal data. |
---|
| 1178 | if (syscalonly) { |
---|
| 1179 | return 0; |
---|
| 1180 | } |
---|
| 1181 | |
---|
| 1182 | } else { |
---|
| 1183 | if (!syscalonly) { |
---|
| 1184 | return 0; |
---|
| 1185 | } |
---|
| 1186 | } |
---|
| 1187 | } |
---|
| 1188 | |
---|
| 1189 | // Last operation was to read header or FG table; now read data. |
---|
| 1190 | break; |
---|
| 1191 | |
---|
| 1192 | case 1: |
---|
| 1193 | // Encountered header while trying to read data; read it. |
---|
| 1194 | EOS = 1; |
---|
| 1195 | jstat = -1; |
---|
| 1196 | break; |
---|
| 1197 | |
---|
| 1198 | case 2: |
---|
| 1199 | // End of scan; read past it. |
---|
| 1200 | jstat = 0; |
---|
| 1201 | break; |
---|
| 1202 | |
---|
| 1203 | case 3: |
---|
| 1204 | // End-of-file; retry applies to real-time mode. |
---|
| 1205 | if (retries++ >= cRetry) { |
---|
| 1206 | return -1; |
---|
| 1207 | } |
---|
| 1208 | |
---|
| 1209 | sleep(10); |
---|
| 1210 | jstat = 0; |
---|
| 1211 | break; |
---|
| 1212 | |
---|
| 1213 | case 4: |
---|
| 1214 | // Encountered FG table while trying to read data; read it. |
---|
| 1215 | jstat = -1; |
---|
| 1216 | break; |
---|
| 1217 | |
---|
| 1218 | case 5: |
---|
| 1219 | // Illegal data at end of block after close/reopen operation; retry. |
---|
| 1220 | jstat = 0; |
---|
| 1221 | break; |
---|
| 1222 | |
---|
| 1223 | default: |
---|
| 1224 | // Shouldn't reach here. |
---|
| 1225 | fprintf(stderr, "Unrecognized RPFITSIN return code: %d (retrying)\n", |
---|
| 1226 | jstat); |
---|
| 1227 | jstat = 0; |
---|
| 1228 | break; |
---|
| 1229 | } |
---|
| 1230 | } |
---|
| 1231 | |
---|
[1372] | 1232 | fprintf(stderr, "ERROR, RPFITS read failed too many times.\n"); |
---|
[1325] | 1233 | return 2; |
---|
| 1234 | } |
---|
| 1235 | |
---|
| 1236 | //-------------------------------------------------------- MBFITSreader::close |
---|
| 1237 | |
---|
| 1238 | // Close the input file. |
---|
| 1239 | |
---|
| 1240 | void MBFITSreader::close(void) |
---|
| 1241 | { |
---|
| 1242 | if (cMBopen) { |
---|
| 1243 | int jstat = 1; |
---|
| 1244 | rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag, |
---|
| 1245 | &cBin, &cIFno, &cSrcNo); |
---|
| 1246 | |
---|
| 1247 | if (cBeams) delete [] cBeams; |
---|
| 1248 | if (cIFs) delete [] cIFs; |
---|
| 1249 | if (cNChan) delete [] cNChan; |
---|
| 1250 | if (cNPol) delete [] cNPol; |
---|
| 1251 | if (cHaveXPol) delete [] cHaveXPol; |
---|
| 1252 | if (cStartChan) delete [] cStartChan; |
---|
| 1253 | if (cEndChan) delete [] cEndChan; |
---|
| 1254 | if (cRefChan) delete [] cRefChan; |
---|
| 1255 | |
---|
| 1256 | if (cVis) delete [] cVis; |
---|
| 1257 | if (cWgt) delete [] cWgt; |
---|
| 1258 | |
---|
| 1259 | if (cBeamSel) delete [] cBeamSel; |
---|
| 1260 | if (cIFSel) delete [] cIFSel; |
---|
| 1261 | if (cChanOff) delete [] cChanOff; |
---|
| 1262 | if (cXpolOff) delete [] cXpolOff; |
---|
| 1263 | if (cBuffer) delete [] cBuffer; |
---|
| 1264 | if (cPosUTC) delete [] cPosUTC; |
---|
| 1265 | |
---|
| 1266 | cMBopen = 0; |
---|
| 1267 | } |
---|
| 1268 | } |
---|