Line data Source code
1 : C-----------------------------------------------------------------------
2 : C SUBROUTINE RPFITSIN
3 : C-----------------------------------------------------------------------
4 : C
5 : C For information on the use of this software, and on the RPFITS
6 : C format, see the file RPFITS.DEFN.
7 : C
8 : C Programmer: Ray Norris
9 : C Date: 25 April 1985
10 : C
11 : C $Id: rpfitsin.f,v 1.48 2018/10/26 03:03:37 wie017 Exp $
12 : C-----------------------------------------------------------------------
13 :
14 22717682 : subroutine RPFITSIN (jstat, vis, weight, baseline, ut, u, v, w,
15 : : flag, bin, if_no, sourceno)
16 :
17 : integer baseline, flag, bin, if_no, sourceno
18 : real weight(*), ut, u, v, w
19 : complex vis(*)
20 :
21 :
22 : include 'rpfits.inc'
23 :
24 : logical async, endhdr, endscan, isopen, new_antenna, open_only,
25 : : starthdr
26 : integer AT_CLOSE, AT_OPEN_READ, AT_READ, AT_SKIP_EOF, AT_UNREAD,
27 : : bufleft, bufleft3, bufptr, grplength, grpptr, i, i1, i2,
28 : : i3, i_buff(640), i_grphdr(11), icard, ierr, illegal, j,
29 : : jstat, k, lun, pcount, SIMPLE
30 : real buffer(640), crpix4, grphdr(11), r1, r2, revis,
31 : : sc_buf(max_sc*max_if*ant_max), pra, pdec
32 : static :: buffer
33 : double precision d2pi
34 : character keyvalue*20, keyword*8, m(32)*80, terr*2
35 :
36 : equivalence (i_buff(1), buffer(1))
37 : equivalence (i_grphdr(1), grphdr(1))
38 : equivalence (sc_buf(1), sc_cal(1,1,1))
39 :
40 : parameter (d2pi = 2d0 * 3.14159265358979323846d0)
41 :
42 : data isopen /.false./
43 : data async /.false./
44 : data new_antenna /.false./
45 : data illegal /32768/
46 :
47 : save
48 :
49 : C-------------------------- DECIDE ON ACTION ---------------------------
50 :
51 1505 : rp_iostat = 0
52 1505 : errmsg = ''
53 :
54 1505 : open_only = jstat.eq.-3
55 :
56 1505 : if (jstat.eq.-3) go to 1000
57 1505 : if (jstat.eq.-2) go to 1000
58 1503 : if (jstat.eq.-1) go to 2000
59 1503 : if (jstat.eq.0) go to 3000
60 2 : if (jstat.eq.1) go to 5000
61 0 : if (jstat.eq.2) go to 6000
62 :
63 0 : write (errmsg, '(a,i3)') 'Illegal value of jstat =', jstat
64 0 : call RPFERR (errmsg)
65 0 : jstat = -1
66 2 : RETURN
67 :
68 : C--------------------------- OPEN FITS FILE ----------------------------
69 :
70 2 : 1000 if (isopen) then
71 0 : call RPFERR ('File is already open.')
72 0 : jstat = -1
73 0 : RETURN
74 : end if
75 :
76 2 : rp_iostat = AT_OPEN_READ (file, async, lun)
77 2 : if (rp_iostat.ne.0) then
78 0 : call RPFERR ('File open error')
79 0 : jstat = -1
80 0 : RETURN
81 : end if
82 2 : isopen = .true.
83 :
84 2 : if (open_only) then
85 0 : jstat = 0
86 0 : RETURN
87 : end if
88 :
89 : C----------------------------- READ HEADER -----------------------------
90 :
91 2 : 2000 if (.not.isopen) then
92 0 : call RPFERR ('File is not open.')
93 0 : jstat = -1
94 0 : RETURN
95 : end if
96 :
97 2 : bufptr = 0
98 2 : n_if = 0
99 2 : icard = 1
100 2 : if (ncard.lt.0) ncard = -1
101 2 : an_found = .false.
102 2 : if_found = .false.
103 2 : su_found = .false.
104 2 : fg_found = .false.
105 2 : nx_found = .false.
106 2 : mt_found = .false.
107 2 : cu_found = .false.
108 2 : pra = 0.0
109 2 : pdec = 0.0
110 :
111 : C Look for start of next header.
112 2 : starthdr = .false.
113 4 : do while (.not.starthdr)
114 2 : rp_iostat = AT_READ (lun, buffer)
115 2 : if (rp_iostat.ne.0) then
116 0 : if (rp_iostat.eq.-1) then
117 0 : jstat = 3
118 0 : RETURN
119 : end if
120 :
121 0 : call RPFERR ('I/O error reading header')
122 0 : jstat = -1
123 0 : RETURN
124 : end if
125 :
126 2 : jstat = SIMPLE (buffer, lun)
127 2 : if (jstat.eq.1) then
128 : C Start of header.
129 2 : starthdr = .true.
130 0 : else if (jstat.eq.3) then
131 : C End-of-file while reading flag table.
132 0 : RETURN
133 0 : else if (jstat.eq.4) then
134 : C Encountered flag table.
135 0 : RETURN
136 0 : else if (jstat.ne.0) then
137 : C Fortran I/O error status.
138 0 : jstat = -1
139 0 : RETURN
140 : end if
141 :
142 1284 : write (m, '(32(20a4,:,/))') (buffer(j),j=1,640)
143 : end do
144 :
145 : C Scan through header, getting the interesting bits.
146 2 : endhdr = .false.
147 10 : do 2500 while (.not.endhdr)
148 8 : if (.not.starthdr) then
149 6 : rp_iostat = AT_READ (lun, buffer)
150 3846 : write (m,'(32(20a4,:,/))') (buffer(j),j=1,640)
151 6 : if (rp_iostat.ne.0) then
152 0 : if (rp_iostat.eq.-1) then
153 0 : jstat = 3
154 0 : RETURN
155 : end if
156 :
157 0 : call RPFERR ('I/O error reading header')
158 0 : jstat = -1
159 0 : RETURN
160 : end if
161 : end if
162 :
163 8 : starthdr = .false.
164 8 : version = ' '
165 214 : do 2400 i = 1, 32
166 : C Parse the RPFITS keyword and keyvalue.
167 208 : keyword = m(i)(1:8)
168 :
169 208 : if (m(i)(11:11).eq.'''') then
170 : C Must be a character value.
171 56 : keyvalue = m(i)(12:31)
172 1176 : do j = 1, 20
173 1176 : if (keyvalue(j:j).eq.'''') then
174 : C Strip off the trailing apostrophe.
175 50 : keyvalue(j:) = ' '
176 : end if
177 : end do
178 : else
179 152 : keyvalue = m(i)(11:30)
180 : end if
181 :
182 : C Lexical chop based on the first letter of the keyword name.
183 208 : if (keyword(:1).le.'C') then
184 : C Keyword names beginning with A to C.
185 70 : if (keyword.eq.'ALTRVAL ') then
186 2 : read (keyvalue, '(g20.12)') vel1
187 68 : else if (keyword.eq.'BUNIT') then
188 2 : bunit = keyvalue(:16)
189 66 : else if (keyword.eq.'CAL') then
190 2 : cal = keyvalue(:16)
191 64 : else if (keyword.eq.'CDELT4') then
192 2 : read (keyvalue, '(g20.12)') dfreq
193 62 : else if (keyword.eq.'CRPIX4') then
194 2 : read (keyvalue, '(g20.12)') crpix4
195 60 : else if (keyword.eq.'CRVAL4') then
196 2 : read (keyvalue, '(g20.12)') freq
197 58 : else if (keyword.eq.'CRVAL5') then
198 2 : read (keyvalue, '(g20.12)') ra
199 2 : if (ra.lt.0d0) ra = ra + d2pi
200 56 : else if (keyword.eq.'CRVAL6') then
201 2 : read (keyvalue, '(g20.12)') dec
202 : end if
203 :
204 138 : else if (keyword(:1).le.'E') then
205 : C Keyword names beginning with D or E.
206 36 : if (keyword.eq.'DATE') then
207 : C Fix old-format dates.
208 2 : call datfit(keyvalue(:10), datwrit, ierr)
209 34 : else if (keyword.eq.'DATE-OBS') then
210 : C Fix old-format dates.
211 2 : call datfit(keyvalue(:10), datobs, ierr)
212 2 : datsys = m(i)(35:36)
213 2 : if (datsys.eq.'UT D') datsys = 'UT'
214 32 : else if (keyword.eq.'DEFEAT ') then
215 2 : read (keyvalue, '(i20)') rp_defeat
216 30 : else if (keyword.eq.'DJMREFP ') then
217 2 : read (keyvalue, '(g20.12)') rp_djmrefp
218 28 : else if (keyword.eq.'DJMREFT ') then
219 2 : read (keyvalue, '(g20.12)') rp_djmreft
220 26 : else if (keyword.eq.'END') then
221 : C END card.
222 0 : endhdr = .true.
223 26 : else if (keyword(1:5).eq.'EPHEM') then
224 24 : read (keyword(6:7), '(i2)') k
225 24 : read (keyvalue, '(g20.12)') rp_c(k)
226 2 : else if (keyword.eq.'EPOCH') then
227 2 : coord = keyvalue(:16)
228 : end if
229 :
230 102 : else if (keyword(:1).le.'N') then
231 : C Keyword names beginning with F to N.
232 30 : if (keyword.eq.'GCOUNT') then
233 2 : read (keyvalue, '(i20)') ncount
234 28 : else if (keyword(1:5).eq.'HUMID') then
235 0 : read (keyword(6:7), '(i2)') k
236 0 : read (keyvalue, '(g20.12)') rp_humid(k)
237 28 : else if (keyword.eq.'INSTRUME') then
238 2 : instrument = keyvalue(:16)
239 26 : else if (keyword.eq.'INTIME') then
240 2 : read (keyvalue, '(i20)') intime
241 24 : else if (keyword.eq.'NAXIS2') then
242 2 : read (keyvalue, '(i20)') data_format
243 2 : write_wt = data_format.eq.3
244 22 : else if (keyword.eq.'NAXIS3') then
245 2 : read (keyvalue, '(i20)') nstok
246 20 : else if (keyword.eq.'NAXIS4') then
247 2 : read (keyvalue, '(i20)') nfreq
248 18 : else if (keyword.eq.'NAXIS7') then
249 : C Note fudge for intermediate format PTI data.
250 0 : read (keyvalue, '(i20)') nstok
251 : end if
252 :
253 72 : else if (keyword(:1).le.'P') then
254 : C Keyword names beginning with M to P.
255 40 : if (keyword.eq.'OBJECT') then
256 2 : object = keyvalue(:16)
257 38 : else if (keyword.eq.'OBSERVER') then
258 2 : rp_observer = keyvalue(:16)
259 36 : else if (keyword.eq.'OBSTYPE') then
260 2 : obstype = keyvalue(:16)
261 34 : else if (keyword.eq.'PCOUNT') then
262 2 : read (keyvalue, '(i20)') pcount
263 32 : else if (keyword.eq.'PMDEC') then
264 2 : read (keyvalue, '(g20.12)') pm_dec
265 30 : else if (keyword.eq.'PMEPOCH') then
266 2 : read (keyvalue, '(g20.12)') pm_epoch
267 28 : else if (keyword.eq.'PMRA') then
268 2 : read (keyvalue, '(g20.12)') pm_ra
269 26 : else if (keyword.eq.'PNTCENTR') then
270 0 : read (m(i)(11:35),'(g12.9,1x,g12.9)') pra,pdec
271 26 : else if (keyword(1:5).eq.'PRESS') then
272 0 : read (keyword(6:7), '(i2)') k
273 0 : read (keyvalue, '(g20.12)') rp_pressure(k)
274 : end if
275 :
276 : else
277 : C Keyword names beginning with Q to Z.
278 32 : if (keyword.eq.'RESTFREQ') then
279 2 : read (keyvalue, '(g20.12)') rfreq
280 30 : else if (keyword.eq.'RPFITS ') then
281 2 : rpfitsversion = keyvalue
282 28 : else if (keyword.eq.'SCANS ') then
283 2 : read (keyvalue, '(i20)') nscan
284 26 : else if (keyword(1:6).eq.'TABLE ') then
285 : C Sort out tables.
286 2 : call RPFITS_READ_TABLE (lun, m, i, endhdr, terr, ierr)
287 2 : if (ierr.ne.0) then
288 0 : if (ierr.eq.1) then
289 0 : jstat = -1
290 : call RPFERR (terr // ' table contains too ' //
291 0 : : 'many entries.')
292 0 : else if (rp_iostat.lt.0) then
293 0 : jstat = 3
294 : else
295 0 : jstat = -1
296 : call RPFERR ('I/O error reading ' // terr //
297 0 : : ' table')
298 : end if
299 :
300 0 : RETURN
301 : end if
302 :
303 24 : else if (keyword(1:5).eq.'TEMPE') then
304 0 : read (keyword(6:7), '(i2)') k
305 0 : read (keyvalue, '(g20.12)') rp_temp(k)
306 24 : else if (keyword.eq.'UTCMTAI ') then
307 2 : read (keyvalue, '(g20.12)') rp_utcmtai
308 22 : else if (keyword.eq.'VELREF ') then
309 2 : read (keyvalue, '(i20)') ivelref
310 20 : else if (keyword.eq.'VERSION ') then
311 2 : version = keyvalue
312 : end if
313 : end if
314 :
315 : C Write into "cards" array if necessary.
316 208 : if (ncard.gt.0) then
317 0 : do j = 1, ncard
318 0 : if (keyword.eq.card(j)(1:8)) then
319 0 : card(j) = m(i)
320 : end if
321 : end do
322 208 : else if (ncard.lt.0) then
323 208 : if (icard.le.max_card .and. .not.endhdr) then
324 412 : card(-ncard) = m(i)
325 206 : icard = icard + 1
326 206 : ncard = ncard - 1
327 : end if
328 : end if
329 :
330 : C Antenna parameters.
331 208 : if (keyword(:7).eq.'ANTENNA') then
332 0 : if (.not.new_antenna) then
333 0 : nant = 0
334 0 : new_antenna = .true.
335 : end if
336 :
337 0 : if (keyword.eq.'ANTENNA') then
338 0 : read (m(i)(11:80), 2200) k, sta(k), x(k), y(k), z(k)
339 : 2200 format (i1,1x,a3,3x,g17.10,3x,g17.10,3x,g17.10)
340 : else
341 : C Old format ('ANTENNA:').
342 0 : read (m(i)(12:71), 2300) k, x(k), y(k), z(k), sta(k)
343 : 2300 format (i1,4x,g13.6,3x,g13.6,3x,g13.6,5x,a3)
344 : end if
345 :
346 0 : nant = nant + 1
347 : end if
348 :
349 208 : if (ENDHDR) go to 2500
350 8 : 2400 continue
351 2 : 2500 continue
352 2 : ncard = ABS(ncard)
353 :
354 : C Set up for reading data.
355 2 : if (data_format.lt.1 .or. data_format.gt.3) then
356 0 : call RPFERR ('NAXIS2 must be 1, 2, or 3.')
357 0 : jstat = -1
358 0 : RETURN
359 : end if
360 :
361 : C Insert default values into table commons if tables weren't found.
362 2 : if (.not.if_found) then
363 0 : n_if = 1
364 0 : if_freq(1) = freq
365 0 : if_invert(1) = 1
366 0 : if_bw(1) = nfreq*dfreq
367 0 : if_nfreq(1) = nfreq
368 0 : if_nstok(1) = nstok
369 0 : if_ref(1) = crpix4
370 0 : do i = 1, 4
371 0 : if_cstok(i,1) = ' '
372 : end do
373 0 : if_simul(1) = 1
374 0 : if_chain(1) = 1
375 : else
376 2 : freq = if_freq(1)
377 2 : nfreq = if_nfreq(1)
378 2 : if (if_nfreq(1).gt.1) then
379 2 : dfreq = if_bw(1)/(if_nfreq(1) - 1)
380 : else
381 0 : dfreq = if_bw(1)/if_nfreq(1)
382 : end if
383 2 : nstok = if_nstok(1)
384 : end if
385 2 : if (.not. su_found) then
386 0 : n_su = 1
387 0 : su_name(1) = object
388 0 : su_ra(1) = ra
389 0 : su_dec(1) = dec
390 : else
391 2 : object = su_name(1)
392 2 : ra = su_ra(1)
393 2 : dec = su_dec(1)
394 : C For single source, record possible pointing centre offset
395 2 : if (n_su.eq.1 .and. (pra.ne.0.0 .or. pdec.ne.0.0)) then
396 0 : su_pra(1) = pra
397 0 : su_pdec(1) = pdec
398 : end if
399 : end if
400 :
401 : C Tidy up.
402 2 : n_if = max(n_if, 1)
403 2 : new_antenna = .false.
404 2 : bufptr = 0
405 :
406 2 : jstat = 0
407 1503 : RETURN
408 :
409 : C----------------------- READ DATA GROUP HEADER ------------------------
410 1501 : 3000 if (.not.isopen) then
411 0 : call RPFERR ('File is not open.')
412 0 : jstat = -1
413 0 : RETURN
414 : end if
415 :
416 : C THE FOLLOWING POINTERS AND COUNTERS ARE USED HERE:
417 : C GRPLENGTH No. of visibilities in group
418 : C GRPPTR Pointer to next visibility in group to be read
419 : C BUFPTR Pointer to next word to be read in current buffer
420 : C BUFLEFT No. of words still to be read from current buffer
421 : C
422 : C Note that data are read in blocks of 5 records = 640 (4byte)
423 : C words.
424 :
425 1501 : grpptr = 1
426 1501 : if_no = 1
427 :
428 1501 : if (bufptr.eq.0 .or. bufptr.eq.641) then
429 1 : rp_iostat = AT_READ (lun, buffer)
430 1 : if (rp_iostat.ne.0) then
431 0 : if (rp_iostat.eq.-1) then
432 0 : jstat = 3
433 0 : RETURN
434 : end if
435 :
436 0 : call RPFERR ('I/O error reading data')
437 0 : jstat = -1
438 0 : RETURN
439 : end if
440 :
441 1 : jstat = SIMPLE (buffer, lun)
442 1 : if (jstat.ne.0) then
443 0 : rp_iostat = AT_UNREAD (lun, buffer)
444 0 : RETURN
445 : end if
446 :
447 1 : bufptr = 1
448 :
449 : end if
450 :
451 :
452 : C READ PARAMETERS FROM FITS FILE
453 : C FORMAT FROM RPFITS IS:
454 : C ------ VIS data ------------- ----------- SYSCAL data ----
455 : C (baseline > 0) (baseline = -1)
456 : C param 1=u in m 0.0
457 : C param 2=v in m 0.0
458 : C param 3=w in m 0.0
459 : C param 4=baseline number -1.0
460 : C param 5=UT in seconds sc_ut: UT in seconds
461 : C param 6= flag (if present) sc_ant
462 : C param 7= bin (if present) sc_if
463 : C param 8=if_no (if present) sc_q
464 : C param 9=sourceno (if present) sc_srcno
465 : C param 10=intbase (if present) intbase (if present)
466 :
467 1501 : 3100 bufleft = 641 - bufptr
468 :
469 : C End of scan?
470 1501 : call VAXI4 (i_buff(bufptr), i1)
471 1501 : endscan = i1.eq.illegal
472 :
473 1501 : if (.not.endscan .and. bufleft.ge.pcount) then
474 : C Old rpfits files may be padded with zeros, so check for u,
475 : C baseline no and UT all zero. Assume that if next vis
476 : C incomplete at end of buffer, next buffer will be all zeros.
477 :
478 1444 : call VAXI4 (i_buff(bufptr+3), i2)
479 1444 : call VAXI4 (i_buff(bufptr+4), i3)
480 1444 : endscan = i1.eq.0 .and. i2.eq.0 .and. i3.eq.0
481 : end if
482 :
483 1501 : if (endscan) then
484 19 : rp_iostat = AT_READ (lun, buffer)
485 19 : if (rp_iostat.ne.0) then
486 1 : if (rp_iostat.eq.-1) then
487 1 : jstat = 3
488 1 : RETURN
489 : end if
490 :
491 0 : call RPFERR ('I/O error reading header')
492 0 : jstat = -1
493 0 : RETURN
494 : end if
495 :
496 18 : jstat = SIMPLE (buffer, lun)
497 18 : if (jstat.ne.0) then
498 0 : rp_iostat = AT_UNREAD (lun, buffer)
499 0 : RETURN
500 : end if
501 :
502 18 : bufptr = 1
503 18 : jstat = 5
504 18 : RETURN
505 : end if
506 :
507 : C ------------NOW READ DATA -------------
508 :
509 1482 : if (bufleft.ge.pcount) then
510 :
511 : C If it will all fit in current buffer, then things are easy.
512 : call GETPARM (jstat, buffer, i_buff, bufptr, bufptr, buffer,
513 : : pcount, u, v, w, baseline, lun, ut, flag, bin, if_no,
514 1444 : : sourceno)
515 1444 : if (jstat.eq.-2) goto 3100
516 1444 : if (jstat.ne.0) RETURN
517 1444 : bufptr = bufptr+pcount
518 :
519 : else
520 : C We can recover only part of the group header. Dispose of what
521 : C we have, then read the remainder from the next batch of data
522 : C (pcount blocks).
523 :
524 247 : do i = 1,bufleft
525 247 : i_grphdr(i) = i_buff(bufptr+i-1)
526 : end do
527 38 : rp_iostat = AT_READ (lun, buffer)
528 38 : if (rp_iostat.ne.0) then
529 0 : if (rp_iostat.eq.-1) then
530 0 : jstat = 3
531 0 : RETURN
532 : end if
533 :
534 0 : call RPFERR ('I/O error reading data')
535 0 : jstat = -1
536 0 : RETURN
537 : end if
538 :
539 38 : jstat = SIMPLE (buffer, lun)
540 38 : if (jstat.ne.0) then
541 0 : rp_iostat = AT_UNREAD (lun, buffer)
542 0 : RETURN
543 : end if
544 :
545 38 : bufptr = pcount-bufleft
546 :
547 : C Extract bufptr items from the next buffer.
548 247 : do i = 1, bufptr
549 247 : i_grphdr(i+bufleft) = i_buff(i)
550 : end do
551 :
552 : call GETPARM (jstat, grphdr, i_grphdr, 1, bufptr, buffer,
553 : : pcount, u, v, w, baseline, lun, ut, flag, bin, if_no,
554 38 : : sourceno)
555 38 : if (jstat.eq.-2) goto 3100
556 38 : if (jstat.ne.0) RETURN
557 :
558 : C Set bufptr to the first visibility in the new buffer.
559 38 : bufptr = bufptr + 1
560 :
561 : end if
562 :
563 :
564 : C Determine GRPLENGTH.
565 1482 : if (baseline.eq.-1) then
566 57 : grplength = sc_q*sc_if*sc_ant
567 1425 : else if (if_no.gt.1) then
568 912 : grplength = if_nfreq(if_no)*if_nstok(if_no)
569 : else
570 513 : grplength = nstok*nfreq
571 : end if
572 :
573 1482 : if (baseline.eq.-1) go to 4000
574 :
575 : C--------------------- READ VISIBILITY DATA GROUP ----------------------
576 :
577 : C The RPFITS data format is determined by the value of NAXIS2:
578 : C
579 : C NAXIS2 word 1 word 2 word 3
580 : C ------ -------- -------- --------
581 : C 1 Real(vis) - -
582 : C 2 Real(vis) Imag(vis) -
583 : C 3 Real(vis) Imag(vis) Weight
584 :
585 1425 : if (data_format.lt.1 .or. data_format.gt.3) then
586 0 : call RPFERR ('NAXIS2 in file must be 1, 2, or 3.')
587 0 : jstat = -1
588 0 : RETURN
589 : end if
590 :
591 24985 : 3500 bufleft = 641 - bufptr
592 24985 : if (bufleft.ge.(data_format*(grplength-grpptr+1))) then
593 : C Entire group can be filled from existing buffer.
594 206302 : do i = grpptr, grplength
595 204877 : if (data_format.eq.1) then
596 0 : call VAXR4 (buffer(bufptr), vis(i))
597 : else
598 204877 : call VAXR4 (buffer(bufptr), r1)
599 204877 : call VAXR4 (buffer(bufptr+1), r2)
600 204877 : vis(i) = CMPLX(r1, r2)
601 :
602 204877 : if (data_format.eq.3) then
603 0 : call VAXR4 (buffer(bufptr+2), weight(i))
604 : end if
605 : end if
606 206302 : bufptr = bufptr + data_format
607 : end do
608 :
609 1425 : jstat = 0
610 1425 : RETURN
611 :
612 : else
613 : C Otherwise things are a bit more complicated, first read
614 : C complete visibilities in old buffer.
615 47120 : bufleft3 = bufleft/data_format
616 7349865 : do i = 1, bufleft3
617 7326305 : if (data_format.eq.1) then
618 0 : call VAXR4 (buffer(bufptr), vis(grpptr+i-1))
619 : else
620 7326305 : call VAXR4 (buffer(bufptr), r1)
621 7326305 : call VAXR4 (buffer(bufptr+1), r2)
622 7326305 : vis(grpptr+i-1) = CMPLX(r1, r2)
623 :
624 7326305 : if (data_format.eq.3) then
625 0 : call VAXR4 (buffer(bufptr+2), weight(grpptr+i-1))
626 : end if
627 : end if
628 7349865 : bufptr = bufptr + data_format
629 : end do
630 23560 : grpptr = grpptr + bufleft3
631 :
632 : C Read the fraction of a visibility left in old buffer.
633 : C Should not happen for data_format = 1.
634 47120 : bufleft = bufleft - data_format*bufleft3
635 23560 : if (bufleft.eq.1) then
636 11286 : call VAXR4 (buffer(640), revis)
637 12274 : else if (bufleft.eq.2 .and. data_format.eq.3) then
638 0 : call VAXR4 (buffer(639), r1)
639 0 : call VAXR4 (buffer(640), r2)
640 0 : vis(grpptr) = CMPLX(r1, r2)
641 : end if
642 :
643 : C Now read in a new buffer.
644 23560 : rp_iostat = AT_READ (lun, buffer)
645 23560 : if (rp_iostat.ne.0) then
646 0 : if (rp_iostat.eq.-1) then
647 0 : jstat = 3
648 0 : RETURN
649 : end if
650 :
651 0 : call RPFERR ('I/O error reading data')
652 0 : jstat = -1
653 0 : RETURN
654 : end if
655 :
656 23560 : jstat = SIMPLE (buffer, lun)
657 23560 : if (jstat.ne.0) then
658 0 : rp_iostat = AT_UNREAD (lun, buffer)
659 0 : RETURN
660 : end if
661 :
662 : C Fill any incomplete visibility (data_format = 2 or 3 only).
663 23560 : if (bufleft.eq.0) then
664 12274 : bufptr = 1
665 :
666 11286 : else if (bufleft.eq.1) then
667 11286 : call VAXR4 (buffer(1), r1)
668 11286 : vis(grpptr) = CMPLX(revis, r1)
669 11286 : if (data_format.eq.3) then
670 0 : call VAXR4 (buffer(2), weight(grpptr))
671 : end if
672 11286 : grpptr = grpptr + 1
673 11286 : bufptr = data_format
674 :
675 0 : else if (bufleft.eq.2 .and. data_format.eq.3) then
676 0 : call VAXR4 (buffer(1), weight(grpptr))
677 0 : grpptr = grpptr + 1
678 0 : bufptr = 2
679 : end if
680 : end if
681 :
682 : C Return to pick up the rest of the group.
683 23617 : go to 3500
684 :
685 : C----------------------- READ SYSCAL DATA GROUP ------------------------
686 :
687 : C Note that in this context GRPLENGTH is in units of words, not
688 : C visibilities.
689 :
690 57 : 4000 bufleft = 641 - bufptr
691 57 : if (bufleft.ge.(grplength-grpptr+1)) then
692 :
693 : C Entire group can be filled from existing buffer.
694 6441 : do i = grpptr, grplength
695 6384 : call VAXR4 (buffer(bufptr), sc_buf(i))
696 6441 : bufptr = bufptr + 1
697 : end do
698 :
699 57 : jstat = 0
700 57 : RETURN
701 :
702 : else
703 : C Otherwise read complete visibilities in old buffer.
704 0 : do i = 1, bufleft
705 0 : call VAXR4 (buffer(bufptr), sc_buf(grpptr+i-1))
706 0 : bufptr = bufptr + 1
707 : end do
708 0 : grpptr = grpptr + bufleft
709 :
710 : C Then read in a new buffer.
711 0 : rp_iostat = AT_READ (lun, buffer)
712 0 : if (rp_iostat.ne.0) then
713 0 : if (rp_iostat.eq.-1) then
714 0 : jstat = 3
715 0 : RETURN
716 : end if
717 :
718 0 : call RPFERR ('I/O error reading data')
719 0 : jstat = -1
720 0 : RETURN
721 : end if
722 :
723 0 : jstat = SIMPLE (buffer, lun)
724 0 : if (jstat.ne.0) then
725 0 : rp_iostat = AT_UNREAD (lun, buffer)
726 0 : RETURN
727 : end if
728 0 : bufptr = 1
729 : end if
730 :
731 : C Go back to pick up the rest of the group.
732 2 : go to 4000
733 :
734 : C--------------------------- CLOSE FITS FILE ---------------------------
735 :
736 2 : 5000 if (isopen) then
737 2 : rp_iostat = AT_CLOSE (lun)
738 2 : if (rp_iostat.ne.0) then
739 0 : call RPFERR ('I/O error closing file')
740 0 : jstat = -1
741 0 : RETURN
742 : end if
743 2 : isopen = .false.
744 : end if
745 :
746 2 : jstat = 0
747 2 : RETURN
748 :
749 : C------------------------- SKIP TO END OF FILE -------------------------
750 :
751 0 : 6000 if (.not.isopen) then
752 0 : call RPFERR ('File is not open.')
753 0 : jstat = -1
754 0 : RETURN
755 : end if
756 :
757 0 : rp_iostat = AT_SKIP_EOF (lun)
758 0 : if (rp_iostat.eq.-1) then
759 0 : jstat = 3
760 : else
761 0 : call RPFERR ('I/O error skipping to EOF')
762 0 : jstat = -1
763 0 : RETURN
764 : end if
765 :
766 0 : return
767 : end
768 :
769 : C-----------------------------------------------------------------------
770 :
771 23619 : integer function SIMPLE (buffer, lun)
772 :
773 : C-----------------------------------------------------------------------
774 : C SIMPLE tests for the start of a new header or FG (flag) table.
775 : C Reads the FG table if encountered.
776 : C-----------------------------------------------------------------------
777 :
778 : include 'rpfits.inc'
779 :
780 : logical endhdr
781 : integer ierr, j, lun
782 : character m(80)*32, terr*2
783 : real buffer(640)
784 :
785 : C Assume not.
786 23619 : SIMPLE = 0
787 :
788 : C Write first 8 characters from buffer into character string.
789 70857 : write (m(1)(1:8),'(2a4)') (buffer(j),j=1,2)
790 :
791 23619 : if (m(1)(1:6).eq.'SIMPLE') then
792 : C Start of header.
793 2 : SIMPLE = 1
794 :
795 23617 : else if (m(1)(1:8).eq.'FG TABLE') then
796 : C Start of FG (flag) table.
797 0 : SIMPLE = 4
798 :
799 0 : write (m, '(32(20a4,:,/))') (buffer(j),j=1,640)
800 0 : call RPFITS_READ_TABLE (lun, m, -1, endhdr, terr, ierr)
801 0 : if (ierr.ne.0) then
802 0 : if (ierr.eq.1) then
803 0 : call RPFERR ('FG table contains too many entries.')
804 0 : SIMPLE = -1
805 0 : else if (rp_iostat.lt.0) then
806 0 : SIMPLE = 3
807 : else
808 0 : SIMPLE = -1
809 0 : call RPFERR ('I/O error reading FG table')
810 : end if
811 : end if
812 : end if
813 :
814 23619 : return
815 : end
816 :
817 : C-----------------------------------------------------------------------
818 :
819 1482 : subroutine GETPARM (jstat, grphdr, i_grphdr, grpptr, bufptr,
820 : : buffer, pcount, u, v, w, baseline, lun, ut,
821 : : flag, bin, if_no, sourceno)
822 :
823 : C-----------------------------------------------------------------------
824 : C Read group header parameters from grphdr and check validity.
825 : C If invalid, scan through the data until valid data are found and
826 : C return the new buffer and bufptr.
827 : C
828 : C jstat is 0 on exit for immediate success, or -2 if success was
829 : C achieved after skipping data, or -1 for a total lack of success.
830 : C-----------------------------------------------------------------------
831 :
832 : include 'rpfits.inc'
833 :
834 : logical ILLPARM
835 : integer baseline, bin, bufptr, flag, grpptr, i_grphdr(11),
836 : : iant, if_no, iif, iq, jstat, lun, pcount, sourceno
837 : real grphdr(11), buffer(640), rbase, u, v, w, ut
838 :
839 : C First 5 parameters are always there - you hope!
840 1482 : call VAXR4 (grphdr(grpptr), u)
841 1482 : call VAXR4 (grphdr(grpptr+1), v)
842 1482 : call VAXR4 (grphdr(grpptr+2), w)
843 1482 : call VAXR4 (grphdr(grpptr+3), rbase)
844 1482 : call VAXR4 (grphdr(grpptr+4), ut)
845 :
846 1482 : if (rbase.lt.0.0) then
847 : C Syscal parameters.
848 57 : call VAXI4 (i_grphdr(grpptr+5), iant)
849 57 : call VAXI4 (i_grphdr(grpptr+6), iif)
850 57 : call VAXI4 (i_grphdr(grpptr+7), iq)
851 : else
852 : C IF number.
853 1425 : call VAXI4 (i_grphdr(grpptr+7), iif)
854 :
855 1425 : if (pcount.ge.11) then
856 : C Otherwise, data_format comes from NAXIS2.
857 1425 : call VAXI4 (i_grphdr(grpptr+10), data_format)
858 : end if
859 : end if
860 :
861 : C Check for illegal parameters.
862 1482 : if (ILLPARM(u, v, w, rbase, ut, iant, iif, iq)) then
863 : C This can be caused by a bad block, so look for more data.
864 0 : call RPFERR ('Corrupted data encountered, skipping...')
865 0 : call SKIPTHRU (jstat, bufptr, buffer, lun, pcount)
866 0 : RETURN
867 : end if
868 :
869 : C Looks ok, pick up remaining parameters.
870 1482 : baseline = NINT(rbase)
871 1482 : if (baseline.eq.-1) then
872 : C Syscal parameters.
873 57 : sc_ut = ut
874 57 : sc_ant = iant
875 57 : sc_if = iif
876 57 : sc_q = iq
877 57 : call VAXI4 (i_grphdr(grpptr+8), sc_srcno)
878 57 : if (pcount.gt.9) then
879 57 : call VAXR4 (REAL(i_grphdr(grpptr+9)), intbase)
880 : else
881 0 : intbase = 0.0
882 : end if
883 :
884 1425 : else if (pcount.gt.5) then
885 1425 : call VAXI4 (i_grphdr(grpptr+5), flag)
886 1425 : call VAXI4 (i_grphdr(grpptr+6), bin)
887 1425 : call VAXI4 (i_grphdr(grpptr+7), if_no)
888 1425 : call VAXI4 (i_grphdr(grpptr+8), sourceno)
889 :
890 1425 : if (pcount.gt.9) then
891 1425 : call VAXR4 (grphdr(grpptr+9), intbase)
892 : else
893 0 : intbase = intime
894 : end if
895 : end if
896 :
897 1482 : jstat = 0
898 1482 : return
899 0 : end
900 :
901 : *-----------------------------------------------------------------------
902 :
903 1482 : logical function ILLPARM (u, v, w, rbase, ut, iant, iif, iq)
904 :
905 : *-----------------------------------------------------------------------
906 : * Check for any illegal parameters; return true if so.
907 : *-----------------------------------------------------------------------
908 :
909 : include 'rpfits.inc'
910 :
911 : integer baseline, iant, iant1, iant2, iif, iq
912 : real u, ut, v, w, rbase
913 :
914 2964 : if (data_format.lt.1 .or. data_format.gt.3) then
915 : * Invalid data format.
916 0 : ILLPARM = .true.
917 :
918 : else if (abs(u).gt.1e10 .or.
919 1482 : : abs(v).gt.1e10 .or.
920 : : abs(w).gt.1e10) then
921 : * Invalid visibility coordinate.
922 0 : ILLPARM = .true.
923 :
924 1482 : else if (rbase.lt.-1.1 .or. rbase.gt.(257*nant+0.1)) then
925 : * Invalid baseline number.
926 0 : ILLPARM = .true.
927 :
928 1482 : else if (ut.lt.0.0 .or. ut.gt.172800.0) then
929 : * Invalid time.
930 0 : ILLPARM = .true.
931 :
932 : else
933 : * Baseline can now safely be converted to integer.
934 1482 : baseline = NINT(rbase)
935 :
936 1482 : if (ABS(rbase - baseline).gt.0.001) then
937 : * This value is not close enough to an integer to be valid.
938 0 : ILLPARM = .true.
939 :
940 : else
941 1482 : if (baseline.eq.-1) then
942 : * Syscal record.
943 : ILLPARM = iant.lt.1 .or. iant.gt.ant_max .or.
944 : : iif.lt.1 .or. iif.gt.max_if .or.
945 57 : : iq.lt.1 .or. iq.gt.100
946 :
947 : else
948 : * Data record.
949 1425 : iant1 = baseline/256
950 1425 : iant2 = MOD(baseline,256)
951 : ILLPARM = iant1.lt.1 .or. iant1.gt.nant .or.
952 : : iant2.lt.1 .or. iant2.gt.nant .or.
953 4275 : : iif.lt.0 .or. iif.gt.max_if
954 : end if
955 : end if
956 : end if
957 :
958 1482 : return
959 : end
960 :
961 : C-----------------------------------------------------------------------
962 :
963 0 : subroutine SKIPTHRU (jstat, bufptr, buffer, lun, pcount)
964 :
965 : C-----------------------------------------------------------------------
966 : C Skip through data looking for recognisable data or header.
967 : C
968 : C Returns jstat = -2 if successful.
969 : C
970 : C rpn 17/11/90
971 : C-----------------------------------------------------------------------
972 :
973 : include 'rpfits.inc'
974 :
975 : logical ILLPARM
976 : integer AT_READ, AT_UNREAD, bufptr, i, iant, iif, iq, j, jstat,
977 : : lun, pcount, SIMPLE
978 : real buffer(640), rbase, u, ut, v, w
979 :
980 0 : do 10 j = 1, 1000
981 : C Read a new block; the remainder of the old one is unlikely to
982 : C contain anything useful (and at most one integration).
983 0 : rp_iostat = AT_READ (lun, buffer)
984 0 : if (rp_iostat.ne.0) then
985 0 : if (rp_iostat.eq.-1) then
986 0 : jstat = 3
987 0 : RETURN
988 : end if
989 :
990 0 : call RPFERR ('Read error')
991 0 : jstat = -1
992 0 : RETURN
993 : end if
994 :
995 : C Check to see if it's a header block.
996 0 : jstat = SIMPLE (buffer, lun)
997 0 : if (jstat.ne.0) then
998 0 : rp_iostat = AT_UNREAD (lun, buffer)
999 0 : RETURN
1000 : end if
1001 0 : bufptr = 1
1002 :
1003 : C Scan through the block looking for something legal.
1004 0 : do i = 1, 640
1005 0 : call VAXR4 (buffer(bufptr), u)
1006 0 : call VAXR4 (buffer(bufptr+1), v)
1007 0 : call VAXR4 (buffer(bufptr+2), w)
1008 0 : call VAXR4 (buffer(bufptr+3), rbase)
1009 0 : call VAXR4 (buffer(bufptr+4), ut)
1010 :
1011 0 : if (rbase.lt.0.0) then
1012 : C Syscal parameters.
1013 0 : call VAXI4 (buffer(bufptr+5), iant)
1014 0 : call VAXI4 (buffer(bufptr+6), iif)
1015 0 : call VAXI4 (buffer(bufptr+7), iq)
1016 : else
1017 : C IF number.
1018 0 : call VAXI4 (buffer(bufptr+7), iif)
1019 :
1020 0 : if (pcount.ge.11) then
1021 : C Otherwise, data_format comes from NAXIS2.
1022 0 : call VAXI4 (buffer(bufptr+10), data_format)
1023 : end if
1024 : end if
1025 :
1026 0 : if (.not.ILLPARM(u, v, w, rbase, ut, iant, iif, iq)) then
1027 0 : goto 999
1028 : end if
1029 :
1030 0 : bufptr = bufptr + 1
1031 0 : if (bufptr.gt.632) goto 10
1032 : end do
1033 0 : 10 continue
1034 :
1035 : C Success!
1036 0 : 999 jstat = -2
1037 0 : return
1038 5814 : end
|
