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

Last change on this file since 3118 was 3089, checked in by Kana Sugimoto, 9 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

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