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 0 : 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 0 : rp_iostat = 0
52 0 : errmsg = ''
53 :
54 0 : open_only = jstat.eq.-3
55 :
56 0 : if (jstat.eq.-3) go to 1000
57 0 : if (jstat.eq.-2) go to 1000
58 0 : if (jstat.eq.-1) go to 2000
59 0 : if (jstat.eq.0) go to 3000
60 0 : 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 0 : RETURN
67 :
68 : C--------------------------- OPEN FITS FILE ----------------------------
69 :
70 0 : 1000 if (isopen) then
71 0 : call RPFERR ('File is already open.')
72 0 : jstat = -1
73 0 : RETURN
74 : end if
75 :
76 0 : rp_iostat = AT_OPEN_READ (file, async, lun)
77 0 : if (rp_iostat.ne.0) then
78 0 : call RPFERR ('File open error')
79 0 : jstat = -1
80 0 : RETURN
81 : end if
82 0 : isopen = .true.
83 :
84 0 : if (open_only) then
85 0 : jstat = 0
86 0 : RETURN
87 : end if
88 :
89 : C----------------------------- READ HEADER -----------------------------
90 :
91 0 : 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 0 : bufptr = 0
98 0 : n_if = 0
99 0 : icard = 1
100 0 : if (ncard.lt.0) ncard = -1
101 0 : an_found = .false.
102 0 : if_found = .false.
103 0 : su_found = .false.
104 0 : fg_found = .false.
105 0 : nx_found = .false.
106 0 : mt_found = .false.
107 0 : cu_found = .false.
108 0 : pra = 0.0
109 0 : pdec = 0.0
110 :
111 : C Look for start of next header.
112 0 : starthdr = .false.
113 0 : do while (.not.starthdr)
114 0 : rp_iostat = AT_READ (lun, buffer)
115 0 : 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 0 : jstat = SIMPLE (buffer, lun)
127 0 : if (jstat.eq.1) then
128 : C Start of header.
129 0 : 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 0 : write (m, '(32(20a4,:,/))') (buffer(j),j=1,640)
143 : end do
144 :
145 : C Scan through header, getting the interesting bits.
146 0 : endhdr = .false.
147 0 : do 2500 while (.not.endhdr)
148 0 : if (.not.starthdr) then
149 0 : rp_iostat = AT_READ (lun, buffer)
150 0 : write (m,'(32(20a4,:,/))') (buffer(j),j=1,640)
151 0 : 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 0 : starthdr = .false.
164 0 : version = ' '
165 0 : do 2400 i = 1, 32
166 : C Parse the RPFITS keyword and keyvalue.
167 0 : keyword = m(i)(1:8)
168 :
169 0 : if (m(i)(11:11).eq.'''') then
170 : C Must be a character value.
171 0 : keyvalue = m(i)(12:31)
172 0 : do j = 1, 20
173 0 : if (keyvalue(j:j).eq.'''') then
174 : C Strip off the trailing apostrophe.
175 0 : keyvalue(j:) = ' '
176 : end if
177 : end do
178 : else
179 0 : keyvalue = m(i)(11:30)
180 : end if
181 :
182 : C Lexical chop based on the first letter of the keyword name.
183 0 : if (keyword(:1).le.'C') then
184 : C Keyword names beginning with A to C.
185 0 : if (keyword.eq.'ALTRVAL ') then
186 0 : read (keyvalue, '(g20.12)') vel1
187 0 : else if (keyword.eq.'BUNIT') then
188 0 : bunit = keyvalue(:16)
189 0 : else if (keyword.eq.'CAL') then
190 0 : cal = keyvalue(:16)
191 0 : else if (keyword.eq.'CDELT4') then
192 0 : read (keyvalue, '(g20.12)') dfreq
193 0 : else if (keyword.eq.'CRPIX4') then
194 0 : read (keyvalue, '(g20.12)') crpix4
195 0 : else if (keyword.eq.'CRVAL4') then
196 0 : read (keyvalue, '(g20.12)') freq
197 0 : else if (keyword.eq.'CRVAL5') then
198 0 : read (keyvalue, '(g20.12)') ra
199 0 : if (ra.lt.0d0) ra = ra + d2pi
200 0 : else if (keyword.eq.'CRVAL6') then
201 0 : read (keyvalue, '(g20.12)') dec
202 : end if
203 :
204 0 : else if (keyword(:1).le.'E') then
205 : C Keyword names beginning with D or E.
206 0 : if (keyword.eq.'DATE') then
207 : C Fix old-format dates.
208 0 : call datfit(keyvalue(:10), datwrit, ierr)
209 0 : else if (keyword.eq.'DATE-OBS') then
210 : C Fix old-format dates.
211 0 : call datfit(keyvalue(:10), datobs, ierr)
212 0 : datsys = m(i)(35:36)
213 0 : if (datsys.eq.'UT D') datsys = 'UT'
214 0 : else if (keyword.eq.'DEFEAT ') then
215 0 : read (keyvalue, '(i20)') rp_defeat
216 0 : else if (keyword.eq.'DJMREFP ') then
217 0 : read (keyvalue, '(g20.12)') rp_djmrefp
218 0 : else if (keyword.eq.'DJMREFT ') then
219 0 : read (keyvalue, '(g20.12)') rp_djmreft
220 0 : else if (keyword.eq.'END') then
221 : C END card.
222 0 : endhdr = .true.
223 0 : else if (keyword(1:5).eq.'EPHEM') then
224 0 : read (keyword(6:7), '(i2)') k
225 0 : read (keyvalue, '(g20.12)') rp_c(k)
226 0 : else if (keyword.eq.'EPOCH') then
227 0 : coord = keyvalue(:16)
228 : end if
229 :
230 0 : else if (keyword(:1).le.'N') then
231 : C Keyword names beginning with F to N.
232 0 : if (keyword.eq.'GCOUNT') then
233 0 : read (keyvalue, '(i20)') ncount
234 0 : 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 0 : else if (keyword.eq.'INSTRUME') then
238 0 : instrument = keyvalue(:16)
239 0 : else if (keyword.eq.'INTIME') then
240 0 : read (keyvalue, '(i20)') intime
241 0 : else if (keyword.eq.'NAXIS2') then
242 0 : read (keyvalue, '(i20)') data_format
243 0 : write_wt = data_format.eq.3
244 0 : else if (keyword.eq.'NAXIS3') then
245 0 : read (keyvalue, '(i20)') nstok
246 0 : else if (keyword.eq.'NAXIS4') then
247 0 : read (keyvalue, '(i20)') nfreq
248 0 : 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 0 : else if (keyword(:1).le.'P') then
254 : C Keyword names beginning with M to P.
255 0 : if (keyword.eq.'OBJECT') then
256 0 : object = keyvalue(:16)
257 0 : else if (keyword.eq.'OBSERVER') then
258 0 : rp_observer = keyvalue(:16)
259 0 : else if (keyword.eq.'OBSTYPE') then
260 0 : obstype = keyvalue(:16)
261 0 : else if (keyword.eq.'PCOUNT') then
262 0 : read (keyvalue, '(i20)') pcount
263 0 : else if (keyword.eq.'PMDEC') then
264 0 : read (keyvalue, '(g20.12)') pm_dec
265 0 : else if (keyword.eq.'PMEPOCH') then
266 0 : read (keyvalue, '(g20.12)') pm_epoch
267 0 : else if (keyword.eq.'PMRA') then
268 0 : read (keyvalue, '(g20.12)') pm_ra
269 0 : else if (keyword.eq.'PNTCENTR') then
270 0 : read (m(i)(11:35),'(g12.9,1x,g12.9)') pra,pdec
271 0 : 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 0 : if (keyword.eq.'RESTFREQ') then
279 0 : read (keyvalue, '(g20.12)') rfreq
280 0 : else if (keyword.eq.'RPFITS ') then
281 0 : rpfitsversion = keyvalue
282 0 : else if (keyword.eq.'SCANS ') then
283 0 : read (keyvalue, '(i20)') nscan
284 0 : else if (keyword(1:6).eq.'TABLE ') then
285 : C Sort out tables.
286 0 : call RPFITS_READ_TABLE (lun, m, i, endhdr, terr, ierr)
287 0 : 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 0 : 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 0 : else if (keyword.eq.'UTCMTAI ') then
307 0 : read (keyvalue, '(g20.12)') rp_utcmtai
308 0 : else if (keyword.eq.'VELREF ') then
309 0 : read (keyvalue, '(i20)') ivelref
310 0 : else if (keyword.eq.'VERSION ') then
311 0 : version = keyvalue
312 : end if
313 : end if
314 :
315 : C Write into "cards" array if necessary.
316 0 : 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 0 : else if (ncard.lt.0) then
323 0 : if (icard.le.max_card .and. .not.endhdr) then
324 0 : card(-ncard) = m(i)
325 0 : icard = icard + 1
326 0 : ncard = ncard - 1
327 : end if
328 : end if
329 :
330 : C Antenna parameters.
331 0 : 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 0 : if (ENDHDR) go to 2500
350 0 : 2400 continue
351 : 2500 continue
352 0 : ncard = ABS(ncard)
353 :
354 : C Set up for reading data.
355 0 : 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 0 : 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 0 : freq = if_freq(1)
377 0 : nfreq = if_nfreq(1)
378 0 : if (if_nfreq(1).gt.1) then
379 0 : dfreq = if_bw(1)/(if_nfreq(1) - 1)
380 : else
381 0 : dfreq = if_bw(1)/if_nfreq(1)
382 : end if
383 0 : nstok = if_nstok(1)
384 : end if
385 0 : 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 0 : object = su_name(1)
392 0 : ra = su_ra(1)
393 0 : dec = su_dec(1)
394 : C For single source, record possible pointing centre offset
395 0 : 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 0 : n_if = max(n_if, 1)
403 0 : new_antenna = .false.
404 0 : bufptr = 0
405 :
406 0 : jstat = 0
407 0 : RETURN
408 :
409 : C----------------------- READ DATA GROUP HEADER ------------------------
410 0 : 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 0 : grpptr = 1
426 0 : if_no = 1
427 :
428 0 : if (bufptr.eq.0 .or. bufptr.eq.641) then
429 0 : rp_iostat = AT_READ (lun, buffer)
430 0 : 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 0 : jstat = SIMPLE (buffer, lun)
442 0 : if (jstat.ne.0) then
443 0 : rp_iostat = AT_UNREAD (lun, buffer)
444 0 : RETURN
445 : end if
446 :
447 0 : 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 0 : 3100 bufleft = 641 - bufptr
468 :
469 : C End of scan?
470 0 : call VAXI4 (i_buff(bufptr), i1)
471 0 : endscan = i1.eq.illegal
472 :
473 0 : 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 0 : call VAXI4 (i_buff(bufptr+3), i2)
479 0 : call VAXI4 (i_buff(bufptr+4), i3)
480 0 : endscan = i1.eq.0 .and. i2.eq.0 .and. i3.eq.0
481 : end if
482 :
483 0 : if (endscan) then
484 0 : rp_iostat = AT_READ (lun, buffer)
485 0 : if (rp_iostat.ne.0) then
486 0 : if (rp_iostat.eq.-1) then
487 0 : jstat = 3
488 0 : 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 0 : jstat = SIMPLE (buffer, lun)
497 0 : if (jstat.ne.0) then
498 0 : rp_iostat = AT_UNREAD (lun, buffer)
499 0 : RETURN
500 : end if
501 :
502 0 : bufptr = 1
503 0 : jstat = 5
504 0 : RETURN
505 : end if
506 :
507 : C ------------NOW READ DATA -------------
508 :
509 0 : 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 0 : : sourceno)
515 0 : if (jstat.eq.-2) goto 3100
516 0 : if (jstat.ne.0) RETURN
517 0 : 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 0 : do i = 1,bufleft
525 0 : i_grphdr(i) = i_buff(bufptr+i-1)
526 : end do
527 0 : rp_iostat = AT_READ (lun, buffer)
528 0 : 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 0 : jstat = SIMPLE (buffer, lun)
540 0 : if (jstat.ne.0) then
541 0 : rp_iostat = AT_UNREAD (lun, buffer)
542 0 : RETURN
543 : end if
544 :
545 0 : bufptr = pcount-bufleft
546 :
547 : C Extract bufptr items from the next buffer.
548 0 : do i = 1, bufptr
549 0 : 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 0 : : sourceno)
555 0 : if (jstat.eq.-2) goto 3100
556 0 : if (jstat.ne.0) RETURN
557 :
558 : C Set bufptr to the first visibility in the new buffer.
559 0 : bufptr = bufptr + 1
560 :
561 : end if
562 :
563 :
564 : C Determine GRPLENGTH.
565 0 : if (baseline.eq.-1) then
566 0 : grplength = sc_q*sc_if*sc_ant
567 0 : else if (if_no.gt.1) then
568 0 : grplength = if_nfreq(if_no)*if_nstok(if_no)
569 : else
570 0 : grplength = nstok*nfreq
571 : end if
572 :
573 0 : 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 0 : 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 0 : 3500 bufleft = 641 - bufptr
592 0 : if (bufleft.ge.(data_format*(grplength-grpptr+1))) then
593 : C Entire group can be filled from existing buffer.
594 0 : do i = grpptr, grplength
595 0 : if (data_format.eq.1) then
596 0 : call VAXR4 (buffer(bufptr), vis(i))
597 : else
598 0 : call VAXR4 (buffer(bufptr), r1)
599 0 : call VAXR4 (buffer(bufptr+1), r2)
600 0 : vis(i) = CMPLX(r1, r2)
601 :
602 0 : if (data_format.eq.3) then
603 0 : call VAXR4 (buffer(bufptr+2), weight(i))
604 : end if
605 : end if
606 0 : bufptr = bufptr + data_format
607 : end do
608 :
609 0 : jstat = 0
610 0 : RETURN
611 :
612 : else
613 : C Otherwise things are a bit more complicated, first read
614 : C complete visibilities in old buffer.
615 0 : bufleft3 = bufleft/data_format
616 0 : do i = 1, bufleft3
617 0 : if (data_format.eq.1) then
618 0 : call VAXR4 (buffer(bufptr), vis(grpptr+i-1))
619 : else
620 0 : call VAXR4 (buffer(bufptr), r1)
621 0 : call VAXR4 (buffer(bufptr+1), r2)
622 0 : vis(grpptr+i-1) = CMPLX(r1, r2)
623 :
624 0 : if (data_format.eq.3) then
625 0 : call VAXR4 (buffer(bufptr+2), weight(grpptr+i-1))
626 : end if
627 : end if
628 0 : bufptr = bufptr + data_format
629 : end do
630 0 : 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 0 : bufleft = bufleft - data_format*bufleft3
635 0 : if (bufleft.eq.1) then
636 0 : call VAXR4 (buffer(640), revis)
637 0 : 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 0 : rp_iostat = AT_READ (lun, buffer)
645 0 : 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 0 : jstat = SIMPLE (buffer, lun)
657 0 : 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 0 : if (bufleft.eq.0) then
664 0 : bufptr = 1
665 :
666 0 : else if (bufleft.eq.1) then
667 0 : call VAXR4 (buffer(1), r1)
668 0 : vis(grpptr) = CMPLX(revis, r1)
669 0 : if (data_format.eq.3) then
670 0 : call VAXR4 (buffer(2), weight(grpptr))
671 : end if
672 0 : grpptr = grpptr + 1
673 0 : 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 0 : 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 0 : 4000 bufleft = 641 - bufptr
691 0 : if (bufleft.ge.(grplength-grpptr+1)) then
692 :
693 : C Entire group can be filled from existing buffer.
694 0 : do i = grpptr, grplength
695 0 : call VAXR4 (buffer(bufptr), sc_buf(i))
696 0 : bufptr = bufptr + 1
697 : end do
698 :
699 0 : jstat = 0
700 0 : 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 0 : go to 4000
733 :
734 : C--------------------------- CLOSE FITS FILE ---------------------------
735 :
736 0 : 5000 if (isopen) then
737 0 : rp_iostat = AT_CLOSE (lun)
738 0 : 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 0 : isopen = .false.
744 : end if
745 :
746 0 : jstat = 0
747 0 : 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 0 : 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 0 : SIMPLE = 0
787 :
788 : C Write first 8 characters from buffer into character string.
789 0 : write (m(1)(1:8),'(2a4)') (buffer(j),j=1,2)
790 :
791 0 : if (m(1)(1:6).eq.'SIMPLE') then
792 : C Start of header.
793 0 : SIMPLE = 1
794 :
795 0 : 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 0 : return
815 : end
816 :
817 : C-----------------------------------------------------------------------
818 :
819 0 : 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 0 : call VAXR4 (grphdr(grpptr), u)
841 0 : call VAXR4 (grphdr(grpptr+1), v)
842 0 : call VAXR4 (grphdr(grpptr+2), w)
843 0 : call VAXR4 (grphdr(grpptr+3), rbase)
844 0 : call VAXR4 (grphdr(grpptr+4), ut)
845 :
846 0 : if (rbase.lt.0.0) then
847 : C Syscal parameters.
848 0 : call VAXI4 (i_grphdr(grpptr+5), iant)
849 0 : call VAXI4 (i_grphdr(grpptr+6), iif)
850 0 : call VAXI4 (i_grphdr(grpptr+7), iq)
851 : else
852 : C IF number.
853 0 : call VAXI4 (i_grphdr(grpptr+7), iif)
854 :
855 0 : if (pcount.ge.11) then
856 : C Otherwise, data_format comes from NAXIS2.
857 0 : call VAXI4 (i_grphdr(grpptr+10), data_format)
858 : end if
859 : end if
860 :
861 : C Check for illegal parameters.
862 0 : 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 0 : baseline = NINT(rbase)
871 0 : if (baseline.eq.-1) then
872 : C Syscal parameters.
873 0 : sc_ut = ut
874 0 : sc_ant = iant
875 0 : sc_if = iif
876 0 : sc_q = iq
877 0 : call VAXI4 (i_grphdr(grpptr+8), sc_srcno)
878 0 : if (pcount.gt.9) then
879 0 : call VAXR4 (REAL(i_grphdr(grpptr+9)), intbase)
880 : else
881 0 : intbase = 0.0
882 : end if
883 :
884 0 : else if (pcount.gt.5) then
885 0 : call VAXI4 (i_grphdr(grpptr+5), flag)
886 0 : call VAXI4 (i_grphdr(grpptr+6), bin)
887 0 : call VAXI4 (i_grphdr(grpptr+7), if_no)
888 0 : call VAXI4 (i_grphdr(grpptr+8), sourceno)
889 :
890 0 : if (pcount.gt.9) then
891 0 : call VAXR4 (grphdr(grpptr+9), intbase)
892 : else
893 0 : intbase = intime
894 : end if
895 : end if
896 :
897 0 : jstat = 0
898 0 : return
899 0 : end
900 :
901 : *-----------------------------------------------------------------------
902 :
903 0 : 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 0 : 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 0 : : abs(v).gt.1e10 .or.
920 : : abs(w).gt.1e10) then
921 : * Invalid visibility coordinate.
922 0 : ILLPARM = .true.
923 :
924 0 : else if (rbase.lt.-1.1 .or. rbase.gt.(257*nant+0.1)) then
925 : * Invalid baseline number.
926 0 : ILLPARM = .true.
927 :
928 0 : 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 0 : baseline = NINT(rbase)
935 :
936 0 : 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 0 : 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 0 : : iq.lt.1 .or. iq.gt.100
946 :
947 : else
948 : * Data record.
949 0 : iant1 = baseline/256
950 0 : iant2 = MOD(baseline,256)
951 : ILLPARM = iant1.lt.1 .or. iant1.gt.nant .or.
952 : : iant2.lt.1 .or. iant2.gt.nant .or.
953 0 : : iif.lt.0 .or. iif.gt.max_if
954 : end if
955 : end if
956 : end if
957 :
958 0 : 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 0 : end
|
