Line data Source code
1 : *=======================================================================
2 : * Copyright (C) 1999,2001,2002
3 : * Associated Universities, Inc. Washington DC, USA.
4 : *
5 : * This library is free software; you can redistribute it and/or
6 : * modify it under the terms of the GNU Library General Public
7 : * License as published by the Free Software Foundation; either
8 : * version 2 of the License, or (at your option) any later version.
9 : *
10 : * This library is distributed in the hope that it will be useful,
11 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 : * GNU Library General Public License for more details.
14 : *
15 : * You should have received a copy of the GNU Library General Public
16 : * License along with this library; if not, write to the Free
17 : * Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,
18 : * MA 02139, USA.
19 : *
20 : * Correspondence concerning AIPS++ should be addressed as follows:
21 : * Internet email: aips2-request@nrao.edu.
22 : * Postal address: AIPS++ Project Office
23 : * National Radio Astronomy Observatory
24 : * 520 Edgemont Road
25 : * Charlottesville, VA 22903-2475 USA
26 : *
27 : * $Id$
28 : *-----------------------------------------------------------------------
29 : C
30 : C Calculate gridded coordinates
31 : C
32 50 : subroutine sgridsdclip (xy, sampling, pos, loc, off)
33 : implicit none
34 : integer sampling
35 : integer loc(2), off(2)
36 : double precision xy(2)
37 : real pos(2)
38 : integer idim
39 :
40 150 : do idim=1,2
41 100 : pos(idim)=xy(idim)+1.0
42 100 : loc(idim)=nint(pos(idim))
43 150 : off(idim)=nint((loc(idim)-pos(idim))*sampling)
44 : end do
45 50 : return
46 : end
47 : C
48 : C Is this on the grid?
49 : C
50 50 : logical function ogridsdclip (nx, ny, loc, support)
51 : implicit none
52 : integer nx, ny, loc(2), support
53 : ogridsdclip=(loc(1)-support.ge.1).and.(loc(1)+support.le.nx).and.
54 50 : $ (loc(2)-support.ge.1).and.(loc(2)+support.le.ny)
55 50 : return
56 : end
57 : C
58 : C Grid a number of visibility records: single dish gridding
59 : C but with complex images including additional process for
60 : C min/max clipping
61 : C
62 : C This subroutine assumes specific clipping process after
63 : C all the data are accumulated
64 : C
65 : C Post-accumulation process is as follows. On each pixel,
66 : C
67 : C - if npoints == 1, grid = gmin, wgrid = wmin, sumwt must
68 : C be updated accordingly
69 : C - if npoints == 2, grid = gmin + gmax, wgrid = wmin + wmax,
70 : C sumwt must be updated accordingly
71 : C - otherwise, leave grid, wgrid, and sumwt as they are
72 : C
73 111 : subroutine ggridsdclip (xy, values, nvispol, nvischan,
74 111 : $ dowt, flag, rflag, weight, nrow, irow,
75 37 : $ grid, wgrid,
76 37 : $ npoints, gmin, wmin, gmax, wmax,
77 : $ nx, ny, npol, nchan,
78 111 : $ support, sampling, convFunc, chanmap, polmap, sumwt)
79 :
80 : implicit none
81 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
82 : complex values(nvispol, nvischan, nrow)
83 : complex grid(nx, ny, npol, nchan)
84 : real wgrid(nx, ny, npol, nchan)
85 : integer npoints(nx, ny, npol, nchan)
86 : complex gmin(nx, ny, npol, nchan)
87 : complex gmax(nx, ny, npol, nchan)
88 : real wmin(nx, ny, npol, nchan)
89 : real wmax(nx, ny, npol, nchan)
90 : double precision xy(2,nrow)
91 : integer flag(nvispol, nvischan, nrow)
92 : integer rflag(nrow)
93 : real weight(nvischan, nrow)
94 : double precision sumwt(npol, nchan)
95 : integer irow
96 : integer support, sampling
97 : integer chanmap(nvischan), polmap(nvispol)
98 : integer dowt
99 :
100 : complex nvalue, thevalue
101 :
102 : real convFunc(*)
103 : real norm
104 : real wt, wtx, wty
105 : real swap, theweight, nweight
106 :
107 : logical ogridsdclip
108 :
109 : real pos(2), rloc(2)
110 : integer loc(2), off(2)
111 : integer rbeg, rend
112 74 : integer irad((2*support+1)**2)
113 : integer ix, iy, ipol, ichan
114 : integer apol, achan
115 : integer ir
116 37 : integer xloc(2*support+1), yloc(2*support+1)
117 : integer ax, ay
118 :
119 37 : if(irow.ge.0) then
120 0 : rbeg=irow+1
121 0 : rend=irow+1
122 : else
123 37 : rbeg=1
124 37 : rend=nrow
125 : end if
126 :
127 87 : do irow=rbeg, rend
128 87 : if(rflag(irow).eq.0) then
129 50 : call sgridsdclip(xy(1,irow), sampling, pos, loc, off)
130 : C if (ogridsdclip(nx, ny, loc, support)) then
131 50 : if (ogridsdclip(nx, ny, loc, 0)) then
132 50 : ir=1
133 50 : norm=-(support+1)*sampling+off(1)
134 50 : rloc(2)=-(support+1)*sampling+off(2)
135 100 : do iy=1,2*support+1
136 50 : rloc(2)=rloc(2)+sampling
137 50 : rloc(1)=norm
138 150 : do ix=1,2*support+1
139 50 : rloc(1)=rloc(1)+sampling
140 50 : irad(ir)=sqrt(rloc(1)**2+rloc(2)**2)+1
141 100 : ir=ir+1
142 : end do
143 : end do
144 50 : xloc(1)=loc(1)-support
145 50 : do ix=2,2*support+1
146 50 : xloc(ix)=xloc(ix-1)+1
147 : end do
148 50 : yloc(1)=loc(2)-support
149 50 : do iy=2,2*support+1
150 50 : yloc(iy)=yloc(iy-1)+1
151 : end do
152 106 : do ichan=1, nvischan
153 56 : achan=chanmap(ichan)+1
154 56 : if((achan.ge.1).and.(achan.le.nchan).and.
155 50 : $ (weight(ichan,irow).gt.0.0)) then
156 112 : do ipol=1, nvispol
157 56 : apol=polmap(ipol)+1
158 : if((flag(ipol,ichan,irow).ne.1).and.
159 112 : $ (apol.ge.1).and.(apol.le.npol)) then
160 56 : if (dowt.eq.1) then
161 0 : thevalue=1.0
162 : else
163 56 : thevalue=conjg(values(ipol,ichan,irow))
164 : end if
165 56 : theweight=weight(ichan,irow)
166 56 : norm=0.0
167 56 : ir=1
168 : C do iy=-support,support
169 112 : do iy=1,2*support+1
170 56 : ay=yloc(iy)
171 112 : if ((ay.ge.1).and.(ay.le.ny)) then
172 : C do ix=-support,support
173 112 : do ix=1,2*support+1
174 56 : ax=xloc(ix)
175 112 : if ((ax.ge.1).and.(ax.le.nx)) then
176 56 : ir = (iy-1)*(2*support+1) + ix
177 56 : wt=convFunc(irad(ir))
178 56 : nweight=theweight*wt
179 56 : nvalue=thevalue*nweight
180 56 : if (npoints(ax,ay,apol,achan).eq.0) then
181 22 : gmin(ax,ay,apol,achan)=thevalue
182 22 : wmin(ax,ay,apol,achan)=nweight
183 : else if
184 34 : $ (npoints(ax,ay,apol,achan).eq.1) then
185 18 : if (real(gmin(ax,ay,apol,achan)).lt.
186 : $ real(thevalue)) then
187 4 : gmax(ax,ay,apol,achan)=thevalue
188 4 : wmax(ax,ay,apol,achan)=nweight
189 : else
190 : gmax(ax,ay,apol,achan)=
191 14 : $ gmin(ax,ay,apol,achan)
192 : wmax(ax,ay,apol,achan)=
193 14 : $ wmin(ax,ay,apol,achan)
194 14 : gmin(ax,ay,apol,achan)=thevalue
195 14 : wmin(ax,ay,apol,achan)=nweight
196 : end if
197 : else
198 16 : if (real(thevalue).le.
199 : $ real(gmin(ax,ay,apol,achan))) then
200 : nvalue=wmin(ax,ay,apol,achan)*
201 14 : $ gmin(ax,ay,apol,achan)
202 14 : gmin(ax,ay,apol,achan)=thevalue
203 14 : swap=nweight
204 14 : nweight=wmin(ax,ay,apol,achan)
205 14 : wmin(ax,ay,apol,achan)=swap
206 2 : else if (real(thevalue).ge.
207 : $ real(gmax(ax,ay,apol,achan))) then
208 : nvalue=wmax(ax,ay,apol,achan)*
209 0 : $ gmax(ax,ay,apol,achan)
210 0 : gmax(ax,ay,apol,achan)=thevalue
211 0 : swap=nweight
212 0 : nweight=wmax(ax,ay,apol,achan)
213 0 : wmax(ax,ay,apol,achan)=swap
214 : end if
215 : grid(ax,ay,apol,achan)=
216 16 : $ grid(ax,ay,apol,achan)+nvalue
217 : wgrid(ax,ay,apol,achan)=
218 16 : $ wgrid(ax,ay,apol,achan)+nweight
219 16 : norm=norm+nweight
220 : end if
221 : npoints(ax,ay,apol,achan)=
222 56 : $ npoints(ax,ay,apol,achan)+1
223 : end if
224 : end do
225 : end if
226 : end do
227 56 : sumwt(apol,achan)=sumwt(apol,achan)+norm
228 : end if
229 : end do
230 : end if
231 : end do
232 : end if
233 : end if
234 : end do
235 37 : return
236 : end
|
