source: branches/alma/external/atnf/PKSIO/MBFITSreader.cc@ 1488

Last change on this file since 1488 was 1453, checked in by TakTsutsumi, 16 years ago

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: many

Test Programs: sd.scantable(), sd.scantable.save()

Put in Release Notes: N/A

Description: copied from current casapy code tree


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