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

Last change on this file since 1583 was 1452, checked in by Malte Marquarding, 16 years ago

update from livedata CVS
File size: 54.9 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.54 2008-11-17 06:51:55 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 simulIF = 0;
591 int maxChan = 0;
592 int maxXpol = 0;
593
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 simulIF = max(simulIF, 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.");
647 close();
648 return 1;
649 }
650
651 // Allocate buffer data storage; the MBrecord constructor zeroes
652 // class members such as cycleNo that are tested in the first pass
653 // below.
654 int nBuff = cNBeamSel * cNBin;
655 cBuffer = new MBrecord[nBuff];
656
657 // Allocate memory for spectral arrays.
658 for (int ibuff = 0; ibuff < nBuff; ibuff++) {
659 cBuffer[ibuff].setNIFs(simulIF);
660 cBuffer[ibuff].allocate(0, maxChan, maxXpol);
661
662 // Signal that this IF in this buffer has been flushed.
663 for (int iIF = 0; iIF < simulIF; iIF++) {
664 cBuffer[ibuff].IFno[iIF] = 0;
665 }
666 }
667
668 cPosUTC = new double[cNBeamSel];
669
670 cFirst = 0;
671 cScanNo = 1;
672 cCycleNo = 0;
673 cPrevUTC = -1.0;
674 }
675
676 // Check for end-of-scan.
677 if (cEOS) {
678 cScanNo++;
679 cCycleNo = 0;
680 cPrevUTC = -1.0;
681 }
682
683 // Check for change-of-day.
684 double cod = 0.0;
685 if ((cUTC + 86400.0) < (cPrevUTC + 600.0)) {
686 // cUTC should continue to increase past 86400 during a single scan.
687 // However, if the RPFITS file contains multiple scans that straddle
688 // midnight then cUTC can jump backwards from the end of one scan to
689 // the start of the next.
690#ifdef PKSIO_DEBUG
691 fprintf(stderr, "Change-of-day on cUTC: %.1f -> %.1f",
692 cPrevUTC, cUTC);
693#endif
694 // Can't change the recorded value of cUTC directly (without also
695 // changing dateobs) so change-of-day must be recorded separately as
696 // an offset to be applied when comparing integration timestamps.
697 cod = 86400.0;
698
699 } else if (cUTC < cPrevUTC - 1.0) {
700 sprintf(cMsg,
701 "WARNING: Cycle %d:%03d-%03d, UTC went backwards from\n"
702 " %.1f to %.1f! Incrementing day number,\n"
703 " positions may be unreliable.", cScanNo, cCycleNo,
704 cCycleNo+1, cPrevUTC, cUTC);
705 logMsg(cMsg);
706 cUTC += 86400.0;
707 }
708
709 if (cNBin > 1) {
710 // Binning mode: correct the time.
711 cUTC += param_.intbase * (cBin - (cNBin + 1)/2.0);
712 }
713
714 // New integration cycle?
715 if ((cUTC+cod) > cPrevUTC) {
716 cCycleNo++;
717 cPrevUTC = cUTC + 0.0001;
718 }
719
720 // Apply beam selection.
721 beamNo = int(cBaseline / 256.0);
722 if (beamNo == 1) {
723 // Store the position of beam 1 for grid convergence corrections.
724 cRA0 = cU;
725 cDec0 = cV;
726 }
727 iBeamSel = cBeamSel[beamNo-1];
728 if (iBeamSel < 0) continue;
729
730 // Sanity check (mainly for MOPS).
731 if (cIFno > cNIF) continue;
732
733 // Apply IF selection.
734 iIFSel = cIFSel[cIFno - 1];
735 if (iIFSel < 0) continue;
736
737 sprintf(cDateObs, "%-10.10s", names_.datobs);
738 cDateObs[10] = '\0';
739
740 // Compute buffer number.
741 iMBuff = cBuffer + iBeamSel;
742 if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
743
744 if (cCycleNo < iMBuff->cycleNo) {
745 // Note that if the first beam and IF are not both selected cEOS
746 // will be cleared by rpget() when the next beam/IF is read.
747 cEOS = 1;
748 }
749
750 // Begin flush cycle?
751 if (cEOS || (iMBuff->nIF && (cUTC+cod) > (iMBuff->utc+0.0001))) {
752 cFlushing = 1;
753 cFlushBin = 0;
754 cFlushIF = 0;
755 }
756
757#ifdef PKSIO_DEBUG
758 char rel = '=';
759 double dt = utcDiff(cUTC, cW);
760 if (dt < 0.0) {
761 rel = '<';
762 } else if (dt > 0.0) {
763 rel = '>';
764 }
765
766 fprintf(stderr, "\n In:%4d%4d%3d%3d %.3f %c %.3f (%+.3fs) - "
767 "%sflushing\n", cScanNo, cCycleNo, beamNo, cIFno, cUTC, rel, cW, dt,
768 cFlushing ? "" : "not ");
769 if (cEOS) {
770 fprintf(stderr, "Start of new scan, flushing previous scan.\n");
771 }
772#endif
773 }
774 }
775
776
777 if (cFlushing) {
778 // Find the oldest integration to flush, noting that the last
779 // integration cycle may be incomplete.
780 beamNo = 0;
781 int cycleNo = 0;
782 for (; cFlushBin < cNBin; cFlushBin++) {
783 for (iBeamSel = 0; iBeamSel < cNBeamSel; iBeamSel++) {
784 iMBuff = cBuffer + iBeamSel + cNBeamSel*cFlushBin;
785
786 // iMBuff->nIF is set to zero (below) to signal that all IFs in
787 // an integration have been flushed.
788 if (iMBuff->nIF) {
789 if (cycleNo == 0 || iMBuff->cycleNo < cycleNo) {
790 beamNo = iMBuff->beamNo;
791 cycleNo = iMBuff->cycleNo;
792 }
793 }
794 }
795
796 if (beamNo) {
797 // Found an integration to flush.
798 break;
799 }
800 }
801
802 if (beamNo) {
803 iBeamSel = cBeamSel[beamNo-1];
804 iMBuff = cBuffer + iBeamSel + cNBeamSel*cFlushBin;
805
806 // Find the IF to flush.
807 for (; cFlushIF < iMBuff->nIF; cFlushIF++) {
808 if (iMBuff->IFno[cFlushIF]) break;
809 }
810
811 } else {
812 // Flush complete.
813 cFlushing = 0;
814 if (cEOF) {
815 return -1;
816 }
817
818 // The last record read must have been the first of a new cycle.
819 beamNo = int(cBaseline / 256.0);
820 iBeamSel = cBeamSel[beamNo-1];
821
822 // Compute buffer number.
823 iMBuff = cBuffer + iBeamSel;
824 if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
825 }
826 }
827
828
829 if (cFlushing && cFlushBin == 0 && cFlushIF == 0 && cInterp) {
830 // Start of flush cycle, interpolate the beam position.
831 //
832 // The position is measured by the control system at a time returned by
833 // RPFITSIN as the 'w' visibility coordinate. The ra and dec, returned
834 // as the 'u' and 'v' visibility coordinates, must be interpolated to
835 // the integration time which RPFITSIN returns as 'cUTC', this usually
836 // being a second or two later. The interpolation method used here is
837 // based on the scan rate.
838 //
839 // "This" RA, Dec, and UTC refers to the position currently stored in
840 // the buffer marked for output (iMBuff). This position is interpolated
841 // to the midpoint of that integration using either
842 // a) the rate currently sitting in iMBuff, which was computed from
843 // the previous integration, otherwise
844 // b) from the position recorded in the "next" integration which is
845 // currently sitting in the RPFITS commons,
846 // so that the position timestamps straddle the midpoint of the
847 // integration and is thereby interpolated rather than extrapolated.
848 //
849 // At the end of a scan, or if the next position has not been updated
850 // or its timestamp does not advance sufficiently, the most recent
851 // determination of the scan rate will be used for extrapolation which
852 // is quantified by the "rate age" measured in seconds beyond the
853 // interval defined by the position timestamps.
854
855 // At this point, iMBuff contains cU, cV, cW, parAngle and focusRot
856 // stored from the previous call to rpget() for this beam (i.e. "this"),
857 // and also raRate, decRate and paRate computed from that integration
858 // and the previous one.
859 double thisRA = iMBuff->ra;
860 double thisDec = iMBuff->dec;
861 double thisUTC = cPosUTC[iBeamSel];
862 double thisPA = iMBuff->parAngle + iMBuff->focusRot;
863
864#ifdef PKSIO_DEBUG
865 fprintf(stderr, "This (%d) ra, dec, UTC: %9.4f %9.4f %10.3f %9.4f\n",
866 iMBuff->cycleNo, thisRA*R2D, thisDec*R2D, thisUTC, thisPA*R2D);
867#endif
868
869 if (cEOF || cEOS) {
870 // Use rates from the last cycle.
871 raRate = iMBuff->raRate;
872 decRate = iMBuff->decRate;
873 paRate = iMBuff->paRate;
874
875 } else {
876 if (cW == thisUTC) {
877 // The control system at Mopra typically does not update the
878 // positions between successive integration cycles at the end of a
879 // scan (nor are they flagged). In this case we use the previously
880 // computed rates, even if from the previous scan since these are
881 // likely to be a better guess than anything else.
882 raRate = iMBuff->raRate;
883 decRate = iMBuff->decRate;
884 paRate = iMBuff->paRate;
885
886 if (cU == thisRA && cV == thisDec) {
887 // Position and timestamp unchanged.
888 pCode = 1;
889
890 } else if (fabs(cU-thisRA) < 0.0001 && fabs(cV-thisDec) < 0.0001) {
891 // Allow small rounding errors (seen infrequently).
892 pCode = 1;
893
894 } else {
895 // (cU,cV) are probably rubbish (not yet seen in practice).
896 pCode = 2;
897 cU = thisRA;
898 cV = thisDec;
899 }
900
901#ifdef PKSIO_DEBUG
902 fprintf(stderr, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
903 "(0.000s)\n", cCycleNo, cU*R2D, cV*R2D, cW);
904#endif
905
906 } else {
907 double nextRA = cU;
908 double nextDec = cV;
909
910 // Check and, if necessary, repair the position timestamp,
911 // remembering that pCode refers to the NEXT cycle.
912 pCode = fixw(cDateObs, cCycleNo, beamNo, cAvRate, thisRA, thisDec,
913 thisUTC, nextRA, nextDec, cW);
914 if (pCode > 0) pCode += 3;
915 double nextUTC = cW;
916
917#ifdef PKSIO_DEBUG
918 fprintf(stderr, "Next (%d) ra, dec, UTC: %9.4f %9.4f %10.3f "
919 "(%+.3fs)\n", cCycleNo, nextRA*R2D, nextDec*R2D, nextUTC,
920 utcDiff(nextUTC, thisUTC));
921#endif
922
923 // Compute the scan rate for this beam.
924 double dUTC = utcDiff(nextUTC, thisUTC);
925 if ((0.0 < dUTC) && (dUTC < 600.0)) {
926 scanRate(cRA0, cDec0, thisRA, thisDec, nextRA, nextDec, dUTC,
927 raRate, decRate);
928
929 // Update the mean scan rate.
930 cAvRate[0] = (cAvRate[0]*cNRate + raRate) / (cNRate + 1);
931 cAvRate[1] = (cAvRate[1]*cNRate + decRate) / (cNRate + 1);
932 cNRate++;
933
934 // Rate of change of position angle.
935 if (sc_.sc_ant <= anten_.nant) {
936 paRate = 0.0;
937 } else {
938 int iOff = sc_.sc_q * (sc_.sc_ant - 1) - 1;
939 double nextPA = sc_.sc_cal[iOff + 4] + sc_.sc_cal[iOff + 7];
940 double paDiff = nextPA - thisPA;
941 if (paDiff > PI) {
942 paDiff -= TWOPI;
943 } else if (paDiff < -PI) {
944 paDiff += TWOPI;
945 }
946 paRate = paDiff / dUTC;
947 }
948
949 if (cInterp == 2) {
950 // Use the same interpolation scheme as the original pksmbfits
951 // client. This incorrectly assumed that (nextUTC - thisUTC) is
952 // equal to the integration time and interpolated by computing a
953 // weighted sum of the positions before and after the required
954 // time.
955
956 double utc = iMBuff->utc;
957 double tw1 = 1.0 - utcDiff(utc, thisUTC) / iMBuff->exposure;
958 double tw2 = 1.0 - utcDiff(nextUTC, utc) / iMBuff->exposure;
959 double gamma = (tw2 / (tw1 + tw2)) * dUTC / (utc - thisUTC);
960
961 // Guard against RA cycling through 24h in either direction.
962 if (fabs(nextRA - thisRA) > PI) {
963 if (nextRA < thisRA) {
964 nextRA += TWOPI;
965 } else {
966 nextRA -= TWOPI;
967 }
968 }
969
970 raRate = gamma * (nextRA - thisRA) / dUTC;
971 decRate = gamma * (nextDec - thisDec) / dUTC;
972 }
973
974 } else {
975 if (cCycleNo == 2 && fabs(utcDiff(cUTC,cW)) < 600.0) {
976 // thisUTC (i.e. cW for the first cycle) is rubbish, and
977 // probably the position as well (extremely rare in practice,
978 // e.g. 97-12-19_1029_235708-18_586e.hpf which actually has the
979 // t/1000 scaling bug in the first cycle).
980 iMBuff->pCode = 3;
981 thisRA = cU;
982 thisDec = cV;
983 thisUTC = cW;
984 raRate = 0.0;
985 decRate = 0.0;
986 paRate = 0.0;
987
988 } else {
989 // cW is rubbish and probably (cU,cV), and possibly the
990 // parallactic angle and everything else as well (rarely seen
991 // in practice, e.g. 97-12-09_0743_235707-58_327c.hpf and
992 // 97-09-01_0034_123717-42_242b.hpf, the latter with bad
993 // parallactic angle).
994 pCode = 3;
995 cU = thisRA;
996 cV = thisDec;
997 cW = thisUTC;
998 raRate = iMBuff->raRate;
999 decRate = iMBuff->decRate;
1000 paRate = iMBuff->paRate;
1001 }
1002 }
1003 }
1004 }
1005
1006
1007 // Choose the closest rate determination.
1008 if (cCycleNo == 1) {
1009 // Scan containing a single integration.
1010 iMBuff->raRate = 0.0;
1011 iMBuff->decRate = 0.0;
1012 iMBuff->paRate = 0.0;
1013
1014 } else {
1015 double dUTC = iMBuff->utc - cPosUTC[iBeamSel];
1016
1017 if (dUTC >= 0.0) {
1018 // In HIPASS/ZOA, the position timestamp, which should always occur
1019 // on the whole second, normally precedes an integration midpoint
1020 // falling on the half-second. Consequently, positive ages are
1021 // always half-integral.
1022 dUTC = utcDiff(iMBuff->utc, cW);
1023 if (dUTC > 0.0) {
1024 iMBuff->rateAge = dUTC;
1025 } else {
1026 iMBuff->rateAge = 0.0f;
1027 }
1028
1029 iMBuff->raRate = raRate;
1030 iMBuff->decRate = decRate;
1031 iMBuff->paRate = paRate;
1032
1033 } else {
1034 // In HIPASS/ZOA, negative ages occur when the integration midpoint,
1035 // occurring on the whole second, precedes the position timestamp.
1036 // Thus negative ages are always an integral number of seconds.
1037 // They have only been seen to occur sporadically in the period
1038 // 1999/05/31 to 1999/11/01, e.g. 1999-07-26_1821_005410-74_007c.hpf
1039 //
1040 // In recent (2008/10/07) Mopra data, small negative ages (~10ms,
1041 // occasionally up to ~300ms) seem to be the norm, with both the
1042 // position timestamp and integration midpoint falling close to but
1043 // not on the integral second.
1044 if (cCycleNo == 2) {
1045 // We have to start with something!
1046 iMBuff->rateAge = dUTC;
1047
1048 } else {
1049 // Although we did not record the relevant position timestamp
1050 // explicitly, it can easily be deduced.
1051 double w = iMBuff->utc - utcDiff(cUTC, iMBuff->utc) -
1052 iMBuff->rateAge;
1053 dUTC = utcDiff(iMBuff->utc, w);
1054
1055 if (dUTC > 0.0) {
1056 iMBuff->rateAge = 0.0f;
1057 } else {
1058 iMBuff->rateAge = dUTC;
1059 }
1060 }
1061
1062 iMBuff->raRate = raRate;
1063 iMBuff->decRate = decRate;
1064 iMBuff->paRate = paRate;
1065 }
1066 }
1067
1068#ifdef PKSIO_DEBUG
1069 double avRate = sqrt(cAvRate[0]*cAvRate[0] + cAvRate[1]*cAvRate[1]);
1070 fprintf(stderr, "RA, Dec, Av & PA rates: %8.4f %8.4f %8.4f %8.4f "
1071 "pCode %d\n", raRate*R2D, decRate*R2D, avRate*R2D, paRate*R2D, pCode);
1072#endif
1073
1074
1075 // Compute the position of this beam for all bins.
1076 for (int idx = 0; idx < cNBin; idx++) {
1077 int jbuff = iBeamSel + cNBeamSel*idx;
1078
1079 cBuffer[jbuff].raRate = iMBuff->raRate;
1080 cBuffer[jbuff].decRate = iMBuff->decRate;
1081 cBuffer[jbuff].paRate = iMBuff->paRate;
1082
1083 double dUTC = utcDiff(cBuffer[jbuff].utc, thisUTC);
1084 if (dUTC > 100.0) {
1085 // Must have cycled through midnight.
1086 dUTC -= 86400.0;
1087 }
1088
1089 applyRate(cRA0, cDec0, thisRA, thisDec,
1090 cBuffer[jbuff].raRate, cBuffer[jbuff].decRate, dUTC,
1091 cBuffer[jbuff].ra, cBuffer[jbuff].dec);
1092
1093#ifdef PKSIO_DEBUG
1094 fprintf(stderr, "Intp (%d) ra, dec, UTC: %9.4f %9.4f %10.3f (pCode, "
1095 "age: %d %.1fs)\n", iMBuff->cycleNo, cBuffer[jbuff].ra*R2D,
1096 cBuffer[jbuff].dec*R2D, cBuffer[jbuff].utc, iMBuff->pCode,
1097 iMBuff->rateAge);
1098#endif
1099 }
1100 }
1101
1102
1103 if (cFlushing) {
1104 // Copy buffer location out one IF at a time.
1105 MBrec.extract(*iMBuff, cFlushIF);
1106 haveData = 1;
1107
1108#ifdef PKSIO_DEBUG
1109 fprintf(stderr, "Out:%4d%4d%3d%3d\n", MBrec.scanNo, MBrec.cycleNo,
1110 MBrec.beamNo, MBrec.IFno[0]);
1111#endif
1112
1113 // Signal that this IF in this buffer location has been flushed.
1114 iMBuff->IFno[cFlushIF] = 0;
1115
1116 if (cFlushIF == iMBuff->nIF - 1) {
1117 // Signal that all IFs in this buffer location have been flushed.
1118 iMBuff->nIF = 0;
1119
1120 // Stop cEOS being set when the next integration is read.
1121 iMBuff->cycleNo = 0;
1122
1123 } else {
1124 // Carry on flushing the other IFs.
1125 continue;
1126 }
1127
1128 // Has the whole buffer been flushed?
1129 if (cFlushBin == cNBin - 1) {
1130 if (cEOS || cEOF) {
1131 // Carry on flushing other buffers.
1132 cFlushIF = 0;
1133 continue;
1134 }
1135
1136 cFlushing = 0;
1137
1138 beamNo = int(cBaseline / 256.0);
1139 iBeamSel = cBeamSel[beamNo-1];
1140
1141 // Compute buffer number.
1142 iMBuff = cBuffer + iBeamSel;
1143 if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
1144 }
1145 }
1146
1147 if (!cFlushing) {
1148 // Buffer this MBrec.
1149 if ((cScanNo > iMBuff->scanNo) && iMBuff->IFno[0]) {
1150 // Sanity check on the number of IFs in the new scan.
1151 if (if_.n_if != cNIF) {
1152 sprintf(cMsg, "WARNING: Scan %d has %d IFs instead of %d, "
1153 "continuing.", cScanNo, if_.n_if, cNIF);
1154 logMsg(cMsg);
1155 }
1156 }
1157
1158 // Sanity check on incomplete integrations within a scan.
1159 if (iMBuff->nIF && (iMBuff->cycleNo != cCycleNo)) {
1160 // Force the incomplete integration to be flushed before proceeding.
1161 cFlushing = 1;
1162 continue;
1163 }
1164
1165#ifdef PKSIO_DEBUG
1166 fprintf(stderr, "Buf:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, cIFno);
1167#endif
1168
1169 // Store IF-independent parameters only for the first IF of a new cycle,
1170 // particularly because this is the only one for which the scan rates
1171 // are computed above.
1172 int firstIF = (iMBuff->nIF == 0);
1173 if (firstIF) {
1174 iMBuff->scanNo = cScanNo;
1175 iMBuff->cycleNo = cCycleNo;
1176
1177 // Times.
1178 strcpy(iMBuff->datobs, cDateObs);
1179 iMBuff->utc = cUTC;
1180 iMBuff->exposure = param_.intbase;
1181
1182 // Source identification.
1183 sprintf(iMBuff->srcName, "%-16.16s",
1184 names_.su_name + (cSrcNo-1)*16);
1185 iMBuff->srcName[16] = '\0';
1186 iMBuff->srcRA = doubles_.su_ra[cSrcNo-1];
1187 iMBuff->srcDec = doubles_.su_dec[cSrcNo-1];
1188
1189 // Rest frequency of the line of interest.
1190 iMBuff->restFreq = doubles_.rfreq;
1191 if (strncmp(names_.instrument, "ATPKSMB", 7) == 0) {
1192 // Fix the HI rest frequency recorded for Parkes multibeam data.
1193 double reffreq = doubles_.freq;
1194 double restfreq = doubles_.rfreq;
1195 if ((restfreq == 0.0 || fabs(restfreq - reffreq) == 0.0) &&
1196 fabs(reffreq - 1420.405752e6) < 100.0) {
1197 iMBuff->restFreq = 1420.405752e6;
1198 }
1199 }
1200
1201 // Observation type.
1202 int j;
1203 for (j = 0; j < 15; j++) {
1204 iMBuff->obsType[j] = names_.card[11+j];
1205 if (iMBuff->obsType[j] == '\'') break;
1206 }
1207 iMBuff->obsType[j] = '\0';
1208
1209 // Beam-dependent parameters.
1210 iMBuff->beamNo = beamNo;
1211
1212 // Beam position at the specified time.
1213 if (cSUpos) {
1214 // Non-ATNF data that does not store the position in (u,v,w).
1215 iMBuff->ra = doubles_.su_ra[cSrcNo-1];
1216 iMBuff->dec = doubles_.su_dec[cSrcNo-1];
1217 } else {
1218 iMBuff->ra = cU;
1219 iMBuff->dec = cV;
1220 }
1221 cPosUTC[iBeamSel] = cW;
1222 iMBuff->pCode = pCode;
1223
1224 // Store rates for next time.
1225 iMBuff->raRate = raRate;
1226 iMBuff->decRate = decRate;
1227 iMBuff->paRate = paRate;
1228 }
1229
1230 // IF-dependent parameters.
1231 int iIF = cIFno - 1;
1232 int startChan = cStartChan[iIF];
1233 int endChan = cEndChan[iIF];
1234 int refChan = cRefChan[iIF];
1235
1236 int nChan = abs(endChan - startChan) + 1;
1237
1238 iIFSel = cIFSel[iIF];
1239 if (iMBuff->IFno[iIFSel] == 0) {
1240 iMBuff->nIF++;
1241 iMBuff->IFno[iIFSel] = cIFno;
1242 } else {
1243 // Integration cycle written to the output file twice (the only known
1244 // example is 1999-05-22_1914_000-031805_03v.hpf).
1245 sprintf(cMsg, "WARNING: Integration cycle %d:%d, beam %2d, \n"
1246 " IF %d was duplicated.", cScanNo, cCycleNo-1,
1247 beamNo, cIFno);
1248 logMsg(cMsg);
1249 }
1250 iMBuff->nChan[iIFSel] = nChan;
1251 iMBuff->nPol[iIFSel] = cNPol[iIF];
1252
1253 iMBuff->fqRefPix[iIFSel] = doubles_.if_ref[iIF];
1254 iMBuff->fqRefVal[iIFSel] = doubles_.if_freq[iIF];
1255 iMBuff->fqDelt[iIFSel] =
1256 if_.if_invert[iIF] * fabs(doubles_.if_bw[iIF] /
1257 (if_.if_nfreq[iIF] - 1));
1258
1259 // Adjust for channel selection.
1260 if (iMBuff->fqRefPix[iIFSel] != refChan) {
1261 iMBuff->fqRefVal[iIFSel] +=
1262 (refChan - iMBuff->fqRefPix[iIFSel]) *
1263 iMBuff->fqDelt[iIFSel];
1264 iMBuff->fqRefPix[iIFSel] = refChan;
1265 }
1266
1267 if (endChan < startChan) {
1268 iMBuff->fqDelt[iIFSel] = -iMBuff->fqDelt[iIFSel];
1269 }
1270
1271
1272 // System temperature.
1273 int iBeam = beamNo - 1;
1274 int scq = sc_.sc_q;
1275 float TsysPol1 = sc_.sc_cal[scq*iBeam + 3];
1276 float TsysPol2 = sc_.sc_cal[scq*iBeam + 4];
1277 iMBuff->tsys[iIFSel][0] = TsysPol1*TsysPol1;
1278 iMBuff->tsys[iIFSel][1] = TsysPol2*TsysPol2;
1279
1280 // Calibration factor; may be changed later if the data is recalibrated.
1281 if (scq > 14) {
1282 // Will only be present for Parkes Multibeam or LBA data.
1283 iMBuff->calfctr[iIFSel][0] = sc_.sc_cal[scq*iBeam + 14];
1284 iMBuff->calfctr[iIFSel][1] = sc_.sc_cal[scq*iBeam + 15];
1285 } else {
1286 iMBuff->calfctr[iIFSel][0] = 0.0f;
1287 iMBuff->calfctr[iIFSel][1] = 0.0f;
1288 }
1289
1290 // Cross-polarization calibration factor (unknown to MBFITS).
1291 for (int j = 0; j < 2; j++) {
1292 iMBuff->xcalfctr[iIFSel][j] = 0.0f;
1293 }
1294
1295 // Baseline parameters (unknown to MBFITS).
1296 iMBuff->haveBase = 0;
1297
1298 // Data (always present in MBFITS).
1299 iMBuff->haveSpectra = 1;
1300
1301 // Flag: bit 0 set if off source.
1302 // bit 1 set if loss of sync in A polarization.
1303 // bit 2 set if loss of sync in B polarization.
1304 unsigned char rpflag =
1305 (unsigned char)(sc_.sc_cal[scq*iBeam + 12] + 0.5f);
1306
1307 // The baseline flag may be set independently.
1308 if (rpflag == 0) rpflag = cFlag;
1309
1310 // Copy and scale data.
1311 int inc = 2 * if_.if_nstok[iIF];
1312 if (endChan < startChan) inc = -inc;
1313
1314 float TsysF;
1315 iMBuff->spectra[iIFSel] = iMBuff->spectra[0] + cChanOff[iIF];
1316 iMBuff->flagged[iIFSel] = iMBuff->flagged[0] + cChanOff[iIF];
1317
1318 float *spectra = iMBuff->spectra[iIFSel];
1319 unsigned char *flagged = iMBuff->flagged[iIFSel];
1320 for (int ipol = 0; ipol < cNPol[iIF]; ipol++) {
1321 if (sc_.sc_cal[scq*iBeam + 3 + ipol] > 0.0f) {
1322 // The correlator has already applied the calibration.
1323 TsysF = 1.0f;
1324 } else {
1325 // The correlator has normalized cVis[k] to a Tsys of 500K.
1326 TsysF = iMBuff->tsys[iIFSel][ipol] / 500.0f;
1327 }
1328
1329 int k = 2 * (if_.if_nstok[iIF]*(startChan - 1) + ipol);
1330 for (int ichan = 0; ichan < nChan; ichan++) {
1331 *(spectra++) = TsysF * cVis[k];
1332 *(flagged++) = rpflag;
1333 k += inc;
1334 }
1335 }
1336
1337 if (cHaveXPol[iIF]) {
1338 int k = 2 * (3*(startChan - 1) + 2);
1339 iMBuff->xpol[iIFSel] = iMBuff->xpol[0] + cXpolOff[iIF];
1340 float *xpol = iMBuff->xpol[iIFSel];
1341 for (int ichan = 0; ichan < nChan; ichan++) {
1342 *(xpol++) = cVis[k];
1343 *(xpol++) = cVis[k+1];
1344 k += inc;
1345 }
1346 }
1347
1348
1349 // Calibration factor applied to the data by the correlator.
1350 if (scq > 14) {
1351 // Will only be present for Parkes Multibeam or LBA data.
1352 iMBuff->tcal[iIFSel][0] = sc_.sc_cal[scq*iBeam + 14];
1353 iMBuff->tcal[iIFSel][1] = sc_.sc_cal[scq*iBeam + 15];
1354 } else {
1355 iMBuff->tcal[iIFSel][0] = 0.0f;
1356 iMBuff->tcal[iIFSel][1] = 0.0f;
1357 }
1358
1359 if (firstIF) {
1360 if (sc_.sc_ant <= anten_.nant) {
1361 // No extra syscal information present.
1362 iMBuff->extraSysCal = 0;
1363 iMBuff->azimuth = 0.0f;
1364 iMBuff->elevation = 0.0f;
1365 iMBuff->parAngle = 0.0f;
1366 iMBuff->focusAxi = 0.0f;
1367 iMBuff->focusTan = 0.0f;
1368 iMBuff->focusRot = 0.0f;
1369 iMBuff->temp = 0.0f;
1370 iMBuff->pressure = 0.0f;
1371 iMBuff->humidity = 0.0f;
1372 iMBuff->windSpeed = 0.0f;
1373 iMBuff->windAz = 0.0f;
1374 strcpy(iMBuff->tcalTime, " ");
1375 iMBuff->refBeam = 0;
1376
1377 } else {
1378 // Additional information for Parkes Multibeam data.
1379 int iOff = scq*(sc_.sc_ant - 1) - 1;
1380 iMBuff->extraSysCal = 1;
1381
1382 iMBuff->azimuth = sc_.sc_cal[iOff + 2];
1383 iMBuff->elevation = sc_.sc_cal[iOff + 3];
1384 iMBuff->parAngle = sc_.sc_cal[iOff + 4];
1385
1386 iMBuff->focusAxi = sc_.sc_cal[iOff + 5] * 1e-3;
1387 iMBuff->focusTan = sc_.sc_cal[iOff + 6] * 1e-3;
1388 iMBuff->focusRot = sc_.sc_cal[iOff + 7];
1389
1390 iMBuff->temp = sc_.sc_cal[iOff + 8];
1391 iMBuff->pressure = sc_.sc_cal[iOff + 9];
1392 iMBuff->humidity = sc_.sc_cal[iOff + 10];
1393 iMBuff->windSpeed = sc_.sc_cal[iOff + 11];
1394 iMBuff->windAz = sc_.sc_cal[iOff + 12];
1395
1396 char *tcalTime = iMBuff->tcalTime;
1397 sprintf(tcalTime, "%-16.16s", (char *)(&sc_.sc_cal[iOff+13]));
1398 tcalTime[16] = '\0';
1399
1400#ifndef AIPS_LITTLE_ENDIAN
1401 // Do byte swapping on the ASCII date string.
1402 for (int j = 0; j < 16; j += 4) {
1403 char ctmp;
1404 ctmp = tcalTime[j];
1405 tcalTime[j] = tcalTime[j+3];
1406 tcalTime[j+3] = ctmp;
1407 ctmp = tcalTime[j+1];
1408 tcalTime[j+1] = tcalTime[j+2];
1409 tcalTime[j+2] = ctmp;
1410 }
1411#endif
1412
1413 // Reference beam number.
1414 float refbeam = sc_.sc_cal[iOff + 17];
1415 if (refbeam > 0.0f || refbeam < 100.0f) {
1416 iMBuff->refBeam = int(refbeam);
1417 } else {
1418 iMBuff->refBeam = 0;
1419 }
1420 }
1421 }
1422 }
1423 }
1424
1425 return 0;
1426}
1427
1428//-------------------------------------------------------- MBFITSreader::rpget
1429
1430// Read the next data record from the RPFITS file.
1431
1432int MBFITSreader::rpget(int syscalonly, int &EOS)
1433{
1434 EOS = 0;
1435
1436 int retries = 0;
1437
1438 // Allow 10 read errors.
1439 int numErr = 0;
1440
1441 int jstat = 0;
1442 while (numErr < 10) {
1443 int lastjstat = jstat;
1444
1445 switch(rpfitsin(jstat)) {
1446 case -1:
1447 // Read failed; retry.
1448 numErr++;
1449 logMsg("WARNING: RPFITS read failed - retrying.");
1450 jstat = 0;
1451 break;
1452
1453 case 0:
1454 // Successful read.
1455 if (lastjstat == 0) {
1456 if (cBaseline == -1) {
1457 // Syscal data.
1458 if (syscalonly) {
1459 return 0;
1460 }
1461
1462 } else {
1463 if (!syscalonly) {
1464 return 0;
1465 }
1466 }
1467 }
1468
1469 // Last operation was to read header or FG table; now read data.
1470 break;
1471
1472 case 1:
1473 // Encountered header while trying to read data; read it.
1474 EOS = 1;
1475 jstat = -1;
1476 break;
1477
1478 case 2:
1479 // End of scan; read past it.
1480 jstat = 0;
1481 break;
1482
1483 case 3:
1484 // End-of-file; retry applies to real-time mode.
1485 if (retries++ >= cRetry) {
1486 return -1;
1487 }
1488
1489 sleep(10);
1490 jstat = 0;
1491 break;
1492
1493 case 4:
1494 // Encountered FG table while trying to read data; read it.
1495 jstat = -1;
1496 break;
1497
1498 case 5:
1499 // Illegal data at end of block after close/reopen operation; retry.
1500 jstat = 0;
1501 break;
1502
1503 default:
1504 // Shouldn't reach here.
1505 sprintf(cMsg, "WARNING: Unrecognized RPFITSIN return code: %d "
1506 "(retrying).", jstat);
1507 logMsg(cMsg);
1508 jstat = 0;
1509 break;
1510 }
1511 }
1512
1513 logMsg("ERROR: RPFITS read failed too many times.");
1514 return 2;
1515}
1516
1517//----------------------------------------------------- MBFITSreader::rpfitsin
1518
1519// Wrapper around RPFITSIN that reports errors. Returned RPFITSIN subroutine
1520// arguments are captured as MBFITSreader member variables.
1521
1522int MBFITSreader::rpfitsin(int &jstat)
1523
1524{
1525 rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
1526 &cBin, &cIFno, &cSrcNo);
1527
1528 // Handle messages from RPFITSIN.
1529 if (names_.errmsg[0] != ' ') {
1530 int i;
1531 for (i = 80; i > 0; i--) {
1532 if (names_.errmsg[i-1] != ' ') break;
1533 }
1534
1535 sprintf(cMsg, "WARNING: Cycle %d:%03d, RPFITSIN reported -\n"
1536 " %.*s", cScanNo, cCycleNo, i, names_.errmsg);
1537 logMsg(cMsg);
1538 }
1539
1540 return jstat;
1541}
1542
1543//------------------------------------------------------- MBFITSreader::fixPos
1544
1545// Check and, if necessary, repair a position timestamp.
1546//
1547// Problems with the position timestamp manifest themselves via the scan rate:
1548//
1549// 1) Zero scan rate pairs, 1997/02/28 to 1998/01/07
1550//
1551// These occur because the position timestamp for the first integration
1552// of the pair is erroneous; the value recorded is t/1000, where t is the
1553// true value.
1554// Earliest known: 97-02-28_1725_132653-42_258a.hpf
1555// Latest known: 98-01-02_1923_095644-50_165c.hpf
1556// (time range chosen to encompass observing runs).
1557//
1558// 2) Slow-fast scan rate pairs (0.013 - 0.020 deg/s),
1559// 1997/03/28 to 1998/01/07.
1560//
1561// The UTC position timestamp is 1.0s later than it should be (never
1562// earlier), almost certainly arising from an error in the telescope
1563// control system.
1564// Earliest known: 97-03-28_0150_010420-74_008d.hpf
1565// Latest known: 98-01-04_1502_065150-02_177c.hpf
1566// (time range chosen to encompass observing runs).
1567//
1568// 3) Slow-fast scan rate pairs (0.015 - 0.018 deg/s),
1569// 1999/05/20 to 2001/07/12 (HIPASS and ZOA),
1570// 2001/09/02 to 2001/12/04 (HIPASS and ZOA),
1571// 2002/03/28 to 2002/05/13 (ZOA only),
1572// 2003/04/26 to 2003/06/09 (ZOA only).
1573// Earliest known: 1999-05-20_1818_175720-50_297e.hpf
1574// Latest known: 2001-12-04_1814_065531p14_173e.hpf (HIPASS)
1575// 2003-06-09_1924_352-085940_-6c.hpf (ZOA)
1576//
1577// Caused by the Linux signalling NaN problem. IEEE "signalling" NaNs
1578// are silently transformed to "quiet" NaNs during assignment by setting
1579// bit 22. This affected RPFITS because of its use of VAX-format
1580// floating-point numbers which, with their permuted bytes, may sometimes
1581// appear as signalling NaNs.
1582//
1583// The problem arose when the linux correlator came online and was
1584// fixed with a workaround to the RPFITS library (repeated episodes
1585// are probably due to use of an older version of the library). It
1586// should not have affected the data significantly because of the
1587// low relative error, which ranges from 0.0000038 to 0.0000076, but
1588// it is important for the computation of scan rates which requires
1589// taking the difference of two large UTC timestamps, one or other
1590// of which will have 0.5s added to it.
1591//
1592// The return value identifies which, if any, of these problems was repaired.
1593
1594int MBFITSreader::fixw(
1595 const char *datobs,
1596 int cycleNo,
1597 int beamNo,
1598 double avRate[2],
1599 double thisRA,
1600 double thisDec,
1601 double thisUTC,
1602 double nextRA,
1603 double nextDec,
1604 float &nextUTC)
1605{
1606 if (strcmp(datobs, "2003-06-09") > 0) {
1607 return 0;
1608
1609 } else if (strcmp(datobs, "1998-01-07") <= 0) {
1610 if (nextUTC < thisUTC && (nextUTC + 86400.0) > (thisUTC + 600.0)) {
1611 // Possible scaling problem.
1612 double diff = nextUTC*1000.0 - thisUTC;
1613 if (0.0 < diff && diff < 600.0) {
1614 nextUTC *= 1000.0;
1615 return 1;
1616 } else {
1617 // Irreparable.
1618 return -1;
1619 }
1620 }
1621
1622 if (cycleNo > 2) {
1623 if (beamNo == 1) {
1624 // This test is only reliable for beam 1.
1625 double dUTC = nextUTC - thisUTC;
1626 if (dUTC < 0.0) dUTC += 86400.0;
1627
1628 // Guard against RA cycling through 24h in either direction.
1629 if (fabs(nextRA - thisRA) > PI) {
1630 if (nextRA < thisRA) {
1631 nextRA += TWOPI;
1632 } else {
1633 nextRA -= TWOPI;
1634 }
1635 }
1636
1637 double dRA = (nextRA - thisRA) * cos(nextDec);
1638 double dDec = nextDec - thisDec;
1639 double arc = sqrt(dRA*dRA + dDec*dDec);
1640
1641 double averate = sqrt(avRate[0]*avRate[0] + avRate[1]*avRate[1]);
1642 double diff1 = fabs(averate - arc/(dUTC-1.0));
1643 double diff2 = fabs(averate - arc/dUTC);
1644 if ((diff1 < diff2) && (diff1 < 0.05*averate)) {
1645 nextUTC -= 1.0;
1646 cCode5 = cycleNo;
1647 return 2;
1648 } else {
1649 cCode5 = 0;
1650 }
1651
1652 } else {
1653 if (cycleNo == cCode5) {
1654 nextUTC -= 1.0;
1655 return 2;
1656 }
1657 }
1658 }
1659
1660 } else if ((strcmp(datobs, "1999-05-20") >= 0 &&
1661 strcmp(datobs, "2001-07-12") <= 0) ||
1662 (strcmp(datobs, "2001-09-02") >= 0 &&
1663 strcmp(datobs, "2001-12-04") <= 0) ||
1664 (strcmp(datobs, "2002-03-28") >= 0 &&
1665 strcmp(datobs, "2002-05-13") <= 0) ||
1666 (strcmp(datobs, "2003-04-26") >= 0 &&
1667 strcmp(datobs, "2003-06-09") <= 0)) {
1668 // Signalling NaN problem, e.g. 1999-07-26_1839_011106-74_009c.hpf.
1669 // Position timestamps should always be an integral number of seconds.
1670 double resid = nextUTC - int(nextUTC);
1671 if (resid == 0.5) {
1672 nextUTC -= 0.5;
1673 return 3;
1674 }
1675 }
1676
1677 return 0;
1678}
1679
1680//-------------------------------------------------------- MBFITSreader::close
1681
1682// Close the input file.
1683
1684void MBFITSreader::close(void)
1685{
1686 if (cMBopen) {
1687 int jstat = 1;
1688 rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
1689 &cBin, &cIFno, &cSrcNo);
1690
1691 if (cBeams) delete [] cBeams;
1692 if (cIFs) delete [] cIFs;
1693 if (cNChan) delete [] cNChan;
1694 if (cNPol) delete [] cNPol;
1695 if (cHaveXPol) delete [] cHaveXPol;
1696 if (cStartChan) delete [] cStartChan;
1697 if (cEndChan) delete [] cEndChan;
1698 if (cRefChan) delete [] cRefChan;
1699
1700 if (cVis) delete [] cVis;
1701 if (cWgt) delete [] cWgt;
1702
1703 if (cBeamSel) delete [] cBeamSel;
1704 if (cIFSel) delete [] cIFSel;
1705 if (cChanOff) delete [] cChanOff;
1706 if (cXpolOff) delete [] cXpolOff;
1707 if (cBuffer) delete [] cBuffer;
1708 if (cPosUTC) delete [] cPosUTC;
1709
1710 cMBopen = 0;
1711 }
1712}
1713
1714//-------------------------------------------------------------------- utcDiff
1715
1716// Subtract two UTCs (s) allowing for any plausible number of cycles through
1717// 86400s, returning a result in the range [-43200, +43200]s.
1718
1719double MBFITSreader::utcDiff(double utc1, double utc2)
1720{
1721 double diff = utc1 - utc2;
1722
1723 if (diff > 43200.0) {
1724 diff -= 86400.0;
1725 while (diff > 43200.0) diff -= 86400.0;
1726 } else if (diff < -43200.0) {
1727 diff += 86400.0;
1728 while (diff < -43200.0) diff += 86400.0;
1729 }
1730
1731 return diff;
1732}
1733
1734//------------------------------------------------------- scanRate & applyRate
1735
1736// Compute and apply the scan rate corrected for grid convergence. (ra0,dec0)
1737// are the coordinates of the central beam, assumed to be the tracking centre.
1738// The rate computed in RA will be a rate of change of angular distance in the
1739// direction of increasing RA at the position of the central beam. Similarly
1740// for declination. Angles in radian, time in s.
1741
1742void MBFITSreader::scanRate(
1743 double ra0,
1744 double dec0,
1745 double ra1,
1746 double dec1,
1747 double ra2,
1748 double dec2,
1749 double dt,
1750 double &raRate,
1751 double &decRate)
1752{
1753 // Transform to a system where the central beam lies on the equator at 12h.
1754 eulerx(ra1, dec1, ra0+HALFPI, -dec0, -HALFPI, ra1, dec1);
1755 eulerx(ra2, dec2, ra0+HALFPI, -dec0, -HALFPI, ra2, dec2);
1756
1757 raRate = (ra2 - ra1) / dt;
1758 decRate = (dec2 - dec1) / dt;
1759}
1760
1761
1762void MBFITSreader::applyRate(
1763 double ra0,
1764 double dec0,
1765 double ra1,
1766 double dec1,
1767 double raRate,
1768 double decRate,
1769 double dt,
1770 double &ra2,
1771 double &dec2)
1772{
1773 // Transform to a system where the central beam lies on the equator at 12h.
1774 eulerx(ra1, dec1, ra0+HALFPI, -dec0, -HALFPI, ra1, dec1);
1775
1776 ra2 = ra1 + (raRate * dt);
1777 dec2 = dec1 + (decRate * dt);
1778
1779 // Transform back.
1780 eulerx(ra2, dec2, -HALFPI, dec0, ra0+HALFPI, ra2, dec2);
1781}
1782
1783//--------------------------------------------------------------------- eulerx
1784
1785void MBFITSreader::eulerx(
1786 double lng0,
1787 double lat0,
1788 double phi0,
1789 double theta,
1790 double phi,
1791 double &lng1,
1792 double &lat1)
1793
1794// Applies the Euler angle based transformation of spherical coordinates.
1795//
1796// phi0 Longitude of the ascending node in the old system, radians. The
1797// ascending node is the point of intersection of the equators of
1798// the two systems such that the equator of the new system crosses
1799// from south to north as viewed in the old system.
1800//
1801// theta Angle between the poles of the two systems, radians. THETA is
1802// positive for a positive rotation about the ascending node.
1803//
1804// phi Longitude of the ascending node in the new system, radians.
1805
1806{
1807 // Compute intermediaries.
1808 double lng0p = lng0 - phi0;
1809 double slng0p = sin(lng0p);
1810 double clng0p = cos(lng0p);
1811 double slat0 = sin(lat0);
1812 double clat0 = cos(lat0);
1813 double ctheta = cos(theta);
1814 double stheta = sin(theta);
1815
1816 double x = clat0*clng0p;
1817 double y = clat0*slng0p*ctheta + slat0*stheta;
1818
1819 // Longitude in the new system.
1820 if (x != 0.0 || y != 0.0) {
1821 lng1 = phi + atan2(y, x);
1822 } else {
1823 // Longitude at the poles in the new system is consistent with that
1824 // specified in the old system.
1825 lng1 = phi + lng0p;
1826 }
1827 lng1 = fmod(lng1, TWOPI);
1828 if (lng1 < 0.0) lng1 += TWOPI;
1829
1830 lat1 = asin(slat0*ctheta - clat0*stheta*slng0p);
1831}
Note: See TracBrowser for help on using the repository browser.