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 | |
---|
49 | using namespace std; |
---|
50 | |
---|
51 | // Numerical constants. |
---|
52 | const double PI = 3.141592653589793238462643; |
---|
53 | const double TWOPI = 2.0 * PI; |
---|
54 | |
---|
55 | //------------------------------------------------- MBFITSreader::MBFITSreader |
---|
56 | |
---|
57 | // Default constructor. |
---|
58 | |
---|
59 | MBFITSreader::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 | |
---|
100 | MBFITSreader::~MBFITSreader() |
---|
101 | { |
---|
102 | close(); |
---|
103 | } |
---|
104 | |
---|
105 | //--------------------------------------------------------- MBFITSreader::open |
---|
106 | |
---|
107 | // Open the RPFITS file for reading. |
---|
108 | |
---|
109 | int 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 | |
---|
300 | int 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 | |
---|
401 | int 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 | |
---|
414 | int 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 | |
---|
429 | int 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 | |
---|
1128 | int 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 | |
---|
1218 | void 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 | } |
---|