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