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 | //# |
---|
29 | //# $Id: MBFITSreader.cc,v 19.34 2007/07/02 06:12:18 cal103 Exp $ |
---|
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 = 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 | |
---|
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 | int jstat = -3; |
---|
130 | rpfitsin_(&jstat, cVis, cWgt, &cBaseline, &cUTC, &cU, &cV, &cW, &cFlag, |
---|
131 | &cBin, &cIFno, &cSrcNo); |
---|
132 | |
---|
133 | if (jstat) { |
---|
134 | fprintf(stderr, "ERROR, failed to open MBFITS file: %s\n", rpname); |
---|
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) { |
---|
154 | fprintf(stderr, "ERROR, failed to read MBFITS header: %s\n", rpname); |
---|
155 | close(); |
---|
156 | return 1; |
---|
157 | } |
---|
158 | |
---|
159 | // Mopra data has some peculiarities. |
---|
160 | cMopra = strncmp(names_.instrument, "ATMOPRA", 7) == 0; |
---|
161 | |
---|
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"); |
---|
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) { |
---|
191 | fprintf(stderr, "ERROR, couldn't determine number of beams.\n"); |
---|
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)) { |
---|
293 | fprintf(stderr, "ERROR, failed to read first syscal record.\n"); |
---|
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 | |
---|
313 | int 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) { |
---|
328 | fprintf(stderr, "ERROR, an MBFITS file has not been opened.\n"); |
---|
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 | |
---|
414 | int 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 | |
---|
427 | int 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 | |
---|
442 | int MBFITSreader::read( |
---|
443 | PKSMBrecord &MBrec) |
---|
444 | { |
---|
445 | int beamNo = -1; |
---|
446 | int haveData, status; |
---|
447 | PKSMBrecord *iMBuff = 0x0; |
---|
448 | |
---|
449 | if (!cMBopen) { |
---|
450 | fprintf(stderr, "ERROR, an MBFITS file has not been opened.\n"); |
---|
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) { |
---|
578 | fprintf(stderr, "ERROR, cannot handle binning mode for multiple " |
---|
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 | |
---|
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 | |
---|
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. |
---|
950 | if (cSUpos) { |
---|
951 | // Non-ATNF data that does not store the position in (u,v,w). |
---|
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 | |
---|
1150 | int 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 | |
---|
1232 | fprintf(stderr, "ERROR, RPFITS read failed too many times.\n"); |
---|
1233 | return 2; |
---|
1234 | } |
---|
1235 | |
---|
1236 | //-------------------------------------------------------- MBFITSreader::close |
---|
1237 | |
---|
1238 | // Close the input file. |
---|
1239 | |
---|
1240 | void 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 | } |
---|