Line data Source code
1 : *=======================================================================
2 : * Copyright (C) 1999-2012
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: fgridft.f 19685 2006-10-05 20:57:29Z rurvashi $
28 : *-----------------------------------------------------------------------
29 : C
30 : C Grid a number of visibility records
31 : C
32 0 : subroutine ggrid (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, sampling, convFunc, chanmap, polmap, sumwt)
36 :
37 : implicit none
38 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
39 : complex values(nvispol, nvischan, nrow)
40 : double complex grid(nx, ny, npol, nchan)
41 : double precision uvw(3, nrow), freq(nvischan), c, scale(2),
42 : $ offset(2)
43 : double precision dphase(nrow), uvdist
44 : complex phasor
45 : integer flag(nvispol, nvischan, nrow)
46 : integer rflag(nrow)
47 : real weight(nvischan, nrow)
48 : double precision sumwt(npol, nchan)
49 : integer rownum
50 : integer support, sampling
51 : integer chanmap(nchan), polmap(npol)
52 : integer dopsf
53 :
54 : double complex nvalue
55 :
56 : double precision convFunc(*)
57 : double precision norm
58 : double precision wt, wtx, wty
59 :
60 : logical ogrid
61 :
62 : double precision pos(2)
63 : integer loc(2), off(2), iloc(2)
64 : integer rbeg, rend
65 : integer ix, iy, ipol, ichan
66 : integer apol, achan, irow
67 :
68 0 : irow=rownum
69 :
70 0 : if(irow.ge.0) then
71 0 : rbeg=irow+1
72 0 : rend=irow+1
73 : else
74 0 : rbeg=1
75 0 : rend=nrow
76 : end if
77 :
78 0 : do irow=rbeg, rend
79 0 : if(rflag(irow).eq.0) then
80 0 : do ichan=1, nvischan
81 0 : achan=chanmap(ichan)+1
82 0 : if((achan.ge.1).and.(achan.le.nchan).and.
83 0 : $ (weight(ichan,irow).ne.0.0)) then
84 : call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c,
85 0 : $ scale, offset, sampling, pos, loc, off, phasor)
86 0 : if (ogrid(nx, ny, loc, support)) then
87 0 : do ipol=1, nvispol
88 0 : apol=polmap(ipol)+1
89 : if((flag(ipol,ichan,irow).ne.1).and.
90 0 : $ (apol.ge.1).and.(apol.le.npol)) then
91 : C If we are making a PSF then we don't want to phase
92 : C rotate but we do want to reproject uvw
93 0 : if(dopsf.eq.1) then
94 0 : nvalue=cmplx(weight(ichan,irow))
95 : else
96 : nvalue=weight(ichan,irow)*
97 0 : $ (values(ipol,ichan,irow)*phasor)
98 : end if
99 0 : norm=0.0
100 0 : do iy=-support,support
101 0 : iloc(2)=abs(sampling*iy+off(2))+1
102 0 : wty=convFunc(iloc(2))
103 0 : do ix=-support,support
104 0 : iloc(1)=abs(sampling*ix+off(1))+1
105 0 : wtx=convFunc(iloc(1))
106 0 : wt=wtx*wty
107 : grid(loc(1)+ix,loc(2)+iy,apol,achan)=
108 : $ grid(loc(1)+ix,loc(2)+iy,apol,achan)+
109 0 : $ nvalue*wt
110 0 : norm=norm+wt
111 : end do
112 : end do
113 : sumwt(apol,achan)=sumwt(apol,achan)+
114 0 : $ weight(ichan,irow)*norm
115 : end if
116 : end do
117 : end if
118 : end if
119 : end do
120 : end if
121 : end do
122 0 : return
123 : end
124 : C
125 : C Single precision gridder
126 0 : subroutine ggrids (uvw, dphase, values, nvispol, nvischan,
127 0 : $ dopsf, flag, rflag, weight, nrow, rownum,
128 0 : $ scale, offset, grid, nx, ny, npol, nchan, freq, c,
129 0 : $ support, sampling, convFunc, chanmap, polmap, sumwt)
130 :
131 : implicit none
132 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
133 : complex values(nvispol, nvischan, nrow)
134 : complex grid(nx, ny, npol, nchan)
135 : double precision uvw(3, nrow), freq(nvischan), c, scale(2),
136 : $ offset(2)
137 : double precision dphase(nrow), uvdist
138 : complex phasor
139 : integer flag(nvispol, nvischan, nrow)
140 : integer rflag(nrow)
141 : real weight(nvischan, nrow)
142 : double precision sumwt(npol, nchan)
143 : integer rownum
144 : integer support, sampling
145 : integer chanmap(nchan), polmap(npol)
146 : integer dopsf
147 :
148 : double complex nvalue
149 :
150 : double precision convFunc(*)
151 : double precision norm
152 : double precision wt, wtx, wty
153 :
154 : logical ogrid
155 :
156 : double precision pos(2)
157 : integer loc(2), off(2), iloc(2)
158 : integer rbeg, rend
159 : integer ix, iy, ipol, ichan
160 : integer apol, achan, irow
161 :
162 0 : irow=rownum
163 :
164 0 : if(irow.ge.0) then
165 0 : rbeg=irow+1
166 0 : rend=irow+1
167 : else
168 0 : rbeg=1
169 0 : rend=nrow
170 : end if
171 :
172 0 : do irow=rbeg, rend
173 0 : if(rflag(irow).eq.0) then
174 0 : do ichan=1, nvischan
175 0 : achan=chanmap(ichan)+1
176 0 : if((achan.ge.1).and.(achan.le.nchan).and.
177 0 : $ (weight(ichan,irow).ne.0.0)) then
178 : call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c,
179 0 : $ scale, offset, sampling, pos, loc, off, phasor)
180 0 : if (ogrid(nx, ny, loc, support)) then
181 0 : do ipol=1, nvispol
182 0 : apol=polmap(ipol)+1
183 : if((flag(ipol,ichan,irow).ne.1).and.
184 0 : $ (apol.ge.1).and.(apol.le.npol)) then
185 : C If we are making a PSF then we don't want to phase
186 : C rotate but we do want to reproject uvw
187 0 : if(dopsf.eq.1) then
188 0 : nvalue=cmplx(weight(ichan,irow))
189 : else
190 : nvalue=weight(ichan,irow)*
191 0 : $ (values(ipol,ichan,irow)*phasor)
192 : end if
193 0 : norm=0.0
194 0 : do iy=-support,support
195 0 : iloc(2)=abs(sampling*iy+off(2))+1
196 0 : wty=convFunc(iloc(2))
197 0 : do ix=-support,support
198 0 : iloc(1)=abs(sampling*ix+off(1))+1
199 0 : wtx=convFunc(iloc(1))
200 0 : wt=wtx*wty
201 : grid(loc(1)+ix,loc(2)+iy,apol,achan)=
202 : $ grid(loc(1)+ix,loc(2)+iy,apol,achan)+
203 0 : $ nvalue*wt
204 0 : norm=norm+wt
205 : end do
206 : end do
207 : sumwt(apol,achan)=sumwt(apol,achan)+
208 0 : $ weight(ichan,irow)*norm
209 : end if
210 : end do
211 : end if
212 : end if
213 : end do
214 : end if
215 : end do
216 0 : return
217 : end
218 :
219 :
220 1052 : subroutine sectggridd(values, nvispol, nvischan,
221 1052 : $ dopsf, flag, rflag, weight, nrow,
222 1052 : $ grid, nx, ny, npol, nchan,
223 1052 : $ support, sampling, convFunc, chanmap, polmap, sumwt, x0,
224 1052 : $ y0, nxsub, nysub, rbeg, rend, loc, off, phasor)
225 : implicit none
226 :
227 : integer, intent(in) :: nx, ny, npol, nchan, nvispol, nvischan,nrow
228 : complex, intent(in) :: values(nvispol, nvischan, nrow)
229 : double complex, intent(inout) :: grid(nx, ny, npol, nchan)
230 : integer, intent(in) :: x0, y0, nxsub, nysub
231 : double precision, intent(in) :: convFunc(*)
232 : integer, intent(in) :: chanmap(nvischan), polmap(nvispol)
233 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
234 : integer, intent(in) :: rflag(nrow)
235 : real, intent(in) :: weight(nvischan, nrow)
236 : double precision, intent(inout) :: sumwt(npol, nchan)
237 : integer, intent(in) :: support, sampling
238 : integer, intent(in) :: dopsf, rbeg, rend
239 :
240 :
241 : integer, intent(in) :: loc(2, nvischan, nrow)
242 : integer, intent(in) :: off(2, nvischan, nrow)
243 : integer :: iloc(2)
244 : complex, intent(in) :: phasor(nvischan, nrow)
245 : double complex :: nvalue
246 :
247 : double precision :: norm
248 : double precision :: wt, wtx, wty
249 :
250 : logical :: onmygrid
251 :
252 : double precision :: pos(2)
253 : integer :: ix, iy, ipol, ichan
254 : integer :: apol, achan, irow
255 : integer :: posx, posy
256 : integer :: msupporty, psupporty
257 : integer :: msupportx, psupportx
258 1052 : double precision :: cfnow(-support:support, -support:support)
259 :
260 :
261 999112 : do irow=rbeg, rend
262 999112 : if(rflag(irow).eq.0) then
263 8737820 : do ichan=1, nvischan
264 7739760 : achan=chanmap(ichan)+1
265 7739760 : if((achan.ge.1).and.(achan.le.nchan).and.
266 998060 : $ (weight(ichan,irow).ne.0.0)) then
267 : C call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c,
268 : C $ scale, offset, sampling, pos, loc, off, phasor)
269 : C write(*,*) loc(1, ichan, irow), loc(2, ichan, irow), irow,ichan
270 7739760 : if (onmygrid(loc(1,ichan, irow), nx, ny, x0, y0, nxsub,
271 : $ nysub, support)) then
272 17346897 : do ipol=1, nvispol
273 11564598 : apol=polmap(ipol)+1
274 : if((flag(ipol,ichan,irow).ne.1).and.
275 17346897 : $ (apol.ge.1).and.(apol.le.npol)) then
276 : C If we are making a PSF then we don't want to phase
277 : C rotate but we do want to reproject uvw
278 11564598 : if(dopsf.eq.1) then
279 11564598 : nvalue=cmplx(weight(ichan,irow))
280 : else
281 : nvalue=weight(ichan,irow)*
282 0 : $ (values(ipol,ichan,irow)*phasor(ichan, irow))
283 : end if
284 11564598 : norm=0.0
285 11564598 : msupporty=-support
286 11564598 : psupporty=support
287 11564598 : psupportx=support
288 11564598 : msupportx=-support
289 92516784 : do iy=-support, support
290 80952186 : iloc(2)=abs(sampling*iy+off(2,ichan,irow))+1
291 80952186 : wty=convFunc(iloc(2))
292 659182086 : do ix=-support, support
293 : iloc(1)=abs(sampling*ix+off(1,ichan,
294 566665302 : $ irow))+1
295 566665302 : wtx=convFunc(iloc(1))
296 566665302 : cfnow(ix, iy)=wtx*wty
297 647617488 : norm=norm+cfnow(ix,iy)
298 : end do
299 : end do
300 92516784 : do iy=-support, support
301 659182086 : do ix=-support, support
302 647617488 : cfnow(ix, iy)=cfnow(ix, iy)/norm
303 : end do
304 : end do
305 11564598 : if((loc(1, ichan, irow)-support) .lt. x0)
306 4276932 : $ msupportx= -(loc(1, ichan, irow)-x0)
307 11564598 : if((loc(1, ichan, irow)+support).ge.(x0+nxsub))
308 4977952 : $ psupportx= x0+nxsub-loc(1, ichan, irow)-1
309 11564598 : if((loc(2, ichan, irow)-support) .lt. y0)
310 4486870 : $ msupporty= -(loc(2, ichan, irow)-y0)
311 11564598 : if((loc(2, ichan, irow)+support).ge.(y0+nysub))
312 4745432 : $ psupporty= y0+nysub-loc(2, ichan, irow)-1
313 : C write(*,*) msupportx, psupportx, msupporty, psupporty
314 11564598 : norm=0.0
315 59298760 : do iy=msupporty,psupporty
316 47734162 : posy=loc(2, ichan, irow)+iy
317 : C if( (posy .lt. (y0+nysub)) .and.
318 : C $ (posy.ge. y0)) then
319 : C iloc(2)=abs(sampling*iy+off(2,ichan,irow))+1
320 : C wty=convFunc(iloc(2))
321 248922880 : do ix=msupportx,psupportx
322 189624120 : posx=loc(1, ichan, irow)+ix
323 : C if( (posx .lt. (x0+nxsub)) .and.
324 : C $ (posx .ge. x0)) then
325 : C write(*,*) posx, posy, loc(1), loc(2), x0, y0, nxsub, nysub
326 : C iloc(1)=abs(sampling*ix+off(1,ichan,
327 : C $ irow))+1
328 : C wtx=convFunc(iloc(1))
329 : C wt=wtx*wty
330 : grid(posx,posy,apol,achan)=
331 : $ grid(posx, posy,apol,achan)+
332 189624120 : $ nvalue*cfnow(ix,iy)
333 237358282 : norm=norm+cfnow(ix,iy)
334 : C write(*,*) iloc(1), iloc(2), posx, posy
335 : C end if
336 : end do
337 : C end if
338 : end do
339 : sumwt(apol,achan)=sumwt(apol,achan)+
340 11564598 : $ weight(ichan,irow)*norm
341 : end if
342 : end do
343 : end if
344 : C if onmygrid
345 : end if
346 : end do
347 : end if
348 : end do
349 1052 : return
350 : end
351 :
352 : C subroutine sectggrids(uvw, dphase, values, nvispol, nvischan,
353 : C $ dopsf, flag, rflag, weight, nrow,
354 : C $ scale, offset, grid, nx, ny, npol, nchan, freq, c,
355 : C $ support, sampling, convFunc, chanmap, polmap, sumwt, x0,
356 : C $ y0, nxsub, nysub, rbeg, rend)
357 :
358 26400 : subroutine sectggrids(values, nvispol, nvischan,
359 26400 : $ dopsf, flag, rflag, weight, nrow,
360 26400 : $ grid, nx, ny, npol, nchan,
361 26400 : $ support, sampling, convFunc, chanmap, polmap, sumwt, x0,
362 26400 : $ y0, nxsub, nysub, rbeg, rend, loc, off, phasor)
363 :
364 :
365 :
366 : implicit none
367 :
368 : integer, intent(in) :: nx,ny,npol,nchan, nvispol, nvischan, nrow
369 : complex, intent(in) :: values(nvispol, nvischan, nrow)
370 : complex, intent(inout) :: grid(nx, ny, npol, nchan)
371 : integer, intent(in) :: x0, y0, nxsub, nysub
372 : double precision, intent(in) :: convFunc(*)
373 : integer, intent(in) :: chanmap(nvischan), polmap(nvispol)
374 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
375 : integer, intent(in) :: rflag(nrow)
376 : real, intent(in) :: weight(nvischan, nrow)
377 : double precision, intent(inout) :: sumwt(npol, nchan)
378 : integer, intent(in) :: support, sampling
379 : integer, intent(in) :: dopsf, rbeg, rend
380 :
381 : integer, intent(in) :: loc(2, nvischan, nrow)
382 : integer, intent(in) :: off(2, nvischan, nrow)
383 : integer :: iloc(2)
384 : complex, intent(in) :: phasor(nvischan, nrow)
385 : double complex :: nvalue
386 :
387 : double precision :: norm
388 : double precision :: wt, wtx, wty
389 :
390 : logical :: onmygrid
391 :
392 : double precision :: pos(2)
393 : integer :: ix, iy, ipol, ichan
394 : integer :: apol, achan, irow
395 : integer :: posx, posy
396 : integer :: msupporty, psupporty
397 : integer :: msupportx, psupportx
398 26400 : double precision :: cfnow(-support:support, -support:support)
399 :
400 :
401 9292800 : do irow=rbeg, rend
402 9292800 : if(rflag(irow).eq.0) then
403 935906400 : do ichan=1, nvischan
404 926640000 : achan=chanmap(ichan)+1
405 926640000 : if((achan.ge.1).and.(achan.le.nchan).and.
406 9266400 : $ (weight(ichan,irow).ne.0.0)) then
407 : C call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c,
408 : C $ scale, offset, sampling, pos, loc, off, phasor)
409 9266400 : if (onmygrid(loc(1,ichan, irow), nx, ny, x0, y0, nxsub,
410 : $ nysub, support)) then
411 3911055 : do ipol=1, nvispol
412 3128844 : apol=polmap(ipol)+1
413 : if((flag(ipol,ichan,irow).ne.1).and.
414 3911055 : $ (apol.ge.1).and.(apol.le.npol)) then
415 : C If we are making a PSF then we don't want to phase
416 : C rotate but we do want to reproject uvw
417 1564422 : if(dopsf.eq.1) then
418 0 : nvalue=cmplx(weight(ichan,irow))
419 : else
420 : nvalue=weight(ichan,irow)*
421 1564422 : $ (values(ipol,ichan,irow)*phasor(ichan, irow))
422 : end if
423 1564422 : norm=0.0
424 12515376 : do iy=-support, support
425 10950954 : iloc(2)=abs(sampling*iy+off(2,ichan,irow))+1
426 10950954 : wty=convFunc(iloc(2))
427 89172054 : do ix=-support, support
428 : iloc(1)=abs(sampling*ix+off(1,ichan,
429 76656678 : $ irow))+1
430 76656678 : wtx=convFunc(iloc(1))
431 76656678 : cfnow(ix, iy)=wtx*wty
432 87607632 : norm=norm+cfnow(ix,iy)
433 : end do
434 : end do
435 12515376 : do iy=-support, support
436 89172054 : do ix=-support, support
437 87607632 : cfnow(ix, iy)=cfnow(ix, iy)/norm
438 : end do
439 : end do
440 1564422 : norm=0.0
441 1564422 : msupporty=-support
442 1564422 : psupporty=support
443 1564422 : psupportx=support
444 1564422 : msupportx=-support
445 1564422 : if((loc(1, ichan, irow)-support) .lt. x0)
446 658284 : $ msupportx= -(loc(1, ichan, irow)-x0)
447 1564422 : if((loc(1, ichan, irow)+support).ge.(x0+nxsub))
448 742560 : $ psupportx= x0+nxsub-loc(1, ichan, irow)-1
449 1564422 : if((loc(2, ichan, irow)-support) .lt. y0)
450 659190 : $ msupporty= -(loc(2, ichan, irow)-y0)
451 1564422 : if((loc(2, ichan, irow)+support).ge.(y0+nysub))
452 786464 : $ psupporty= y0+nysub-loc(2, ichan, irow)-1
453 7010128 : do iy=msupporty,psupporty
454 5445706 : posy=loc(2, ichan, irow)+iy
455 27648928 : do ix=msupportx,psupportx
456 20638800 : posx=loc(1, ichan, irow)+ix
457 : grid(posx,posy,apol,achan)=
458 : $ grid(posx, posy,apol,achan)+
459 20638800 : $ nvalue*cfnow(ix, iy)
460 26084506 : norm=norm+cfnow(ix,iy)
461 : end do
462 : end do
463 : sumwt(apol,achan)=sumwt(apol,achan)+
464 1564422 : $ weight(ichan,irow)*norm
465 : end if
466 : end do
467 : end if
468 : end if
469 : end do
470 : end if
471 : end do
472 26400 : return
473 : end
474 :
475 :
476 :
477 : C
478 : C Degrid a number of visibility records
479 : C
480 0 : subroutine dgrid (uvw, dphase, values, nvispol, nvischan,
481 0 : $ flag, rflag,
482 0 : $ nrow, rownum, scale, offset, grid, nx, ny, npol, nchan, freq,
483 : $ c, support, sampling, convFunc, chanmap, polmap)
484 :
485 : implicit none
486 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
487 : complex values(nvispol, nvischan, nrow)
488 : complex grid(nx, ny, npol, nchan)
489 : double precision uvw(3, nrow), freq(nvischan), c, scale(2),
490 : $ offset(2)
491 : double precision dphase(nrow), uvdist
492 : complex phasor
493 : integer flag(nvispol, nvischan, nrow)
494 : integer rflag(nrow)
495 : integer rownum
496 : integer support, sampling
497 : integer chanmap(*), polmap(*)
498 :
499 : complex nvalue
500 :
501 : double precision convFunc(*)
502 : real norm
503 :
504 : logical ogrid
505 :
506 : double precision pos(2)
507 : integer loc(2), off(2), iloc(2)
508 : integer rbeg, rend
509 : integer ix, iy, ipol, ichan
510 : integer apol, achan, irow
511 : real wt, wtx, wty
512 :
513 0 : irow=rownum
514 :
515 0 : if(irow.ge.0) then
516 0 : rbeg=irow+1
517 0 : rend=irow+1
518 : else
519 0 : rbeg=1
520 0 : rend=nrow
521 : end if
522 :
523 0 : do irow=rbeg, rend
524 0 : if(rflag(irow).eq.0) then
525 0 : do ichan=1, nvischan
526 0 : achan=chanmap(ichan)+1
527 0 : if((achan.ge.1).and.(achan.le.nchan)) then
528 : call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c,
529 0 : $ scale, offset, sampling, pos, loc, off, phasor)
530 0 : if (ogrid(nx, ny, loc, support)) then
531 0 : do ipol=1, nvispol
532 0 : apol=polmap(ipol)+1
533 : if((flag(ipol,ichan,irow).ne.1).and.
534 0 : $ (apol.ge.1).and.(apol.le.npol)) then
535 0 : nvalue=0.0
536 0 : norm=0.0
537 0 : do iy=-support,support
538 0 : iloc(2)=abs(sampling*iy+off(2))+1
539 0 : wty=convFunc(iloc(2))
540 0 : do ix=-support,support
541 0 : iloc(1)=abs(sampling*ix+off(1))+1
542 0 : wtx=convFunc(iloc(1))
543 0 : wt=wtx*wty
544 0 : norm=norm+wt
545 : nvalue=nvalue+wt*
546 0 : $ grid(loc(1)+ix,loc(2)+iy,apol,achan)
547 : end do
548 : end do
549 : values(ipol,ichan,irow)=(nvalue*conjg(phasor))
550 0 : $ /norm
551 : end if
552 : end do
553 : end if
554 : end if
555 : end do
556 : end if
557 : end do
558 0 : return
559 : end
560 :
561 :
562 6560 : subroutine sectdgrid (values, nvispol, nvischan,
563 6560 : $ flag, rflag,
564 6560 : $ nrow, grid, nx, ny, npol, nchan,
565 : $ support, sampling, convFunc, chanmap, polmap, rbeg, rend,
566 6560 : $ loc,off,phasor)
567 :
568 : implicit none
569 : integer, intent(in) :: nx,ny,npol,nchan,nvispol,nvischan,nrow
570 : complex, intent(inout) :: values(nvispol, nvischan, nrow)
571 : complex, intent(in) :: grid(nx, ny, npol, nchan)
572 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
573 : integer, intent(in) :: rflag(nrow)
574 : integer, intent(in) :: support, sampling
575 : integer, intent(in) :: chanmap(*), polmap(*)
576 :
577 : complex :: nvalue
578 :
579 : double precision, intent(in) :: convFunc(*)
580 : real norm
581 :
582 : logical ogrid
583 : integer, intent(in) :: loc(2, nvischan, nrow)
584 : integer, intent(in) :: off(2, nvischan, nrow)
585 : complex, intent(in) :: phasor(nvischan, nrow)
586 : integer :: iloc(2)
587 : integer, intent(in) :: rbeg, rend
588 : integer ix, iy, ipol, ichan
589 : integer apol, achan, irow
590 : real wt, wtx, wty
591 :
592 :
593 117500 : do irow=rbeg, rend
594 117500 : if(rflag(irow).eq.0) then
595 7171680 : do ichan=1, nvischan
596 7060740 : achan=chanmap(ichan)+1
597 7171680 : if((achan.ge.1).and.(achan.le.nchan)) then
598 : C call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c,
599 : C $ scale, offset, sampling, pos, loc, off, phasor)
600 110940 : if (ogrid(nx, ny, loc(1,ichan, irow) , support)) then
601 473220 : do ipol=1, nvispol
602 362280 : apol=polmap(ipol)+1
603 : if((flag(ipol,ichan,irow).ne.1).and.
604 473220 : $ (apol.ge.1).and.(apol.le.npol)) then
605 221880 : nvalue=0.0
606 221880 : norm=0.0
607 1775040 : do iy=-support,support
608 : iloc(2)=abs(sampling*iy+off(2,ichan,
609 1553160 : $ irow))+1
610 1553160 : wty=convFunc(iloc(2))
611 12647160 : do ix=-support,support
612 : iloc(1)=abs(sampling*ix+off(1,ichan,
613 10872120 : $ irow))+1
614 10872120 : wtx=convFunc(iloc(1))
615 10872120 : wt=wtx*wty
616 10872120 : norm=norm+wt
617 : nvalue=nvalue+wt*
618 : $ grid(loc(1, ichan, irow)+ix,
619 12425280 : $ loc(2, ichan, irow)+iy,apol,achan)
620 : end do
621 : end do
622 : values(ipol,ichan,irow)=(nvalue*conjg(
623 : $ phasor(ichan, irow)))
624 221880 : $ /norm
625 : end if
626 : end do
627 : end if
628 : end if
629 : end do
630 : end if
631 : end do
632 6560 : return
633 : end
634 :
635 800 : subroutine sectdgridjy (values, nvispol, nvischan,
636 800 : $ flag, rflag,
637 800 : $ nrow, grid, nx, ny, npol, nchan,
638 : $ support, sampling, convFunc, chanmap, polmap, rbeg, rend,
639 800 : $ loc,off,phasor, scalechan)
640 :
641 : implicit none
642 : integer, intent(in) :: nx,ny,npol,nchan,nvispol,nvischan,nrow
643 : complex, intent(inout) :: values(nvispol, nvischan, nrow)
644 : complex, intent(in) :: grid(nx, ny, npol, nchan)
645 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
646 : integer, intent(in) :: rflag(nrow)
647 : integer, intent(in) :: support, sampling
648 : integer, intent(in) :: chanmap(*), polmap(*)
649 : complex :: nvalue
650 :
651 : double precision, intent(in) :: convFunc(*), scalechan(*)
652 : real norm
653 :
654 : logical ogrid
655 : integer, intent(in) :: loc(2, nvischan, nrow)
656 : integer, intent(in) :: off(2, nvischan, nrow)
657 : complex, intent(in) :: phasor(nvischan, nrow)
658 : integer :: iloc(2)
659 : integer, intent(in) :: rbeg, rend
660 : integer ix, iy, ipol, ichan
661 : integer apol, achan, irow
662 : real wt, wtx, wty
663 :
664 :
665 71000 : do irow=rbeg, rend
666 71000 : if(rflag(irow).eq.0) then
667 7090200 : do ichan=1, nvischan
668 7020000 : achan=chanmap(ichan)+1
669 7090200 : if((achan.ge.1).and.(achan.le.nchan)) then
670 : C call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c,
671 : C $ scale, offset, sampling, pos, loc, off, phasor)
672 70200 : if (ogrid(nx, ny, loc(1,ichan, irow) , support)) then
673 351000 : do ipol=1, nvispol
674 280800 : apol=polmap(ipol)+1
675 : if((flag(ipol,ichan,irow).ne.1).and.
676 351000 : $ (apol.ge.1).and.(apol.le.npol)) then
677 140400 : nvalue=0.0
678 140400 : norm=0.0
679 1123200 : do iy=-support,support
680 : iloc(2)=abs(sampling*iy+off(2,ichan,
681 982800 : $ irow))+1
682 982800 : wty=convFunc(iloc(2))
683 8002800 : do ix=-support,support
684 : iloc(1)=abs(sampling*ix+off(1,ichan,
685 6879600 : $ irow))+1
686 6879600 : wtx=convFunc(iloc(1))
687 6879600 : wt=wtx*wty
688 6879600 : norm=norm+wt
689 : nvalue=nvalue+wt*scalechan(ichan)*
690 : $ grid(loc(1, ichan, irow)+ix,
691 7862400 : $ loc(2, ichan, irow)+iy,apol,achan)
692 : end do
693 : end do
694 : values(ipol,ichan,irow)=(nvalue*conjg(
695 : $ phasor(ichan, irow)))
696 140400 : $ /norm
697 : end if
698 : end do
699 : end if
700 : end if
701 : end do
702 : end if
703 : end do
704 800 : return
705 : end
706 :
707 :
708 :
709 : C
710 : C Calculate gridded coordinates and the phasor needed for
711 : C phase rotation.
712 : C
713 0 : subroutine sgrid (uvw, dphase, freq, c, scale, offset, sampling,
714 : $ pos, loc, off, phasor)
715 : implicit none
716 : integer sampling
717 : integer loc(2), off(2)
718 : double precision uvw(3), freq, c, scale(2), offset(2)
719 : double precision pos(2)
720 : double precision dphase, phase
721 : complex phasor
722 : integer idim
723 : double precision pi
724 : data pi/3.14159265358979323846/
725 :
726 0 : do idim=1,2
727 0 : pos(idim)=scale(idim)*uvw(idim)*freq/c+(offset(idim)+1.0)
728 0 : loc(idim)=nint(pos(idim))
729 0 : off(idim)=nint((loc(idim)-pos(idim))*sampling)
730 : end do
731 0 : phase=-2.0D0*pi*dphase*freq/c
732 0 : phasor=cmplx(cos(phase), sin(phase))
733 0 : return
734 : end
735 : C
736 : C Is this on the grid?
737 : C
738 181140 : logical function ogrid (nx, ny, loc, support)
739 : implicit none
740 : integer nx, ny, loc(2), support
741 : ogrid=(loc(1)-support.ge.1).and.(loc(1)+support.le.nx).and.
742 181140 : $ (loc(2)-support.ge.1).and.(loc(2)+support.le.ny)
743 181140 : return
744 : end
745 :
746 17006160 : logical function onmygrid(loc, nx, ny, nx0, ny0, nxsub, nysub,
747 : $ support)
748 : implicit none
749 : integer nx, ny, nx0, ny0, nxsub, nysub, loc(2), support
750 : logical ogrid
751 : C logical onmygrid
752 : C onmygrid=ogrid(nx, ny, loc, support)
753 : C if(.not. onmygrid) then
754 : C return
755 : C end if
756 : C onmygrid=(loc(1)-nx0 .ge. (-support)) .and. ((loc(1)-nx0-nxsub)
757 : C $ .le. (support)) .and.((loc(2)-ny0) .ge. (-support)) .and.
758 : C $ ((loc(2)-ny0-nysub) .le. (support))
759 : C----------------------------
760 : integer loc1sub, loc1plus, loc2sub, loc2plus
761 17006160 : loc1sub=loc(1)-support
762 17006160 : loc1plus=loc(1)+support
763 17006160 : loc2sub=loc(2)-support
764 17006160 : loc2plus=loc(2)+support
765 : C----------------
766 : C onmygrid=(loc1plus .ge. nx0) .and. (loc1sub
767 : C $ .le. (nx0+nxsub)) .and.(loc2plus .ge. ny0) .and.
768 : C $ (loc2sub .le. (ny0+nysub))
769 : C--------------
770 : onmygrid=(loc1plus .ge. nx0) .and. (loc1sub
771 : $ .le. (nx0+nxsub)) .and.(loc2plus .ge. ny0) .and.
772 : $ (loc2sub .le. (ny0+nysub)) .and. (loc1sub.ge.1)
773 : $ .and.(loc1plus.le.nx).and.
774 17006160 : $ (loc2sub.ge.1).and.(loc2plus.le.ny)
775 :
776 17006160 : return
777 : end
778 711455 : subroutine locuvw(uvw, dphase, freq, nvchan, scale, offset,
779 711455 : $ sampling, loc, off, phasor, irow, doW, Cinv)
780 : implicit none
781 : double precision, intent(in) :: uvw(3, *), dphase(*), freq(*)
782 : integer, intent(in) :: nvchan, sampling, irow, doW
783 : double precision, intent(in) :: scale(3), offset(3), Cinv
784 : integer, intent(inout) :: loc(2+doW, nvchan, *), off(2+doW,
785 : $ nvchan, *)
786 : complex, intent(inout) :: phasor(nvchan, *)
787 : integer :: f, k, row
788 : double precision :: phase, pos
789 : double precision :: pi
790 : data pi/3.14159265358979323846/
791 :
792 711455 : row=irow+1
793 :
794 44807135 : do f=1, nvchan
795 132287040 : do k=1,2
796 88191360 : pos=scale(k)*uvw(k,row)*freq(f)*Cinv+(offset(k)+1.0)
797 88191360 : loc(k, f, row)=nint(pos)
798 132287040 : off(k, f, row)=nint((loc(k, f, row)-pos)*sampling)
799 : end do
800 44095680 : phase=-2.0D0*pi*dphase(row)*freq(f)*Cinv
801 44095680 : phasor(f,row)=cmplx(cos(phase), sin(phase))
802 44807135 : if(doW .eq. 1) then
803 : pos=sqrt(abs(scale(3)*uvw(3, row)*freq(f)*Cinv))+offset(3)
804 14040000 : $ +1.0
805 14040000 : loc(3,f, row)=nint(pos)
806 14040000 : off(3, f, row)=0
807 : end if
808 : end do
809 711455 : return
810 : end
|
