Line data Source code
1 : C-----------------------------------------------------------------------
2 : C
3 : C SUBROUTINE RPFITSOUT
4 : C
5 : C $Id: rpfitsout.f,v 1.35 2015/01/19 05:13:00 rey052 Exp $
6 : C-----------------------------------------------------------------------
7 0 : subroutine RPFITSOUT (jstat, vis, weight, baseline, ut, u, v,
8 : + w, flag, bin, if_no, sourceno)
9 :
10 : C This routine is for writing an RPFITS file.
11 : C Its function when called depends primarily on the value of
12 : C JSTAT, as follows:
13 : C JSTAT=-3......Open file
14 : C JSTAT=-2......Open file and write a header
15 : C JSTAT=-1......Write a header
16 : C JSTAT=0.......Write a data group to the file
17 : C JSTAT=1.......Close the file.
18 : C JSTAT=2.......Write an FG table to the end of the data
19 : C JSTAT=3.......Flush buffer, close file, reopen file
20 : C When JSTAT=1 or -1, all the other arguments are dummy, and are
21 : C left unchanged.
22 : C When JSTAT=0, all the parameters and the contents of common blocks
23 : C should be set before entry. They are all unchanged on exit, except
24 : C DATOBS which, if blank on entry, will contain the current UT
25 : C date on exit.
26 : C
27 : C On exit, JSTAT indicates the success of the operation:
28 : C JSTAT=-1.........Operation unsuccessful
29 : C JSTAT=0..........Operation successful
30 : C
31 : C-----------------DUMMY ARGUMENTS---------------------------------------
32 :
33 : integer baseline, flag, bin, if_no, sourceno
34 : real weight(*), ut, u, v, w
35 : complex vis(*)
36 :
37 : C-----------------OTHER BITS & PIECES-----------------------------------
38 :
39 : include 'rpfits.inc'
40 :
41 : double precision D2PI
42 : parameter (D2PI = 2d0 * 3.14159265358979323846d0)
43 :
44 : logical async, open, scan_head
45 : integer AT_CLOSE, AT_CREATE, AT_REOPEN_WRITE, AT_WRITE, bufleft,
46 : + bufleft3, bufptr, cdelt3, crval3, grplength, grpptr, i,
47 : + i_buff(640), i_grphdr(11), icard, illegal, instat,
48 : + jstat, k, length, lun, n_words, newdim, pcount
49 : real buffer(640), crpix4, grphdr(11),
50 : + sc_buf(max_sc*max_if*ant_max)
51 : character key*8, utdate*12
52 :
53 : equivalence (i_buff(1), buffer(1))
54 : equivalence (i_grphdr(1), grphdr(1))
55 : equivalence (sc_buf(1), sc_cal(1,1,1))
56 :
57 : common /filelun/ lun
58 :
59 : save buffer, bufptr
60 :
61 : C-------------------TEMPLATE FOR RPFITS HEADERS-------------------------
62 :
63 : integer ACTUALDIM, MAX_HEADER
64 : parameter (ACTUALDIM = 70)
65 : parameter (MAX_HEADER = 2240)
66 :
67 : character m(MAX_HEADER)*80
68 : integer ichr(640), j
69 : character mout(32)*80
70 : real chr(640)
71 : equivalence (ichr, chr, mout)
72 :
73 : data illegal /32768/
74 : data bufptr /1/
75 : data async /.false./
76 : data length /2560/
77 : data open /.false./
78 : data pcount /11/
79 : data rpfitsversion /'2.12 '/
80 :
81 : data (m(i), i = 1,15) /
82 : + 'SIMPLE = F / NONCONFORMIST',
83 : + 'FORMAT = ''RPFITS'' / RPFITS',
84 : + 'SCANS = -1 / No. of scans in file',
85 : + 'BITPIX = -32 / Values are real',
86 : + 'NAXIS = 6 /',
87 : + 'NAXIS1 = 0 / Required for grouped data',
88 : + 'NAXIS2 = 3 / Complex=real,imag.,weight',
89 : + 'NAXIS3 = ? / No. of Stokes parameters',
90 : + 'NAXIS4 = ? / No. of frequencies',
91 : + 'NAXIS5 = 1 / Right Ascension (EPOCH)',
92 : + 'NAXIS6 = 1 / Declination (EPOCH)',
93 : + 'RPFITS = ''0.0 '' / RPFITS version number',
94 : + 'GROUPS = T / Data structured in groups',
95 : + 'PCOUNT = 11 / No. of random parameters',
96 : + 'GCOUNT = ? / No. of groups'/
97 :
98 : data (m(i), i = 16,25) /
99 : + 'BUNIT = ''JY '' / Unit of flux',
100 : + 'BLANK = -32768 / Value of blank pixel',
101 : + 'OBJECT = '' '' / Source Name',
102 : + 'INSTRUME= ''DUMMY '' / Instrument',
103 : + 'CAL = ''NONE '' / Cal applied',
104 : + 'EPOCH = ''J2000 '' / Epoch of RA & Dec.',
105 : + 'OBSERVER= '' '' /',
106 : + 'DATE = '' '' / UT Date data written',
107 : + 'DATE-OBS= '' '' / UT Date data generated',
108 : + 'HISTORY AT'' '' /'/
109 :
110 :
111 : data (m(i), i = 26,37) /
112 : + 'CTYPE2 = ''COMPLEX'' / Fringe visibilities',
113 : + 'CRPIX2 = 1.0 /',
114 : + 'CRVAL2 = 1.0 /',
115 : + 'CDELT2 = 1.0 /',
116 :
117 : + 'CTYPE3 = ''STOKES'' /',
118 : + 'CRPIX3 = 1 /',
119 : + 'CRVAL3 = 1 /',
120 : + 'CDELT3 = 1 /',
121 :
122 : + 'CTYPE4 = ''FREQ '' / Frequency in Hz',
123 : + 'CRPIX4 = ? / Ref. pixel= centre channel',
124 : + 'CRVAL4 = ? / Freq. (Hz) of ref. pixel',
125 : + 'CDELT4 = 1.0 / Frequency interval (Hz)'/
126 :
127 : data (m(i), i = 38,45) /
128 : + 'CTYPE5 = ''RA'' / Right Ascension in radians',
129 : + 'CRPIX5 = 1 /',
130 : + 'CRVAL5 = 0.0 /',
131 : + 'CDELT5 = 1.0 /',
132 :
133 : + 'CTYPE6 = ''DEC'' / Declination in radians',
134 : + 'CRPIX6 = 1 /',
135 : + 'CRVAL6 = ? /',
136 : + 'CDELT6 = 1.0 /'/
137 :
138 : data (m(i), i = 46,56) /
139 : + 'PTYPE1 = ''UU '' / U coordinate in metres',
140 : + 'PTYPE2 = ''VV '' / V coordinate in metres',
141 : + 'PTYPE3 = ''WW '' / W coordinate in metres',
142 : + 'PTYPE4 = ''BASELINE'' / =aerial(j)*256+aerial(i)',
143 : + 'PTYPE5 = ''UT '' / UT in seconds',
144 : + 'PTYPE6 = ''FLAG '' / Data flag',
145 : + 'PTYPE7 = ''BIN '' / Pulsar bin no.',
146 : + 'PTYPE8 = ''IF_NO'' / IF no (index to IF table)',
147 : + 'PTYPE9 = ''SOURCENO'' / SRC no (index to SU table)',
148 : + 'PTYPE10 = ''INTBAS'' / Integration time (secs)',
149 : + 'VELREF = 0 / Velocity reference system'/
150 :
151 : data (m(i), i = 57,66) /
152 : + 'RESTFREQ= 0.0 / Line rest frequency in Hz',
153 : + 'ALTRVAL = 0.0 / Alternate reference value',
154 : + 'ALTRPIX = 1.0 / Alternate reference pixel',
155 : + 'INTIME = 0 / Integration time in seconds',
156 : + 'HISTORY = ''? '' /',
157 : + 'TABLES = 1 / May be tables',
158 : + 'DEFEAT = 0 / Ephemeris defeat switch',
159 : + 'UTCMTAI = 0 / UTC-TAI in seconds',
160 : + 'DJMREFP = 0 / Param from USNO circular',
161 : + 'DJMREFT = 0 / Param from USNO circular'/
162 : data (m(i), i = 67,70) /
163 : + 'VERSION = '' '' / Set by calling program ',
164 : + 'PMRA = 0 / Proper motion RA (sect/day)',
165 : + 'PMDEC = 0 / Proper motion DEC (asec/day)',
166 : + 'PMEPOCH = 0 / Ref epoch for proper motion'/
167 :
168 : C---------------------DECIDE ON FUNCTION--------------------------------
169 :
170 0 : instat=jstat
171 0 : if (jstat.eq.-3) go to 100
172 0 : if (jstat.eq.-2) go to 100
173 0 : if (jstat.eq.-1) go to 100
174 0 : if (jstat.eq.0) go to 2000
175 0 : if (jstat.eq.1) go to 6000
176 0 : if (jstat.eq.2) go to 1000
177 0 : if (jstat.eq.3) go to 6000
178 : write (6, *)
179 0 : + ' Illegal value of jstat in routine RPFITSOUT is ', jstat
180 0 : RETURN
181 :
182 : C------------------------WRITE HEADER ----------------------------------
183 :
184 : C Only this first bit differs for file or scan headers.
185 0 : 100 if (jstat.eq.-1) then
186 0 : scan_head = .true.
187 0 : jstat = 0
188 :
189 :
190 : C If an incomplete buffer remains fill it with reserved operands
191 : C for easy detection then write it out.
192 0 : if (bufptr.gt.1) then
193 0 : do i = bufptr, 640
194 0 : call I4VAX (illegal, i_buff(i))
195 : end do
196 :
197 0 : rp_iostat = AT_WRITE (lun, buffer, length)
198 0 : if (rp_iostat.ne.0) then
199 0 : jstat = -1
200 0 : write (6, *) 'Cannot empty buffer'
201 0 : RETURN
202 : end if
203 : end if
204 : else
205 0 : if (open) then
206 0 : write (6, *) 'File is already open'
207 0 : jstat = -1
208 0 : RETURN
209 : end if
210 0 : scan_head = .false.
211 0 : jstat = 0
212 0 : rp_iostat = AT_CREATE (file, async, 0, lun)
213 0 : if (rp_iostat.ne.0) then
214 0 : jstat = -1
215 0 : write (6, *) 'Cannot open file'
216 0 : RETURN
217 : end if
218 :
219 0 : nscan = 0
220 : end if
221 :
222 0 : if (instat.eq.-3) RETURN
223 :
224 : C------------------- SET UP STOKES PARAMETERS --------------------------
225 :
226 : C First see how many words are to be written.
227 0 : if (data_format .eq. 1) then
228 0 : n_words = 1
229 0 : write_wt = .false.
230 0 : else if (data_format .eq. 2 ) then
231 0 : n_words = 2
232 0 : write_wt = .false.
233 0 : else if (data_format .eq. 3 ) then
234 0 : n_words = 3
235 0 : write_wt = .true.
236 : else
237 0 : write (6, *) 'rpfitsout:data_format must be 1, 2 or 3'
238 0 : jstat = -1
239 0 : RETURN
240 : end if
241 :
242 : c This section really needs to be re-done sometime to handle data
243 : c with non-standard stokes order: JER.
244 0 : if (feed_type(1,1)(1:1).eq.'I' .or. feed_type(1,1).eq.' ') then
245 0 : crval3 = 1
246 0 : cdelt3 = 1
247 0 : else if (feed_type(1,1).eq.'R' .or. feed_type(1,1).eq.'L') then
248 0 : crval3 = -1
249 0 : cdelt3 = -1
250 0 : else if (feed_type(1,1).eq.'X' .or. feed_type(1,1).eq.'Y') then
251 0 : crval3 = -5
252 0 : cdelt3 = -1
253 : else
254 0 : crval3 = 5
255 0 : cdelt3 = 1
256 : end if
257 :
258 :
259 : C------------------ SET UP HEADER PARAMETERS ---------------------------
260 :
261 : C Get UTC date.
262 0 : call datfit (' ', utdate, jstat)
263 :
264 : C Create current UTC date in the form yyyy-mm-dd
265 0 : if (datobs(1:2).eq.' ') then
266 0 : if (utdate.eq.' ') then
267 0 : jstat = -1
268 0 : datobs = ' '
269 0 : write (6, *) 'Failed to get current date for header.'
270 0 : RETURN
271 : else
272 0 : datobs = utdate
273 : endif
274 : end if
275 :
276 : C------------ SORT OUT OBSOLETE HEADER QUANTITIES IF NECESSARY ---------
277 0 : n_if = MAX (1, n_if)
278 0 : if (if_freq(1).eq.0.) if_freq(1) = freq
279 0 : if (if_bw(1).eq.0.) if_bw(1) = nfreq*dfreq
280 0 : if (if_nfreq(1).eq.0) if_nfreq(1) = nfreq
281 0 : if (if_nstok(1).eq.0) if_nstok(1) = nstok
282 0 : if (if_ref(1).eq.0) if_ref(1) = (nfreq+1)/2
283 0 : if (if_simul(1) .eq.0) if_simul(1) = 1
284 0 : if (if_chain(1) .eq.0) if_chain(1) = 1
285 0 : freq = if_freq(1)
286 0 : nfreq = if_nfreq(1)
287 0 : if (dfreq.eq.0.) dfreq = if_bw(1)/if_nfreq(1)
288 0 : nstok = if_nstok(1)
289 0 : crpix4 = if_ref(1)
290 :
291 0 : n_su = MAX (1, n_su)
292 0 : if (object.eq.' ') then
293 0 : object = su_name(1)
294 0 : else if (su_name(1).eq.' ') then
295 0 : su_name(1) = object
296 : end if
297 :
298 0 : if (ra.eq.0d0 .and. dec.eq.0d0) then
299 0 : ra = su_ra(1)
300 0 : dec = su_dec(1)
301 0 : else if (su_ra(1).eq.0d0 .and. su_dec(1).eq.0d0) then
302 0 : su_ra(1) = ra
303 0 : su_dec(1) = dec
304 : end if
305 :
306 : C-------------------- WRITE FILE HEADER --------------------------------
307 :
308 :
309 0 : do i = 1, ACTUALDIM
310 0 : key = m(i)(1:8)
311 0 : if (key.EQ.'SCANS ') then
312 0 : if (scan_head) then
313 0 : write (m(i)(11:30), '(i20)') 0
314 : else
315 0 : write (m(i)(11:30), '(i20)') -1
316 : end if
317 : end if
318 :
319 0 : if (key(1:5).eq.'NAXIS') then
320 0 : if (key.EQ.'NAXIS2') then
321 0 : write (m(i)(11:30), '(i20)') n_words
322 0 : else if (key.EQ.'NAXIS3') then
323 0 : write (m(i)(11:30), '(i20)') nstok
324 0 : else if (key.EQ.'NAXIS4') then
325 0 : write (m(i)(11:30), '(i20)') nfreq
326 : end if
327 0 : else if (key(1:5).eq.'CRVAL') then
328 0 : if (key.EQ.'CRVAL3') then
329 0 : write (m(i)(11:30), '(i20)') crval3
330 0 : else if (key.EQ.'CRVAL4') then
331 0 : write (m(i)(11:30), '(g20.12)') freq
332 0 : call RJUSTY(m(i)(11:30))
333 0 : else if (key.EQ.'CRVAL5') then
334 0 : if (ra.lt.0d0) ra = ra + D2PI
335 0 : write (m(i)(11:30), '(g20.12)') ra
336 0 : call RJUSTY(m(i)(11:30))
337 0 : else if (key.EQ.'CRVAL6') then
338 0 : write (m(i)(11:30), '(g20.12)') dec
339 0 : call RJUSTY(m(i)(11:30))
340 : end if
341 0 : else if (key(1:5).eq.'CDELT') then
342 0 : if (key.EQ.'CDELT3') then
343 0 : write (m(i)(11:30), '(i20)') cdelt3
344 0 : else if (key.EQ.'CDELT4') then
345 0 : write (m(i)(11:30), '(g20.12)') dfreq
346 0 : call RJUSTY(m(i)(11:30))
347 : end if
348 0 : else if (key.EQ.'CRPIX4') then
349 0 : write (m(i)(11:30), '(g20.5)') crpix4
350 0 : else if (key.EQ.'GCOUNT') then
351 0 : write (m(i)(11:30), '(i20)') ncount
352 0 : else if (key.EQ.'INTIME') then
353 0 : write (m(i)(11:30), '(i20)') intime
354 0 : else if (key.EQ.'BUNIT ') then
355 0 : call LJUSTY(bunit)
356 0 : write(m(i)(12:29),'(a16,2x)') bunit
357 0 : else if (key.EQ.'OBJECT ') then
358 0 : call LJUSTY(object)
359 0 : write(m(i)(12:29),'(a8,10x)') object
360 0 : else if (key.EQ.'INSTRUME') then
361 0 : call LJUSTY(instrument)
362 0 : write(m(i)(12:29),'(a16,2x)') instrument
363 0 : else if (key.EQ.'CAL ') then
364 0 : call LJUSTY(cal)
365 0 : write(m(i)(12:29),'(a16,2x)') cal
366 0 : else if (key.EQ.'OBSERVER') then
367 0 : call LJUSTY(rp_observer)
368 0 : write(m(i)(12:29),'(a16,2x)') rp_observer
369 0 : else if (key.EQ.'RPFITS ') then
370 0 : call LJUSTY(rpfitsversion)
371 0 : write(m(i)(12:31),'(a20)') rpfitsversion
372 0 : else if (key.EQ.'VERSION ') then
373 0 : call LJUSTY(version)
374 0 : write(m(i)(12:31),'(a20)') version
375 0 : else if (key.EQ.'DATE ') then
376 0 : m(i)(12:23) = utdate
377 0 : else if (key.EQ.'DATE-OBS') then
378 0 : m(i)(12:23) = datobs
379 0 : else if (m(i)(1:9).EQ.'HISTORY =') then
380 0 : m(i)(12:31) = 'RPFITSOUT ' // datobs(:10)
381 0 : else if (key.EQ.'VELREF ') then
382 0 : write (m(i)(11:30), '(i20)') ivelref
383 0 : else if (key.EQ.'RESTFREQ') then
384 0 : write (m(i)(11:30), '(g20.12)') rfreq
385 0 : call RJUSTY(m(i)(11:30))
386 0 : else if (key.EQ.'ALTRVAL ') then
387 0 : write (m(i)(11:30), '(g20.12)') vel1
388 0 : call RJUSTY(m(i)(11:30))
389 0 : else if (key.EQ.'EPOCH ') then
390 0 : m(i)(12:19) = coord
391 0 : else if (key.EQ.'DEFEAT ') then
392 0 : write (m(i)(11:30), '(i20)') rp_defeat
393 0 : else if (key.EQ.'UTCMTAI ') then
394 0 : write (m(i)(11:30), '(g20.12)') rp_utcmtai
395 0 : call RJUSTY(m(i)(11:30))
396 0 : else if (key.EQ.'DJMREFP ') then
397 0 : write (m(i)(11:30), '(g20.12)') rp_djmrefp
398 0 : call RJUSTY(m(i)(11:30))
399 0 : else if (key.EQ.'DJMREFT ') then
400 0 : write (m(i)(11:30), '(g20.12)') rp_djmreft
401 0 : call RJUSTY(m(i)(11:30))
402 0 : else if (key.EQ.'PMRA ') then
403 0 : write (m(i)(11:30), '(g20.12)') pm_ra
404 0 : call RJUSTY(m(i)(11:30))
405 0 : else if (key.EQ.'PMDEC ') then
406 0 : write (m(i)(11:30), '(g20.12)') pm_dec
407 0 : call RJUSTY(m(i)(11:30))
408 0 : else if (key.EQ.'PMEPOCH ') then
409 0 : write (m(i)(11:30), '(g20.12)') pm_epoch
410 0 : call RJUSTY(m(i)(11:30))
411 : end if
412 : end do
413 0 : i = ACTUALDIM
414 :
415 : C Write antenna cards (coordinates in metres), and met stuff
416 : C for PTI data.
417 0 : if ((Index(instrument,'PTI').gt.0) .and. (x(1).ne.0.0)) then
418 0 : do k = 1, nant
419 0 : i = 4*(k-1) + 1 + ACTUALDIM
420 0 : write (m(i), 900) k, sta(k)(1:3), x(k), y(k), z(k)
421 : 900 format ('ANTENNA N=',I1,1x,a3,
422 : + ' X=',g17.10,' Y=',g17.10,' Z=',g17.10)
423 0 : write (m(i+1), 910) 'PRESS', k, ' = ', rp_pressure(k)
424 0 : write (m(i+2), 910) 'TEMPE', k, ' = ', rp_temp(k)
425 0 : write (m(i+3), 910) 'HUMID', k, ' = ', rp_humid(k)
426 : 910 format (a5,i2,a3,g20.6)
427 : end do
428 0 : i = ACTUALDIM + 4*nant
429 : end if
430 :
431 : C Write ephemeris parameters.
432 0 : do k = 1, 12
433 0 : write (m(i+k), 910) 'EPHEM', k, ' = ', rp_c(k)
434 0 : call RJUSTY(m(i+k)(11:30))
435 : end do
436 0 : i = i + 12
437 :
438 : C Write user-defined cards.
439 0 : do icard = 1, ncard
440 0 : write (m(i+icard), '(a)') card(icard)
441 : end do
442 0 : i = i + ncard
443 :
444 : C Write out tables.
445 0 : if (n_if.ge.1) call WRITE_IF_TABLE (i, m)
446 0 : if (n_su.ge.1) call WRITE_SU_TABLE (i, m)
447 0 : if (n_fg.ge.1) call WRITE_FG_TABLE (i, m)
448 0 : if (n_mt.ge.1) call WRITE_MT_TABLE (i, m)
449 0 : if (n_cu.ge.1) call WRITE_CU_TABLE (i, m)
450 0 : if (nant.ge.1) call WRITE_AN_TABLE (i, m)
451 :
452 : C Write it all out.
453 0 : newdim = i + 1
454 0 : write (m(newdim), '(a)') 'END'
455 :
456 0 : do i = 1, (newdim-1)/32+1
457 0 : do j = 1, 32
458 0 : mout(j) = m(j + 32*(i-1))
459 : end do
460 : C 'ichr', 'chr' and 'mout' are equivalenced
461 : C read (mout, '(32(20a4,:,/))') (ichr(j), j=1,640)
462 0 : rp_iostat = AT_WRITE (lun, chr, length)
463 0 : if (rp_iostat.ne.0) then
464 0 : jstat = -1
465 0 : write (6, *) 'Cannot write header'
466 0 : RETURN
467 : end if
468 : end do
469 :
470 : C Tidy up ready for next time.
471 0 : do icard = ACTUALDIM + 1, newdim
472 0 : m(icard) = ' '
473 : end do
474 :
475 :
476 0 : jstat = 0
477 0 : bufptr = 1
478 0 : nscan = nscan + 1
479 0 : n_if = MAX (n_if, 1)
480 0 : RETURN
481 :
482 : C--------------------------- WRITE FG TABLE ----------------------------
483 :
484 : 1000 continue
485 :
486 : C Fill the remainder of the buffer with reserved operands for easy
487 : C detection.
488 0 : if (bufptr.gt.1) then
489 0 : do i = bufptr, 640
490 0 : call I4VAX (illegal, i_buff(i))
491 : end do
492 :
493 : C Flush buffer.
494 0 : rp_iostat = AT_WRITE (lun, buffer, length)
495 0 : if (rp_iostat.ne.0) then
496 0 : jstat = -1
497 0 : write (6, *) 'Cannot empty buffer'
498 0 : RETURN
499 : end if
500 : end if
501 0 : bufptr = 1
502 :
503 : C Write table into buffer.
504 0 : newdim = 0
505 0 : if (n_fg.ge.1) call write_fg_table (newdim, m)
506 :
507 : C Flush buffer.
508 0 : do i = 1, (newdim-1)/32+1
509 0 : do j = 1, 32
510 0 : mout(j) = m(j + 32*(i-1))
511 : end do
512 : C 'ichr', 'chr' and 'mout' are equivalenced
513 0 : rp_iostat = AT_WRITE (lun, chr, length)
514 0 : if (rp_iostat.ne.0) then
515 0 : jstat = -1
516 0 : write (6, *) 'Cannot write FG table'
517 0 : RETURN
518 : end if
519 : end do
520 :
521 0 : jstat = 0
522 0 : RETURN
523 :
524 : C------------------------WRITE DATA TO FITS FILE------------------------
525 :
526 : 2000 continue
527 :
528 : C THE FOLLOWING POINTERS AND COUNTERS ARE USED HERE:
529 : C GRPLENGTH No. of visibilities in group
530 : C GRPPTR Pointer to next visibility in group to be written
531 : C BUFPTR Pointer to next word to be written to current buffer
532 : C BUFLEFT No. of words still to be written to current buffer
533 :
534 : C-----------------------------------------------------------------------
535 : C Note: data are written in blocks of 5 records = 640(4byte) words
536 :
537 : C Determine grplength.
538 0 : if (baseline.eq.-1) then
539 0 : grplength = sc_ant*sc_if*sc_q
540 0 : else if (if_no.gt.1) then
541 0 : grplength = if_nfreq(if_no)*if_nstok(if_no)
542 : else
543 0 : grplength = nstok*nfreq
544 : end if
545 :
546 0 : grpptr = 1
547 :
548 :
549 : C WRITE PARAMETERS TO FILE, FORMAT OF RPFITS IS:
550 : C----------- VIS data ------------- ----------- SYSCAL data -------
551 : C (baseline > 0) (baseline = -1)
552 : C param 1=u in m 0.0
553 : C param 2=v in m 0.0
554 : C param 3=w in m 0.0
555 : C param 4=baseline number -1.0
556 : C param 5=UT in seconds sc_ut: UT in seconds
557 : C param 6= flag (if present) sc_ant
558 : C param 7= bin (if present) sc_if
559 : C param 8=if_no (if present) sc_q
560 : C param 9=sourceno (if present) sc_srcno
561 : C param 10=intbase integration time 0.0
562 : C param 11=data_format 0
563 :
564 0 : bufleft = 641-bufptr
565 0 : if (baseline.eq.-1) then
566 0 : if (bufleft.ge.pcount) then
567 0 : call R4VAX (0.0, buffer(bufptr))
568 0 : call R4VAX (0.0, buffer(bufptr+1))
569 0 : call R4VAX (0.0, buffer(bufptr+2))
570 0 : call R4VAX (FLOAT(baseline), buffer(bufptr+3))
571 0 : call R4VAX (sc_ut, buffer(bufptr+4))
572 0 : call I4VAX (sc_ant, i_buff(bufptr+5))
573 0 : call I4VAX (sc_if, i_buff(bufptr+6))
574 0 : call I4VAX (sc_q, i_buff(bufptr+7))
575 0 : call I4VAX (sc_srcno, i_buff(bufptr+8))
576 0 : call R4VAX (0.0, buffer(bufptr+9))
577 0 : call I4VAX (0, i_buff(bufptr+10))
578 0 : bufptr = bufptr + pcount
579 0 : ut = sc_ut
580 : else
581 0 : call R4VAX (0.0, grphdr(1))
582 0 : call R4VAX (0.0, grphdr(2))
583 0 : call R4VAX (0.0, grphdr(3))
584 0 : call R4VAX (FLOAT(baseline), grphdr(4))
585 0 : call R4VAX (sc_ut, grphdr(5))
586 0 : call I4VAX (sc_ant, i_grphdr(6))
587 0 : call I4VAX (sc_if, i_grphdr(7))
588 0 : call I4VAX (sc_q, i_grphdr(8))
589 0 : call I4VAX (sc_srcno, i_grphdr(9))
590 0 : call R4VAX (0.0, grphdr(10))
591 0 : call I4VAX (0, i_grphdr(11))
592 0 : ut = sc_ut
593 :
594 0 : do i = 1, bufleft
595 0 : i_buff(bufptr+i-1) = i_grphdr(i)
596 : end do
597 :
598 0 : rp_iostat = AT_WRITE (lun, buffer, length)
599 0 : if (rp_iostat.ne.0) then
600 0 : jstat = -1
601 0 : write (6, *) 'Cannot write data (1)'
602 0 : RETURN
603 : end if
604 0 : bufptr = pcount - bufleft
605 0 : do i = 1, bufptr
606 0 : i_buff(i) = i_grphdr(i+bufleft)
607 : end do
608 0 : bufptr = bufptr+1
609 : end if
610 : else
611 0 : if (bufleft.ge.pcount) then
612 0 : call R4VAX (u, buffer(bufptr))
613 0 : call R4VAX (v, buffer(bufptr+1))
614 0 : call R4VAX (w, buffer(bufptr+2))
615 0 : call R4VAX (FLOAT(baseline), buffer(bufptr+3))
616 0 : call R4VAX (ut, buffer(bufptr+4))
617 0 : call I4VAX (flag, i_buff(bufptr+5))
618 0 : call I4VAX (bin, i_buff(bufptr+6))
619 0 : call I4VAX (if_no, i_buff(bufptr+7))
620 0 : call I4VAX (sourceno, i_buff(bufptr+8))
621 0 : call R4VAX (intbase, buffer(bufptr+9))
622 0 : call I4VAX (data_format, i_buff(bufptr+10))
623 :
624 0 : bufptr = bufptr + pcount
625 : else
626 0 : call R4VAX (u, grphdr(1))
627 0 : call R4VAX (v, grphdr(2))
628 0 : call R4VAX (w, grphdr(3))
629 0 : call R4VAX (FLOAT(baseline), grphdr(4))
630 0 : call R4VAX (ut, grphdr(5))
631 0 : call I4VAX (flag, i_grphdr(6))
632 0 : call I4VAX (bin, i_grphdr(7))
633 0 : call I4VAX (if_no, i_grphdr(8))
634 0 : call I4VAX (sourceno, i_grphdr(9))
635 0 : call R4VAX (intbase, grphdr(10))
636 0 : call I4VAX (data_format, i_grphdr(11))
637 :
638 0 : do i = 1, bufleft
639 0 : i_buff(bufptr+i-1) = i_grphdr(i)
640 : end do
641 0 : rp_iostat = AT_WRITE (lun, buffer, length)
642 0 : if (rp_iostat.ne.0) then
643 0 : jstat = -1
644 0 : write (6, *) 'Cannot write data (1)'
645 0 : RETURN
646 : end if
647 0 : bufptr = pcount - bufleft
648 0 : do i = 1, bufptr
649 0 : i_buff(i) = i_grphdr(i+bufleft)
650 : end do
651 0 : bufptr = bufptr + 1
652 : end if
653 : end if
654 :
655 0 : if (bufptr.eq.0.or.bufptr.eq.641) then
656 0 : rp_iostat = AT_WRITE (lun, buffer, length)
657 0 : if (rp_iostat.ne.0) then
658 0 : jstat = -1
659 0 : write (6, *) 'Cannot write data (2)'
660 0 : RETURN
661 : end if
662 0 : bufptr = 1
663 : end if
664 :
665 : C ---------------------- WRITE VIS DATA TO FILE ------------------------
666 0 : if (baseline.eq.-1) go to 4000
667 : C FORMAT FROM RPFITS IS:
668 : C data_format 3 2 1
669 : C word 1 = Re(vis) Re(vis) Real(vis)
670 : C word 2 = Imag(vis) Imag(vis) -
671 : C word 3 = weight - -
672 :
673 : 3500 continue
674 : C First see how many words are to be written.
675 0 : if (data_format .eq. 1) then
676 0 : n_words = 1
677 0 : write_wt = .false.
678 0 : else if (data_format .eq. 2 ) then
679 0 : n_words = 2
680 0 : write_wt = .false.
681 0 : else if (data_format .eq. 3 ) then
682 0 : n_words = 3
683 0 : write_wt = .true.
684 : else
685 0 : write (6, *) 'rpfitsout:data_format must be 1, 2 or 3'
686 0 : jstat = -1
687 0 : RETURN
688 : end if
689 :
690 : C type *,'data_format is ',data_format
691 : C type *,'n_words is ',n_words
692 : C type *,'write_wt is ',write_wt
693 :
694 : C If entire group can be put in existing buffer then do so.
695 0 : bufleft = 641 - bufptr
696 0 : if (bufleft.ge.(n_words*(grplength-grpptr+1))) then
697 0 : do i = grpptr, grplength
698 0 : if (n_words.eq.1) then
699 0 : call R4VAX (REAL(vis(i)), buffer(bufptr))
700 : else
701 0 : call R4VAX (REAL(vis(i)), buffer(bufptr))
702 0 : call R4VAX (AIMAG(vis(i)), buffer(bufptr+1))
703 : end if
704 0 : if (write_wt) call R4VAX (weight(i), buffer(bufptr+2))
705 0 : bufptr = bufptr+n_words
706 : end do
707 0 : jstat = 0
708 0 : RETURN
709 :
710 : else
711 : C Otherwise things are a bit more complicated, first write
712 : C complete visibilities to old buffer.
713 0 : bufleft3 = bufleft/n_words
714 0 : do i = 1, bufleft3
715 0 : call R4VAX (REAL(vis(grpptr+i-1)), buffer(bufptr))
716 0 : if (n_words.gt.1)
717 0 : + call R4VAX (AIMAG(vis(grpptr+i-1)), buffer(bufptr+1))
718 :
719 0 : if (write_wt)
720 0 : + call R4VAX (weight(grpptr+i-1), buffer(bufptr+2))
721 0 : bufptr = bufptr+n_words
722 : end do
723 0 : grpptr = grpptr+bufleft3
724 :
725 : C Now write the fraction of a visibility left into old buffer.
726 0 : bufleft = bufleft-n_words*bufleft3
727 :
728 0 : if (bufleft.eq.1) then
729 0 : call R4VAX (REAL(vis(grpptr)), buffer(640))
730 0 : else if (bufleft.eq.2 .and. n_words.eq.3) then
731 0 : call R4VAX (REAL(vis(grpptr)), buffer(639))
732 0 : call R4VAX (AIMAG(vis(grpptr)), buffer(640))
733 : end if
734 :
735 : C Start a new buffer.
736 0 : rp_iostat = AT_WRITE (lun, buffer, length)
737 0 : if (rp_iostat.ne.0) then
738 0 : jstat = -1
739 0 : write (6, *) 'Cannot write data (3)'
740 0 : RETURN
741 : end if
742 :
743 : C Fill any incomplete visibility, and reset BUFPTR.
744 0 : if (bufleft.eq.0) then
745 0 : bufptr = 1
746 0 : else if (bufleft.eq.1) then
747 0 : call R4VAX (AIMAG(vis(grpptr)), buffer(1))
748 0 : if (write_wt) call R4VAX (weight(grpptr), buffer(2))
749 0 : grpptr = grpptr+1
750 0 : bufptr = n_words
751 0 : else if (bufleft.eq.2 .and. write_wt) then
752 0 : call R4VAX (weight(grpptr), buffer(1))
753 0 : grpptr = grpptr+1
754 0 : bufptr = n_words-1
755 : end if
756 : end if
757 :
758 : C Loop to write out the rest of the group.
759 0 : go to 3500
760 :
761 :
762 : C---------------------- WRITE SYSCAL DATA TO FILE ----------------------
763 : 4000 continue
764 :
765 0 : bufleft = 641 - bufptr
766 0 : if (bufleft.ge.(grplength-grpptr+1)) then
767 : C Easy, the group fits into buffer.
768 0 : do i = grpptr, grplength
769 0 : call R4VAX (sc_buf(i), buffer(bufptr))
770 0 : bufptr = bufptr + 1
771 : end do
772 :
773 0 : jstat = 0
774 0 : RETURN
775 :
776 : else
777 : C Fill the buffer.
778 0 : do i = 1, bufleft
779 0 : call R4VAX (sc_buf(grpptr+i-1), buffer(bufptr))
780 0 : bufptr = bufptr + 1
781 : end do
782 :
783 0 : grpptr = grpptr + bufleft
784 :
785 : C Write the buffer out.
786 0 : rp_iostat = AT_WRITE (lun, buffer, length)
787 0 : if (rp_iostat.ne.0) then
788 0 : jstat = -1
789 0 : write (6, *) 'Cannot write data (3)'
790 0 : RETURN
791 : end if
792 :
793 : C Start a new buffer.
794 0 : bufptr = 1
795 : end if
796 :
797 : C Return to write out the rest of the group.
798 0 : go to 4000
799 :
800 :
801 : C--------------------------- CLOSE FILE --------------------------------
802 :
803 : 6000 continue
804 :
805 0 : if (bufptr.gt.1) then
806 : C Fill buffer with reserved operands.
807 0 : do i = bufptr, 640
808 0 : call I4VAX (illegal, i_buff(i))
809 : end do
810 :
811 : C Write out buffer.
812 0 : rp_iostat = AT_WRITE (lun, buffer, length)
813 0 : if (rp_iostat.ne.0) then
814 0 : jstat = -1
815 0 : write (6, *) 'Cannot empty buffer'
816 0 : RETURN
817 : end if
818 : end if
819 :
820 : C Reset buffer pointer.
821 0 : bufptr = 1
822 :
823 : C Close file.
824 0 : rp_iostat = AT_CLOSE (lun)
825 0 : if (rp_iostat.ne.0) then
826 0 : jstat = -1
827 0 : write (6, *) 'Cannot close file'
828 0 : RETURN
829 : end if
830 :
831 0 : if (jstat.eq.1) then
832 0 : jstat = 0
833 0 : open = .false.
834 : else
835 : C REOPEN file.
836 0 : rp_iostat = AT_REOPEN_WRITE (file, lun)
837 0 : if (rp_iostat.ne.0) then
838 0 : jstat = -1
839 0 : write (6, *) 'Cannot reopen file'
840 : else
841 0 : jstat = 0
842 : end if
843 : end if
844 :
845 0 : RETURN
846 :
847 : end
|
