source: branches/asap-3.x/external/atnf/PKSIO/MBFITSreader.cc@ 2441

Last change on this file since 2441 was 1720, checked in by Malte Marquarding, 15 years ago

Update from livedata CVS repository
File size: 55.7 KB
Line 
1//#---------------------------------------------------------------------------
2//# MBFITSreader.cc: ATNF single-dish RPFITS reader.
3//#---------------------------------------------------------------------------
4//# livedata - processing pipeline for single-dish, multibeam spectral data.
5//# Copyright (C) 2000-2009, Australia Telescope National Facility, CSIRO
6//#
7//# This file is part of livedata.
8//#
9//# livedata is free software: you can redistribute it and/or modify it under
10//# the terms of the GNU General Public License as published by the Free
11//# Software Foundation, either version 3 of the License, or (at your option)
12//# any later version.
13//#
14//# livedata is distributed in the hope that it will be useful, but WITHOUT
15//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
17//# more details.
18//#
19//# You should have received a copy of the GNU General Public License along
20//# with livedata. If not, see <http://www.gnu.org/licenses/>.
21//#
22//# Correspondence concerning livedata may be directed to:
23//# Internet email: mcalabre@atnf.csiro.au
24//# Postal address: Dr. Mark Calabretta
25//# Australia Telescope National Facility, CSIRO
26//# PO Box 76
27//# Epping NSW 1710
28//# AUSTRALIA
29//#
30//# http://www.atnf.csiro.au/computing/software/livedata.html
31//# $Id: MBFITSreader.cc,v 19.57 2009-10-30 06:34:36 cal103 Exp $
32//#---------------------------------------------------------------------------
33//# The MBFITSreader class reads single dish RPFITS files (such as Parkes
34//# Multibeam MBFITS files).
35//#
36//# Original: 2000/07/28 Mark Calabretta
37//#---------------------------------------------------------------------------
38
39#include <atnf/pks/pks_maths.h>
40#include <atnf/PKSIO/MBFITSreader.h>
41#include <atnf/PKSIO/MBrecord.h>
42
43#include <casa/math.h>
44#include <casa/iostream.h>
45#include <casa/stdio.h>
46#include <casa/stdlib.h>
47#include <casa/string.h>
48#include <unistd.h>
49
50#include <RPFITS.h>
51
52using namespace std;
53
54// Numerical constants.
55const double PI = 3.141592653589793238462643;
56const double TWOPI = 2.0 * PI;
57const double HALFPI = PI / 2.0;
58const double R2D = 180.0 / PI;
59
60//------------------------------------------------- MBFITSreader::MBFITSreader
61
62// Default constructor.
63
64MBFITSreader::MBFITSreader(
65 const int retry,
66 const int interpolate)
67{
68 cRetry = retry;
69 if (cRetry > 10) {
70 cRetry = 10;
71 }
72
73 cInterp = interpolate;
74 if (cInterp < 0 || cInterp > 2) {
75 cInterp = 1;
76 }
77
78 // Initialize pointers.
79 cBeams = 0x0;
80 cIFs = 0x0;
81 cNChan = 0x0;
82 cNPol = 0x0;
83 cHaveXPol = 0x0;
84 cStartChan = 0x0;
85 cEndChan = 0x0;
86 cRefChan = 0x0;
87
88 cVis = 0x0;
89 cWgt = 0x0;
90
91 cBeamSel = 0x0;
92 cIFSel = 0x0;
93 cChanOff = 0x0;
94 cXpolOff = 0x0;
95 cBuffer = 0x0;
96 cPosUTC = 0x0;
97
98 cMBopen = 0;
99
100 // Tell RPFITSIN not to report errors directly.
101 iostat_.errlun = -1;
102
103 // By default, messages are written to stderr.
104 initMsg();
105}
106
107//------------------------------------------------ MBFITSreader::~MBFITSreader
108
109// Destructor.
110
111MBFITSreader::~MBFITSreader()
112{
113 close();
114}
115
116//--------------------------------------------------------- MBFITSreader::open
117
118// Open the RPFITS file for reading.
119
120int MBFITSreader::open(
121 char *rpname,
122 int &nBeam,
123 int* &beams,
124 int &nIF,
125 int* &IFs,
126 int* &nChan,
127 int* &nPol,
128 int* &haveXPol,
129 int &haveBase,
130 int &haveSpectra,
131 int &extraSysCal)
132{
133 // Clear the message stack.
134 clearMsg();
135
136 if (cMBopen) {
137 close();
138 }
139
140 strcpy(names_.file, rpname);
141
142 // Open the RPFITS file.
143 int jstat = -3;
144 if (rpfitsin(jstat)) {
145 sprintf(cMsg, "ERROR: Failed to open MBFITS file\n %s", rpname);
146 logMsg(cMsg);
147 return 1;
148 }
149
150 cMBopen = 1;
151
152 // Tell RPFITSIN that we want the OBSTYPE card.
153 int j;
154 param_.ncard = 1;
155 for (j = 0; j < 80; j++) {
156 names_.card[j] = ' ';
157 }
158 strncpy(names_.card, "OBSTYPE", 7);
159
160 // Read the first header.
161 jstat = -1;
162 if (rpfitsin(jstat)) {
163 sprintf(cMsg, "ERROR: Failed to read MBFITS header in file\n"
164 " %s", rpname);
165 logMsg(cMsg);
166 close();
167 return 1;
168 }
169
170 // Mopra data has some peculiarities.
171 cMopra = strncmp(names_.instrument, "ATMOPRA", 7) == 0;
172
173 // Non-ATNF data may not store the position in (u,v,w).
174 if (strncmp(names_.sta, "tid", 3) == 0) {
175 sprintf(cMsg, "WARNING: Found Tidbinbilla data");
176 cSUpos = 1;
177 } else if (strncmp(names_.sta, "HOB", 3) == 0) {
178 sprintf(cMsg, "WARNING: Found Hobart data");
179 cSUpos = 1;
180 } else if (strncmp(names_.sta, "CED", 3) == 0) {
181 sprintf(cMsg, "WARNING: Found Ceduna data");
182 cSUpos = 1;
183 } else {
184 cSUpos = 0;
185 }
186
187 if (cSUpos) {
188 strcat(cMsg, ", using telescope position\n from SU table.");
189 logMsg(cMsg);
190 cInterp = 0;
191 }
192
193 // Mean scan rate (for timestamp repairs).
194 cNRate = 0;
195 cAvRate[0] = 0.0;
196 cAvRate[1] = 0.0;
197 cCode5 = 0;
198
199
200 // Find the maximum beam number.
201 cNBeam = 0;
202 for (int iBeam = 0; iBeam < anten_.nant; iBeam++) {
203 if (anten_.ant_num[iBeam] > cNBeam) {
204 cNBeam = anten_.ant_num[iBeam];
205 }
206 }
207
208 if (cNBeam <= 0) {
209 logMsg("ERROR: Couldn't determine number of beams.");
210 close();
211 return 1;
212 }
213
214 // Construct the beam mask.
215 cBeams = new int[cNBeam];
216 for (int iBeam = 0; iBeam < cNBeam; iBeam++) {
217 cBeams[iBeam] = 0;
218 }
219
220 // ...beams present in the data.
221 for (int iBeam = 0; iBeam < anten_.nant; iBeam++) {
222 // Guard against dubious beam numbers, e.g. zeroes in
223 // 1999-09-29_1632_024848p14_071b.hpf and the four scans following.
224 // Note that the actual beam number is decoded from the 'baseline' random
225 // parameter for each spectrum and is only used for beam selection.
226 int beamNo = anten_.ant_num[iBeam];
227 if (beamNo != iBeam+1) {
228 char sta[8];
229 strncpy(sta, names_.sta+(8*iBeam), 8);
230 char *cp = sta + 7;
231 while (*cp == ' ') *(cp--) = '\0';
232
233 sprintf(cMsg,
234 "WARNING: RPFITSIN returned beam number %2d for AN table\n"
235 " entry %2d with name '%.8s'", beamNo, iBeam+1, sta);
236
237 char text[8];
238 sprintf(text, "MB%2.2d", iBeam+1);
239 cp = cMsg + strlen(cMsg);
240 if (strncmp(sta, text, 8) == 0) {
241 beamNo = iBeam + 1;
242 sprintf(cp, "; using beam number %2d.", beamNo);
243 } else {
244 sprintf(cp, ".");
245 }
246
247 logMsg(cMsg);
248 }
249
250 if (0 < beamNo && beamNo <= cNBeam) {
251 cBeams[beamNo-1] = 1;
252 }
253 }
254
255 // Passing back the address of the array allows PKSFITSreader::select() to
256 // modify its elements directly.
257 nBeam = cNBeam;
258 beams = cBeams;
259
260
261 // Number of IFs.
262 cNIF = if_.n_if;
263 cIFs = new int[cNIF];
264 for (int iIF = 0; iIF < cNIF; iIF++) {
265 cIFs[iIF] = 1;
266 }
267
268 // Passing back the address of the array allows PKSFITSreader::select() to
269 // modify its elements directly.
270 nIF = cNIF;
271 IFs = cIFs;
272
273
274 // Number of channels and polarizations.
275 cNChan = new int[cNIF];
276 cNPol = new int[cNIF];
277 cHaveXPol = new int[cNIF];
278 cGetXPol = 0;
279
280 int maxProd = 0;
281 for (int iIF = 0; iIF < cNIF; iIF++) {
282 cNChan[iIF] = if_.if_nfreq[iIF];
283 cNPol[iIF] = if_.if_nstok[iIF];
284 cNChan[iIF] -= cNChan[iIF]%2;
285
286 // Do we have cross-polarization data?
287 if ((cHaveXPol[iIF] = cNPol[iIF] > 2)) {
288 // Cross-polarization data is handled separately.
289 cNPol[iIF] = 2;
290
291 // Default is to get it if we have it.
292 cGetXPol = 1;
293 }
294
295 // Maximum number of spectral products in any IF.
296 int nProd = if_.if_nfreq[iIF] * if_.if_nstok[iIF];
297 if (maxProd < nProd) maxProd = nProd;
298 }
299
300 // Allocate memory for RPFITSIN subroutine arguments.
301 if (cVis) delete [] cVis;
302 if (cWgt) delete [] cWgt;
303 cVis = new float[2*maxProd];
304 cWgt = new float[maxProd];
305
306 nChan = cNChan;
307 nPol = cNPol;
308 haveXPol = cHaveXPol;
309
310
311 // Default channel range selection.
312 cStartChan = new int[cNIF];
313 cEndChan = new int[cNIF];
314 cRefChan = new int[cNIF];
315
316 for (int iIF = 0; iIF < cNIF; iIF++) {
317 cStartChan[iIF] = 1;
318 cEndChan[iIF] = cNChan[iIF];
319 cRefChan[iIF] = cNChan[iIF]/2 + 1;
320 }
321
322 cGetSpectra = 1;
323
324
325 // No baseline parameters in MBFITS.
326 haveBase = 0;
327
328 // Always have spectra in MBFITS.
329 haveSpectra = cHaveSpectra = 1;
330
331
332 // Integration cycle time (s).
333 cIntTime = param_.intime;
334
335 // Can't deduce binning mode till later.
336 cNBin = 0;
337
338
339 // Read the first syscal record.
340 if (rpget(1, cEOS)) {
341 logMsg("ERROR: Failed to read first syscal record.");
342 close();
343 return 1;
344 }
345
346 // Additional information for Parkes Multibeam data?
347 extraSysCal = (sc_.sc_ant > anten_.nant);
348
349
350 cFirst = 1;
351 cEOF = 0;
352 cFlushing = 0;
353
354 return 0;
355}
356
357//---------------------------------------------------- MBFITSreader::getHeader
358
359// Get parameters describing the data.
360
361int MBFITSreader::getHeader(
362 char observer[32],
363 char project[32],
364 char telescope[32],
365 double antPos[3],
366 char obsType[32],
367 char bunit[32],
368 float &equinox,
369 char radecsys[32],
370 char dopplerFrame[32],
371 char datobs[32],
372 double &utc,
373 double &refFreq,
374 double &bandwidth)
375{
376 if (!cMBopen) {
377 logMsg("ERROR: An MBFITS file has not been opened.");
378 return 1;
379 }
380
381 sprintf(observer, "%-16.16s", names_.rp_observer);
382 sprintf(project, "%-16.16s", names_.object);
383 sprintf(telescope, "%-16.16s", names_.instrument);
384
385 // Observatory coordinates (ITRF), in m.
386 antPos[0] = doubles_.x[0];
387 antPos[1] = doubles_.y[0];
388 antPos[2] = doubles_.z[0];
389
390 // This is the only sure way to identify the telescope, maybe.
391 if (strncmp(names_.sta, "MB0", 3) == 0) {
392 // Parkes Multibeam.
393 sprintf(telescope, "%-16.16s", "ATPKSMB");
394 antPos[0] = -4554232.087;
395 antPos[1] = 2816759.046;
396 antPos[2] = -3454035.950;
397
398 } else if (strncmp(names_.sta, "HOH", 3) == 0) {
399 // Parkes HOH receiver.
400 sprintf(telescope, "%-16.16s", "ATPKSHOH");
401 antPos[0] = -4554232.087;
402 antPos[1] = 2816759.046;
403 antPos[2] = -3454035.950;
404
405 } else if (strncmp(names_.sta, "CA0", 3) == 0) {
406 // An ATCA antenna, use the array centre position.
407 sprintf(telescope, "%-16.16s", "ATCA");
408 antPos[0] = -4750915.837;
409 antPos[1] = 2792906.182;
410 antPos[2] = -3200483.747;
411
412 // ATCA-104. Updated position at epoch 2007/06/24 from Chris Phillips.
413 // antPos[0] = -4751640.182; // ± 0.008
414 // antPos[1] = 2791700.322; // ± 0.006
415 // antPos[2] = -3200490.668; // ± 0.007
416 //
417 } else if (strncmp(names_.sta, "MOP", 3) == 0) {
418 // Mopra. Updated position at epoch 2007/06/24 from Chris Phillips.
419 sprintf(telescope, "%-16.16s", "ATMOPRA");
420 antPos[0] = -4682769.444; // ± 0.009
421 antPos[1] = 2802618.963; // ± 0.006
422 antPos[2] = -3291758.864; // ± 0.008
423
424 } else if (strncmp(names_.sta, "HOB", 3) == 0) {
425 // Hobart.
426 sprintf(telescope, "%-16.16s", "HOBART");
427 antPos[0] = -3950236.735;
428 antPos[1] = 2522347.567;
429 antPos[2] = -4311562.569;
430
431 } else if (strncmp(names_.sta, "CED", 3) == 0) {
432 // Ceduna. Updated position at epoch 2007/06/24 from Chris Phillips.
433 sprintf(telescope, "%-16.16s", "CEDUNA");
434 antPos[0] = -3753443.168; // ± 0.017
435 antPos[1] = 3912709.794; // ± 0.017
436 antPos[2] = -3348067.060; // ± 0.016
437
438 } else if (strncmp(names_.sta, "tid", 3) == 0) {
439 // DSS.
440 sprintf(telescope, "%-16.16s", "DSS-43");
441 antPos[0] = -4460894.727;
442 antPos[1] = 2682361.530;
443 antPos[2] = -3674748.424;
444 }
445
446 // Observation type.
447 int j;
448 for (j = 0; j < 31; j++) {
449 obsType[j] = names_.card[11+j];
450 if (obsType[j] == '\'') break;
451 }
452 obsType[j] = '\0';
453
454 // Brightness unit.
455 sprintf(bunit, "%-16.16s", names_.bunit);
456 if (strcmp(bunit, "JY") == 0) {
457 bunit[1] = 'y';
458 } else if (strcmp(bunit, "JY/BEAM") == 0) {
459 strcpy(bunit, "Jy/beam");
460 }
461
462 // Coordinate frames.
463 equinox = 2000.0f;
464 strcpy(radecsys, "FK5");
465 strcpy(dopplerFrame, "TOPOCENT");
466
467 // Time at start of observation.
468 sprintf(datobs, "%-10.10s", names_.datobs);
469 utc = cUTC;
470
471 // Spectral parameters.
472 refFreq = doubles_.if_freq[0];
473 bandwidth = doubles_.if_bw[0];
474
475 return 0;
476}
477
478//-------------------------------------------------- MBFITSreader::getFreqInfo
479
480// Get frequency parameters for each IF.
481
482int MBFITSreader::getFreqInfo(
483 int &nIF,
484 double* &startFreq,
485 double* &endFreq)
486{
487 // This is RPFITS - can't do it!
488 return 1;
489}
490
491//---------------------------------------------------- MBFITSreader::findRange
492
493// Find the range of the data selected in time and position.
494
495int MBFITSreader::findRange(
496 int &nRow,
497 int &nSel,
498 char dateSpan[2][32],
499 double utcSpan[2],
500 double* &positions)
501{
502 // This is RPFITS - can't do it!
503 return 1;
504}
505
506//--------------------------------------------------------- MBFITSreader::read
507
508// Read the next data record (if you're feeling lucky).
509
510int MBFITSreader::read(
511 MBrecord &MBrec)
512{
513 int beamNo = -1;
514 int haveData, pCode = 0, status;
515 double raRate = 0.0, decRate = 0.0, paRate = 0.0;
516 MBrecord *iMBuff = 0x0;
517
518 if (!cMBopen) {
519 logMsg("ERROR: An MBFITS file has not been opened.");
520 return 1;
521 }
522
523 // Positions recorded in the input records usually do not coincide with the
524 // midpoint of the integration and hence the input must be buffered so that
525 // true positions may be interpolated.
526 //
527 // On the first call nBeamSel buffers of length nBin, are allocated and
528 // filled, where nBin is the number of time bins.
529 //
530 // The input records for binned, single beam data with multiple simultaneous
531 // IFs are ordered by IF within each integration rather than by bin number
532 // and hence are not in time order. No multibeam data exists with
533 // nBin > 1 but the likelihood that the input records would be in beam/IF
534 // order and the requirement that output records be in time order would
535 // force an elaborate double-buffering system and we do not support it.
536 //
537 // Once all buffers are filled, the next record for each beam pertains to
538 // the next integration and should contain new position information allowing
539 // the proper position for each spectrum in the buffer to be interpolated.
540 // The buffers are then flushed in time order. For single beam data there
541 // is only one buffer and reads from the MBFITS file are suspended while the
542 // flush is in progress. For multibeam data each buffer is of unit length
543 // so the flush completes immediately and the new record takes its place.
544
545 haveData = 0;
546 while (!haveData) {
547 int iBeamSel = -1, iIFSel = -1;
548
549 if (!cFlushing) {
550 if (cEOF) {
551 return -1;
552 }
553
554 // Read the next record.
555 pCode = 0;
556 if ((status = rpget(0, cEOS)) == -1) {
557 // EOF.
558 cEOF = 1;
559 cFlushing = 1;
560 cFlushBin = 0;
561 cFlushIF = 0;
562
563#ifdef PKSIO_DEBUG
564 fprintf(stderr, "\nEnd-of-file detected, flushing last cycle.\n");
565#endif
566
567 } else if (status) {
568 // IO error.
569 return 1;
570
571 } else {
572 if (cFirst) {
573 // First data; cBeamSel[] stores the buffer index for each beam.
574 cNBeamSel = 0;
575 cBeamSel = new int[cNBeam];
576
577 for (int iBeam = 0; iBeam < cNBeam; iBeam++) {
578 if (cBeams[iBeam]) {
579 // Buffer offset for this beam.
580 cBeamSel[iBeam] = cNBeamSel++;
581 } else {
582 // Signal that the beam is not selected.
583 cBeamSel[iBeam] = -1;
584 }
585 }
586
587 // Set up bookkeeping arrays for IFs.
588 cIFSel = new int[cNIF];
589 cChanOff = new int[cNIF];
590 cXpolOff = new int[cNIF];
591
592 int maxChan = 0;
593 int maxXpol = 0;
594
595 cSimulIF = 0;
596 for (int iIF = 0; iIF < cNIF; iIF++) {
597 if (cIFs[iIF]) {
598 // Buffer index for each IF within each simultaneous set.
599 cIFSel[iIF] = 0;
600
601 // Array offsets for each IF within each simultaneous set.
602 cChanOff[iIF] = 0;
603 cXpolOff[iIF] = 0;
604
605 // Look for earlier IFs in the same simultaneous set.
606 for (int jIF = 0; jIF < iIF; jIF++) {
607 if (!cIFs[jIF]) continue;
608
609 if (if_.if_simul[jIF] == if_.if_simul[iIF]) {
610 // Got one, increment indices.
611 cIFSel[iIF]++;
612
613 cChanOff[iIF] += cNChan[jIF] * cNPol[jIF];
614 if (cHaveXPol[jIF]) {
615 cXpolOff[iIF] += 2 * cNChan[jIF];
616 }
617 }
618 }
619
620 // Maximum number of selected IFs in any simultaneous set.
621 cSimulIF = max(cSimulIF, cIFSel[iIF]+1);
622
623 // Maximum memory required for any simultaneous set.
624 maxChan = max(maxChan, cChanOff[iIF] + cNChan[iIF]*cNPol[iIF]);
625 if (cHaveXPol[iIF]) {
626 maxXpol = max(maxXpol, cXpolOff[iIF] + 2*cNChan[iIF]);
627 }
628
629 } else {
630 // Signal that the IF is not selected.
631 cIFSel[iIF] = -1;
632 }
633 }
634
635 // Check for binning mode observations.
636 if (param_.intbase > 0.0f) {
637 cNBin = int((cIntTime / param_.intbase) + 0.5);
638
639 // intbase sometimes contains rubbish.
640 if (cNBin == 0) {
641 cNBin = 1;
642 }
643 } else {
644 cNBin = 1;
645 }
646
647 if (cNBin > 1 && cNBeamSel > 1) {
648 logMsg("ERROR: Cannot handle binning mode for multiple beams.\n"
649 " Select a single beam for input.");
650 close();
651 return 1;
652 }
653
654 // Allocate buffer data storage; the MBrecord constructor zeroes
655 // class members such as cycleNo that are tested in the first pass
656 // below.
657 int nBuff = cNBeamSel * cNBin;
658 cBuffer = new MBrecord[nBuff];
659
660 // Allocate memory for spectral arrays.
661 for (int ibuff = 0; ibuff < nBuff; ibuff++) {
662 cBuffer[ibuff].setNIFs(cSimulIF);
663 cBuffer[ibuff].allocate(0, maxChan, maxXpol);
664
665 // Signal that this IF in this buffer has been flushed.
666 for (int iIF = 0; iIF < cSimulIF; iIF++) {
667 cBuffer[ibuff].IFno[iIF] = 0;
668 }
669 }
670
671 cPosUTC = new double[cNBeamSel];
672
673 cFirst = 0;
674 cScanNo = 1;
675 cCycleNo = 0;
676 cPrevUTC = -1.0;
677 }
678
679 // Check for end-of-scan.
680 if (cEOS) {
681 cScanNo++;
682 cCycleNo = 0;
683 cPrevUTC = -1.0;
684 }
685
686 // Apply beam and IF selection before the change-of-day test to allow
687 // a single selected beam and IF to be handled in binning-mode.
688 beamNo = int(cBaseline / 256.0);
689 if (beamNo == 1) {
690 // Store the position of beam 1 for grid convergence corrections.
691 cRA0 = cU;
692 cDec0 = cV;
693 }
694 iBeamSel = cBeamSel[beamNo-1];
695 if (iBeamSel < 0) continue;
696
697 // Sanity check (mainly for MOPS).
698 if (cIFno > cNIF) continue;
699
700 // Apply IF selection; iIFSel == 0 for the first selected IF, == 1
701 // for the second, etc.
702 iIFSel = cIFSel[cIFno - 1];
703 if (iIFSel < 0) continue;
704
705
706 if (cNBin > 1) {
707 // Binning mode: correct the time.
708 cUTC += param_.intbase * (cBin - (cNBin + 1)/2.0);
709 }
710
711 // Check for change-of-day.
712 double cod = 0.0;
713 if ((cUTC + 86400.0) < (cPrevUTC + 600.0)) {
714 // cUTC should continue to increase past 86400 during a single scan.
715 // However, if the RPFITS file contains multiple scans that straddle
716 // midnight then cUTC can jump backwards from the end of one scan to
717 // the start of the next.
718#ifdef PKSIO_DEBUG
719 fprintf(stderr, "Change-of-day on cUTC: %.1f -> %.1f\n",
720 cPrevUTC, cUTC);
721#endif
722 // Can't change the recorded value of cUTC directly (without also
723 // changing dateobs) so change-of-day must be recorded separately as
724 // an offset to be applied when comparing integration timestamps.
725 cod = 86400.0;
726 }
727
728 if ((cUTC+cod) < cPrevUTC - 1.0) {
729 if (cBin == 1 && iIFSel) {
730 // Multiple-IF, binning-mode data is only partially time ordered.
731#ifdef PKSIO_DEBUG
732 fprintf(stderr, "New IF in multiple-IF, binning-mode data.\n");
733#endif
734 cCycleNo -= cNBin;
735 cPrevUTC = -1.0;
736
737 } else {
738 // All other data should be fully time ordered.
739 sprintf(cMsg,
740 "WARNING: Cycle %d:%03d-%03d, UTC went backwards from\n"
741 " %.1f to %.1f! Incrementing day number,\n"
742 " positions may be unreliable.", cScanNo, cCycleNo,
743 cCycleNo+1, cPrevUTC, cUTC);
744 logMsg(cMsg);
745 cUTC += 86400.0;
746 }
747 }
748
749 // New integration cycle?
750 if ((cUTC+cod) > cPrevUTC) {
751 cCycleNo++;
752 cPrevUTC = cUTC + 0.0001;
753 }
754
755 sprintf(cDateObs, "%-10.10s", names_.datobs);
756 cDateObs[10] = '\0';
757
758 // Compute buffer number.
759 iMBuff = cBuffer + iBeamSel;
760 if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
761
762 if (cCycleNo < iMBuff->cycleNo) {
763 // Note that if the first beam and IF are not both selected cEOS
764 // will be cleared by rpget() when the next beam/IF is read.
765 cEOS = 1;
766 }
767
768 // Begin flush cycle?
769 if (cEOS || (iMBuff->nIF && (cUTC+cod) > (iMBuff->utc+0.0001))) {
770 cFlushing = 1;
771 cFlushBin = 0;
772 cFlushIF = 0;
773 }
774
775#ifdef PKSIO_DEBUG
776 char rel = '=';
777 double dt = utcDiff(cUTC, cW);
778 if (dt < 0.0) {
779 rel = '<';
780 } else if (dt > 0.0) {
781 rel = '>';
782 }
783
784 fprintf(stderr, "\n In:%4d%4d%3d%3d %.3f %c %.3f (%+.3fs) - "
785 "%sflushing\n", cScanNo, cCycleNo, beamNo, cIFno, cUTC, rel, cW, dt,
786 cFlushing ? "" : "not ");
787 if (cEOS) {
788 fprintf(stderr, "Start of new scan, flushing previous scan.\n");
789 }
790#endif
791 }
792 }
793
794
795 if (cFlushing) {
796 // Find the oldest integration to flush, noting that the last
797 // integration cycle may be incomplete.
798 beamNo = 0;
799 int cycleNo = 0;
800 for (; cFlushBin < cNBin; cFlushBin++) {
801 for (iBeamSel = 0; iBeamSel < cNBeamSel; iBeamSel++) {
802 iMBuff = cBuffer + iBeamSel + cNBeamSel*cFlushBin;
803
804 // iMBuff->nIF is decremented (below) and if zero signals that all
805 // IFs in an integration have been flushed.
806 if (iMBuff->nIF) {
807 if (cycleNo == 0 || iMBuff->cycleNo < cycleNo) {
808 beamNo = iMBuff->beamNo;
809 cycleNo = iMBuff->cycleNo;
810 }
811 }
812 }
813
814 if (beamNo) {
815 // Found an integration to flush.
816 break;
817 }
818
819 // Start with the first IF in the next bin.
820 cFlushIF = 0;
821 }
822
823 if (beamNo) {
824 iBeamSel = cBeamSel[beamNo-1];
825 iMBuff = cBuffer + iBeamSel + cNBeamSel*cFlushBin;
826
827 // Find the IF to flush.
828 for (; cFlushIF < cSimulIF; cFlushIF++) {
829 if (iMBuff->IFno[cFlushIF]) break;
830 }
831
832 } else {
833 // Flush complete.
834 cFlushing = 0;
835 if (cEOF) {
836 return -1;
837 }
838
839 // The last record read must have been the first of a new cycle.
840 beamNo = int(cBaseline / 256.0);
841 iBeamSel = cBeamSel[beamNo-1];
842
843 // Compute buffer number.
844 iMBuff = cBuffer + iBeamSel;
845 if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
846 }
847 }
848
849
850 if (cInterp && cFlushing == 1) {
851 // Start of flush cycle, interpolate the beam position.
852 //
853 // The position is measured by the control system at a time returned by
854 // RPFITSIN as the 'w' visibility coordinate. The ra and dec, returned
855 // as the 'u' and 'v' visibility coordinates, must be interpolated to
856 // the integration time which RPFITSIN returns as 'cUTC', this usually
857 // being a second or two later. The interpolation method used here is
858 // based on the scan rate.
859 //
860 // "This" RA, Dec, and UTC refers to the position currently stored in
861 // the buffer marked for output (iMBuff). This position is interpolated
862 // to the midpoint of that integration using either
863 // a) the rate currently sitting in iMBuff, which was computed from
864 // the previous integration, otherwise
865 // b) from the position recorded in the "next" integration which is
866 // currently sitting in the RPFITS commons,
867 // so that the position timestamps straddle the midpoint of the
868 // integration and is thereby interpolated rather than extrapolated.
869 //
870 // At the end of a scan, or if the next position has not been updated
871 // or its timestamp does not advance sufficiently, the most recent
872 // determination of the scan rate will be used for extrapolation which
873 // is quantified by the "rate age" measured in seconds beyond the
874 // interval defined by the position timestamps.
875
876 // At this point, iMBuff contains cU, cV, cW, parAngle and focusRot
877 // stored from the previous call to rpget() for this beam (i.e. "this"),
878 // and also raRate, decRate and paRate computed from that integration
879 // and the previous one.
880 double thisRA = iMBuff->ra;
881 double thisDec = iMBuff->dec;
882 double thisUTC = cPosUTC[iBeamSel];
883 double thisPA = iMBuff->parAngle + iMBuff->focusRot;
884
885#ifdef PKSIO_DEBUG
886 fprintf(stderr, "This (%d) ra, dec, UTC: %9.4f %9.4f %10.3f %9.4f\n",
887 iMBuff->cycleNo, thisRA*R2D, thisDec*R2D, thisUTC, thisPA*R2D);
888#endif
889
890 if (cEOF || cEOS) {
891 // Use rates from the last cycle.
892 raRate = iMBuff->raRate;
893 decRate = iMBuff->decRate;
894 paRate = iMBuff->paRate;
895
896 } else {
897 if (cW == thisUTC) {
898 // The control system at Mopra typically does not update the
899 // positions between successive integration cycles at the end of a
900 // scan (nor are they flagged). In this case we use the previously
901 // computed rates, even if from the previous scan since these are
902 // likely to be a better guess than anything else.
903 raRate = iMBuff->raRate;
904 decRate = iMBuff->decRate;
905 paRate = iMBuff->paRate;
906
907 if (cU == thisRA && cV == thisDec) {
908 // Position and timestamp unchanged.
909 pCode = 1;
910
911 } else if (fabs(cU-thisRA) < 0.0001 && fabs(cV-thisDec) < 0.0001) {
912 // Allow small rounding errors (seen infrequently).
913 pCode = 1;
914
915 } else {
916 // (cU,cV) are probably rubbish (not yet seen in practice).
917 pCode = 2;
918 cU = thisRA;
919 cV = thisDec;
920 }
921
922#ifdef PKSIO_DEBUG
923 fprintf(stderr, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
924 "(0.000s)\n", cCycleNo, cU*R2D, cV*R2D, cW);
925#endif
926
927 } else {
928 double nextRA = cU;
929 double nextDec = cV;
930
931 // Check and, if necessary, repair the position timestamp,
932 // remembering that pCode refers to the NEXT cycle.
933 pCode = fixw(cDateObs, cCycleNo, beamNo, cAvRate, thisRA, thisDec,
934 thisUTC, nextRA, nextDec, cW);
935 if (pCode > 0) pCode += 3;
936 double nextUTC = cW;
937
938#ifdef PKSIO_DEBUG
939 fprintf(stderr, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
940 "(%+.3fs)\n", cCycleNo, nextRA*R2D, nextDec*R2D, nextUTC,
941 utcDiff(nextUTC, thisUTC));
942#endif
943
944 // Compute the scan rate for this beam.
945 double dUTC = utcDiff(nextUTC, thisUTC);
946 if ((0.0 < dUTC) && (dUTC < 600.0)) {
947 scanRate(cRA0, cDec0, thisRA, thisDec, nextRA, nextDec, dUTC,
948 raRate, decRate);
949
950 // Update the mean scan rate.
951 cAvRate[0] = (cAvRate[0]*cNRate + raRate) / (cNRate + 1);
952 cAvRate[1] = (cAvRate[1]*cNRate + decRate) / (cNRate + 1);
953 cNRate++;
954
955 // Rate of change of position angle.
956 if (sc_.sc_ant <= anten_.nant) {
957 paRate = 0.0;
958 } else {
959 int iOff = sc_.sc_q * (sc_.sc_ant - 1) - 1;
960 double nextPA = sc_.sc_cal[iOff + 4] + sc_.sc_cal[iOff + 7];
961 double paDiff = nextPA - thisPA;
962 if (paDiff > PI) {
963 paDiff -= TWOPI;
964 } else if (paDiff < -PI) {
965 paDiff += TWOPI;
966 }
967 paRate = paDiff / dUTC;
968 }
969
970 if (cInterp == 2) {
971 // Use the same interpolation scheme as the original pksmbfits
972 // client. This incorrectly assumed that (nextUTC - thisUTC) is
973 // equal to the integration time and interpolated by computing a
974 // weighted sum of the positions before and after the required
975 // time.
976
977 double utc = iMBuff->utc;
978 double tw1 = 1.0 - utcDiff(utc, thisUTC) / iMBuff->exposure;
979 double tw2 = 1.0 - utcDiff(nextUTC, utc) / iMBuff->exposure;
980 double gamma = (tw2 / (tw1 + tw2)) * dUTC / (utc - thisUTC);
981
982 // Guard against RA cycling through 24h in either direction.
983 if (fabs(nextRA - thisRA) > PI) {
984 if (nextRA < thisRA) {
985 nextRA += TWOPI;
986 } else {
987 nextRA -= TWOPI;
988 }
989 }
990
991 raRate = gamma * (nextRA - thisRA) / dUTC;
992 decRate = gamma * (nextDec - thisDec) / dUTC;
993 }
994
995 } else {
996 if (cCycleNo == 2 && fabs(utcDiff(cUTC,cW)) < 600.0) {
997 // thisUTC (i.e. cW for the first cycle) is rubbish, and
998 // probably the position as well (extremely rare in practice,
999 // e.g. 97-12-19_1029_235708-18_586e.hpf which actually has the
1000 // t/1000 scaling bug in the first cycle).
1001 iMBuff->pCode = 3;
1002 thisRA = cU;
1003 thisDec = cV;
1004 thisUTC = cW;
1005 raRate = 0.0;
1006 decRate = 0.0;
1007 paRate = 0.0;
1008
1009 } else {
1010 // cW is rubbish and probably (cU,cV), and possibly the
1011 // parallactic angle and everything else as well (rarely seen
1012 // in practice, e.g. 97-12-09_0743_235707-58_327c.hpf and
1013 // 97-09-01_0034_123717-42_242b.hpf, the latter with bad
1014 // parallactic angle).
1015 pCode = 3;
1016 cU = thisRA;
1017 cV = thisDec;
1018 cW = thisUTC;
1019 raRate = iMBuff->raRate;
1020 decRate = iMBuff->decRate;
1021 paRate = iMBuff->paRate;
1022 }
1023 }
1024 }
1025 }
1026
1027
1028 // Choose the closest rate determination.
1029 if (cCycleNo == 1) {
1030 // Scan containing a single integration.
1031 iMBuff->raRate = 0.0;
1032 iMBuff->decRate = 0.0;
1033 iMBuff->paRate = 0.0;
1034
1035 } else {
1036 double dUTC = iMBuff->utc - cPosUTC[iBeamSel];
1037
1038 if (dUTC >= 0.0) {
1039 // In HIPASS/ZOA, the position timestamp, which should always occur
1040 // on the whole second, normally precedes an integration midpoint
1041 // falling on the half-second. Consequently, positive ages are
1042 // always half-integral.
1043 dUTC = utcDiff(iMBuff->utc, cW);
1044 if (dUTC > 0.0) {
1045 iMBuff->rateAge = dUTC;
1046 } else {
1047 iMBuff->rateAge = 0.0f;
1048 }
1049
1050 iMBuff->raRate = raRate;
1051 iMBuff->decRate = decRate;
1052 iMBuff->paRate = paRate;
1053
1054 } else {
1055 // In HIPASS/ZOA, negative ages occur when the integration midpoint,
1056 // occurring on the whole second, precedes the position timestamp.
1057 // Thus negative ages are always an integral number of seconds.
1058 // They have only been seen to occur sporadically in the period
1059 // 1999/05/31 to 1999/11/01, e.g. 1999-07-26_1821_005410-74_007c.hpf
1060 //
1061 // In recent (2008/10/07) Mopra data, small negative ages (~10ms,
1062 // occasionally up to ~300ms) seem to be the norm, with both the
1063 // position timestamp and integration midpoint falling close to but
1064 // not on the integral second.
1065 if (cCycleNo == 2) {
1066 // We have to start with something!
1067 iMBuff->rateAge = dUTC;
1068
1069 } else {
1070 // Although we did not record the relevant position timestamp
1071 // explicitly, it can easily be deduced.
1072 double w = iMBuff->utc - utcDiff(cUTC, iMBuff->utc) -
1073 iMBuff->rateAge;
1074 dUTC = utcDiff(iMBuff->utc, w);
1075
1076 if (dUTC > 0.0) {
1077 iMBuff->rateAge = 0.0f;
1078 } else {
1079 iMBuff->rateAge = dUTC;
1080 }
1081 }
1082
1083 iMBuff->raRate = raRate;
1084 iMBuff->decRate = decRate;
1085 iMBuff->paRate = paRate;
1086 }
1087 }
1088
1089#ifdef PKSIO_DEBUG
1090 double avRate = sqrt(cAvRate[0]*cAvRate[0] + cAvRate[1]*cAvRate[1]);
1091 fprintf(stderr, "RA, Dec, Av & PA rates: %8.4f %8.4f %8.4f %8.4f "
1092 "pCode %d\n", raRate*R2D, decRate*R2D, avRate*R2D, paRate*R2D, pCode);
1093#endif
1094
1095
1096 // Compute the position of this beam for all bins.
1097 for (int idx = 0; idx < cNBin; idx++) {
1098 int jbuff = iBeamSel + cNBeamSel*idx;
1099
1100 cBuffer[jbuff].raRate = iMBuff->raRate;
1101 cBuffer[jbuff].decRate = iMBuff->decRate;
1102 cBuffer[jbuff].paRate = iMBuff->paRate;
1103
1104 double dUTC = utcDiff(cBuffer[jbuff].utc, thisUTC);
1105 if (dUTC > 100.0) {
1106 // Must have cycled through midnight.
1107 dUTC -= 86400.0;
1108 }
1109
1110 applyRate(cRA0, cDec0, thisRA, thisDec,
1111 cBuffer[jbuff].raRate, cBuffer[jbuff].decRate, dUTC,
1112 cBuffer[jbuff].ra, cBuffer[jbuff].dec);
1113
1114#ifdef PKSIO_DEBUG
1115 fprintf(stderr, "Intp (%d) ra, dec, UTC: %9.4f %9.4f %10.3f (pCode, "
1116 "age: %d %.1fs)\n", iMBuff->cycleNo, cBuffer[jbuff].ra*R2D,
1117 cBuffer[jbuff].dec*R2D, cBuffer[jbuff].utc, iMBuff->pCode,
1118 iMBuff->rateAge);
1119#endif
1120 }
1121
1122 cFlushing = 2;
1123 }
1124
1125
1126 if (cFlushing) {
1127 // Copy buffer location out one IF at a time.
1128 MBrec.extract(*iMBuff, cFlushIF);
1129 haveData = 1;
1130
1131#ifdef PKSIO_DEBUG
1132 fprintf(stderr, "Out:%4d%4d%3d%3d\n", MBrec.scanNo, MBrec.cycleNo,
1133 MBrec.beamNo, MBrec.IFno[0]);
1134#endif
1135
1136 // Signal that this IF in this buffer location has been flushed.
1137 iMBuff->IFno[cFlushIF] = 0;
1138
1139 iMBuff->nIF--;
1140 if (iMBuff->nIF == 0) {
1141 // All IFs in this buffer location have been flushed. Stop cEOS
1142 // being set when the next integration is read.
1143 iMBuff->cycleNo = 0;
1144
1145 } else {
1146 // Carry on flushing the other IFs.
1147 continue;
1148 }
1149
1150 // Has the whole buffer been flushed?
1151 if (cFlushBin == cNBin - 1) {
1152 if (cEOS || cEOF) {
1153 // Carry on flushing other buffers.
1154 cFlushIF = 0;
1155 continue;
1156 }
1157
1158 cFlushing = 0;
1159
1160 beamNo = int(cBaseline / 256.0);
1161 iBeamSel = cBeamSel[beamNo-1];
1162
1163 // Compute buffer number.
1164 iMBuff = cBuffer + iBeamSel;
1165 if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
1166 }
1167 }
1168
1169 if (!cFlushing) {
1170 // Buffer this MBrec.
1171 if ((cScanNo > iMBuff->scanNo) && iMBuff->IFno[0]) {
1172 // Sanity check on the number of IFs in the new scan.
1173 if (if_.n_if != cNIF) {
1174 sprintf(cMsg, "WARNING: Scan %d has %d IFs instead of %d, "
1175 "continuing.", cScanNo, if_.n_if, cNIF);
1176 logMsg(cMsg);
1177 }
1178 }
1179
1180 // Sanity check on incomplete integrations within a scan.
1181 if (iMBuff->nIF && (iMBuff->cycleNo != cCycleNo)) {
1182 // Force the incomplete integration to be flushed before proceeding.
1183 cFlushing = 1;
1184 continue;
1185 }
1186
1187#ifdef PKSIO_DEBUG
1188 fprintf(stderr, "Buf:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, cIFno);
1189#endif
1190
1191 // Store IF-independent parameters only for the first IF of a new cycle,
1192 // particularly because this is the only one for which the scan rates
1193 // are computed above.
1194 int firstIF = (iMBuff->nIF == 0);
1195 if (firstIF) {
1196 iMBuff->scanNo = cScanNo;
1197 iMBuff->cycleNo = cCycleNo;
1198
1199 // Times.
1200 strcpy(iMBuff->datobs, cDateObs);
1201 iMBuff->utc = cUTC;
1202 iMBuff->exposure = param_.intbase;
1203
1204 // Source identification.
1205 sprintf(iMBuff->srcName, "%-16.16s",
1206 names_.su_name + (cSrcNo-1)*16);
1207 iMBuff->srcName[16] = '\0';
1208 iMBuff->srcRA = doubles_.su_ra[cSrcNo-1];
1209 iMBuff->srcDec = doubles_.su_dec[cSrcNo-1];
1210
1211 // Rest frequency of the line of interest.
1212 iMBuff->restFreq = doubles_.rfreq;
1213 if (strncmp(names_.instrument, "ATPKSMB", 7) == 0) {
1214 // Fix the HI rest frequency recorded for Parkes multibeam data.
1215 double reffreq = doubles_.freq;
1216 double restfreq = doubles_.rfreq;
1217 if ((restfreq == 0.0 || fabs(restfreq - reffreq) == 0.0) &&
1218 fabs(reffreq - 1420.405752e6) < 100.0) {
1219 iMBuff->restFreq = 1420.405752e6;
1220 }
1221 }
1222
1223 // Observation type.
1224 int j;
1225 for (j = 0; j < 15; j++) {
1226 iMBuff->obsType[j] = names_.card[11+j];
1227 if (iMBuff->obsType[j] == '\'') break;
1228 }
1229 iMBuff->obsType[j] = '\0';
1230
1231 // Beam-dependent parameters.
1232 iMBuff->beamNo = beamNo;
1233
1234 // Beam position at the specified time.
1235 if (cSUpos) {
1236 // Non-ATNF data that does not store the position in (u,v,w).
1237 iMBuff->ra = doubles_.su_ra[cSrcNo-1];
1238 iMBuff->dec = doubles_.su_dec[cSrcNo-1];
1239 } else {
1240 iMBuff->ra = cU;
1241 iMBuff->dec = cV;
1242 }
1243 cPosUTC[iBeamSel] = cW;
1244 iMBuff->pCode = pCode;
1245
1246 // Store rates for next time.
1247 iMBuff->raRate = raRate;
1248 iMBuff->decRate = decRate;
1249 iMBuff->paRate = paRate;
1250 }
1251
1252 // IF-dependent parameters.
1253 int iIF = cIFno - 1;
1254 int startChan = cStartChan[iIF];
1255 int endChan = cEndChan[iIF];
1256 int refChan = cRefChan[iIF];
1257
1258 int nChan = abs(endChan - startChan) + 1;
1259
1260 iIFSel = cIFSel[iIF];
1261 if (iMBuff->IFno[iIFSel] == 0) {
1262 iMBuff->nIF++;
1263 iMBuff->IFno[iIFSel] = cIFno;
1264 } else {
1265 // Integration cycle written to the output file twice (the only known
1266 // example is 1999-05-22_1914_000-031805_03v.hpf).
1267 sprintf(cMsg, "WARNING: Integration cycle %d:%d, beam %2d, \n"
1268 " IF %d was duplicated.", cScanNo, cCycleNo-1,
1269 beamNo, cIFno);
1270 logMsg(cMsg);
1271 }
1272 iMBuff->nChan[iIFSel] = nChan;
1273 iMBuff->nPol[iIFSel] = cNPol[iIF];
1274
1275 iMBuff->fqRefPix[iIFSel] = doubles_.if_ref[iIF];
1276 iMBuff->fqRefVal[iIFSel] = doubles_.if_freq[iIF];
1277 iMBuff->fqDelt[iIFSel] =
1278 if_.if_invert[iIF] * fabs(doubles_.if_bw[iIF] /
1279 (if_.if_nfreq[iIF] - 1));
1280
1281 // Adjust for channel selection.
1282 if (iMBuff->fqRefPix[iIFSel] != refChan) {
1283 iMBuff->fqRefVal[iIFSel] +=
1284 (refChan - iMBuff->fqRefPix[iIFSel]) *
1285 iMBuff->fqDelt[iIFSel];
1286 iMBuff->fqRefPix[iIFSel] = refChan;
1287 }
1288
1289 if (endChan < startChan) {
1290 iMBuff->fqDelt[iIFSel] = -iMBuff->fqDelt[iIFSel];
1291 }
1292
1293
1294 // System temperature.
1295 int iBeam = beamNo - 1;
1296 int scq = sc_.sc_q;
1297 float TsysPol1 = sc_.sc_cal[scq*iBeam + 3];
1298 float TsysPol2 = sc_.sc_cal[scq*iBeam + 4];
1299 iMBuff->tsys[iIFSel][0] = TsysPol1*TsysPol1;
1300 iMBuff->tsys[iIFSel][1] = TsysPol2*TsysPol2;
1301
1302 // Calibration factor; may be changed later if the data is recalibrated.
1303 if (scq > 14) {
1304 // Will only be present for Parkes Multibeam or LBA data.
1305 iMBuff->calfctr[iIFSel][0] = sc_.sc_cal[scq*iBeam + 14];
1306 iMBuff->calfctr[iIFSel][1] = sc_.sc_cal[scq*iBeam + 15];
1307 } else {
1308 iMBuff->calfctr[iIFSel][0] = 0.0f;
1309 iMBuff->calfctr[iIFSel][1] = 0.0f;
1310 }
1311
1312 // Cross-polarization calibration factor (unknown to MBFITS).
1313 for (int j = 0; j < 2; j++) {
1314 iMBuff->xcalfctr[iIFSel][j] = 0.0f;
1315 }
1316
1317 // Baseline parameters (unknown to MBFITS).
1318 iMBuff->haveBase = 0;
1319
1320 // Data (always present in MBFITS).
1321 iMBuff->haveSpectra = 1;
1322
1323 // Flag: bit 0 set if off source.
1324 // bit 1 set if loss of sync in A polarization.
1325 // bit 2 set if loss of sync in B polarization.
1326 unsigned char rpflag =
1327 (unsigned char)(sc_.sc_cal[scq*iBeam + 12] + 0.5f);
1328
1329 // The baseline flag may be set independently.
1330 if (rpflag == 0) rpflag = cFlag;
1331
1332 // Copy and scale data.
1333 int inc = 2 * if_.if_nstok[iIF];
1334 if (endChan < startChan) inc = -inc;
1335
1336 float TsysF;
1337 iMBuff->spectra[iIFSel] = iMBuff->spectra[0] + cChanOff[iIF];
1338 iMBuff->flagged[iIFSel] = iMBuff->flagged[0] + cChanOff[iIF];
1339
1340 float *spectra = iMBuff->spectra[iIFSel];
1341 unsigned char *flagged = iMBuff->flagged[iIFSel];
1342 for (int ipol = 0; ipol < cNPol[iIF]; ipol++) {
1343 if (sc_.sc_cal[scq*iBeam + 3 + ipol] > 0.0f) {
1344 // The correlator has already applied the calibration.
1345 TsysF = 1.0f;
1346 } else {
1347 // The correlator has normalized cVis[k] to a Tsys of 500K.
1348 TsysF = iMBuff->tsys[iIFSel][ipol] / 500.0f;
1349 }
1350
1351 int k = 2 * (if_.if_nstok[iIF]*(startChan - 1) + ipol);
1352 for (int ichan = 0; ichan < nChan; ichan++) {
1353 *(spectra++) = TsysF * cVis[k];
1354 *(flagged++) = rpflag;
1355 k += inc;
1356 }
1357 }
1358
1359 if (cHaveXPol[iIF]) {
1360 int k = 2 * (3*(startChan - 1) + 2);
1361 iMBuff->xpol[iIFSel] = iMBuff->xpol[0] + cXpolOff[iIF];
1362 float *xpol = iMBuff->xpol[iIFSel];
1363 for (int ichan = 0; ichan < nChan; ichan++) {
1364 *(xpol++) = cVis[k];
1365 *(xpol++) = cVis[k+1];
1366 k += inc;
1367 }
1368 }
1369
1370
1371 // Calibration factor applied to the data by the correlator.
1372 if (scq > 14) {
1373 // Will only be present for Parkes Multibeam or LBA data.
1374 iMBuff->tcal[iIFSel][0] = sc_.sc_cal[scq*iBeam + 14];
1375 iMBuff->tcal[iIFSel][1] = sc_.sc_cal[scq*iBeam + 15];
1376 } else {
1377 iMBuff->tcal[iIFSel][0] = 0.0f;
1378 iMBuff->tcal[iIFSel][1] = 0.0f;
1379 }
1380
1381 if (firstIF) {
1382 if (sc_.sc_ant <= anten_.nant) {
1383 // No extra syscal information present.
1384 iMBuff->extraSysCal = 0;
1385 iMBuff->azimuth = 0.0f;
1386 iMBuff->elevation = 0.0f;
1387 iMBuff->parAngle = 0.0f;
1388 iMBuff->focusAxi = 0.0f;
1389 iMBuff->focusTan = 0.0f;
1390 iMBuff->focusRot = 0.0f;
1391 iMBuff->temp = 0.0f;
1392 iMBuff->pressure = 0.0f;
1393 iMBuff->humidity = 0.0f;
1394 iMBuff->windSpeed = 0.0f;
1395 iMBuff->windAz = 0.0f;
1396 strcpy(iMBuff->tcalTime, " ");
1397 iMBuff->refBeam = 0;
1398
1399 } else {
1400 // Additional information for Parkes Multibeam data.
1401 int iOff = scq*(sc_.sc_ant - 1) - 1;
1402 iMBuff->extraSysCal = 1;
1403
1404 iMBuff->azimuth = sc_.sc_cal[iOff + 2];
1405 iMBuff->elevation = sc_.sc_cal[iOff + 3];
1406 iMBuff->parAngle = sc_.sc_cal[iOff + 4];
1407
1408 iMBuff->focusAxi = sc_.sc_cal[iOff + 5] * 1e-3;
1409 iMBuff->focusTan = sc_.sc_cal[iOff + 6] * 1e-3;
1410 iMBuff->focusRot = sc_.sc_cal[iOff + 7];
1411
1412 iMBuff->temp = sc_.sc_cal[iOff + 8];
1413 iMBuff->pressure = sc_.sc_cal[iOff + 9];
1414 iMBuff->humidity = sc_.sc_cal[iOff + 10];
1415 iMBuff->windSpeed = sc_.sc_cal[iOff + 11];
1416 iMBuff->windAz = sc_.sc_cal[iOff + 12];
1417
1418 char *tcalTime = iMBuff->tcalTime;
1419 sprintf(tcalTime, "%-16.16s", (char *)(&sc_.sc_cal[iOff+13]));
1420 tcalTime[16] = '\0';
1421
1422#ifndef AIPS_LITTLE_ENDIAN
1423 // Do byte swapping on the ASCII date string.
1424 for (int j = 0; j < 16; j += 4) {
1425 char ctmp;
1426 ctmp = tcalTime[j];
1427 tcalTime[j] = tcalTime[j+3];
1428 tcalTime[j+3] = ctmp;
1429 ctmp = tcalTime[j+1];
1430 tcalTime[j+1] = tcalTime[j+2];
1431 tcalTime[j+2] = ctmp;
1432 }
1433#endif
1434
1435 // Reference beam number.
1436 float refbeam = sc_.sc_cal[iOff + 17];
1437 if (refbeam > 0.0f || refbeam < 100.0f) {
1438 iMBuff->refBeam = int(refbeam);
1439 } else {
1440 iMBuff->refBeam = 0;
1441 }
1442 }
1443 }
1444 }
1445 }
1446
1447 return 0;
1448}
1449
1450//-------------------------------------------------------- MBFITSreader::rpget
1451
1452// Read the next data record from the RPFITS file.
1453
1454int MBFITSreader::rpget(int syscalonly, int &EOS)
1455{
1456 EOS = 0;
1457
1458 int retries = 0;
1459
1460 // Allow 10 read errors.
1461 int numErr = 0;
1462
1463 int jstat = 0;
1464 while (numErr < 10) {
1465 int lastjstat = jstat;
1466
1467 switch(rpfitsin(jstat)) {
1468 case -1:
1469 // Read failed; retry.
1470 numErr++;
1471 logMsg("WARNING: RPFITS read failed - retrying.");
1472 jstat = 0;
1473 break;
1474
1475 case 0:
1476 // Successful read.
1477 if (lastjstat == 0) {
1478 if (cBaseline == -1) {
1479 // Syscal data.
1480 if (syscalonly) {
1481 return 0;
1482 }
1483
1484 } else {
1485 if (!syscalonly) {
1486 return 0;
1487 }
1488 }
1489 }
1490
1491 // Last operation was to read header or FG table; now read data.
1492 break;
1493
1494 case 1:
1495 // Encountered header while trying to read data; read it.
1496 EOS = 1;
1497 jstat = -1;
1498 break;
1499
1500 case 2:
1501 // End of scan; read past it.
1502 jstat = 0;
1503 break;
1504
1505 case 3:
1506 // End-of-file; retry applies to real-time mode.
1507 if (retries++ >= cRetry) {
1508 return -1;
1509 }
1510
1511 sleep(10);
1512 jstat = 0;
1513 break;
1514
1515 case 4:
1516 // Encountered FG table while trying to read data; read it.
1517 jstat = -1;
1518 break;
1519
1520 case 5:
1521 // Illegal data at end of block after close/reopen operation; retry.
1522 jstat = 0;
1523 break;
1524
1525 default:
1526 // Shouldn't reach here.
1527 sprintf(cMsg, "WARNING: Unrecognized RPFITSIN return code: %d "
1528 "(retrying).", jstat);
1529 logMsg(cMsg);
1530 jstat = 0;
1531 break;
1532 }
1533 }
1534
1535 logMsg("ERROR: RPFITS read failed too many times.");
1536 return 2;
1537}
1538
1539//----------------------------------------------------- MBFITSreader::rpfitsin
1540
1541// Wrapper around RPFITSIN that reports errors. Returned RPFITSIN subroutine
1542// arguments are captured as MBFITSreader member variables.
1543
1544int MBFITSreader::rpfitsin(int &jstat)
1545
1546{
1547 rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
1548 &cBin, &cIFno, &cSrcNo);
1549
1550 // Handle messages from RPFITSIN.
1551 if (names_.errmsg[0] != ' ') {
1552 int i;
1553 for (i = 80; i > 0; i--) {
1554 if (names_.errmsg[i-1] != ' ') break;
1555 }
1556
1557 sprintf(cMsg, "WARNING: Cycle %d:%03d, RPFITSIN reported -\n"
1558 " %.*s", cScanNo, cCycleNo, i, names_.errmsg);
1559 logMsg(cMsg);
1560 }
1561
1562 return jstat;
1563}
1564
1565//------------------------------------------------------- MBFITSreader::fixPos
1566
1567// Check and, if necessary, repair a position timestamp.
1568//
1569// Problems with the position timestamp manifest themselves via the scan rate:
1570//
1571// 1) Zero scan rate pairs, 1997/02/28 to 1998/01/07
1572//
1573// These occur because the position timestamp for the first integration
1574// of the pair is erroneous; the value recorded is t/1000, where t is the
1575// true value.
1576// Earliest known: 97-02-28_1725_132653-42_258a.hpf
1577// Latest known: 98-01-02_1923_095644-50_165c.hpf
1578// (time range chosen to encompass observing runs).
1579//
1580// 2) Slow-fast scan rate pairs (0.013 - 0.020 deg/s),
1581// 1997/03/28 to 1998/01/07.
1582//
1583// The UTC position timestamp is 1.0s later than it should be (never
1584// earlier), almost certainly arising from an error in the telescope
1585// control system.
1586// Earliest known: 97-03-28_0150_010420-74_008d.hpf
1587// Latest known: 98-01-04_1502_065150-02_177c.hpf
1588// (time range chosen to encompass observing runs).
1589//
1590// 3) Slow-fast scan rate pairs (0.015 - 0.018 deg/s),
1591// 1999/05/20 to 2001/07/12 (HIPASS and ZOA),
1592// 2001/09/02 to 2001/12/04 (HIPASS and ZOA),
1593// 2002/03/28 to 2002/05/13 (ZOA only),
1594// 2003/04/26 to 2003/06/09 (ZOA only).
1595// Earliest known: 1999-05-20_1818_175720-50_297e.hpf
1596// Latest known: 2001-12-04_1814_065531p14_173e.hpf (HIPASS)
1597// 2003-06-09_1924_352-085940_-6c.hpf (ZOA)
1598//
1599// Caused by the Linux signalling NaN problem. IEEE "signalling" NaNs
1600// are silently transformed to "quiet" NaNs during assignment by setting
1601// bit 22. This affected RPFITS because of its use of VAX-format
1602// floating-point numbers which, with their permuted bytes, may sometimes
1603// appear as signalling NaNs.
1604//
1605// The problem arose when the linux correlator came online and was
1606// fixed with a workaround to the RPFITS library (repeated episodes
1607// are probably due to use of an older version of the library). It
1608// should not have affected the data significantly because of the
1609// low relative error, which ranges from 0.0000038 to 0.0000076, but
1610// it is important for the computation of scan rates which requires
1611// taking the difference of two large UTC timestamps, one or other
1612// of which will have 0.5s added to it.
1613//
1614// The return value identifies which, if any, of these problems was repaired.
1615
1616int MBFITSreader::fixw(
1617 const char *datobs,
1618 int cycleNo,
1619 int beamNo,
1620 double avRate[2],
1621 double thisRA,
1622 double thisDec,
1623 double thisUTC,
1624 double nextRA,
1625 double nextDec,
1626 float &nextUTC)
1627{
1628 if (strcmp(datobs, "2003-06-09") > 0) {
1629 return 0;
1630
1631 } else if (strcmp(datobs, "1998-01-07") <= 0) {
1632 if (nextUTC < thisUTC && (nextUTC + 86400.0) > (thisUTC + 600.0)) {
1633 // Possible scaling problem.
1634 double diff = nextUTC*1000.0 - thisUTC;
1635 if (0.0 < diff && diff < 600.0) {
1636 nextUTC *= 1000.0;
1637 return 1;
1638 } else {
1639 // Irreparable.
1640 return -1;
1641 }
1642 }
1643
1644 if (cycleNo > 2) {
1645 if (beamNo == 1) {
1646 // This test is only reliable for beam 1.
1647 double dUTC = nextUTC - thisUTC;
1648 if (dUTC < 0.0) dUTC += 86400.0;
1649
1650 // Guard against RA cycling through 24h in either direction.
1651 if (fabs(nextRA - thisRA) > PI) {
1652 if (nextRA < thisRA) {
1653 nextRA += TWOPI;
1654 } else {
1655 nextRA -= TWOPI;
1656 }
1657 }
1658
1659 double dRA = (nextRA - thisRA) * cos(nextDec);
1660 double dDec = nextDec - thisDec;
1661 double arc = sqrt(dRA*dRA + dDec*dDec);
1662
1663 double averate = sqrt(avRate[0]*avRate[0] + avRate[1]*avRate[1]);
1664 double diff1 = fabs(averate - arc/(dUTC-1.0));
1665 double diff2 = fabs(averate - arc/dUTC);
1666 if ((diff1 < diff2) && (diff1 < 0.05*averate)) {
1667 nextUTC -= 1.0;
1668 cCode5 = cycleNo;
1669 return 2;
1670 } else {
1671 cCode5 = 0;
1672 }
1673
1674 } else {
1675 if (cycleNo == cCode5) {
1676 nextUTC -= 1.0;
1677 return 2;
1678 }
1679 }
1680 }
1681
1682 } else if ((strcmp(datobs, "1999-05-20") >= 0 &&
1683 strcmp(datobs, "2001-07-12") <= 0) ||
1684 (strcmp(datobs, "2001-09-02") >= 0 &&
1685 strcmp(datobs, "2001-12-04") <= 0) ||
1686 (strcmp(datobs, "2002-03-28") >= 0 &&
1687 strcmp(datobs, "2002-05-13") <= 0) ||
1688 (strcmp(datobs, "2003-04-26") >= 0 &&
1689 strcmp(datobs, "2003-06-09") <= 0)) {
1690 // Signalling NaN problem, e.g. 1999-07-26_1839_011106-74_009c.hpf.
1691 // Position timestamps should always be an integral number of seconds.
1692 double resid = nextUTC - int(nextUTC);
1693 if (resid == 0.5) {
1694 nextUTC -= 0.5;
1695 return 3;
1696 }
1697 }
1698
1699 return 0;
1700}
1701
1702//-------------------------------------------------------- MBFITSreader::close
1703
1704// Close the input file.
1705
1706void MBFITSreader::close(void)
1707{
1708 if (cMBopen) {
1709 int jstat = 1;
1710 rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
1711 &cBin, &cIFno, &cSrcNo);
1712
1713 if (cBeams) delete [] cBeams;
1714 if (cIFs) delete [] cIFs;
1715 if (cNChan) delete [] cNChan;
1716 if (cNPol) delete [] cNPol;
1717 if (cHaveXPol) delete [] cHaveXPol;
1718 if (cStartChan) delete [] cStartChan;
1719 if (cEndChan) delete [] cEndChan;
1720 if (cRefChan) delete [] cRefChan;
1721
1722 if (cVis) delete [] cVis;
1723 if (cWgt) delete [] cWgt;
1724
1725 if (cBeamSel) delete [] cBeamSel;
1726 if (cIFSel) delete [] cIFSel;
1727 if (cChanOff) delete [] cChanOff;
1728 if (cXpolOff) delete [] cXpolOff;
1729 if (cBuffer) delete [] cBuffer;
1730 if (cPosUTC) delete [] cPosUTC;
1731
1732 cMBopen = 0;
1733 }
1734}
1735
1736//-------------------------------------------------------------------- utcDiff
1737
1738// Subtract two UTCs (s) allowing for any plausible number of cycles through
1739// 86400s, returning a result in the range [-43200, +43200]s.
1740
1741double MBFITSreader::utcDiff(double utc1, double utc2)
1742{
1743 double diff = utc1 - utc2;
1744
1745 if (diff > 43200.0) {
1746 diff -= 86400.0;
1747 while (diff > 43200.0) diff -= 86400.0;
1748 } else if (diff < -43200.0) {
1749 diff += 86400.0;
1750 while (diff < -43200.0) diff += 86400.0;
1751 }
1752
1753 return diff;
1754}
1755
1756//------------------------------------------------------- scanRate & applyRate
1757
1758// Compute and apply the scan rate corrected for grid convergence. (ra0,dec0)
1759// are the coordinates of the central beam, assumed to be the tracking centre.
1760// The rate computed in RA will be a rate of change of angular distance in the
1761// direction of increasing RA at the position of the central beam. Similarly
1762// for declination. Angles in radian, time in s.
1763
1764void MBFITSreader::scanRate(
1765 double ra0,
1766 double dec0,
1767 double ra1,
1768 double dec1,
1769 double ra2,
1770 double dec2,
1771 double dt,
1772 double &raRate,
1773 double &decRate)
1774{
1775 // Transform to a system where the central beam lies on the equator at 12h.
1776 eulerx(ra1, dec1, ra0+HALFPI, -dec0, -HALFPI, ra1, dec1);
1777 eulerx(ra2, dec2, ra0+HALFPI, -dec0, -HALFPI, ra2, dec2);
1778
1779 raRate = (ra2 - ra1) / dt;
1780 decRate = (dec2 - dec1) / dt;
1781}
1782
1783
1784void MBFITSreader::applyRate(
1785 double ra0,
1786 double dec0,
1787 double ra1,
1788 double dec1,
1789 double raRate,
1790 double decRate,
1791 double dt,
1792 double &ra2,
1793 double &dec2)
1794{
1795 // Transform to a system where the central beam lies on the equator at 12h.
1796 eulerx(ra1, dec1, ra0+HALFPI, -dec0, -HALFPI, ra1, dec1);
1797
1798 ra2 = ra1 + (raRate * dt);
1799 dec2 = dec1 + (decRate * dt);
1800
1801 // Transform back.
1802 eulerx(ra2, dec2, -HALFPI, dec0, ra0+HALFPI, ra2, dec2);
1803}
1804
1805//--------------------------------------------------------------------- eulerx
1806
1807void MBFITSreader::eulerx(
1808 double lng0,
1809 double lat0,
1810 double phi0,
1811 double theta,
1812 double phi,
1813 double &lng1,
1814 double &lat1)
1815
1816// Applies the Euler angle based transformation of spherical coordinates.
1817//
1818// phi0 Longitude of the ascending node in the old system, radians. The
1819// ascending node is the point of intersection of the equators of
1820// the two systems such that the equator of the new system crosses
1821// from south to north as viewed in the old system.
1822//
1823// theta Angle between the poles of the two systems, radians. THETA is
1824// positive for a positive rotation about the ascending node.
1825//
1826// phi Longitude of the ascending node in the new system, radians.
1827
1828{
1829 // Compute intermediaries.
1830 double lng0p = lng0 - phi0;
1831 double slng0p = sin(lng0p);
1832 double clng0p = cos(lng0p);
1833 double slat0 = sin(lat0);
1834 double clat0 = cos(lat0);
1835 double ctheta = cos(theta);
1836 double stheta = sin(theta);
1837
1838 double x = clat0*clng0p;
1839 double y = clat0*slng0p*ctheta + slat0*stheta;
1840
1841 // Longitude in the new system.
1842 if (x != 0.0 || y != 0.0) {
1843 lng1 = phi + atan2(y, x);
1844 } else {
1845 // Longitude at the poles in the new system is consistent with that
1846 // specified in the old system.
1847 lng1 = phi + lng0p;
1848 }
1849 lng1 = fmod(lng1, TWOPI);
1850 if (lng1 < 0.0) lng1 += TWOPI;
1851
1852 lat1 = asin(slat0*ctheta - clat0*stheta*slng0p);
1853}
Note: See TracBrowser for help on using the repository browser.