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: fwproj.f 17791 2004-08-25 02:28:46Z cvsmgr $
28 : *-----------------------------------------------------------------------
29 : C
30 : C Grid a number of visibility records
31 : C
32 0 : subroutine gwgrid (uvw, dphase, values, nvispol, nvischan,
33 0 : $ dopsf, flag, rflag, weight, nrow, rownum,
34 0 : $ scale, offset, grid, nx, ny, npol, nchan, freq, c,
35 0 : $ support, convsize, sampling, wconvsize, convfunc,
36 0 : $ chanmap, polmap,
37 0 : $ sumwt)
38 :
39 : implicit none
40 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
41 : complex values(nvispol, nvischan, nrow)
42 : double complex grid(nx, ny, npol, nchan)
43 : double precision uvw(3, nrow), freq(nvischan), c, scale(3),
44 : $ offset(3)
45 : double precision dphase(nrow), uvdist
46 : complex phasor
47 : integer flag(nvispol, nvischan, nrow)
48 : integer rflag(nrow)
49 : real weight(nvischan, nrow), phase
50 : double precision sumwt(npol, nchan)
51 : integer rownum
52 : integer support(*), rsupport
53 : integer chanmap(nchan), polmap(npol)
54 : integer dopsf
55 :
56 : double complex nvalue
57 :
58 : integer convsize, sampling, wconvsize
59 : complex convfunc(convsize/2-1, convsize/2-1, wconvsize),
60 : $ cwt
61 :
62 : real norm
63 : real wt
64 :
65 : logical owp
66 :
67 : double precision pos(3)
68 : integer loc(3), off(3), iloc(3)
69 : integer rbeg, rend
70 : integer ix, iy, ipol, ichan
71 : integer apol, achan, irow
72 : double precision pi
73 : data pi/3.14159265358979323846/
74 :
75 0 : irow=rownum
76 :
77 0 : if(irow.ge.0) then
78 0 : rbeg=irow+1
79 0 : rend=irow+1
80 : else
81 0 : rbeg=1
82 0 : rend=nrow
83 : end if
84 :
85 0 : do irow=rbeg, rend
86 0 : if(rflag(irow).eq.0) then
87 0 : do ichan=1, nvischan
88 0 : achan=chanmap(ichan)+1
89 0 : if((achan.ge.1).and.(achan.le.nchan).and.
90 0 : $ (weight(ichan,irow).ne.0.0)) then
91 : call swp(uvw(1,irow), dphase(irow), freq(ichan), c,
92 0 : $ scale, offset, sampling, pos, loc, off, phasor)
93 0 : iloc(3)=max(1, min(wconvsize, loc(3)))
94 0 : rsupport=support(iloc(3))
95 0 : if (owp(nx, ny, wconvsize, loc, rsupport)) then
96 0 : do ipol=1, nvispol
97 0 : apol=polmap(ipol)+1
98 : if((flag(ipol,ichan,irow).ne.1).and.
99 0 : $ (apol.ge.1).and.(apol.le.npol)) then
100 : C If we are making a PSF then we don't want to phase
101 : C rotate but we do want to reproject uvw
102 0 : if(dopsf.eq.1) then
103 0 : nvalue=cmplx(weight(ichan,irow))
104 : else
105 : nvalue=weight(ichan,irow)*
106 0 : $ (values(ipol,ichan,irow)*phasor)
107 : end if
108 : C norm will be the value we would get for the peak
109 : C at the phase center. We will want to normalize
110 : C the final image by this term.
111 0 : norm=0.0
112 0 : do iy=-rsupport,rsupport
113 0 : iloc(2)=1+abs(iy*sampling+off(2))
114 : C !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(ix)
115 : C !SOMP DO
116 0 : do ix=-rsupport,rsupport
117 0 : iloc(1)=1+abs(ix*sampling+off(1))
118 0 : if(uvw(3,irow).gt.0.0) then
119 : cwt=conjg(convfunc(iloc(1),
120 0 : $ iloc(2), iloc(3)))
121 : else
122 : cwt=convfunc(iloc(1),
123 0 : $ iloc(2), iloc(3))
124 : end if
125 : grid(loc(1)+ix,
126 : $ loc(2)+iy,apol,achan)=
127 : $ grid(loc(1)+ix,
128 : $ loc(2)+iy,apol,achan)+
129 0 : $ nvalue*cwt
130 0 : norm=norm+real(cwt)
131 : end do
132 : C !$OMP END DO
133 : C !$OMP END PARALLEL
134 : end do
135 : sumwt(apol,achan)=sumwt(apol,achan)+
136 0 : $ weight(ichan,irow)*norm
137 : end if
138 : end do
139 : else
140 : C write(*,*) uvw(3,irow), pos(1), pos(2), pos(3),
141 : C $ loc(1), loc(2), loc(3)
142 : end if
143 : end if
144 : end do
145 : end if
146 : end do
147 0 : return
148 : end
149 :
150 :
151 201544 : subroutine sectgwgridd (uvw, values, nvispol, nvischan,
152 201544 : $ dopsf, flag, rflag, weight, nrow,
153 201544 : $ grid, nx, ny, npol, nchan,
154 201544 : $ support, convsize, sampling, wconvsize, convfunc,
155 201544 : $ chanmap, polmap,
156 201544 : $ sumwt, x0,
157 201544 : $ y0, nxsub, nysub, rbeg, rend, loc, off, phasor)
158 :
159 : implicit none
160 :
161 : integer, intent(in) :: nx,ny, npol, nchan, nvispol, nvischan, nrow
162 : double precision, intent(in) :: uvw(3, nrow)
163 : complex, intent(in) :: values(nvispol, nvischan, nrow)
164 : double complex, intent(inout) :: grid(nx, ny, npol, nchan)
165 : integer, intent(in) :: x0, y0, nxsub, nysub
166 : complex, intent(in) :: phasor(nvischan, nrow)
167 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
168 : integer, intent(in) :: rflag(nrow)
169 : real, intent(in) :: weight(nvischan, nrow)
170 : double precision, intent(inout) :: sumwt(npol, nchan)
171 : integer, intent(in) :: support(*)
172 : integer :: rsupport
173 : integer, intent(in) :: chanmap(nchan), polmap(npol)
174 : integer, intent(in) :: dopsf
175 :
176 : double complex :: nvalue
177 :
178 : integer, intent(in) :: convsize, sampling, wconvsize
179 : complex, intent(in) :: convfunc(convsize/2-1,convsize/2-1,
180 : $ wconvsize)
181 : complex :: cwt
182 :
183 : real :: norm
184 : real :: wt
185 :
186 : logical :: onmywgrid
187 :
188 : integer, intent(in) :: loc(3, nvischan, nrow)
189 : integer, intent(in) :: off(3, nvischan, nrow)
190 : integer :: iloc(3)
191 : integer, intent(in) :: rbeg, rend
192 : integer :: ix, iy, ipol, ichan
193 : integer :: apol, achan, irow
194 : integer :: posx, posy, msupportx, msupporty, psupportx, psupporty
195 : C complex :: cfunc(convsize/2-1, convsize/2-1)
196 : C integer :: npoint
197 :
198 70943488 : do irow=rbeg, rend
199 70943488 : if(rflag(irow).eq.0) then
200 196483056 : do ichan=1, nvischan
201 129226216 : achan=chanmap(ichan)+1
202 129226216 : if((achan.ge.1).and.(achan.le.nchan).and.
203 67256840 : $ (weight(ichan,irow).ne.0.0)) then
204 90943976 : iloc(3)=max(1, min(wconvsize, loc(3, ichan, irow)))
205 90943976 : rsupport=support(iloc(3))
206 : C write(*,*) irow, ichan, iloc(3), rsupport
207 90943976 : if (onmywgrid(loc(1, ichan, irow), nx, ny, wconvsize,
208 : $ x0, y0, nxsub, nysub, rsupport, msupportx, msupporty
209 : $ , psupportx, psupporty)) then
210 : CCC I thought removing the if in the inner loop will be faster
211 : CCC but not really as i have to calculate more conjg at this stage
212 : C npoint=rsupport*sampling+1
213 : C if(uvw(3,irow).gt.0.0) then
214 : C cfunc(1:npoint, 1:npoint)=conjg(
215 : C $ convfunc(1:npoint,1:npoint,iloc(3)))
216 : C else
217 : C cfunc(1:npoint, 1:npoint)=
218 : C $ convfunc(1:npoint,1:npoint,iloc(3))
219 : C end if
220 : CCCC
221 12070848 : do ipol=1, nvispol
222 8047232 : apol=polmap(ipol)+1
223 : if((flag(ipol,ichan,irow).ne.1).and.
224 12070848 : $ (apol.ge.1).and.(apol.le.npol)) then
225 : C If we are making a PSF then we don't want to phase
226 : C rotate but we do want to reproject uvw
227 8047232 : if(dopsf.eq.1) then
228 3147572 : nvalue=cmplx(weight(ichan,irow))
229 : else
230 : nvalue=weight(ichan,irow)*
231 4899660 : $ (values(ipol,ichan,irow)*phasor(ichan,irow))
232 : end if
233 : C norm will be the value we would get for the peak
234 : C at the phase center. We will want to normalize
235 : C the final image by this term.
236 8047232 : norm=0.0
237 : C msupporty=-rsupport
238 : C psupporty=rsupport
239 : C psupportx=rsupport
240 : C msupportx=-rsupport
241 : C if((loc(1, ichan, irow)-rsupport) .lt. x0)
242 : C $ msupportx= -(loc(1, ichan, irow)-x0)
243 : C if((loc(1, ichan, irow)+rsupport).ge.(x0+nxsub))
244 : C $ psupportx= x0+nxsub-loc(1, ichan, irow)-1
245 : C if((loc(2, ichan, irow)-rsupport) .lt. y0)
246 : C $ msupporty= -(loc(2, ichan, irow)-y0)
247 : C if((loc(2, ichan, irow)+rsupport).ge.(y0+nysub))
248 : C $ psupporty= y0+nysub-loc(2, ichan, irow)-1
249 75397748 : do iy=msupporty,psupporty
250 67350516 : posy=loc(2, ichan, irow)+iy
251 67350516 : iloc(2)=1+abs(iy*sampling+off(2, ichan,irow))
252 729683614 : do ix=msupportx,psupportx
253 654285866 : posx=loc(1, ichan, irow)+ix
254 : iloc(1)=1+abs(ix*sampling+
255 654285866 : $ off(1,ichan,irow))
256 : cwt=convfunc(iloc(1),
257 654285866 : $ iloc(2), iloc(3))
258 654285866 : if(uvw(3,irow).gt.0.0) cwt=conjg(cwt)
259 : C else
260 : C cwt=convfunc(iloc(1),
261 : C $ iloc(2), iloc(3))
262 : C end if
263 : grid(posx,posy,apol,achan)=
264 : $ grid(posx,posy,apol,achan)+
265 654285866 : $ nvalue*cwt
266 721636382 : norm=norm+real(cwt)
267 : end do
268 : end do
269 : sumwt(apol,achan)=sumwt(apol,achan)+
270 8047232 : $ weight(ichan,irow)*norm
271 : end if
272 : end do
273 : else
274 : C write(*,*) uvw(3,irow), pos(1), pos(2), pos(3),
275 : C $ loc(1), loc(2), loc(3)
276 : end if
277 : end if
278 : end do
279 : end if
280 : end do
281 201544 : return
282 : end
283 : C
284 : C single precision gridder
285 12800 : subroutine sectgwgrids (uvw, values, nvispol, nvischan,
286 12800 : $ dopsf, flag, rflag, weight, nrow,
287 12800 : $ grid, nx, ny, npol, nchan,
288 12800 : $ support, convsize, sampling, wconvsize, convfunc,
289 12800 : $ chanmap, polmap,
290 12800 : $ sumwt, x0,
291 12800 : $ y0, nxsub, nysub, rbeg, rend, loc, off, phasor)
292 :
293 : implicit none
294 :
295 : integer, intent(in) :: nx,ny, npol, nchan, nvispol, nvischan, nrow
296 : double precision, intent(in) :: uvw(3, nrow)
297 : complex, intent(in) :: values(nvispol, nvischan, nrow)
298 : complex, intent(inout) :: grid(nx, ny, npol, nchan)
299 : integer, intent(in) :: x0, y0, nxsub, nysub
300 : complex, intent(in) :: phasor(nvischan, nrow)
301 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
302 : integer, intent(in) :: rflag(nrow)
303 : real, intent(in) :: weight(nvischan, nrow)
304 : double precision, intent(inout) :: sumwt(npol, nchan)
305 : integer, intent(in) :: support(*)
306 : integer :: rsupport
307 : integer, intent(in) :: chanmap(nchan), polmap(npol)
308 : integer, intent(in) :: dopsf
309 :
310 : double complex :: nvalue
311 :
312 : integer, intent(in) :: convsize, sampling, wconvsize
313 : complex, intent(in) :: convfunc(convsize/2-1,convsize/2-1,
314 : $ wconvsize)
315 : complex :: cwt
316 :
317 : real :: norm
318 : real :: wt
319 :
320 : logical :: onmywgrid
321 :
322 : integer, intent(in) :: loc(3, nvischan, nrow)
323 : integer, intent(in) :: off(3, nvischan, nrow)
324 : integer :: iloc(3)
325 : integer, intent(in) :: rbeg, rend
326 : integer :: ix, iy, ipol, ichan
327 : integer :: apol, achan, irow
328 : integer :: posx, posy, msupportx, msupporty, psupportx, psupporty
329 : C complex :: cfunc(convsize/2-1, convsize/2-1)
330 : C integer :: npoint
331 :
332 4505600 : do irow=rbeg, rend
333 4505600 : if(rflag(irow).eq.0) then
334 453772800 : do ichan=1, nvischan
335 449280000 : achan=chanmap(ichan)+1
336 449280000 : if((achan.ge.1).and.(achan.le.nchan).and.
337 4492800 : $ (weight(ichan,irow).ne.0.0)) then
338 4492800 : iloc(3)=max(1, min(wconvsize, loc(3, ichan, irow)))
339 4492800 : rsupport=support(iloc(3))
340 4492800 : if (onmywgrid(loc(1, ichan, irow), nx, ny, wconvsize,
341 : $ x0, y0, nxsub, nysub, rsupport, msupportx,
342 : $ msupporty, psupportx, psupporty)) then
343 : CCC I thought removing the if in the inner loop will be faster
344 : CCC but not really as i have to calculate more conjg at this stage
345 : C npoint=rsupport*sampling+1
346 : C if(uvw(3,irow).gt.0.0) then
347 : C cfunc(1:npoint, 1:npoint)=conjg(
348 : C $ convfunc(1:npoint,1:npoint,iloc(3)))
349 : C else
350 : C cfunc(1:npoint, 1:npoint)=
351 : C $ convfunc(1:npoint,1:npoint,iloc(3))
352 : C end if
353 : CCCC
354 2090520 : do ipol=1, nvispol
355 1672416 : apol=polmap(ipol)+1
356 : if((flag(ipol,ichan,irow).ne.1).and.
357 2090520 : $ (apol.ge.1).and.(apol.le.npol)) then
358 : C If we are making a PSF then we don't want to phase
359 : C rotate but we do want to reproject uvw
360 836208 : if(dopsf.eq.1) then
361 0 : nvalue=cmplx(weight(ichan,irow))
362 : else
363 : nvalue=weight(ichan,irow)*
364 836208 : $ (values(ipol,ichan,irow)*phasor(ichan,irow))
365 : end if
366 : C norm will be the value we would get for the peak
367 : C at the phase center. We will want to normalize
368 : C the final image by this term.
369 836208 : norm=0.0
370 : C msupporty=-rsupport
371 : C psupporty=rsupport
372 : C psupportx=rsupport
373 : C msupportx=-rsupport
374 : C if((loc(1, ichan, irow)-rsupport) .lt. x0)
375 : C $ msupportx= -(loc(1, ichan, irow)-x0)
376 : C if((loc(1, ichan, irow)+rsupport).ge.(x0+nxsub))
377 : C $ psupportx= x0+nxsub-loc(1, ichan, irow)-1
378 : C if((loc(2, ichan, irow)-rsupport) .lt. y0)
379 : C $ msupporty= -(loc(2, ichan, irow)-y0)
380 : C if((loc(2, ichan, irow)+rsupport).ge.(y0+nysub))
381 : C $ psupporty= y0+nysub-loc(2, ichan, irow)-1
382 3834342 : do iy=msupporty,psupporty
383 2998134 : posy=loc(2, ichan, irow)+iy
384 2998134 : iloc(2)=1+abs(iy*sampling+off(2, ichan,irow))
385 15206742 : do ix=msupportx,psupportx
386 11372400 : posx=loc(1, ichan, irow)+ix
387 : iloc(1)=1+abs(ix*sampling+
388 11372400 : $ off(1,ichan,irow))
389 : cwt=convfunc(iloc(1),
390 11372400 : $ iloc(2), iloc(3))
391 11372400 : if(uvw(3,irow).gt.0.0) cwt=conjg(cwt)
392 : C else
393 : C cwt=convfunc(iloc(1),
394 : C $ iloc(2), iloc(3))
395 : C end if
396 : grid(posx,posy,apol,achan)=
397 : $ grid(posx,posy,apol,achan)+
398 11372400 : $ nvalue*cwt
399 14370534 : norm=norm+real(cwt)
400 : end do
401 : end do
402 : sumwt(apol,achan)=sumwt(apol,achan)+
403 836208 : $ weight(ichan,irow)*norm
404 : end if
405 : end do
406 : else
407 : C write(*,*) uvw(3,irow), pos(1), pos(2), pos(3),
408 : C $ loc(1), loc(2), loc(3)
409 : end if
410 : end if
411 : end do
412 : end if
413 : end do
414 12800 : return
415 : end
416 :
417 :
418 : C
419 : C Degrid a number of visibility records
420 : C
421 0 : subroutine dwgrid (uvw, dphase, values, nvispol, nvischan,
422 0 : $ flag, rflag,
423 0 : $ nrow, rownum, scale, offset, grid, nx, ny, npol, nchan, freq,
424 0 : $ c, support, convsize, sampling, wconvsize, convfunc,
425 : $ chanmap, polmap)
426 :
427 : implicit none
428 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
429 : complex values(nvispol, nvischan, nrow)
430 : complex grid(nx, ny, npol, nchan)
431 : double precision uvw(3, nrow), freq(nvischan), c, scale(3),
432 : $ offset(3)
433 : double precision dphase(nrow), uvdist
434 : complex phasor
435 : integer flag(nvispol, nvischan, nrow)
436 : integer rflag(nrow)
437 : integer rownum
438 : integer support(*), rsupport
439 : integer chanmap(*), polmap(*)
440 :
441 : complex nvalue
442 :
443 : integer convsize, wconvsize, sampling
444 : complex convfunc(convsize/2-1, convsize/2-1, wconvsize),
445 : $ cwt
446 :
447 : real norm, phase
448 :
449 : logical owp
450 :
451 : double precision pos(3)
452 : integer loc(3), off(3), iloc(3)
453 : integer rbeg, rend
454 : integer ix, iy, ipol, ichan
455 : integer apol, achan, irow
456 : real wt, wtx, wty
457 : double precision pi
458 : data pi/3.14159265358979323846/
459 :
460 0 : irow=rownum
461 :
462 0 : if(irow.ge.0) then
463 0 : rbeg=irow+1
464 0 : rend=irow+1
465 : else
466 0 : rbeg=1
467 0 : rend=nrow
468 : end if
469 : C
470 :
471 0 : do irow=rbeg, rend
472 0 : if(rflag(irow).eq.0) then
473 0 : do ichan=1, nvischan
474 0 : achan=chanmap(ichan)+1
475 0 : if((achan.ge.1).and.(achan.le.nchan)) then
476 : call swp(uvw(1,irow), dphase(irow), freq(ichan), c,
477 0 : $ scale, offset, sampling, pos, loc, off, phasor)
478 0 : iloc(3)=max(1, min(wconvsize, loc(3)))
479 0 : rsupport=support(iloc(3))
480 0 : if (owp(nx, ny, wconvsize, loc, rsupport)) then
481 0 : do ipol=1, nvispol
482 0 : apol=polmap(ipol)+1
483 : if((flag(ipol,ichan,irow).ne.1).and.
484 0 : $ (apol.ge.1).and.(apol.le.npol)) then
485 :
486 0 : nvalue=0.0
487 0 : do iy=-rsupport,rsupport
488 0 : iloc(2)=1+abs(iy*sampling+off(2))
489 0 : do ix=-rsupport,rsupport
490 0 : iloc(1)=1+abs(ix*sampling+off(1))
491 0 : if(uvw(3,irow).gt.0.0) then
492 : cwt=conjg(convfunc(iloc(1),
493 0 : $ iloc(2), iloc(3)))
494 : else
495 : cwt=convfunc(iloc(1),
496 0 : $ iloc(2), iloc(3))
497 : end if
498 : nvalue=nvalue+conjg(cwt)*
499 0 : $ grid(loc(1)+ix,loc(2)+iy,apol,achan)
500 : end do
501 : end do
502 0 : values(ipol,ichan,irow)=nvalue*conjg(phasor)
503 : end if
504 : end do
505 : end if
506 : end if
507 : end do
508 : end if
509 : end do
510 0 : return
511 : end
512 : C
513 :
514 14896 : subroutine sectdwgrid (uvw, values, nvispol, nvischan,
515 14896 : $ flag, rflag,
516 14896 : $ nrow, grid, nx, ny, npol, nchan,
517 14896 : $ support, convsize, sampling, wconvsize, convfunc,
518 14896 : $ chanmap, polmap, rbeg, rend, loc,off,phasor)
519 :
520 : implicit none
521 : integer, intent(in) :: nrow,nx,ny,npol, nchan, nvispol, nvischan
522 : complex, intent(inout) :: values(nvispol, nvischan, nrow)
523 : complex, intent(in) :: grid(nx, ny, npol, nchan)
524 : double precision, intent(in) :: uvw(3, nrow)
525 : complex , intent(in) :: phasor(nvischan, nrow)
526 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
527 : integer, intent(in) :: rflag(nrow)
528 : integer, intent(in) :: support(*), sampling
529 : integer, intent(in) :: chanmap(*), polmap(*)
530 : integer, intent(in) :: rbeg, rend
531 : integer, intent(in) :: loc(3, nvischan, nrow)
532 : integer, intent(in) :: off(3, nvischan, nrow)
533 : complex :: nvalue
534 :
535 : integer, intent(in) :: convsize, wconvsize
536 : complex, intent(in) :: convfunc(convsize/2-1, convsize/2-1,
537 : $ wconvsize)
538 : complex :: cwt
539 :
540 : real :: norm, phase
541 :
542 : logical :: owp
543 :
544 : integer :: iloc(3), rsupport
545 : integer :: ix, iy, ipol, ichan
546 : integer :: apol, achan, irow
547 : real :: wt, wtx, wty
548 :
549 :
550 622126 : do irow=rbeg, rend
551 622126 : if(rflag(irow).eq.0) then
552 8535876 : do ichan=1, nvischan
553 7952194 : achan=chanmap(ichan)+1
554 8535876 : if((achan.ge.1).and.(achan.le.nchan)) then
555 643498 : iloc(3)=max(1, min(wconvsize, loc(3, ichan, irow)))
556 643498 : rsupport=support(iloc(3))
557 643498 : if (owp(nx, ny, wconvsize, loc(1,ichan,irow), rsupport))
558 : $ then
559 2067981 : do ipol=1, nvispol
560 1425454 : apol=polmap(ipol)+1
561 : if((flag(ipol,ichan,irow).ne.1).and.
562 2067981 : $ (apol.ge.1).and.(apol.le.npol)) then
563 :
564 1285054 : nvalue=0.0
565 15976868 : do iy=-rsupport,rsupport
566 14691814 : iloc(2)=1+abs(iy*sampling+off(2,ichan,irow))
567 195129378 : do ix=-rsupport,rsupport
568 : iloc(1)=1+abs(ix*sampling+
569 179152510 : $ off(1,ichan,irow))
570 :
571 179152510 : if(uvw(3,irow).gt.0.0) then
572 : cwt=conjg(convfunc(iloc(1),
573 88636362 : $ iloc(2), iloc(3)))
574 : else
575 : cwt=convfunc(iloc(1),
576 90516148 : $ iloc(2), iloc(3))
577 : end if
578 : nvalue=nvalue+conjg(cwt)*
579 : $ grid(loc(1,ichan,irow)+ix,loc(2,ichan,
580 193844324 : $ irow)+iy,apol,achan)
581 : end do
582 : end do
583 : values(ipol,ichan,irow)=nvalue*conjg(
584 1285054 : $ phasor(ichan, irow))
585 : end if
586 : end do
587 : end if
588 : end if
589 : end do
590 : end if
591 : end do
592 14896 : return
593 : end
594 : C
595 :
596 :
597 :
598 : C Calculate gridded coordinates and the phasor needed for
599 : C phase rotation.
600 : C
601 0 : subroutine swp (uvw, dphase, freq, c, scale, offset,
602 : $ sampling, pos, loc, off, phasor)
603 : implicit none
604 : integer loc(3), off(3), sampling
605 : double precision uvw(3), freq, c, scale(3), offset(3)
606 : double precision pos(3)
607 : double precision dphase, phase
608 : complex phasor
609 : integer idim
610 : double precision pi
611 : data pi/3.14159265358979323846/
612 :
613 : C pos(3)=(scale(3)*uvw(3)*freq/c)*(scale(3)*uvw(3)*freq/c)
614 : C $ +offset(3)+1.0;
615 : C pos(3)=(scale(3)*uvw(3)*freq/c)+offset(3)+1.0;
616 0 : pos(3)=sqrt(abs(scale(3)*uvw(3)*freq/c))+offset(3)+1.0
617 0 : loc(3)=nint(pos(3))
618 0 : off(3)=0
619 :
620 0 : do idim=1,2
621 : pos(idim)=scale(idim)*uvw(idim)*freq/c+
622 0 : $ (offset(idim)+1.0)
623 0 : loc(idim)=nint(pos(idim))
624 0 : off(idim)=nint((loc(idim)-pos(idim))*sampling)
625 : end do
626 :
627 0 : phase=-2.0D0*pi*dphase*freq/c
628 0 : phasor=cmplx(cos(phase), sin(phase))
629 0 : return
630 : end
631 643498 : logical function owp (nx, ny, nw, loc, support)
632 : implicit none
633 : integer nx, ny, nw, loc(3), support
634 : owp=(support.gt.0).and.
635 : $ (loc(1)-support.ge.1).and.(loc(1)+support.le.nx).and.
636 : $ (loc(2)-support.ge.1).and.(loc(2)+support.le.ny).and.
637 643498 : $ (loc(3).ge.1).and.(loc(3).le.nw)
638 643498 : return
639 : end
640 95436776 : logical function onmywgrid(loc, nx, ny, nw, nx0, ny0, nxsub,nysub,
641 : $ support, msuppx, msuppy, psuppx, psuppy)
642 : implicit none
643 : integer, intent(in) :: nx, ny, nw, nx0, ny0, nxsub, nysub, loc(3),
644 : $ support
645 : integer, intent(out) :: msuppx, msuppy, psuppx, psuppy
646 : integer :: loc1sub, loc1plus, loc2sub, loc2plus
647 95436776 : msuppx=merge(-support, nx0-loc(1), loc(1)-support >= nx0)
648 95436776 : msuppy=merge(-support, ny0-loc(2), loc(2)-support >= ny0)
649 : psuppx=merge(support, nx0+nxsub-loc(1)-1 , (loc(1)+support)
650 95436776 : $ < (nx0+nxsub))
651 : psuppy=merge(support, ny0+nysub-loc(2)-1 , (loc(2)+support)
652 95436776 : $ < (ny0+nysub))
653 95436776 : loc1sub=loc(1)-support
654 95436776 : loc1plus=loc(1)+support
655 95436776 : loc2sub=loc(2)-support
656 95436776 : loc2plus=loc(2)+support
657 :
658 : onmywgrid=(support.gt.0).and. (loc(3).ge.1) .and. (loc(3).le.nw)
659 : $ .and. (loc1plus .ge. nx0) .and. (loc1sub
660 : $ .le. (nx0+nxsub)) .and.(loc2plus .ge. ny0) .and.
661 95436776 : $ (loc2sub .le. (ny0+nysub))
662 95436776 : return
663 : end
|
