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: casa-feedback@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 0 : 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 0 : do idim=1,2
41 0 : pos(idim)=xy(idim)+1.0
42 0 : loc(idim)=nint(pos(idim))
43 0 : off(idim)=nint((loc(idim)-pos(idim))*sampling)
44 : end do
45 0 : return
46 : end
47 : C
48 : C Is this on the grid?
49 : C
50 0 : 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 0 : $ (loc(2)-support.ge.1).and.(loc(2)+support.le.ny)
55 0 : 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 0 : subroutine ggridsdclip (xy, values, nvispol, nvischan,
74 0 : $ dowt, flag, rflag, weight, nrow, irow,
75 0 : $ grid, wgrid,
76 0 : $ npoints, gmin, wmin, gmax, wmax,
77 : $ nx, ny, npol, nchan,
78 0 : $ 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 0 : integer irad((2*support+1)**2)
113 : integer ix, iy, ipol, ichan
114 : integer apol, achan
115 : integer ir
116 0 : integer xloc(2*support+1), yloc(2*support+1)
117 : integer ax, ay
118 :
119 0 : if(irow.ge.0) then
120 0 : rbeg=irow+1
121 0 : rend=irow+1
122 : else
123 0 : rbeg=1
124 0 : rend=nrow
125 : end if
126 :
127 0 : do irow=rbeg, rend
128 0 : if(rflag(irow).eq.0) then
129 0 : call sgridsdclip(xy(1,irow), sampling, pos, loc, off)
130 : C if (ogridsdclip(nx, ny, loc, support)) then
131 0 : if (ogridsdclip(nx, ny, loc, 0)) then
132 0 : ir=1
133 0 : norm=-(support+1)*sampling+off(1)
134 0 : rloc(2)=-(support+1)*sampling+off(2)
135 0 : do iy=1,2*support+1
136 0 : rloc(2)=rloc(2)+sampling
137 0 : rloc(1)=norm
138 0 : do ix=1,2*support+1
139 0 : rloc(1)=rloc(1)+sampling
140 0 : irad(ir)=sqrt(rloc(1)**2+rloc(2)**2)+1
141 0 : ir=ir+1
142 : end do
143 : end do
144 0 : xloc(1)=loc(1)-support
145 0 : do ix=2,2*support+1
146 0 : xloc(ix)=xloc(ix-1)+1
147 : end do
148 0 : yloc(1)=loc(2)-support
149 0 : do iy=2,2*support+1
150 0 : yloc(iy)=yloc(iy-1)+1
151 : end do
152 0 : do ichan=1, nvischan
153 0 : achan=chanmap(ichan)+1
154 0 : if((achan.ge.1).and.(achan.le.nchan).and.
155 0 : $ (weight(ichan,irow).gt.0.0)) then
156 0 : do ipol=1, nvispol
157 0 : apol=polmap(ipol)+1
158 : if((flag(ipol,ichan,irow).ne.1).and.
159 0 : $ (apol.ge.1).and.(apol.le.npol)) then
160 0 : if (dowt.eq.1) then
161 0 : thevalue=1.0
162 : else
163 0 : thevalue=conjg(values(ipol,ichan,irow))
164 : end if
165 0 : theweight=weight(ichan,irow)
166 0 : norm=0.0
167 0 : ir=1
168 : C do iy=-support,support
169 0 : do iy=1,2*support+1
170 0 : ay=yloc(iy)
171 0 : if ((ay.ge.1).and.(ay.le.ny)) then
172 : C do ix=-support,support
173 0 : do ix=1,2*support+1
174 0 : ax=xloc(ix)
175 0 : if ((ax.ge.1).and.(ax.le.nx)) then
176 0 : ir = (iy-1)*(2*support+1) + ix
177 0 : wt=convFunc(irad(ir))
178 0 : nweight=theweight*wt
179 0 : nvalue=thevalue*nweight
180 0 : if (npoints(ax,ay,apol,achan).eq.0) then
181 0 : gmin(ax,ay,apol,achan)=thevalue
182 0 : wmin(ax,ay,apol,achan)=nweight
183 : else if
184 0 : $ (npoints(ax,ay,apol,achan).eq.1) then
185 0 : if (real(gmin(ax,ay,apol,achan)).lt.
186 : $ real(thevalue)) then
187 0 : gmax(ax,ay,apol,achan)=thevalue
188 0 : wmax(ax,ay,apol,achan)=nweight
189 : else
190 : gmax(ax,ay,apol,achan)=
191 0 : $ gmin(ax,ay,apol,achan)
192 : wmax(ax,ay,apol,achan)=
193 0 : $ wmin(ax,ay,apol,achan)
194 0 : gmin(ax,ay,apol,achan)=thevalue
195 0 : wmin(ax,ay,apol,achan)=nweight
196 : end if
197 : else
198 0 : if (real(thevalue).le.
199 : $ real(gmin(ax,ay,apol,achan))) then
200 : nvalue=wmin(ax,ay,apol,achan)*
201 0 : $ gmin(ax,ay,apol,achan)
202 0 : gmin(ax,ay,apol,achan)=thevalue
203 0 : swap=nweight
204 0 : nweight=wmin(ax,ay,apol,achan)
205 0 : wmin(ax,ay,apol,achan)=swap
206 0 : 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 0 : $ grid(ax,ay,apol,achan)+nvalue
217 : wgrid(ax,ay,apol,achan)=
218 0 : $ wgrid(ax,ay,apol,achan)+nweight
219 0 : norm=norm+nweight
220 : end if
221 : npoints(ax,ay,apol,achan)=
222 0 : $ npoints(ax,ay,apol,achan)+1
223 : end if
224 : end do
225 : end if
226 : end do
227 0 : 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 0 : return
236 : end
|
