source: trunk/external/atnf/PKSIO/MBFITSreader.cc@ 1705

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

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