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

Last change on this file since 1375 was 1372, checked in by mar637, 17 years ago

latest update from livedata CVS. This has the Hobart position fix
File size: 35.8 KB
RevLine 
[1325]1//#---------------------------------------------------------------------------
2//# MBFITSreader.cc: ATNF single-dish RPFITS reader.
3//#---------------------------------------------------------------------------
4//# Copyright (C) 2000-2007
5//# Mark Calabretta, ATNF
6//#
7//# This library is free software; you can redistribute it and/or modify it
8//# under the terms of the GNU Library General Public License as published by
9//# the Free Software Foundation; either version 2 of the License, or (at your
10//# option) any later version.
11//#
12//# This library is distributed in the hope that it will be useful, but WITHOUT
13//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15//# License for more details.
16//#
17//# You should have received a copy of the GNU Library General Public License
18//# along with this library; if not, write to the Free Software Foundation,
19//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20//#
21//# Correspondence concerning this software should be addressed as follows:
22//# Internet email: mcalabre@atnf.csiro.au.
23//# Postal address: Dr. Mark Calabretta,
24//# Australia Telescope National Facility,
25//# P.O. Box 76,
26//# Epping, NSW, 2121,
27//# AUSTRALIA
28//#
[1372]29//# $Id: MBFITSreader.cc,v 19.34 2007/07/02 06:12:18 cal103 Exp $
[1325]30//#---------------------------------------------------------------------------
31//# The MBFITSreader class reads single dish RPFITS files (such as Parkes
32//# Multibeam MBFITS files).
33//#
34//# Original: 2000/07/28 Mark Calabretta
35//#---------------------------------------------------------------------------
36
37#include <atnf/PKSIO/MBFITSreader.h>
38#include <atnf/PKSIO/PKSMBrecord.h>
39
40#include <RPFITS.h>
41
42#include <casa/math.h>
43#include <casa/iostream.h>
44#include <casa/stdio.h>
45#include <casa/stdlib.h>
46#include <casa/string.h>
47#include <unistd.h>
48
49using namespace std;
50
51// Numerical constants.
52const double PI = 3.141592653589793238462643;
53const double TWOPI = 2.0 * PI;
54
55//------------------------------------------------- MBFITSreader::MBFITSreader
56
57// Default constructor.
58
59MBFITSreader::MBFITSreader(
60 const int retry,
61 const int interpolate)
62{
63 cRetry = retry;
64 if (cRetry > 10) {
65 cRetry = 10;
66 }
67
68 cInterp = interpolate;
69 if (cInterp < 0 || cInterp > 2) {
70 cInterp = 1;
71 }
72
73 // Initialize pointers.
74 cBeams = 0x0;
75 cIFs = 0x0;
76 cNChan = 0x0;
77 cNPol = 0x0;
78 cHaveXPol = 0x0;
79 cStartChan = 0x0;
80 cEndChan = 0x0;
81 cRefChan = 0x0;
82
83 cVis = 0x0;
84 cWgt = 0x0;
85
86 cBeamSel = 0x0;
87 cIFSel = 0x0;
88 cChanOff = 0x0;
89 cXpolOff = 0x0;
90 cBuffer = 0x0;
91 cPosUTC = 0x0;
92
93 cMBopen = 0;
94}
95
96//------------------------------------------------ MBFITSreader::~MBFITSreader
97
98// Destructor.
99
100MBFITSreader::~MBFITSreader()
101{
102 close();
103}
104
105//--------------------------------------------------------- MBFITSreader::open
106
107// Open the RPFITS file for reading.
108
109int MBFITSreader::open(
110 char *rpname,
111 int &nBeam,
112 int* &beams,
113 int &nIF,
114 int* &IFs,
115 int* &nChan,
116 int* &nPol,
117 int* &haveXPol,
118 int &haveBase,
119 int &haveSpectra,
120 int &extraSysCal)
121{
122 if (cMBopen) {
123 close();
124 }
125
126 strcpy(names_.file, rpname);
127
128 // Open the RPFITS file.
129 int jstat = -3;
130 rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
131 &cBin, &cIFno, &cSrcNo);
132
133 if (jstat) {
[1372]134 fprintf(stderr, "ERROR, failed to open MBFITS file: %s\n", rpname);
[1325]135 return 1;
136 }
137
138 cMBopen = 1;
139
140 // Tell RPFITSIN that we want the OBSTYPE card.
141 int j;
142 param_.ncard = 1;
143 for (j = 0; j < 80; j++) {
144 names_.card[j] = ' ';
145 }
146 strncpy(names_.card, "OBSTYPE", 7);
147
148 // Read the first header.
149 jstat = -1;
150 rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
151 &cBin, &cIFno, &cSrcNo);
152
153 if (jstat) {
[1372]154 fprintf(stderr, "ERROR, failed to read MBFITS header: %s\n", rpname);
[1325]155 close();
156 return 1;
157 }
158
159 // Mopra data has some peculiarities.
160 cMopra = strncmp(names_.instrument, "ATMOPRA", 7) == 0;
161
[1372]162 // Non-ATNF data may not store the position in (u,v,w).
163 if (strncmp(names_.sta, "tid", 3) == 0) {
164 fprintf(stderr, "WARNING, found Tidbinbilla data");
165 cSUpos = 1;
166 } else if (strncmp(names_.sta, "HOB", 3) == 0) {
167 fprintf(stderr, "WARNING, found Hobart data");
168 cSUpos = 1;
169 } else if (strncmp(names_.sta, "CED", 3) == 0) {
170 fprintf(stderr, "WARNING, found Ceduna data");
171 cSUpos = 1;
172 } else {
173 cSUpos = 0;
174 }
175
176 if (cSUpos) {
177 fprintf(stderr, ", using telescope position from SU table.\n");
[1325]178 cInterp = 0;
179 }
180
181
182 // Find the maximum beam number.
183 cNBeam = 0;
184 for (int iBeam = 0; iBeam < anten_.nant; iBeam++) {
185 if (anten_.ant_num[iBeam] > cNBeam) {
186 cNBeam = anten_.ant_num[iBeam];
187 }
188 }
189
190 if (cNBeam <= 0) {
[1372]191 fprintf(stderr, "ERROR, couldn't determine number of beams.\n");
[1325]192 close();
193 return 1;
194 }
195
196 // Construct the beam mask.
197 cBeams = new int[cNBeam];
198 for (int iBeam = 0; iBeam < cNBeam; iBeam++) {
199 cBeams[iBeam] = 0;
200 }
201
202 // ...beams present in the data.
203 for (int iBeam = 0; iBeam < anten_.nant; iBeam++) {
204 cBeams[anten_.ant_num[iBeam] - 1] = 1;
205 }
206
207 // Passing back the address of the array allows PKSFITSreader::select() to
208 // modify its elements directly.
209 nBeam = cNBeam;
210 beams = cBeams;
211
212
213 // Number of IFs.
214 cNIF = if_.n_if;
215 cIFs = new int[cNIF];
216 for (int iIF = 0; iIF < cNIF; iIF++) {
217 cIFs[iIF] = 1;
218 }
219
220 // Passing back the address of the array allows PKSFITSreader::select() to
221 // modify its elements directly.
222 nIF = cNIF;
223 IFs = cIFs;
224
225
226 // Number of channels and polarizations.
227 cNChan = new int[cNIF];
228 cNPol = new int[cNIF];
229 cHaveXPol = new int[cNIF];
230 cGetXPol = 0;
231
232 int maxProd = 0;
233 for (int iIF = 0; iIF < cNIF; iIF++) {
234 cNChan[iIF] = if_.if_nfreq[iIF];
235 cNPol[iIF] = if_.if_nstok[iIF];
236 cNChan[iIF] -= cNChan[iIF]%2;
237
238 // Do we have cross-polarization data?
239 if ((cHaveXPol[iIF] = cNPol[iIF] > 2)) {
240 // Cross-polarization data is handled separately.
241 cNPol[iIF] = 2;
242
243 // Default is to get it if we have it.
244 cGetXPol = 1;
245 }
246
247 // Maximum number of spectral products in any IF.
248 int nProd = if_.if_nfreq[iIF] * if_.if_nstok[iIF];
249 if (maxProd < nProd) maxProd = nProd;
250 }
251
252 // Allocate memory for RPFITSIN subroutine arguments.
253 if (cVis) delete [] cVis;
254 if (cWgt) delete [] cWgt;
255 cVis = new float[2*maxProd];
256 cWgt = new float[maxProd];
257
258 nChan = cNChan;
259 nPol = cNPol;
260 haveXPol = cHaveXPol;
261
262
263 // Default channel range selection.
264 cStartChan = new int[cNIF];
265 cEndChan = new int[cNIF];
266 cRefChan = new int[cNIF];
267
268 for (int iIF = 0; iIF < cNIF; iIF++) {
269 cStartChan[iIF] = 1;
270 cEndChan[iIF] = cNChan[iIF];
271 cRefChan[iIF] = cNChan[iIF]/2 + 1;
272 }
273
274 cGetSpectra = 1;
275
276
277 // No baseline parameters in MBFITS.
278 haveBase = 0;
279
280 // Always have spectra in MBFITS.
281 haveSpectra = cHaveSpectra = 1;
282
283
284 // Integration cycle time (s).
285 cIntTime = param_.intime;
286
287 // Can't deduce binning mode till later.
288 cNBin = 0;
289
290
291 // Read the first syscal record.
292 if (rpget(1, cEOS)) {
[1372]293 fprintf(stderr, "ERROR, failed to read first syscal record.\n");
[1325]294 close();
295 return 1;
296 }
297
298 // Additional information for Parkes Multibeam data?
299 extraSysCal = (sc_.sc_ant > anten_.nant);
300
301
302 cFirst = 1;
303 cEOF = 0;
304 cFlushing = 0;
305
306 return 0;
307}
308
309//---------------------------------------------------- MBFITSreader::getHeader
310
311// Get parameters describing the data.
312
313int MBFITSreader::getHeader(
314 char observer[32],
315 char project[32],
316 char telescope[32],
317 double antPos[3],
318 char obsType[32],
319 float &equinox,
320 char radecsys[32],
321 char dopplerFrame[32],
322 char datobs[32],
323 double &utc,
324 double &refFreq,
325 double &bandwidth)
326{
327 if (!cMBopen) {
[1372]328 fprintf(stderr, "ERROR, an MBFITS file has not been opened.\n");
[1325]329 return 1;
330 }
331
332 sprintf(observer, "%-16.16s", names_.rp_observer);
333 sprintf(project, "%-16.16s", names_.object);
334 sprintf(telescope, "%-16.16s", names_.instrument);
335
336 // Observatory coordinates (ITRF), in m.
337 antPos[0] = doubles_.x[0];
338 antPos[1] = doubles_.y[0];
339 antPos[2] = doubles_.z[0];
340
341 // This is the only sure way to identify the telescope, maybe.
342 if (strncmp(names_.sta, "MB0", 3) == 0) {
343 // Parkes Multibeam.
344 sprintf(telescope, "%-16.16s", "ATPKSMB");
345 antPos[0] = -4554232.087;
346 antPos[1] = 2816759.046;
347 antPos[2] = -3454035.950;
348 } else if (strncmp(names_.sta, "HOH", 3) == 0) {
349 // Parkes HOH receiver.
350 sprintf(telescope, "%-16.16s", "ATPKSHOH");
351 antPos[0] = -4554232.087;
352 antPos[1] = 2816759.046;
353 antPos[2] = -3454035.950;
354 } else if (strncmp(names_.sta, "CA0", 3) == 0) {
355 // An ATCA antenna, use the array centre position.
356 sprintf(telescope, "%-16.16s", "ATCA");
357 antPos[0] = -4750915.837;
358 antPos[1] = 2792906.182;
359 antPos[2] = -3200483.747;
360 } else if (strncmp(names_.sta, "MOP", 3) == 0) {
361 // Mopra.
362 sprintf(telescope, "%-16.16s", "ATMOPRA");
363 antPos[0] = -4682768.630;
364 antPos[1] = 2802619.060;
365 antPos[2] = -3291759.900;
366 } else if (strncmp(names_.sta, "HOB", 3) == 0) {
367 // Hobart.
368 sprintf(telescope, "%-16.16s", "HOBART");
369 antPos[0] = -3950236.735;
370 antPos[1] = 2522347.567;
371 antPos[2] = -4311562.569;
372 } else if (strncmp(names_.sta, "CED", 3) == 0) {
373 // Ceduna.
374 sprintf(telescope, "%-16.16s", "CEDUNA");
375 antPos[0] = -3749943.657;
376 antPos[1] = 3909017.709;
377 antPos[2] = -3367518.309;
378 } else if (strncmp(names_.sta, "tid", 3) == 0) {
379 // DSS.
380 sprintf(telescope, "%-16.16s", "DSS-43");
381 antPos[0] = -4460894.727;
382 antPos[1] = 2682361.530;
383 antPos[2] = -3674748.424;
384 }
385
386 // Observation type.
387 int j;
388 for (j = 0; j < 31; j++) {
389 obsType[j] = names_.card[11+j];
390 if (obsType[j] == '\'') break;
391 }
392 obsType[j] = '\0';
393
394 // Coordinate frames.
395 equinox = 2000.0f;
396 strcpy(radecsys, "FK5");
397 strcpy(dopplerFrame, "TOPOCENT");
398
399 // Time at start of observation.
400 sprintf(datobs, "%-10.10s", names_.datobs);
401 utc = cUTC;
402
403 // Spectral parameters.
404 refFreq = doubles_.if_freq[0];
405 bandwidth = doubles_.if_bw[0];
406
407 return 0;
408}
409
410//-------------------------------------------------- MBFITSreader::getFreqInfo
411
412// Get frequency parameters for each IF.
413
414int MBFITSreader::getFreqInfo(
415 int &nIF,
416 double* &startFreq,
417 double* &endFreq)
418{
419 // This is RPFITS - can't do it!
420 return 1;
421}
422
423//---------------------------------------------------- MBFITSreader::findRange
424
425// Find the range of the data selected in time and position.
426
427int MBFITSreader::findRange(
428 int &nRow,
429 int &nSel,
430 char dateSpan[2][32],
431 double utcSpan[2],
432 double* &positions)
433{
434 // This is RPFITS - can't do it!
435 return 1;
436}
437
438//--------------------------------------------------------- MBFITSreader::read
439
440// Read the next data record.
441
442int MBFITSreader::read(
443 PKSMBrecord &MBrec)
444{
445 int beamNo = -1;
446 int haveData, status;
447 PKSMBrecord *iMBuff = 0x0;
448
449 if (!cMBopen) {
[1372]450 fprintf(stderr, "ERROR, an MBFITS file has not been opened.\n");
[1325]451 return 1;
452 }
453
454 // Positions recorded in the input records do not coincide with the midpoint
455 // of the integration and hence the input must be buffered so that true
456 // positions may be interpolated.
457 //
458 // On the first call nBeamSel buffers of length nBin, are allocated and
459 // filled, where nBin is the number of time bins.
460 //
461 // The input records for binned, single beam data with multiple simultaneous
462 // IFs are ordered by IF within each integration rather than by bin number
463 // and hence are not in time order. No multibeam data exists with
464 // nBin > 1 but the likelihood that the input records would be in beam/IF
465 // order and the requirement that output records be in time order would
466 // force an elaborate double-buffering system and we do not support it.
467 //
468 // Once all buffers are filled, the next record for each beam pertains to
469 // the next integration and should contain new position information allowing
470 // the proper position for each spectrum in the buffer to be interpolated.
471 // The buffers are then flushed in time order. For single beam data there
472 // is only one buffer and reads from the MBFITS file are suspended while the
473 // flush is in progress. For multibeam data each buffer is of unit length
474 // so the flush completes immediately and the new record takes its place.
475
476 haveData = 0;
477 while (!haveData) {
478 int iBeamSel = -1, iIFSel = -1;
479
480 if (!cFlushing) {
481 if (cEOF) {
482 return -1;
483 }
484
485 // Read the next record.
486 if ((status = rpget(0, cEOS)) == -1) {
487 // EOF.
488 cEOF = 1;
489 cFlushing = 1;
490 cFlushBin = 0;
491 cFlushIF = 0;
492
493#ifdef PKSIO_DEBUG
494 printf("End-of-file detected, flushing last scan.\n");
495#endif
496
497 } else if (status) {
498 // IO error.
499 return 1;
500
501 } else {
502 if (cFirst) {
503 // First data; cBeamSel[] stores the buffer index for each beam.
504 cNBeamSel = 0;
505 cBeamSel = new int[cNBeam];
506
507 for (int iBeam = 0; iBeam < cNBeam; iBeam++) {
508 if (cBeams[iBeam]) {
509 // Buffer offset for this beam.
510 cBeamSel[iBeam] = cNBeamSel++;
511 } else {
512 // Signal that the beam is not selected.
513 cBeamSel[iBeam] = -1;
514 }
515 }
516
517 // Set up bookkeeping arrays for IFs.
518 cIFSel = new int[cNIF];
519 cChanOff = new int[cNIF];
520 cXpolOff = new int[cNIF];
521
522 int simulIF = 0;
523 int maxChan = 0;
524 int maxXpol = 0;
525
526 for (int iIF = 0; iIF < cNIF; iIF++) {
527 if (cIFs[iIF]) {
528 // Buffer index for each IF within each simultaneous set.
529 cIFSel[iIF] = 0;
530
531 // Array offsets for each IF within each simultaneous set.
532 cChanOff[iIF] = 0;
533 cXpolOff[iIF] = 0;
534
535 // Look for earlier IFs in the same simultaneous set.
536 for (int jIF = 0; jIF < iIF; jIF++) {
537 if (!cIFs[jIF]) continue;
538
539 if (if_.if_simul[jIF] == if_.if_simul[iIF]) {
540 // Got one, increment indices.
541 cIFSel[iIF]++;
542
543 cChanOff[iIF] += cNChan[jIF] * cNPol[jIF];
544 if (cHaveXPol[jIF]) {
545 cXpolOff[iIF] += 2 * cNChan[jIF];
546 }
547 }
548 }
549
550 // Maximum number of selected IFs in any simultaneous set.
551 simulIF = max(simulIF, cIFSel[iIF]+1);
552
553 // Maximum memory required for any simultaneous set.
554 maxChan = max(maxChan, cChanOff[iIF] + cNChan[iIF]*cNPol[iIF]);
555 if (cHaveXPol[iIF]) {
556 maxXpol = max(maxXpol, cXpolOff[iIF] + 2*cNChan[iIF]);
557 }
558
559 } else {
560 // Signal that the IF is not selected.
561 cIFSel[iIF] = -1;
562 }
563 }
564
565 // Check for binning mode observations.
566 if (param_.intbase > 0.0f) {
567 cNBin = int((cIntTime / param_.intbase) + 0.5);
568
569 // intbase sometimes contains rubbish.
570 if (cNBin == 0) {
571 cNBin = 1;
572 }
573 } else {
574 cNBin = 1;
575 }
576
577 if (cNBin > 1 && cNBeamSel > 1) {
[1372]578 fprintf(stderr, "ERROR, cannot handle binning mode for multiple "
[1325]579 "beams.\n");
580 close();
581 return 1;
582 }
583
584 // Allocate buffer data storage; the PKSMBrecord constructor zeroes
585 // class members such as cycleNo that are tested in the first pass
586 // below.
587 int nBuff = cNBeamSel * cNBin;
588 cBuffer = new PKSMBrecord[nBuff];
589
590 // Allocate memory for spectral arrays.
591 for (int ibuff = 0; ibuff < nBuff; ibuff++) {
592 cBuffer[ibuff].setNIFs(simulIF);
593 cBuffer[ibuff].allocate(0, maxChan, maxXpol);
594 }
595
596 cPosUTC = new double[cNBeamSel];
597
598 cFirst = 0;
599 cScanNo = 1;
600 cCycleNo = 0;
601 cPrevUTC = 0.0;
602 cStaleness = new int[cNBeamSel];
603 for (int iBeamSel = 0; iBeamSel < cNBeamSel; iBeamSel++) {
604 cStaleness[iBeamSel] = 0;
605 }
606 }
607
608 // Check for end-of-scan.
609 if (cEOS) {
610 cScanNo++;
611 cCycleNo = 0;
612 cPrevUTC = 0.0;
613 }
614
615 // Apply beam selection.
616 beamNo = int(cBaseline / 256.0);
617 iBeamSel = cBeamSel[beamNo-1];
618 if (iBeamSel < 0) continue;
619
620 // Sanity check (mainly for MOPS).
621 if (cIFno > cNIF) continue;
622
623 // Apply IF selection.
624 iIFSel = cIFSel[cIFno - 1];
625 if (iIFSel < 0) continue;
626
627 sprintf(cDateObs, "%-10.10s", names_.datobs);
628
629 // Check for change-of-day.
630 if (cUTC < cPrevUTC - 85800.0) {
631 cUTC += 86400.0;
632 }
633
634 if (cNBin > 1) {
635 // Binning mode: correct the time.
636 cUTC += param_.intbase * (cBin - (cNBin + 1)/2.0);
637 }
638
639 // New integration cycle?
640 if (cUTC > cPrevUTC) {
641 cCycleNo++;
642 cPrevUTC = cUTC + 0.0001;
643 }
644
645 // Compute buffer number.
646 iMBuff = cBuffer + iBeamSel;
647 if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
648
649 if (cCycleNo < iMBuff->cycleNo) {
650 // Note that if the first beam and IF are not both selected cEOS
651 // will be cleared by rpget() when the next beam/IF is read.
652 cEOS = 1;
653 }
654
655 // Begin flush cycle?
656 if (cEOS || (iMBuff->nIF && cUTC > iMBuff->utc + 0.0001)) {
657 cFlushing = 1;
658 cFlushBin = 0;
659 cFlushIF = 0;
660 }
661
662#ifdef PKSIO_DEBUG
663 printf(" In:%4d%4d%3d%3d\n", cScanNo, cCycleNo, beamNo, cIFno);
664 if (cEOS) printf("Start of new scan, flushing previous scan.\n");
665#endif
666 }
667 }
668
669
670 if (cFlushing) {
671 // Find the oldest integration to flush, noting that the last
672 // integration cycle may be incomplete.
673 beamNo = 0;
674 int cycleNo = 0;
675 for (; cFlushBin < cNBin; cFlushBin++) {
676 for (iBeamSel = 0; iBeamSel < cNBeamSel; iBeamSel++) {
677 iMBuff = cBuffer + iBeamSel + cNBeamSel*cFlushBin;
678
679 // iMBuff->nIF is set to zero (below) to signal that all IFs in
680 // an integration have been flushed.
681 if (iMBuff->nIF) {
682 if (cycleNo == 0 || iMBuff->cycleNo < cycleNo) {
683 beamNo = iMBuff->beamNo;
684 cycleNo = iMBuff->cycleNo;
685 }
686 }
687 }
688
689 if (beamNo) {
690 // Found an integration to flush.
691 break;
692 }
693 }
694
695 if (beamNo) {
696 iBeamSel = cBeamSel[beamNo-1];
697 iMBuff = cBuffer + iBeamSel + cNBeamSel*cFlushBin;
698
699 // Find the IF to flush.
700 for (; cFlushIF < iMBuff->nIF; cFlushIF++) {
701 if (iMBuff->IFno[cFlushIF]) break;
702 }
703
704 } else {
705 // Flush complete.
706 cFlushing = 0;
707 if (cEOF) {
708 return -1;
709 }
710
711 // The last record read must have been the first of a new cycle.
712 beamNo = int(cBaseline / 256.0);
713 iBeamSel = cBeamSel[beamNo-1];
714
715 // Compute buffer number.
716 iMBuff = cBuffer + iBeamSel;
717 if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
718 }
719 }
720
721
722 if (cFlushing && cFlushBin == 0 && cFlushIF == 0 && cInterp) {
723 // Interpolate the beam position at the start of the flush cycle.
724#ifdef PKSIO_DEBUG
725 printf("Doing position interpolation for beam %d.\n", iMBuff->beamNo);
726#endif
727
728 double prevRA = iMBuff->ra;
729 double prevDec = iMBuff->dec;
730 double prevUTC = cPosUTC[iBeamSel];
731
732 if (!cEOF && !cEOS) {
733 // The position is measured by the control system at a time returned
734 // by RPFITSIN as the 'w' visibility coordinate. The ra and dec,
735 // returned as the 'u' and 'v' visibility coordinates, must be
736 // interpolated to the integration time which RPFITSIN returns as
737 // 'cUTC', this usually being a second or two later.
738 //
739 // Note that the time recorded as the 'w' visibility coordinate
740 // cycles through 86400 back to 0 at midnight, whereas that in 'cUTC'
741 // continues to increase past 86400.
742
743 double thisRA = cU;
744 double thisDec = cV;
745 double thisUTC = cW;
746
747 if (thisUTC < prevUTC) {
748 // Must have cycled through midnight.
749 thisUTC += 86400.0;
750 }
751
752 // Guard against RA cycling through 24h in either direction.
753 if (fabs(thisRA - prevRA) > PI) {
754 if (thisRA < prevRA) {
755 thisRA += TWOPI;
756 } else {
757 thisRA -= TWOPI;
758 }
759 }
760
761 // The control system at Mopra typically does not update the
762 // positions between successive integration cycles at the end of a
763 // scan (nor are they flagged). In this case we use the previously
764 // computed rates, even if from the previous scan since these are
765 // likely to be a better guess than anything else.
766
767 double dUTC = thisUTC - prevUTC;
768
769 // Scan rate for this beam.
770 if (dUTC > 0.0) {
771 iMBuff->raRate = (thisRA - prevRA) / dUTC;
772 iMBuff->decRate = (thisDec - prevDec) / dUTC;
773
774 if (cInterp == 2) {
775 // Use the same interpolation scheme as the original pksmbfits
776 // client. This incorrectly assumed that (thisUTC - prevUTC) is
777 // equal to the integration time and interpolated by computing a
778 // weighted sum of the positions before and after the required
779 // time.
780
781 double utc = iMBuff->utc;
782 if (utc - prevUTC > 100.0) {
783 // Must have cycled through midnight.
784 utc -= 86400.0;
785 }
786
787 double tw1 = 1.0 - (utc - prevUTC) / iMBuff->exposure;
788 double tw2 = 1.0 - (thisUTC - utc) / iMBuff->exposure;
789 double gamma = (tw2 / (tw1 + tw2)) * dUTC / (utc - prevUTC);
790
791 iMBuff->raRate *= gamma;
792 iMBuff->decRate *= gamma;
793 }
794
795 cStaleness[iBeamSel] = 0;
796
797 } else {
798 // Issue warnings.
799 int nch = 0;
800 fprintf(stderr, "WARNING, scan %d,%n cycle %d: Position ",
801 iMBuff->scanNo, &nch, iMBuff->cycleNo);
802
803 if (dUTC < 0.0) {
804 fprintf(stderr, "timestamp went backwards!\n");
805 } else {
806 if (thisRA != prevRA || thisDec != prevDec) {
807 fprintf(stderr, "changed but timestamp unchanged!\n");
808 } else {
809 fprintf(stderr, "and timestamp unchanged!\n");
810 }
811 }
812
813 cStaleness[iBeamSel]++;
814 fprintf(stderr, "%-*s Using stale scan rate, staleness = %d "
815 "cycle%s.\n", nch, "WARNING,", cStaleness[iBeamSel],
816 (cStaleness[iBeamSel] == 1) ? "" : "s");
817
818 if (thisRA != prevRA || thisDec != prevDec) {
819 if (iMBuff->raRate == 0.0 && iMBuff->decRate == 0.0) {
820 fprintf(stderr, "%-*s But the previous rate was zero! "
821 "Position will be inaccurate.\n", nch, "WARNING,");
822 }
823 }
824 }
825 }
826
827 // Compute the position of this beam for all bins.
828 for (int idx = 0; idx < cNBin; idx++) {
829 int jbuff = iBeamSel + cNBeamSel*idx;
830
831 cBuffer[jbuff].raRate = iMBuff->raRate;
832 cBuffer[jbuff].decRate = iMBuff->decRate;
833
834 double dutc = cBuffer[jbuff].utc - prevUTC;
835 if (dutc > 100.0) {
836 // Must have cycled through midnight.
837 dutc -= 86400.0;
838 }
839
840 cBuffer[jbuff].ra = prevRA + cBuffer[jbuff].raRate * dutc;
841 cBuffer[jbuff].dec = prevDec + cBuffer[jbuff].decRate * dutc;
842 if (cBuffer[jbuff].ra < 0.0) {
843 cBuffer[jbuff].ra += TWOPI;
844 } else if (cBuffer[jbuff].ra > TWOPI) {
845 cBuffer[jbuff].ra -= TWOPI;
846 }
847 }
848 }
849
850
851 if (cFlushing) {
852 // Copy buffer location out one IF at a time.
853 MBrec.extract(*iMBuff, cFlushIF);
854 haveData = 1;
855
856#ifdef PKSIO_DEBUG
857 printf("Out:%4d%4d%3d%3d\n", MBrec.scanNo, MBrec.cycleNo, MBrec.beamNo,
858 MBrec.IFno[0]);
859#endif
860
861 // Signal that this IF in this buffer location has been flushed.
862 iMBuff->IFno[cFlushIF] = 0;
863
864 if (cFlushIF == iMBuff->nIF - 1) {
865 // Signal that all IFs in this buffer location have been flushed.
866 iMBuff->nIF = 0;
867
868 // Stop cEOS being set when the next integration is read.
869 iMBuff->cycleNo = 0;
870
871 } else {
872 // Carry on flushing the other IFs.
873 continue;
874 }
875
876 // Has the whole buffer been flushed?
877 if (cFlushBin == cNBin - 1) {
878 if (cEOS || cEOF) {
879 // Carry on flushing other buffers.
880 cFlushIF = 0;
881 continue;
882 }
883
884 cFlushing = 0;
885
886 beamNo = int(cBaseline / 256.0);
887 iBeamSel = cBeamSel[beamNo-1];
888
889 // Compute buffer number.
890 iMBuff = cBuffer + iBeamSel;
891 if (cNBin > 1) iMBuff += cNBeamSel*(cBin-1);
892 }
893 }
894
895 if (!cFlushing) {
896 // Buffer this MBrec.
897 if (cCycleNo == 1 && iMBuff->IFno[0]) {
898 // Sanity check on the number of IFs in the new scan.
899 if (if_.n_if != cNIF) {
900 fprintf(stderr, "WARNING, scan %d has %d IFs instead of %d, "
901 "continuing.\n", cScanNo, if_.n_if, cNIF);
902 }
903 }
904
[1372]905 // Sanity check on incomplete integrations within a scan.
906 if (iMBuff->nIF && (iMBuff->cycleNo != cCycleNo)) {
907 // Force the incomplete integration to be flushed before proceeding.
908 cFlushing = 1;
909 continue;
910 }
911
[1325]912 iMBuff->scanNo = cScanNo;
913 iMBuff->cycleNo = cCycleNo;
914
915 // Times.
916 strncpy(iMBuff->datobs, cDateObs, 10);
917 iMBuff->utc = cUTC;
918 iMBuff->exposure = param_.intbase;
919
920 // Source identification.
921 sprintf(iMBuff->srcName, "%-16.16s",
922 names_.su_name + (cSrcNo-1)*16);
923 iMBuff->srcRA = doubles_.su_ra[cSrcNo-1];
924 iMBuff->srcDec = doubles_.su_dec[cSrcNo-1];
925
926 // Rest frequency of the line of interest.
927 iMBuff->restFreq = doubles_.rfreq;
928 if (strncmp(names_.instrument, "ATPKSMB", 7) == 0) {
929 // Fix the HI rest frequency recorded for Parkes multibeam data.
930 double reffreq = doubles_.freq;
931 double restfreq = doubles_.rfreq;
932 if ((restfreq == 0.0 || fabs(restfreq - reffreq) == 0.0) &&
933 fabs(reffreq - 1420.40575e6) < 100.0) {
934 iMBuff->restFreq = 1420.40575e6;
935 }
936 }
937
938 // Observation type.
939 int j;
940 for (j = 0; j < 15; j++) {
941 iMBuff->obsType[j] = names_.card[11+j];
942 if (iMBuff->obsType[j] == '\'') break;
943 }
944 iMBuff->obsType[j] = '\0';
945
946 // Beam-dependent parameters.
947 iMBuff->beamNo = beamNo;
948
949 // Beam position at the specified time.
[1372]950 if (cSUpos) {
951 // Non-ATNF data that does not store the position in (u,v,w).
[1325]952 iMBuff->ra = doubles_.su_ra[cSrcNo-1];
953 iMBuff->dec = doubles_.su_dec[cSrcNo-1];
954 } else {
955 iMBuff->ra = cU;
956 iMBuff->dec = cV;
957 }
958 cPosUTC[iBeamSel] = cW;
959
960 // IF-dependent parameters.
961 int iIF = cIFno - 1;
962 int startChan = cStartChan[iIF];
963 int endChan = cEndChan[iIF];
964 int refChan = cRefChan[iIF];
965
966 int nChan = abs(endChan - startChan) + 1;
967
968 iIFSel = cIFSel[iIF];
969 iMBuff->nIF++;
970 iMBuff->IFno[iIFSel] = cIFno;
971 iMBuff->nChan[iIFSel] = nChan;
972 iMBuff->nPol[iIFSel] = cNPol[iIF];
973
974 iMBuff->fqRefPix[iIFSel] = doubles_.if_ref[iIF];
975 iMBuff->fqRefVal[iIFSel] = doubles_.if_freq[iIF];
976 iMBuff->fqDelt[iIFSel] =
977 if_.if_invert[iIF] * fabs(doubles_.if_bw[iIF] /
978 (if_.if_nfreq[iIF] - 1));
979
980 // Adjust for channel selection.
981 if (iMBuff->fqRefPix[iIFSel] != refChan) {
982 iMBuff->fqRefVal[iIFSel] +=
983 (refChan - iMBuff->fqRefPix[iIFSel]) *
984 iMBuff->fqDelt[iIFSel];
985 iMBuff->fqRefPix[iIFSel] = refChan;
986 }
987
988 if (endChan < startChan) {
989 iMBuff->fqDelt[iIFSel] = -iMBuff->fqDelt[iIFSel];
990 }
991
992
993 // System temperature.
994 int iBeam = beamNo - 1;
995 int scq = sc_.sc_q;
996 float TsysPol1 = sc_.sc_cal[scq*iBeam + 3];
997 float TsysPol2 = sc_.sc_cal[scq*iBeam + 4];
998 iMBuff->tsys[iIFSel][0] = TsysPol1*TsysPol1;
999 iMBuff->tsys[iIFSel][1] = TsysPol2*TsysPol2;
1000
1001 // Calibration factor; may be changed later if the data is recalibrated.
1002 if (scq > 14) {
1003 // Will only be present for Parkes Multibeam or LBA data.
1004 iMBuff->calfctr[iIFSel][0] = sc_.sc_cal[scq*iBeam + 14];
1005 iMBuff->calfctr[iIFSel][1] = sc_.sc_cal[scq*iBeam + 15];
1006 } else {
1007 iMBuff->calfctr[iIFSel][0] = 0.0f;
1008 iMBuff->calfctr[iIFSel][1] = 0.0f;
1009 }
1010
1011 // Cross-polarization calibration factor (unknown to MBFITS).
1012 for (int j = 0; j < 2; j++) {
1013 iMBuff->xcalfctr[iIFSel][j] = 0.0f;
1014 }
1015
1016 // Baseline parameters (unknown to MBFITS).
1017 iMBuff->haveBase = 0;
1018
1019 // Data (always present in MBFITS).
1020 iMBuff->haveSpectra = 1;
1021
1022 // Flag: bit 0 set if off source.
1023 // bit 1 set if loss of sync in A polarization.
1024 // bit 2 set if loss of sync in B polarization.
1025 unsigned char rpflag =
1026 (unsigned char)(sc_.sc_cal[scq*iBeam + 12] + 0.5f);
1027
1028 // The baseline flag may be set independently.
1029 if (rpflag == 0) rpflag = cFlag;
1030
1031 // Copy and scale data.
1032 int inc = 2 * if_.if_nstok[iIF];
1033 if (endChan < startChan) inc = -inc;
1034
1035 float TsysF;
1036 iMBuff->spectra[iIFSel] = iMBuff->spectra[0] + cChanOff[iIF];
1037 iMBuff->flagged[iIFSel] = iMBuff->flagged[0] + cChanOff[iIF];
1038
1039 float *spectra = iMBuff->spectra[iIFSel];
1040 unsigned char *flagged = iMBuff->flagged[iIFSel];
1041 for (int ipol = 0; ipol < cNPol[iIF]; ipol++) {
1042 if (sc_.sc_cal[scq*iBeam + 3 + ipol] > 0.0f) {
1043 // The correlator has already applied the calibration.
1044 TsysF = 1.0f;
1045 } else {
1046 // The correlator has normalized cVis[k] to a Tsys of 500K.
1047 TsysF = iMBuff->tsys[iIFSel][ipol] / 500.0f;
1048 }
1049
1050 int k = 2 * (if_.if_nstok[iIF]*(startChan - 1) + ipol);
1051 for (int ichan = 0; ichan < nChan; ichan++) {
1052 *(spectra++) = TsysF * cVis[k];
1053 *(flagged++) = rpflag;
1054 k += inc;
1055 }
1056 }
1057
1058 if (cHaveXPol[iIF]) {
1059 int k = 2 * (3*(startChan - 1) + 2);
1060 iMBuff->xpol[iIFSel] = iMBuff->xpol[0] + cXpolOff[iIF];
1061 float *xpol = iMBuff->xpol[iIFSel];
1062 for (int ichan = 0; ichan < nChan; ichan++) {
1063 *(xpol++) = cVis[k];
1064 *(xpol++) = cVis[k+1];
1065 k += inc;
1066 }
1067 }
1068
1069
1070 // Parallactic angle.
1071 iMBuff->parAngle = sc_.sc_cal[scq*iBeam + 11];
1072
1073 // Calibration factor applied to the data by the correlator.
1074 if (scq > 14) {
1075 // Will only be present for Parkes Multibeam or LBA data.
1076 iMBuff->tcal[iIFSel][0] = sc_.sc_cal[scq*iBeam + 14];
1077 iMBuff->tcal[iIFSel][1] = sc_.sc_cal[scq*iBeam + 15];
1078 } else {
1079 iMBuff->tcal[iIFSel][0] = 0.0f;
1080 iMBuff->tcal[iIFSel][1] = 0.0f;
1081 }
1082
1083 if (sc_.sc_ant <= anten_.nant) {
1084 // No extra syscal information present.
1085 iMBuff->extraSysCal = 0;
1086 iMBuff->azimuth = 0.0f;
1087 iMBuff->elevation = 0.0f;
1088 iMBuff->parAngle = 0.0f;
1089 iMBuff->focusAxi = 0.0f;
1090 iMBuff->focusTan = 0.0f;
1091 iMBuff->focusRot = 0.0f;
1092 iMBuff->temp = 0.0f;
1093 iMBuff->pressure = 0.0f;
1094 iMBuff->humidity = 0.0f;
1095 iMBuff->windSpeed = 0.0f;
1096 iMBuff->windAz = 0.0f;
1097 strcpy(iMBuff->tcalTime, " ");
1098 iMBuff->refBeam = 0;
1099
1100 } else {
1101 // Additional information for Parkes Multibeam data.
1102 int iOff = scq*(sc_.sc_ant - 1) - 1;
1103 iMBuff->extraSysCal = 1;
1104 iMBuff->azimuth = sc_.sc_cal[iOff + 2];
1105 iMBuff->elevation = sc_.sc_cal[iOff + 3];
1106 iMBuff->parAngle = sc_.sc_cal[iOff + 4];
1107 iMBuff->focusAxi = sc_.sc_cal[iOff + 5] * 1e-3;
1108 iMBuff->focusTan = sc_.sc_cal[iOff + 6] * 1e-3;
1109 iMBuff->focusRot = sc_.sc_cal[iOff + 7];
1110 iMBuff->temp = sc_.sc_cal[iOff + 8];
1111 iMBuff->pressure = sc_.sc_cal[iOff + 9];
1112 iMBuff->humidity = sc_.sc_cal[iOff + 10];
1113 iMBuff->windSpeed = sc_.sc_cal[iOff + 11];
1114 iMBuff->windAz = sc_.sc_cal[iOff + 12];
1115
1116 char *tcalTime = iMBuff->tcalTime;
1117 sprintf(tcalTime, "%-16.16s", (char *)(&sc_.sc_cal[iOff+13]));
1118
1119#ifndef AIPS_LITTLE_ENDIAN
1120 // Do byte swapping on the ASCII date string.
1121 for (int j = 0; j < 16; j += 4) {
1122 char ctmp;
1123 ctmp = tcalTime[j];
1124 tcalTime[j] = tcalTime[j+3];
1125 tcalTime[j+3] = ctmp;
1126 ctmp = tcalTime[j+1];
1127 tcalTime[j+1] = tcalTime[j+2];
1128 tcalTime[j+2] = ctmp;
1129 }
1130#endif
1131
1132 // Reference beam number.
1133 float refbeam = sc_.sc_cal[iOff + 17];
1134 if (refbeam > 0.0f || refbeam < 100.0f) {
1135 iMBuff->refBeam = int(refbeam);
1136 } else {
1137 iMBuff->refBeam = 0;
1138 }
1139 }
1140 }
1141 }
1142
1143 return 0;
1144}
1145
1146//-------------------------------------------------------- MBFITSreader::rpget
1147
1148// Read the next data record from the RPFITS file.
1149
1150int MBFITSreader::rpget(int syscalonly, int &EOS)
1151{
1152 EOS = 0;
1153
1154 int retries = 0;
1155
1156 // Allow 10 read errors.
1157 int numErr = 0;
1158
1159 int jstat = 0;
1160 while (numErr < 10) {
1161 int lastjstat = jstat;
1162 rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
1163 &cBin, &cIFno, &cSrcNo);
1164
1165 switch(jstat) {
1166 case -1:
1167 // Read failed; retry.
1168 numErr++;
1169 fprintf(stderr, "RPFITS read failed - retrying.\n");
1170 jstat = 0;
1171 break;
1172
1173 case 0:
1174 // Successful read.
1175 if (lastjstat == 0) {
1176 if (cBaseline == -1) {
1177 // Syscal data.
1178 if (syscalonly) {
1179 return 0;
1180 }
1181
1182 } else {
1183 if (!syscalonly) {
1184 return 0;
1185 }
1186 }
1187 }
1188
1189 // Last operation was to read header or FG table; now read data.
1190 break;
1191
1192 case 1:
1193 // Encountered header while trying to read data; read it.
1194 EOS = 1;
1195 jstat = -1;
1196 break;
1197
1198 case 2:
1199 // End of scan; read past it.
1200 jstat = 0;
1201 break;
1202
1203 case 3:
1204 // End-of-file; retry applies to real-time mode.
1205 if (retries++ >= cRetry) {
1206 return -1;
1207 }
1208
1209 sleep(10);
1210 jstat = 0;
1211 break;
1212
1213 case 4:
1214 // Encountered FG table while trying to read data; read it.
1215 jstat = -1;
1216 break;
1217
1218 case 5:
1219 // Illegal data at end of block after close/reopen operation; retry.
1220 jstat = 0;
1221 break;
1222
1223 default:
1224 // Shouldn't reach here.
1225 fprintf(stderr, "Unrecognized RPFITSIN return code: %d (retrying)\n",
1226 jstat);
1227 jstat = 0;
1228 break;
1229 }
1230 }
1231
[1372]1232 fprintf(stderr, "ERROR, RPFITS read failed too many times.\n");
[1325]1233 return 2;
1234}
1235
1236//-------------------------------------------------------- MBFITSreader::close
1237
1238// Close the input file.
1239
1240void MBFITSreader::close(void)
1241{
1242 if (cMBopen) {
1243 int jstat = 1;
1244 rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag,
1245 &cBin, &cIFno, &cSrcNo);
1246
1247 if (cBeams) delete [] cBeams;
1248 if (cIFs) delete [] cIFs;
1249 if (cNChan) delete [] cNChan;
1250 if (cNPol) delete [] cNPol;
1251 if (cHaveXPol) delete [] cHaveXPol;
1252 if (cStartChan) delete [] cStartChan;
1253 if (cEndChan) delete [] cEndChan;
1254 if (cRefChan) delete [] cRefChan;
1255
1256 if (cVis) delete [] cVis;
1257 if (cWgt) delete [] cWgt;
1258
1259 if (cBeamSel) delete [] cBeamSel;
1260 if (cIFSel) delete [] cIFSel;
1261 if (cChanOff) delete [] cChanOff;
1262 if (cXpolOff) delete [] cXpolOff;
1263 if (cBuffer) delete [] cBuffer;
1264 if (cPosUTC) delete [] cPosUTC;
1265
1266 cMBopen = 0;
1267 }
1268}
Note: See TracBrowser for help on using the repository browser.