source: trunk/external-alma/atnf/PKSIO/MBFITSreader.cc@ 2965

Last change on this file since 2965 was 1868, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): atnf

Description: Describe your changes here...

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