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 16617992 : subroutine sectggridd(values, nvispol, nvischan,
221 16617992 : $ dopsf, flag, rflag, weight, nrow,
222 16617992 : $ grid, nx, ny, npol, nchan,
223 16617992 : $ support, sampling, convFunc, chanmap, polmap, sumwt, x0,
224 16617992 : $ 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 16617992 : double precision :: cfnow(-support:support, -support:support)
259 :
260 :
261 5716392680 : do irow=rbeg, rend
262 5716392680 : if(rflag(irow).eq.0) then
263 48196014660 : do ichan=1, nvischan
264 42651879956 : achan=chanmap(ichan)+1
265 42651879956 : if((achan.ge.1).and.(achan.le.nchan).and.
266 5544134704 : $ (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 40683206292 : if (onmygrid(loc(1,ichan, irow), nx, ny, x0, y0, nxsub,
271 : $ nysub, support)) then
272 17249060517 : do ipol=1, nvispol
273 11499998520 : apol=polmap(ipol)+1
274 : if((flag(ipol,ichan,irow).ne.1).and.
275 17249060517 : $ (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 10991852776 : if(dopsf.eq.1) then
279 4584949468 : nvalue=cmplx(weight(ichan,irow))
280 : else
281 : nvalue=weight(ichan,irow)*
282 6406903308 : $ (values(ipol,ichan,irow)*phasor(ichan, irow))
283 : end if
284 10991852776 : norm=0.0
285 10991852776 : msupporty=-support
286 10991852776 : psupporty=support
287 10991852776 : psupportx=support
288 10991852776 : msupportx=-support
289 87934822208 : do iy=-support, support
290 76942969432 : iloc(2)=abs(sampling*iy+off(2,ichan,irow))+1
291 76942969432 : wty=convFunc(iloc(2))
292 >62653*10^7 : do ix=-support, support
293 : iloc(1)=abs(sampling*ix+off(1,ichan,
294 >53860*10^7 : $ irow))+1
295 >53860*10^7 : wtx=convFunc(iloc(1))
296 >53860*10^7 : cfnow(ix, iy)=wtx*wty
297 >61554*10^7 : norm=norm+cfnow(ix,iy)
298 : end do
299 : end do
300 87934822208 : do iy=-support, support
301 >62653*10^7 : do ix=-support, support
302 >61554*10^7 : cfnow(ix, iy)=cfnow(ix, iy)/norm
303 : end do
304 : end do
305 10991852776 : if((loc(1, ichan, irow)-support) .lt. x0)
306 5309909806 : $ msupportx= -(loc(1, ichan, irow)-x0)
307 10991852776 : if((loc(1, ichan, irow)+support).ge.(x0+nxsub))
308 6195845206 : $ psupportx= x0+nxsub-loc(1, ichan, irow)-1
309 10991852776 : if((loc(2, ichan, irow)-support) .lt. y0)
310 5230516868 : $ msupporty= -(loc(2, ichan, irow)-y0)
311 10991852776 : if((loc(2, ichan, irow)+support).ge.(y0+nysub))
312 6162964310 : $ psupporty= y0+nysub-loc(2, ichan, irow)-1
313 : C write(*,*) msupportx, psupportx, msupporty, psupporty
314 10991852776 : norm=0.0
315 44863344750 : do iy=msupporty,psupporty
316 33871491974 : 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 >17072*10^7 : do ix=msupportx,psupportx
322 >12586*10^7 : 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 >12586*10^7 : $ nvalue*cfnow(ix,iy)
333 >15973*10^7 : 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 10991852776 : $ 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 16617992 : 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 77336 : subroutine sectggrids(values, nvispol, nvischan,
359 77336 : $ dopsf, flag, rflag, weight, nrow,
360 77336 : $ grid, nx, ny, npol, nchan,
361 77336 : $ support, sampling, convFunc, chanmap, polmap, sumwt, x0,
362 77336 : $ 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 77336 : double precision :: cfnow(-support:support, -support:support)
399 :
400 :
401 25091312 : do irow=rbeg, rend
402 25091312 : if(rflag(irow).eq.0) then
403 1070526672 : do ichan=1, nvischan
404 1046038456 : achan=chanmap(ichan)+1
405 1046038456 : if((achan.ge.1).and.(achan.le.nchan).and.
406 24488216 : $ (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 123928760 : if (onmygrid(loc(1,ichan, irow), nx, ny, x0, y0, nxsub,
410 : $ nysub, support)) then
411 97195632 : do ipol=1, nvispol
412 65318562 : apol=polmap(ipol)+1
413 : if((flag(ipol,ichan,irow).ne.1).and.
414 97195632 : $ (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 51146388 : if(dopsf.eq.1) then
418 49581966 : 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 51146388 : norm=0.0
424 111679308 : do iy=-support, support
425 60532920 : iloc(2)=abs(sampling*iy+off(2,ichan,irow))+1
426 60532920 : wty=convFunc(iloc(2))
427 237917952 : do ix=-support, support
428 : iloc(1)=abs(sampling*ix+off(1,ichan,
429 126238644 : $ irow))+1
430 126238644 : wtx=convFunc(iloc(1))
431 126238644 : cfnow(ix, iy)=wtx*wty
432 186771564 : norm=norm+cfnow(ix,iy)
433 : end do
434 : end do
435 111679308 : do iy=-support, support
436 237917952 : do ix=-support, support
437 186771564 : cfnow(ix, iy)=cfnow(ix, iy)/norm
438 : end do
439 : end do
440 51146388 : norm=0.0
441 51146388 : msupporty=-support
442 51146388 : psupporty=support
443 51146388 : psupportx=support
444 51146388 : msupportx=-support
445 51146388 : if((loc(1, ichan, irow)-support) .lt. x0)
446 658284 : $ msupportx= -(loc(1, ichan, irow)-x0)
447 51146388 : if((loc(1, ichan, irow)+support).ge.(x0+nxsub))
448 4879990 : $ psupportx= x0+nxsub-loc(1, ichan, irow)-1
449 51146388 : if((loc(2, ichan, irow)-support) .lt. y0)
450 659190 : $ msupporty= -(loc(2, ichan, irow)-y0)
451 51146388 : if((loc(2, ichan, irow)+support).ge.(y0+nysub))
452 4693580 : $ psupporty= y0+nysub-loc(2, ichan, irow)-1
453 102266944 : do iy=msupporty,psupporty
454 51120556 : posy=loc(2, ichan, irow)+iy
455 164891000 : do ix=msupportx,psupportx
456 62624056 : posx=loc(1, ichan, irow)+ix
457 : grid(posx,posy,apol,achan)=
458 : $ grid(posx, posy,apol,achan)+
459 62624056 : $ nvalue*cfnow(ix, iy)
460 113744612 : norm=norm+cfnow(ix,iy)
461 : end do
462 : end do
463 : sumwt(apol,achan)=sumwt(apol,achan)+
464 51146388 : $ 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 77336 : 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 965194 : subroutine sectdgrid (values, nvispol, nvischan,
563 965194 : $ flag, rflag,
564 965194 : $ nrow, grid, nx, ny, npol, nchan,
565 : $ support, sampling, convFunc, chanmap, polmap, rbeg, rend,
566 965194 : $ 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 56068336 : do irow=rbeg, rend
594 56068336 : if(rflag(irow).eq.0) then
595 521365933 : do ichan=1, nvischan
596 467600784 : achan=chanmap(ichan)+1
597 521365933 : 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 460314024 : if (ogrid(nx, ny, loc(1,ichan, irow) , support)) then
601 1380731073 : do ipol=1, nvispol
602 920615582 : apol=polmap(ipol)+1
603 : if((flag(ipol,ichan,irow).ne.1).and.
604 1380731073 : $ (apol.ge.1).and.(apol.le.npol)) then
605 895248406 : nvalue=0.0
606 895248406 : norm=0.0
607 7161987248 : do iy=-support,support
608 : iloc(2)=abs(sampling*iy+off(2,ichan,
609 6266738842 : $ irow))+1
610 6266738842 : wty=convFunc(iloc(2))
611 51029159142 : do ix=-support,support
612 : iloc(1)=abs(sampling*ix+off(1,ichan,
613 43867171894 : $ irow))+1
614 43867171894 : wtx=convFunc(iloc(1))
615 43867171894 : wt=wtx*wty
616 43867171894 : norm=norm+wt
617 : nvalue=nvalue+wt*
618 : $ grid(loc(1, ichan, irow)+ix,
619 50133910736 : $ 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 895248406 : $ /norm
625 : end if
626 : end do
627 : end if
628 : end if
629 : end do
630 : end if
631 : end do
632 965194 : return
633 : end
634 :
635 1088 : subroutine sectdgridjy (values, nvispol, nvischan,
636 1088 : $ flag, rflag,
637 1088 : $ nrow, grid, nx, ny, npol, nchan,
638 : $ support, sampling, convFunc, chanmap, polmap, rbeg, rend,
639 1088 : $ 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 96560 : do irow=rbeg, rend
666 96560 : if(rflag(irow).eq.0) then
667 7368192 : do ichan=1, nvischan
668 7272720 : achan=chanmap(ichan)+1
669 7368192 : 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 322920 : if (ogrid(nx, ny, loc(1,ichan, irow) , support)) then
673 1109160 : do ipol=1, nvispol
674 786240 : apol=polmap(ipol)+1
675 : if((flag(ipol,ichan,irow).ne.1).and.
676 1109160 : $ (apol.ge.1).and.(apol.le.npol)) then
677 645840 : nvalue=0.0
678 645840 : norm=0.0
679 5166720 : do iy=-support,support
680 : iloc(2)=abs(sampling*iy+off(2,ichan,
681 4520880 : $ irow))+1
682 4520880 : wty=convFunc(iloc(2))
683 36812880 : do ix=-support,support
684 : iloc(1)=abs(sampling*ix+off(1,ichan,
685 31646160 : $ irow))+1
686 31646160 : wtx=convFunc(iloc(1))
687 31646160 : wt=wtx*wty
688 31646160 : norm=norm+wt
689 : nvalue=nvalue+wt*scalechan(ichan)*
690 : $ grid(loc(1, ichan, irow)+ix,
691 36167040 : $ 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 645840 : $ /norm
697 : end if
698 : end do
699 : end if
700 : end if
701 : end do
702 : end if
703 : end do
704 1088 : 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 460636944 : 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 460636944 : $ (loc(2)-support.ge.1).and.(loc(2)+support.le.ny)
743 460636944 : return
744 : end
745 :
746 40807135052 : 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 40807135052 : loc1sub=loc(1)-support
762 40807135052 : loc1plus=loc(1)+support
763 40807135052 : loc2sub=loc(2)-support
764 40807135052 : 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 40807135052 : $ (loc2sub.ge.1).and.(loc2plus.le.ny)
775 :
776 40807135052 : return
777 : end
778 250105170 : subroutine locuvw(uvw, dphase, freq, nvchan, scale, offset,
779 250105170 : $ 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 250105170 : row=irow+1
793 :
794 2280094027 : do f=1, nvchan
795 6089966571 : do k=1,2
796 4059977714 : pos=scale(k)*uvw(k,row)*freq(f)*Cinv+(offset(k)+1.0)
797 4059977714 : loc(k, f, row)=nint(pos)
798 6089966571 : off(k, f, row)=nint((loc(k, f, row)-pos)*sampling)
799 : end do
800 2029988857 : phase=-2.0D0*pi*dphase(row)*freq(f)*Cinv
801 2029988857 : phasor(f,row)=cmplx(cos(phase), sin(phase))
802 2280094027 : if(doW .eq. 1) then
803 : pos=sqrt(abs(scale(3)*uvw(3, row)*freq(f)*Cinv))+offset(3)
804 18268146 : $ +1.0
805 18268146 : loc(3,f, row)=nint(pos)
806 18268146 : off(3, f, row)=0
807 : end if
808 : end do
809 250105170 : return
810 : end
|
