Line data Source code
1 :
2 : *=======================================================================
3 : * Copyright (C) 1999-2013
4 : * Associated Universities, Inc. Washington DC, USA.
5 : *
6 : * This library is free software; you can redistribute it and/or
7 : * modify it under the terms of the GNU Library General Public
8 : * License as published by the Free Software Foundation; either
9 : * version 2 of the License, or (at your option) any later version.
10 : *
11 : * This library is distributed in the hope that it will be useful,
12 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : * GNU Library General Public License for more details.
15 : *
16 : * You should have received a copy of the GNU Library General Public
17 : * License along with this library; if not, write to the Free
18 : * Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,
19 : * MA 02139, USA.
20 : *
21 : * Correspondence concerning AIPS++ should be addressed as follows:
22 : * Internet email: casa-feedback@nrao.edu.
23 : * Postal address: AIPS++ Project Office
24 : * National Radio Astronomy Observatory
25 : * 520 Edgemont Road
26 : * Charlottesville, VA 22903-2475 USA
27 : *
28 : * $Id$
29 : *-----------------------------------------------------------------------
30 : C
31 : C Grid a number of visibility records
32 : C
33 0 : subroutine gmosd (uvw, dphase, values, nvispol, nvischan,
34 0 : $ dopsf, flag, rflag, weight, nrow, rownum,
35 0 : $ scale, offset, grid, nx, ny, npol, nchan, freq, c,
36 0 : $ support, convsize, sampling, convfunc,
37 0 : $ chanmap, polmap,
38 0 : $ sumwt, weightgrid, convweight, doweightgrid, convplanemap,
39 : $ nconvplane)
40 :
41 : implicit none
42 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
43 : complex values(nvispol, nvischan, nrow)
44 : double complex grid(nx, ny, npol, nchan)
45 :
46 : double precision uvw(3, nrow), freq(nvischan), c, scale(2),
47 : $ offset(2)
48 : double precision dphase(nrow), uvdist
49 : double precision xlast, ylast
50 : complex phasor
51 : integer flag(nvispol, nvischan, nrow)
52 : integer rflag(nrow)
53 : real weight(nvischan, nrow), phase
54 : double precision sumwt(npol, nchan)
55 : integer rownum
56 : integer support
57 : integer chanmap(nchan), polmap(npol)
58 : integer dopsf
59 : double complex weightgrid(nx, ny, npol, nchan)
60 : integer doweightgrid
61 :
62 : double complex nvalue
63 : double complex nweight
64 : integer convsize, sampling
65 : integer nconvplane
66 : integer convplanemap(nrow)
67 : complex convfunc(convsize, convsize, nconvplane), cwt, crot
68 : complex convweight(convsize, convsize, nconvplane)
69 : integer sampsupp
70 :
71 :
72 : complex sconv(-sampling*(support+1):sampling*(support+1),
73 0 : $ -sampling*(support+1):sampling*(support+1), nconvplane)
74 : complex sconv2(-sampling*(support+1):sampling*(support+1),
75 0 : $ -sampling*(support+1):sampling*(support+1), nconvplane)
76 : real sumsconv
77 : real sumsconv2
78 : real ratioofbeam
79 :
80 : real norm
81 : real wt
82 :
83 : logical omos, doshift
84 :
85 : real pos(3)
86 : integer loc(2), off(2), iloc(2)
87 : integer iiloc(2)
88 : integer rbeg, rend
89 : integer ix, iy, iz, ipol, ichan
90 : integer apol, achan, aconvplane, irow
91 : double precision pi
92 : data pi/3.14159265358979323846/
93 : real maxsconv2, minsconv2
94 0 : sampsupp=(support+1)*sampling
95 0 : irow=rownum
96 :
97 0 : sumsconv=0
98 0 : sumsconv2=0
99 0 : ratioofbeam=1.0
100 :
101 : CCCCCCCCCCCCCCCCCCCCCCCC
102 : C minsconv2=1e17
103 : C maxsconv2=0.0
104 : CCCCCCCCCCCCCCCCCCCCCCCC
105 : C write(*,*) scale, offset
106 0 : do iz=1, nconvplane
107 0 : do iy=-sampsupp,sampsupp
108 0 : iloc(2)=iy+(convsize)/2+1
109 0 : do ix=-sampsupp,sampsupp
110 0 : iloc(1)=ix+(convsize)/2+1
111 0 : sconv(ix,iy,iz)=(convfunc(iloc(1), iloc(2),iz))
112 0 : sconv2(ix,iy,iz)=convweight(iloc(1), iloc(2),iz)
113 : C if(maxsconv2 .lt. abs(sconv2(ix, iy, iz))) then
114 : C maxsconv2=abs(sconv2(ix, iy, iz))
115 : C end if
116 : C if(minsconv2 .gt. abs(sconv2(ix, iy, iz))) then
117 : C minsconv2=abs(sconv2(ix, iy, iz))
118 : C end if
119 : end do
120 : end do
121 : end do
122 :
123 0 : doshift=.FALSE.
124 :
125 0 : if(irow.ge.0) then
126 0 : rbeg=irow+1
127 0 : rend=irow+1
128 : else
129 0 : rbeg=1
130 0 : rend=nrow
131 : end if
132 :
133 : C write(*,*) 'max of sconvs ', maxsconv2, minsconv2, sampsupp,
134 : C $ convsize
135 :
136 :
137 :
138 0 : xlast=0.0
139 0 : ylast=0.0
140 0 : do irow=rbeg, rend
141 0 : aconvplane=convplanemap(irow)+1
142 0 : if(rflag(irow).eq.0) then
143 0 : do ichan=1, nvischan
144 0 : achan=chanmap(ichan)+1
145 0 : if((achan.ge.1).and.(achan.le.nchan).and.
146 0 : $ (weight(ichan,irow).ne.0.0)) then
147 : call smos(uvw(1,irow), dphase(irow), freq(ichan), c,
148 0 : $ scale, offset, sampling, pos, loc, off, phasor)
149 0 : if (omos(nx, ny, loc, support)) then
150 0 : do ipol=1, nvispol
151 0 : apol=polmap(ipol)+1
152 : if((flag(ipol,ichan,irow).ne.1).and.
153 0 : $ (apol.ge.1).and.(apol.le.npol)) then
154 : C If we are making a PSF then we don't want to phase
155 : C rotate but we do want to reproject uvw
156 0 : if(dopsf.eq.1) then
157 0 : nvalue=cmplx(weight(ichan,irow))
158 : else
159 : nvalue=weight(ichan,irow)*
160 0 : $ (values(ipol,ichan,irow)*phasor)
161 : end if
162 0 : if(doweightgrid .gt. 0) then
163 0 : nweight=cmplx(weight(ichan,irow))
164 : end if
165 :
166 : C norm will be the value we would get for the peak
167 : C at the phase center. We will want to normalize
168 : C the final image by this term.
169 0 : norm=0.0
170 0 : if(sampling.eq.1) then
171 0 : do iy=-support,support
172 0 : do ix=-support,support
173 : grid(loc(1)+ix,
174 : $ loc(2)+iy,apol,achan)=
175 : $ grid(loc(1)+ix,
176 : $ loc(2)+iy,apol,achan)+
177 0 : $ nvalue*sconv(ix,iy, aconvplane)
178 :
179 0 : if(doweightgrid .gt. 0) then
180 0 : iloc(1)=nx/2+1+ix
181 0 : iloc(2)=ny/2+1+iy
182 : weightgrid(iloc(1),iloc(2),
183 : $ apol,achan)= weightgrid(
184 : $ iloc(1),iloc(2),apol,achan)
185 0 : $ + nweight*sconv2(ix,iy,aconvplane)
186 :
187 : end if
188 : end do
189 : end do
190 : else
191 : C write(*,*)off
192 0 : do iy=-support, support
193 0 : iloc(2)=(sampling*iy)+off(2)
194 0 : do ix=-support, support
195 0 : iloc(1)=(sampling*ix)+off(1)
196 : cwt=sconv(iloc(1),
197 0 : $ iloc(2),aconvplane)
198 : C write(*,*) support, iloc
199 : grid(loc(1)+ix,
200 : $ loc(2)+iy,apol,achan)=
201 : $ grid(loc(1)+ix,
202 : $ loc(2)+iy,apol,achan)+
203 0 : $ nvalue*cwt
204 0 : if(doweightgrid .gt. 0) then
205 : cwt=sconv2(sampling*ix,
206 : $ sampling*iy,
207 0 : $ aconvplane)
208 0 : iiloc(1)=nx/2+1+ix
209 0 : iiloc(2)=ny/2+1+iy
210 : weightgrid(iiloc(1),iiloc(2),
211 : $ apol,achan)= weightgrid(
212 : $ iiloc(1),iiloc(2),apol,achan)
213 0 : $ + nweight*cwt
214 :
215 : end if
216 : end do
217 : end do
218 : end if
219 : sumwt(apol, achan)= sumwt(apol,achan)+
220 0 : $ weight(ichan,irow)
221 : end if
222 : end do
223 : end if
224 : end if
225 : end do
226 : end if
227 : end do
228 0 : return
229 : end
230 : C
231 :
232 0 : subroutine gmosd2 (uvw, dphase, values, nvispol, nvischan,
233 0 : $ dopsf, flag, rflag, weight, nrow, rownum,
234 0 : $ scale, offset, grid, nx, ny, npol, nchan, freq, c,
235 0 : $ support, convsize, sampling, convfunc,
236 0 : $ chanmap, polmap,
237 0 : $ sumwt, weightgrid, convweight, doweightgrid, convplanemap,
238 0 : $ convchanmap, convpolmap,
239 : $ nconvplane, nconvchan, nconvpol)
240 :
241 : implicit none
242 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
243 : complex values(nvispol, nvischan, nrow)
244 : double complex grid(nx, ny, npol, nchan)
245 :
246 : double precision uvw(3, nrow), freq(nvischan), c, scale(2),
247 : $ offset(2)
248 : double precision dphase(nrow), uvdist
249 : double precision xlast, ylast
250 : complex phasor
251 : integer flag(nvispol, nvischan, nrow)
252 : integer rflag(nrow)
253 : real weight(nvischan, nrow), phase
254 : double precision sumwt(npol, nchan)
255 : integer rownum
256 : integer support
257 : integer chanmap(nchan), polmap(npol)
258 : integer dopsf
259 : double complex weightgrid(nx, ny, npol, nchan)
260 : integer doweightgrid
261 :
262 : double complex nvalue
263 : double complex nweight
264 : integer convsize, sampling
265 : integer nconvplane, nconvchan, nconvpol
266 : integer convplanemap(nrow)
267 : integer convchanmap(nvischan)
268 : integer convpolmap(nvispol)
269 : complex convfunc(convsize, convsize, nconvpol, nconvchan,
270 : $ nconvplane), cwt, crot
271 : complex convweight(convsize, convsize, nconvpol, nconvchan,
272 : $ nconvplane)
273 : integer sampsupp
274 :
275 :
276 : complex sconv(-sampling*(support+1):sampling*(support+1),
277 : $ -sampling*(support+1):sampling*(support+1), nconvplane)
278 : complex sconv2(-sampling*(support+1):sampling*(support+1),
279 : $ -sampling*(support+1):sampling*(support+1), nconvplane)
280 : real sumsconv
281 : real sumsconv2
282 : real ratioofbeam
283 :
284 : real norm
285 : real wt
286 :
287 : logical omos, doshift
288 :
289 : real pos(3)
290 : integer loc(2), off(2), iloc(2)
291 : integer iiloc(2)
292 : integer rbeg, rend
293 : integer ix, iy, iz, ipol, ichan, xind, yind
294 : integer apol, achan, aconvplane, irow
295 : integer aconvpol, aconvchan, xind2, yind2
296 : double precision pi
297 : data pi/3.14159265358979323846/
298 : real maxsconv2, minsconv2
299 0 : sampsupp=(support+1)*sampling
300 0 : irow=rownum
301 :
302 0 : sumsconv=0
303 0 : sumsconv2=0
304 0 : ratioofbeam=1.0
305 :
306 : CCCCCCCCCCCCCCCCCCCCCCCC
307 : C minsconv2=1e17
308 : C maxsconv2=0.0
309 : CCCCCCCCCCCCCCCCCCCCCCCC
310 : C write(*,*) scale, offset
311 : C do iz=1, nconvplane
312 : C do ichan=1, nconvchan
313 : C do ipol=1, nconvpol
314 : C do iy=-sampsupp,sampsupp
315 : C iloc(2)=iy+(convsize)/2+1
316 : C do ix=-sampsupp,sampsupp
317 : C iloc(1)=ix+(convsize)/2+1
318 : C sconv(ix,iy,iz)=(convfunc(iloc(1), iloc(2),ipol,
319 : C $ ichan, iz))
320 : C sconv2(ix,iy,iz)=convweight(iloc(1), iloc(2),iz)
321 : CC if(maxsconv2 .lt. abs(sconv2(ix, iy, iz))) then
322 : CC maxsconv2=abs(sconv2(ix, iy, iz))
323 : CC end if
324 : CC if(minsconv2 .gt. abs(sconv2(ix, iy, iz))) then
325 : CC minsconv2=abs(sconv2(ix, iy, iz))
326 : CC end if
327 : C end do
328 : C end do
329 : C end do
330 :
331 0 : doshift=.FALSE.
332 :
333 0 : if(irow.ge.0) then
334 0 : rbeg=irow+1
335 0 : rend=irow+1
336 : else
337 0 : rbeg=1
338 0 : rend=nrow
339 : end if
340 :
341 : C write(*,*) 'max of sconvs ', maxsconv2, minsconv2, sampsupp,
342 : C $ convsize
343 :
344 :
345 :
346 0 : xlast=0.0
347 0 : ylast=0.0
348 0 : do irow=rbeg, rend
349 0 : aconvplane=convplanemap(irow)+1
350 0 : if(rflag(irow).eq.0) then
351 0 : do ichan=1, nvischan
352 0 : achan=chanmap(ichan)+1
353 0 : aconvchan=convchanmap(ichan)+1
354 0 : if((achan.ge.1).and.(achan.le.nchan).and.
355 0 : $ (weight(ichan,irow).ne.0.0)) then
356 : call smos(uvw(1,irow), dphase(irow), freq(ichan), c,
357 0 : $ scale, offset, sampling, pos, loc, off, phasor)
358 0 : if (omos(nx, ny, loc, support)) then
359 0 : do ipol=1, nvispol
360 0 : apol=polmap(ipol)+1
361 0 : aconvpol=convpolmap(ipol)+1
362 : if((flag(ipol,ichan,irow).ne.1).and.
363 0 : $ (apol.ge.1).and.(apol.le.npol)) then
364 : C If we are making a PSF then we don't want to phase
365 : C rotate but we do want to reproject uvw
366 0 : if(dopsf.eq.1) then
367 0 : nvalue=cmplx(weight(ichan,irow))
368 : else
369 : nvalue=weight(ichan,irow)*
370 0 : $ (values(ipol,ichan,irow)*phasor)
371 : end if
372 0 : if(doweightgrid .gt. 0) then
373 0 : nweight=cmplx(weight(ichan,irow))
374 : end if
375 :
376 : C norm will be the value we would get for the peak
377 : C at the phase center. We will want to normalize
378 : C the final image by this term.
379 0 : norm=0.0
380 0 : if(sampling.eq.1) then
381 0 : do iy=-support,support
382 0 : yind=iy+(convsize)/2+1
383 0 : do ix=-support,support
384 0 : xind=ix+(convsize)/2+1
385 : grid(loc(1)+ix,
386 : $ loc(2)+iy,apol,achan)=
387 : $ grid(loc(1)+ix,
388 : $ loc(2)+iy,apol,achan)+
389 : $ nvalue*convfunc(xind, yind,
390 0 : $ aconvpol, aconvchan, aconvplane)
391 :
392 0 : if(doweightgrid .gt. 0) then
393 0 : iloc(1)=nx/2+1+ix
394 0 : iloc(2)=ny/2+1+iy
395 : weightgrid(iloc(1),iloc(2),
396 : $ apol,achan)= weightgrid(
397 : $ iloc(1),iloc(2),apol,achan)
398 : $ + nweight*convweight(xind, yind,
399 0 : $ aconvpol, aconvchan, aconvplane)
400 :
401 : end if
402 : end do
403 : end do
404 : else
405 : C write(*,*)off
406 0 : do iy=-support, support
407 0 : iloc(2)=(sampling*iy)+off(2)
408 0 : yind=iloc(2)+(convsize)/2+1
409 0 : do ix=-support, support
410 0 : iloc(1)=(sampling*ix)+off(1)
411 0 : xind=iloc(1)+(convsize)/2+1
412 : cwt=convfunc(xind, yind,
413 0 : $ aconvpol, aconvchan, aconvplane)
414 : C write(*,*) support, iloc
415 : grid(loc(1)+ix,
416 : $ loc(2)+iy,apol,achan)=
417 : $ grid(loc(1)+ix,
418 : $ loc(2)+iy,apol,achan)+
419 0 : $ nvalue*cwt
420 0 : if(doweightgrid .gt. 0) then
421 0 : xind2=sampling*ix+(convsize)/2+1
422 0 : yind2=sampling*iy+(convsize)/2+1
423 : cwt=convweight(xind2,
424 0 : $ yind2, aconvpol, aconvchan, aconvplane)
425 0 : iiloc(1)=nx/2+1+ix
426 0 : iiloc(2)=ny/2+1+iy
427 : weightgrid(iiloc(1),iiloc(2),
428 : $ apol,achan)= weightgrid(
429 : $ iiloc(1),iiloc(2),apol,achan)
430 0 : $ + nweight*cwt
431 :
432 : end if
433 : end do
434 : end do
435 : end if
436 : sumwt(apol, achan)= sumwt(apol,achan)+
437 0 : $ weight(ichan,irow)
438 : end if
439 : end do
440 : end if
441 : end if
442 : end do
443 : end if
444 : end do
445 0 : return
446 : end
447 : C
448 :
449 44379 : subroutine gmoswgtd (nvispol, nvischan,
450 44379 : $ flag, rflag, weight, nrow,
451 : $ nx, ny, npol, nchan,
452 : $ support, convsize, sampling,
453 44379 : $ chanmap, polmap,
454 44379 : $ weightgrid, sumwt, convweight, convplanemap,
455 44379 : $ convchanmap, convpolmap,
456 : $ nconvplane, nconvchan, nconvpol, rbeg,
457 44379 : $ rend, loc, off, phasor)
458 :
459 : implicit none
460 : integer, intent(in) :: nx,ny,npol,nchan, nvispol, nvischan, nrow
461 :
462 :
463 : integer, intent(in) :: loc(2, nvischan, nrow)
464 : integer, intent(in) :: off(2, nvischan, nrow)
465 : complex, intent(in) :: phasor(nvischan, nrow)
466 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
467 : integer, intent(in) :: rflag(nrow)
468 : real, intent(in) :: weight(nvischan, nrow)
469 : integer, intent(in) :: support
470 : integer, intent(in) :: chanmap(nchan), polmap(npol)
471 : double complex, intent(inout) :: weightgrid(nx, ny, npol, nchan)
472 : double precision, intent(inout) :: sumwt(npol, nchan)
473 : double complex :: nweight
474 : integer, intent(in) :: convsize, sampling
475 : integer, intent(in) :: nconvplane, nconvchan, nconvpol
476 : integer, intent(in) :: convplanemap(nrow)
477 : integer, intent(in) :: convchanmap(nvischan)
478 : integer, intent(in) :: convpolmap(nvispol)
479 : complex :: cwt
480 : complex, intent(in) :: convweight(convsize, convsize, nconvpol,
481 : $ nconvchan, nconvplane)
482 :
483 :
484 : real :: norm
485 : real :: wt
486 :
487 : logical :: onmosgrid
488 : logical :: doconj
489 :
490 : integer :: iloc(2)
491 : integer :: iiloc(2)
492 : integer, intent(in) :: rbeg, rend
493 : integer :: ix, iy, iz, ipol, ichan, xind, yind
494 : integer :: apol, achan, aconvplane, irow
495 : integer :: aconvpol, aconvchan, xind2, yind2
496 : integer :: posx, posy, msupportx, msupporty, psupportx, psupporty
497 : logical :: centin
498 :
499 :
500 :
501 :
502 12474441 : do irow=rbeg, rend
503 12430062 : aconvplane=abs(convplanemap(irow))+1
504 12430062 : doconj = (convplanemap(irow) < 0)
505 12474441 : if(rflag(irow).eq.0) then
506 43679840 : do ichan=1, nvischan
507 31382192 : achan=chanmap(ichan)+1
508 31382192 : aconvchan=convchanmap(ichan)+1
509 31382192 : if((achan.ge.1).and.(achan.le.nchan).and.
510 12297648 : $ (weight(ichan,irow).ne.0.0)) then
511 28159683 : if (onmosgrid(loc(1, ichan, irow), nx, ny, 1, 1,
512 : $ nx, ny, support, msupportx, msupporty,
513 : $ psupportx, psupporty, centin)) then
514 84452544 : do ipol=1, nvispol
515 56301696 : apol=polmap(ipol)+1
516 56301696 : aconvpol=convpolmap(ipol)+1
517 : if((flag(ipol,ichan,irow).ne.1).and.
518 84452544 : $ (apol.ge.1).and.(apol.le.npol)) then
519 : C If we are making a PSF then we don't want to phase
520 : C rotate but we do want to reproject uvw
521 :
522 :
523 52446888 : nweight=cmplx(weight(ichan,irow))
524 : sumwt(apol, achan)=sumwt(apol, achan)+
525 52446888 : $ weight(ichan, irow)
526 :
527 : C norm will be the value we would get for the peak
528 : C at the phase center. We will want to normalize
529 : C the final image by this term.
530 52446888 : norm=0.0
531 : C write(*,*)off
532 1494058308 : do iy=msupporty, psupporty
533 46370621850 : do ix=msupportx, psupportx
534 :
535 44876563542 : xind2=sampling*ix+(convsize)/2+1
536 44876563542 : yind2=sampling*iy+(convsize)/2+1
537 : cwt=convweight(xind2,
538 : $ yind2, aconvpol,
539 44876563542 : $ aconvchan, aconvplane)
540 :
541 44876563542 : iiloc(1)=nx/2+1+ix
542 44876563542 : iiloc(2)=ny/2+1+iy
543 : weightgrid(iiloc(1),iiloc(2),
544 : $ apol,achan)= weightgrid(
545 : $ iiloc(1),iiloc(2),apol,achan)
546 46318174962 : $ + nweight*cwt
547 :
548 :
549 : end do
550 : end do
551 :
552 : end if
553 : end do
554 : end if
555 : end if
556 : end do
557 : end if
558 : end do
559 44379 : return
560 : end
561 :
562 :
563 : C same as gmoswgtd except with varying support
564 0 : subroutine gmoswgtd2 (nvispol, nvischan,
565 0 : $ flag, rflag, weight, nrow,
566 : $ nx, ny, npol, nchan,
567 0 : $ supports, convsize, sampling,
568 0 : $ chanmap, polmap,
569 0 : $ weightgrid, sumwt, convweight, convplanemap,
570 0 : $ convchanmap, convpolmap,
571 : $ nconvplane, nconvchan, nconvpol, rbeg,
572 0 : $ rend, loc, off, phasor)
573 :
574 : implicit none
575 : integer, intent(in) :: nx,ny,npol,nchan, nvispol, nvischan, nrow
576 :
577 :
578 : integer, intent(in) :: loc(2, nvischan, nrow)
579 : integer, intent(in) :: off(2, nvischan, nrow)
580 : complex, intent(in) :: phasor(nvischan, nrow)
581 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
582 : integer, intent(in) :: rflag(nrow)
583 : real, intent(in) :: weight(nvischan, nrow)
584 :
585 : integer, intent(in) :: chanmap(nchan), polmap(npol)
586 : double complex, intent(inout) :: weightgrid(nx, ny, npol, nchan)
587 : double precision, intent(inout) :: sumwt(npol, nchan)
588 : double complex :: nweight
589 : integer, intent(in) :: convsize, sampling
590 : integer, intent(in) :: nconvplane, nconvchan, nconvpol
591 : integer, intent(in) :: convplanemap(nrow)
592 : integer, intent(in) :: convchanmap(nvischan)
593 : integer, intent(in) :: convpolmap(nvispol)
594 : complex :: cwt
595 : complex, intent(in) :: convweight(convsize, convsize, nconvpol,
596 : $ nconvchan, nconvplane)
597 : integer, intent(in) :: supports(nconvplane)
598 :
599 : real :: norm
600 : real :: wt
601 0 : complex :: cfunc(convsize, convsize)
602 : logical :: onmosgrid
603 :
604 : integer :: iloc(2)
605 : integer :: iiloc(2)
606 : integer, intent(in) :: rbeg, rend
607 : integer :: ix, iy, iz, ipol, ichan, xind, yind
608 : integer :: apol, achan, aconvplane, irow
609 : integer :: aconvpol, aconvchan, xind2, yind2
610 : integer :: posx, posy, msupportx, msupporty, psupportx, psupporty
611 : logical :: centin
612 : integer :: support
613 :
614 :
615 :
616 0 : do irow=rbeg, rend
617 0 : aconvplane=abs(convplanemap(irow))+1
618 0 : support=supports(aconvplane)
619 0 : if(rflag(irow).eq.0) then
620 0 : do ichan=1, nvischan
621 0 : achan=chanmap(ichan)+1
622 0 : aconvchan=convchanmap(ichan)+1
623 0 : if((achan.ge.1).and.(achan.le.nchan).and.
624 0 : $ (weight(ichan,irow).ne.0.0)) then
625 0 : if (onmosgrid(loc(1, ichan, irow), nx, ny, 1, 1,
626 : $ nx, ny, support, msupportx, msupporty,
627 : $ psupportx, psupporty, centin)) then
628 0 : do ipol=1, nvispol
629 0 : apol=polmap(ipol)+1
630 0 : aconvpol=convpolmap(ipol)+1
631 : if((flag(ipol,ichan,irow).ne.1).and.
632 0 : $ (apol.ge.1).and.(apol.le.npol)) then
633 : C If we are making a PSF then we don't want to phase
634 : C rotate but we do want to reproject uvw
635 :
636 :
637 0 : nweight=cmplx(weight(ichan,irow))
638 : sumwt(apol, achan)=sumwt(apol, achan)+
639 0 : $ weight(ichan, irow)
640 :
641 : C norm will be the value we would get for the peak
642 : C at the phase center. We will want to normalize
643 : C the final image by this term.
644 0 : norm=0.0
645 : C write(*,*)off
646 0 : cfunc=convweight(:,:,aconvpol, aconvchan, aconvplane)
647 0 : do iy=msupporty, psupporty
648 0 : do ix=msupportx, psupportx
649 :
650 0 : xind2=sampling*ix+(convsize)/2+1
651 0 : yind2=sampling*iy+(convsize)/2+1
652 : cwt=cfunc(xind2,
653 0 : $ yind2)
654 :
655 0 : iiloc(1)=nx/2+1+ix
656 0 : iiloc(2)=ny/2+1+iy
657 : weightgrid(iiloc(1),iiloc(2),
658 : $ apol,achan)= weightgrid(
659 : $ iiloc(1),iiloc(2),apol,achan)
660 0 : $ + nweight*cwt
661 :
662 :
663 : end do
664 : end do
665 :
666 : end if
667 : end do
668 : end if
669 : end if
670 : end do
671 : end if
672 : end do
673 0 : return
674 : end
675 :
676 :
677 :
678 :
679 :
680 : C Single precision weight grid image...Damn you fortran...no templates
681 0 : subroutine gmoswgts (nvispol, nvischan,
682 0 : $ flag, rflag, weight, nrow,
683 : $ nx, ny, npol, nchan,
684 : $ support, convsize, sampling,
685 0 : $ chanmap, polmap,
686 0 : $ weightgrid, sumwt, convweight, convplanemap,
687 0 : $ convchanmap, convpolmap,
688 : $ nconvplane, nconvchan, nconvpol, rbeg,
689 0 : $ rend, loc, off, phasor)
690 :
691 : implicit none
692 : integer, intent(in) :: nx,ny,npol,nchan, nvispol, nvischan, nrow
693 :
694 :
695 : integer, intent(in) :: loc(2, nvischan, nrow)
696 : integer, intent(in) :: off(2, nvischan, nrow)
697 : complex, intent(in) :: phasor(nvischan, nrow)
698 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
699 : integer, intent(in) :: rflag(nrow)
700 : real, intent(in) :: weight(nvischan, nrow)
701 : integer, intent(in) :: support
702 : integer, intent(in) :: chanmap(nchan), polmap(npol)
703 : complex, intent(inout) :: weightgrid(nx, ny, npol, nchan)
704 : double precision, intent(inout) :: sumwt(npol, nchan)
705 : complex :: nweight
706 : integer, intent(in) :: convsize, sampling
707 : integer, intent(in) :: nconvplane, nconvchan, nconvpol
708 : integer, intent(in) :: convplanemap(nrow)
709 : integer, intent(in) :: convchanmap(nvischan)
710 : integer, intent(in) :: convpolmap(nvispol)
711 : complex :: cwt
712 : complex, intent(in) :: convweight(convsize, convsize, nconvpol,
713 : $ nconvchan, nconvplane)
714 :
715 :
716 : real :: norm
717 : real :: wt
718 :
719 : logical :: onmosgrid
720 :
721 : integer :: iloc(2)
722 : integer :: iiloc(2)
723 : integer, intent(in) :: rbeg, rend
724 : integer :: ix, iy, iz, ipol, ichan, xind, yind
725 : integer :: apol, achan, aconvplane, irow
726 : integer :: aconvpol, aconvchan, xind2, yind2
727 : integer :: posx, posy, msupportx, msupporty, psupportx, psupporty
728 : logical :: centin
729 :
730 :
731 :
732 :
733 0 : do irow=rbeg, rend
734 0 : aconvplane=convplanemap(irow)+1
735 0 : if(rflag(irow).eq.0) then
736 0 : do ichan=1, nvischan
737 0 : achan=chanmap(ichan)+1
738 0 : aconvchan=convchanmap(ichan)+1
739 0 : if((achan.ge.1).and.(achan.le.nchan).and.
740 0 : $ (weight(ichan,irow).ne.0.0)) then
741 0 : if (onmosgrid(loc(1, ichan, irow), nx, ny, 1, 1,
742 : $ nx, ny, support, msupportx, msupporty,
743 : $ psupportx, psupporty, centin)) then
744 0 : do ipol=1, nvispol
745 0 : apol=polmap(ipol)+1
746 0 : aconvpol=convpolmap(ipol)+1
747 : if((flag(ipol,ichan,irow).ne.1).and.
748 0 : $ (apol.ge.1).and.(apol.le.npol)) then
749 : C If we are making a PSF then we don't want to phase
750 : C rotate but we do want to reproject uvw
751 :
752 :
753 0 : nweight=cmplx(weight(ichan,irow))
754 : sumwt(apol, achan)=sumwt(apol, achan)+
755 0 : $ weight(ichan, irow)
756 :
757 : C norm will be the value we would get for the peak
758 : C at the phase center. We will want to normalize
759 : C the final image by this term.
760 0 : norm=0.0
761 : C write(*,*)off
762 0 : do iy=msupporty, psupporty
763 0 : do ix=msupportx, psupportx
764 :
765 0 : xind2=sampling*ix+(convsize)/2+1
766 0 : yind2=sampling*iy+(convsize)/2+1
767 : cwt=convweight(xind2,
768 0 : $ yind2, aconvpol, aconvchan, aconvplane)
769 0 : iiloc(1)=nx/2+1+ix
770 0 : iiloc(2)=ny/2+1+iy
771 : weightgrid(iiloc(1),iiloc(2),
772 : $ apol,achan)= weightgrid(
773 : $ iiloc(1),iiloc(2),apol,achan)
774 0 : $ + nweight*cwt
775 :
776 :
777 : end do
778 : end do
779 :
780 : end if
781 : end do
782 : end if
783 : end if
784 : end do
785 : end if
786 : end do
787 0 : return
788 : end
789 :
790 :
791 4262464 : subroutine sectgmosd2 (values, nvispol, nvischan,
792 4262464 : $ dopsf, flag, rflag, weight, nrow,
793 4262464 : $ grid, nx, ny, npol, nchan,
794 4262464 : $ support, convsize, sampling, convfunc,
795 4262464 : $ chanmap, polmap,
796 4262464 : $ sumwt, convplanemap,
797 4262464 : $ convchanmap, convpolmap,
798 : $ nconvplane, nconvchan, nconvpol, x0, y0, nxsub, nysub, rbeg,
799 4262464 : $ rend, loc, off, phasor)
800 :
801 : implicit none
802 : integer, intent(in) :: nx,ny,npol,nchan, nvispol, nvischan, nrow
803 : complex, intent(in) :: values(nvispol, nvischan, nrow)
804 : double complex, intent(inout) :: grid(nx, ny, npol, nchan)
805 :
806 : integer, intent(in) :: x0, y0, nxsub, nysub
807 : integer, intent(in) :: loc(2, nvischan, nrow)
808 : integer, intent(in) :: off(2, nvischan, nrow)
809 : complex, intent(in) :: phasor(nvischan, nrow)
810 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
811 : integer, intent(in) :: rflag(nrow)
812 : real, intent(in) :: weight(nvischan, nrow)
813 : double precision, intent(inout) :: sumwt(npol, nchan)
814 : integer, intent(in) :: support
815 : integer, intent(in) :: chanmap(nchan), polmap(npol)
816 : integer, intent(in) :: dopsf
817 :
818 : double complex :: nvalue
819 : integer, intent(in) :: convsize, sampling
820 : integer, intent(in) :: nconvplane, nconvchan, nconvpol
821 : integer, intent(in) :: convplanemap(nrow)
822 : integer, intent(in) :: convchanmap(nvischan)
823 : integer, intent(in) :: convpolmap(nvispol)
824 : complex, intent(in) :: convfunc(convsize, convsize, nconvpol,
825 : $ nconvchan, nconvplane)
826 : complex :: cwt
827 :
828 :
829 : real :: norm
830 : real :: wt
831 :
832 : logical :: onmosgrid
833 :
834 : integer :: iloc(2)
835 : integer :: iiloc(2)
836 : integer, intent(in) :: rbeg, rend
837 : integer :: ix, iy, iz, ipol, ichan, xind, yind
838 : integer :: apol, achan, aconvplane, irow
839 : integer :: aconvpol, aconvchan, xind2, yind2
840 : integer :: posx, posy, msupportx, msupporty, psupportx, psupporty
841 : logical :: centin
842 : logical :: doconj
843 :
844 :
845 :
846 1009468480 : do irow=rbeg, rend
847 : C sign of convplanemap determines if to use conjg
848 1005206016 : aconvplane=abs(convplanemap(irow))+1
849 1005206016 : doconj = (convplanemap(irow) < 0)
850 1009468480 : if(rflag(irow).eq.0) then
851 3354341632 : do ichan=1, nvischan
852 2358637824 : achan=chanmap(ichan)+1
853 2358637824 : aconvchan=convchanmap(ichan)+1
854 2358637824 : if((achan.ge.1).and.(achan.le.nchan).and.
855 995703808 : $ (weight(ichan,irow).ne.0.0)) then
856 2151514176 : if (onmosgrid(loc(1, ichan, irow), nx, ny, x0, y0,
857 : $ nxsub, nysub, support, msupportx, msupporty,
858 : $ psupportx, psupporty, centin)) then
859 352536435 : do ipol=1, nvispol
860 235024290 : apol=polmap(ipol)+1
861 235024290 : aconvpol=convpolmap(ipol)+1
862 : if((flag(ipol,ichan,irow).ne.1).and.
863 352536435 : $ (apol.ge.1).and.(apol.le.npol)) then
864 : C If we are making a PSF then we don't want to phase
865 : C rotate but we do want to reproject uvw
866 196732914 : if(dopsf.eq.1) then
867 102549672 : nvalue=cmplx(weight(ichan,irow))
868 : else
869 : nvalue=weight(ichan,irow)*
870 94183242 : $ (values(ipol,ichan,irow)*phasor(ichan, irow))
871 : end if
872 :
873 :
874 : C norm will be the value we would get for the peak
875 : C at the phase center. We will want to normalize
876 : C the final image by this term.
877 196732914 : norm=0.0
878 3019512640 : do iy=msupporty, psupporty
879 : iloc(2)=(sampling*iy)+
880 2822779726 : $ off(2, ichan, irow)
881 2822779726 : yind=iloc(2)+(convsize)/2+1
882 57634006122 : do ix=msupportx, psupportx
883 : iloc(1)=(sampling*ix)+
884 54614493482 : $ off(1, ichan, irow)
885 54614493482 : xind=iloc(1)+(convsize)/2+1
886 : cwt=convfunc(xind, yind,
887 54614493482 : $ aconvpol, aconvchan, aconvplane)
888 54614493482 : if(doconj) cwt=conjg(cwt)
889 : C write(*,*) support, iloc
890 : C write(*,*) loc(1, ichan, irow)+ix,loc(2, ichan, irow)+iy,xind,yind
891 : grid(loc(1, ichan, irow)+ix,
892 : $ loc(2, ichan, irow)+iy,apol,achan)=
893 : $ grid(loc(1, ichan, irow)+ix,
894 : $ loc(2, ichan, irow)+iy,apol,achan)+
895 57437273208 : $ nvalue*cwt
896 : end do
897 : end do
898 196732914 : if(centin) then
899 : sumwt(apol, achan)= sumwt(apol,achan)+
900 63348304 : $ weight(ichan,irow)
901 : endif
902 : end if
903 : end do
904 : C if onmos
905 : end if
906 : end if
907 : end do
908 : end if
909 : end do
910 4262464 : return
911 : end
912 : C Single Precision version
913 0 : subroutine sectgmoss2 (values, nvispol, nvischan,
914 0 : $ dopsf, flag, rflag, weight, nrow,
915 0 : $ grid, nx, ny, npol, nchan,
916 0 : $ support, convsize, sampling, convfunc,
917 0 : $ chanmap, polmap,
918 0 : $ sumwt, convplanemap,
919 0 : $ convchanmap, convpolmap,
920 : $ nconvplane, nconvchan, nconvpol, x0, y0, nxsub, nysub, rbeg,
921 0 : $ rend, loc, off, phasor)
922 :
923 : implicit none
924 : integer, intent(in) :: nx,ny,npol,nchan, nvispol, nvischan, nrow
925 : complex, intent(in) :: values(nvispol, nvischan, nrow)
926 : complex, intent(inout) :: grid(nx, ny, npol, nchan)
927 :
928 : integer, intent(in) :: x0, y0, nxsub, nysub
929 : integer, intent(in) :: loc(2, nvischan, nrow)
930 : integer, intent(in) :: off(2, nvischan, nrow)
931 : complex, intent(in) :: phasor(nvischan, nrow)
932 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
933 : integer, intent(in) :: rflag(nrow)
934 : real, intent(in) :: weight(nvischan, nrow)
935 : double precision, intent(inout) :: sumwt(npol, nchan)
936 : integer, intent(in) :: support
937 : integer, intent(in) :: chanmap(nchan), polmap(npol)
938 : integer, intent(in) :: dopsf
939 :
940 : complex :: nvalue
941 : integer, intent(in) :: convsize, sampling
942 : integer, intent(in) :: nconvplane, nconvchan, nconvpol
943 : integer, intent(in) :: convplanemap(nrow)
944 : integer, intent(in) :: convchanmap(nvischan)
945 : integer, intent(in) :: convpolmap(nvispol)
946 : complex, intent(in) :: convfunc(convsize, convsize, nconvpol,
947 : $ nconvchan, nconvplane)
948 : complex :: cwt
949 :
950 :
951 : real :: norm
952 : real :: wt
953 :
954 : logical :: onmosgrid
955 :
956 : integer :: iloc(2)
957 : integer :: iiloc(2)
958 : integer, intent(in) :: rbeg, rend
959 : integer :: ix, iy, iz, ipol, ichan, xind, yind
960 : integer :: apol, achan, aconvplane, irow
961 : integer :: aconvpol, aconvchan, xind2, yind2
962 : integer :: posx, posy, msupportx, msupporty, psupportx, psupporty
963 : logical :: centin
964 :
965 :
966 :
967 :
968 0 : do irow=rbeg, rend
969 0 : aconvplane=convplanemap(irow)+1
970 0 : if(rflag(irow).eq.0) then
971 0 : do ichan=1, nvischan
972 0 : achan=chanmap(ichan)+1
973 0 : aconvchan=convchanmap(ichan)+1
974 0 : if((achan.ge.1).and.(achan.le.nchan).and.
975 0 : $ (weight(ichan,irow).ne.0.0)) then
976 0 : if (onmosgrid(loc(1, ichan, irow), nx, ny, x0, y0,
977 : $ nxsub, nysub, support, msupportx, msupporty,
978 : $ psupportx, psupporty, centin)) then
979 0 : do ipol=1, nvispol
980 0 : apol=polmap(ipol)+1
981 0 : aconvpol=convpolmap(ipol)+1
982 : if((flag(ipol,ichan,irow).ne.1).and.
983 0 : $ (apol.ge.1).and.(apol.le.npol)) then
984 : C If we are making a PSF then we don't want to phase
985 : C rotate but we do want to reproject uvw
986 0 : if(dopsf.eq.1) then
987 0 : nvalue=cmplx(weight(ichan,irow))
988 : else
989 : nvalue=weight(ichan,irow)*
990 0 : $ (values(ipol,ichan,irow)*phasor(ichan, irow))
991 : end if
992 :
993 :
994 : C norm will be the value we would get for the peak
995 : C at the phase center. We will want to normalize
996 : C the final image by this term.
997 0 : norm=0.0
998 0 : do iy=msupporty, psupporty
999 : iloc(2)=(sampling*iy)+
1000 0 : $ off(2, ichan, irow)
1001 0 : yind=iloc(2)+(convsize)/2+1
1002 0 : do ix=msupportx, psupportx
1003 : iloc(1)=(sampling*ix)+
1004 0 : $ off(1, ichan, irow)
1005 0 : xind=iloc(1)+(convsize)/2+1
1006 : cwt=convfunc(xind, yind,
1007 0 : $ aconvpol, aconvchan, aconvplane)
1008 : C write(*,*) support, iloc
1009 : grid(loc(1, ichan, irow)+ix,
1010 : $ loc(2, ichan, irow)+iy,apol,achan)=
1011 : $ grid(loc(1, ichan, irow)+ix,
1012 : $ loc(2, ichan, irow)+iy,apol,achan)+
1013 0 : $ nvalue*cwt
1014 : end do
1015 : end do
1016 0 : if(centin) then
1017 : sumwt(apol, achan)= sumwt(apol,achan)+
1018 0 : $ weight(ichan,irow)
1019 : end if
1020 : end if
1021 : end do
1022 : end if
1023 : end if
1024 : end do
1025 : end if
1026 : end do
1027 0 : return
1028 : end
1029 : C
1030 :
1031 : C Same as sectgmosd2 except for varrying support across rows
1032 0 : subroutine sectgmosd3 (values, nvispol, nvischan,
1033 0 : $ dopsf, flag, rflag, weight, nrow,
1034 0 : $ grid, nx, ny, npol, nchan,
1035 0 : $ supports, convsize, sampling, convfunc,
1036 0 : $ chanmap, polmap,
1037 0 : $ sumwt, convplanemap,
1038 0 : $ convchanmap, convpolmap,
1039 : $ nconvplane, nconvchan, nconvpol, x0, y0, nxsub, nysub, rbeg,
1040 0 : $ rend, loc, off, phasor)
1041 :
1042 : implicit none
1043 : integer, intent(in) :: nx,ny,npol,nchan, nvispol, nvischan, nrow
1044 : complex, intent(in) :: values(nvispol, nvischan, nrow)
1045 : double complex, intent(inout) :: grid(nx, ny, npol, nchan)
1046 :
1047 : integer, intent(in) :: x0, y0, nxsub, nysub
1048 : integer, intent(in) :: loc(2, nvischan, nrow)
1049 : integer, intent(in) :: off(2, nvischan, nrow)
1050 : complex, intent(in) :: phasor(nvischan, nrow)
1051 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
1052 : integer, intent(in) :: rflag(nrow)
1053 : real, intent(in) :: weight(nvischan, nrow)
1054 : double precision, intent(inout) :: sumwt(npol, nchan)
1055 : integer, intent(in) :: chanmap(nchan), polmap(npol)
1056 : integer, intent(in) :: dopsf
1057 :
1058 : double complex :: nvalue
1059 : integer, intent(in) :: convsize, sampling
1060 : integer, intent(in) :: nconvplane, nconvchan, nconvpol
1061 : integer, intent(in) :: convplanemap(nrow)
1062 : integer, intent(in) :: convchanmap(nvischan)
1063 : integer, intent(in) :: convpolmap(nvispol)
1064 : complex, intent(in) :: convfunc(convsize, convsize, nconvpol,
1065 : $ nconvchan, nconvplane)
1066 : integer, intent(in) :: supports(nconvplane)
1067 : complex :: cwt
1068 0 : complex :: cfunc(convsize, convsize)
1069 : integer :: support
1070 : real :: norm
1071 : real :: wt
1072 :
1073 : logical :: onmosgrid
1074 :
1075 : integer :: iloc(2)
1076 : integer :: iiloc(2)
1077 : integer, intent(in) :: rbeg, rend
1078 : integer :: ix, iy, iz, ipol, ichan, xind, yind
1079 : integer :: apol, achan, aconvplane, irow
1080 : integer :: aconvpol, aconvchan, xind2, yind2
1081 : integer :: posx, posy, msupportx, msupporty, psupportx, psupporty
1082 : logical :: centin
1083 : logical :: doconj
1084 :
1085 :
1086 :
1087 0 : do irow=rbeg, rend
1088 : C sign of convplanemap determines if to use conjg
1089 0 : aconvplane=abs(convplanemap(irow))+1
1090 0 : support=supports(aconvplane)
1091 : C write(*,*) 'support', support, 'acpl', aconvplane, convsize
1092 0 : doconj = (convplanemap(irow) < 0)
1093 0 : if(rflag(irow).eq.0) then
1094 0 : do ichan=1, nvischan
1095 0 : achan=chanmap(ichan)+1
1096 0 : aconvchan=convchanmap(ichan)+1
1097 0 : if((achan.ge.1).and.(achan.le.nchan).and.
1098 0 : $ (weight(ichan,irow).ne.0.0)) then
1099 0 : if (onmosgrid(loc(1, ichan, irow), nx, ny, x0, y0,
1100 : $ nxsub, nysub, support, msupportx, msupporty,
1101 : $ psupportx, psupporty, centin)) then
1102 0 : do ipol=1, nvispol
1103 0 : apol=polmap(ipol)+1
1104 0 : aconvpol=convpolmap(ipol)+1
1105 : if((flag(ipol,ichan,irow).ne.1).and.
1106 0 : $ (apol.ge.1).and.(apol.le.npol)) then
1107 : C If we are making a PSF then we don't want to phase
1108 : C rotate but we do want to reproject uvw
1109 0 : if(dopsf.eq.1) then
1110 0 : nvalue=cmplx(weight(ichan,irow))
1111 : else
1112 : nvalue=weight(ichan,irow)*
1113 0 : $ (values(ipol,ichan,irow)*phasor(ichan, irow))
1114 : end if
1115 :
1116 0 : cfunc=convfunc(:,:, aconvpol, aconvchan, aconvplane)
1117 : C norm will be the value we would get for the peak
1118 : C at the phase center. We will want to normalize
1119 : C the final image by this term.
1120 0 : norm=0.0
1121 0 : do iy=msupporty, psupporty
1122 : iloc(2)=(sampling*iy)+
1123 0 : $ off(2, ichan, irow)
1124 0 : yind=iloc(2)+(convsize)/2+1
1125 0 : do ix=msupportx, psupportx
1126 : iloc(1)=(sampling*ix)+
1127 0 : $ off(1, ichan, irow)
1128 0 : xind=iloc(1)+(convsize)/2+1
1129 0 : cwt=cfunc(xind, yind)
1130 0 : if(doconj) cwt=conjg(cwt)
1131 : C write(*,*) support, iloc
1132 : C write(*,*) loc(1, ichan, irow)+ix,loc(2, ichan, irow)+iy,xind,yind
1133 : grid(loc(1, ichan, irow)+ix,
1134 : $ loc(2, ichan, irow)+iy,apol,achan)=
1135 : $ grid(loc(1, ichan, irow)+ix,
1136 : $ loc(2, ichan, irow)+iy,apol,achan)+
1137 0 : $ nvalue*cwt
1138 : end do
1139 : end do
1140 0 : if(centin) then
1141 : sumwt(apol, achan)= sumwt(apol,achan)+
1142 0 : $ weight(ichan,irow)
1143 : endif
1144 : end if
1145 : end do
1146 : C if onmos
1147 : end if
1148 : end if
1149 : end do
1150 : end if
1151 : end do
1152 0 : return
1153 : end
1154 :
1155 :
1156 :
1157 0 : subroutine gmoss (uvw, dphase, values, nvispol, nvischan,
1158 0 : $ dopsf, flag, rflag, weight, nrow, rownum,
1159 0 : $ scale, offset, grid, nx, ny, npol, nchan, freq, c,
1160 0 : $ support, convsize, sampling, convfunc,
1161 0 : $ chanmap, polmap,
1162 0 : $ sumwt, weightgrid, convweight, doweightgrid, convplanemap,
1163 : $ nconvplane)
1164 :
1165 : implicit none
1166 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
1167 : complex values(nvispol, nvischan, nrow)
1168 : complex grid(nx, ny, npol, nchan)
1169 :
1170 : double precision uvw(3, nrow), freq(nvischan), c, scale(2),
1171 : $ offset(2)
1172 : double precision dphase(nrow), uvdist
1173 : double precision xlast, ylast
1174 : complex phasor
1175 : integer flag(nvispol, nvischan, nrow)
1176 : integer rflag(nrow)
1177 : real weight(nvischan, nrow), phase
1178 : double precision sumwt(npol, nchan)
1179 : integer rownum
1180 : integer support
1181 : integer chanmap(nchan), polmap(npol)
1182 : integer dopsf
1183 : complex weightgrid(nx, ny, npol, nchan)
1184 : integer doweightgrid
1185 :
1186 : complex nvalue
1187 : complex nweight
1188 : integer convsize, sampling
1189 : integer nconvplane
1190 : integer convplanemap(nrow)
1191 : complex convfunc(convsize, convsize, nconvplane), cwt, crot
1192 : complex convweight(convsize, convsize, nconvplane)
1193 :
1194 :
1195 : complex sconv(-(support+1)*sampling:(support+1)*sampling,
1196 0 : $ -(support+1)*sampling:(support+1)*sampling, nconvplane)
1197 : complex sconv2(-(support+1)*sampling:(support+1)*sampling,
1198 0 : $ -(support+1)*sampling:(support+1)*sampling, nconvplane)
1199 : real sumsconv
1200 : real sumsconv2
1201 : real ratioofbeam
1202 :
1203 : real norm
1204 : real wt
1205 :
1206 : logical omos, doshift
1207 :
1208 : real pos(3)
1209 : integer loc(2), off(2), iloc(2)
1210 : integer iiloc(2)
1211 : integer rbeg, rend
1212 : integer ix, iy, iz, ipol, ichan
1213 : integer apol, achan, aconvplane, irow
1214 : double precision pi
1215 : data pi/3.14159265358979323846/
1216 : integer sampsupp
1217 0 : sampsupp=(support+1)*sampling
1218 0 : irow=rownum
1219 :
1220 0 : sumsconv=0
1221 0 : sumsconv2=0
1222 0 : ratioofbeam=1.0
1223 0 : do iz=1, nconvplane
1224 0 : do iy=-sampsupp,sampsupp
1225 0 : iloc(2)=iy+convsize/2+1
1226 0 : do ix=-sampsupp,sampsupp
1227 0 : iloc(1)=ix+convsize/2+1
1228 0 : sconv(ix,iy,iz)=(convfunc(iloc(1), iloc(2),iz))
1229 0 : sconv2(ix,iy,iz)=convweight(iloc(1), iloc(2),iz)
1230 : end do
1231 : end do
1232 : end do
1233 0 : doshift=.FALSE.
1234 :
1235 0 : if(irow.ge.0) then
1236 0 : rbeg=irow+1
1237 0 : rend=irow+1
1238 : else
1239 0 : rbeg=1
1240 0 : rend=nrow
1241 : end if
1242 :
1243 0 : xlast=0.0
1244 0 : ylast=0.0
1245 0 : do irow=rbeg, rend
1246 0 : aconvplane=convplanemap(irow)+1
1247 0 : if(rflag(irow).eq.0) then
1248 0 : do ichan=1, nvischan
1249 0 : achan=chanmap(ichan)+1
1250 0 : if((achan.ge.1).and.(achan.le.nchan).and.
1251 0 : $ (weight(ichan,irow).ne.0.0)) then
1252 : call smos(uvw(1,irow), dphase(irow), freq(ichan), c,
1253 0 : $ scale, offset, sampling, pos, loc, off, phasor)
1254 0 : if (omos(nx, ny, loc, support)) then
1255 0 : do ipol=1, nvispol
1256 0 : apol=polmap(ipol)+1
1257 : if((flag(ipol,ichan,irow).ne.1).and.
1258 0 : $ (apol.ge.1).and.(apol.le.npol)) then
1259 : C If we are making a PSF then we don't want to phase
1260 : C rotate but we do want to reproject uvw
1261 0 : if(dopsf.eq.1) then
1262 0 : nvalue=cmplx(weight(ichan,irow))
1263 : else
1264 : nvalue=weight(ichan,irow)*
1265 0 : $ (values(ipol,ichan,irow)*phasor)
1266 : end if
1267 0 : if(doweightgrid .gt. 0) then
1268 0 : nweight=cmplx(weight(ichan,irow))
1269 : end if
1270 :
1271 : C norm will be the value we would get for the peak
1272 : C at the phase center. We will want to normalize
1273 : C the final image by this term.
1274 0 : norm=0.0
1275 0 : if(sampling.eq.1) then
1276 0 : do iy=-support,support
1277 0 : do ix=-support,support
1278 : grid(loc(1)+ix,
1279 : $ loc(2)+iy,apol,achan)=
1280 : $ grid(loc(1)+ix,
1281 : $ loc(2)+iy,apol,achan)+
1282 0 : $ nvalue*sconv(ix,iy, aconvplane)
1283 :
1284 0 : if(doweightgrid .gt. 0) then
1285 0 : iloc(1)=nx/2+1+ix
1286 0 : iloc(2)=ny/2+1+iy
1287 : weightgrid(iloc(1),iloc(2),
1288 : $ apol,achan)= weightgrid(
1289 : $ iloc(1),iloc(2),apol,achan)
1290 0 : $ + nweight*sconv2(ix,iy,aconvplane)
1291 :
1292 : end if
1293 : end do
1294 : end do
1295 : else
1296 0 : do iy=-support,support
1297 0 : iloc(2)=(sampling*iy+off(2))+1
1298 0 : do ix=-support, support
1299 0 : iloc(1)=(sampling*ix+off(1))+1
1300 : cwt=sconv(iloc(1),
1301 0 : $ iloc(2),aconvplane)
1302 : grid(loc(1)+ix,
1303 : $ loc(2)+iy,apol,achan)=
1304 : $ grid(loc(1)+ix,
1305 : $ loc(2)+iy,apol,achan)+
1306 0 : $ nvalue*cwt
1307 0 : if(doweightgrid .gt. 0) then
1308 : cwt=sconv2(iloc(1), iloc(2),
1309 0 : $ aconvplane)
1310 0 : iiloc(1)=nx/2+1+ix
1311 0 : iiloc(2)=ny/2+1+iy
1312 : weightgrid(iiloc(1),iiloc(2),
1313 : $ apol,achan)= weightgrid(
1314 : $ iiloc(1),iiloc(2),apol,achan)
1315 0 : $ + nweight*cwt
1316 :
1317 : end if
1318 : end do
1319 : end do
1320 : end if
1321 : sumwt(apol, achan)= sumwt(apol,achan)+
1322 0 : $ weight(ichan,irow)
1323 : end if
1324 : end do
1325 : end if
1326 : end if
1327 : end do
1328 : end if
1329 : end do
1330 0 : return
1331 : end
1332 : C
1333 0 : subroutine gmoss2 (uvw, dphase, values, nvispol, nvischan,
1334 0 : $ dopsf, flag, rflag, weight, nrow, rownum,
1335 0 : $ scale, offset, grid, nx, ny, npol, nchan, freq, c,
1336 0 : $ support, convsize, sampling, convfunc,
1337 0 : $ chanmap, polmap,
1338 0 : $ sumwt, weightgrid, convweight, doweightgrid, convplanemap,
1339 0 : $ convchanmap, convpolmap,
1340 : $ nconvplane, nconvchan, nconvpol)
1341 :
1342 : implicit none
1343 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
1344 : complex values(nvispol, nvischan, nrow)
1345 : complex grid(nx, ny, npol, nchan)
1346 :
1347 : double precision uvw(3, nrow), freq(nvischan), c, scale(2),
1348 : $ offset(2)
1349 : double precision dphase(nrow), uvdist
1350 : double precision xlast, ylast
1351 : complex phasor
1352 : integer flag(nvispol, nvischan, nrow)
1353 : integer rflag(nrow)
1354 : real weight(nvischan, nrow), phase
1355 : double precision sumwt(npol, nchan)
1356 : integer rownum
1357 : integer support
1358 : integer chanmap(nchan), polmap(npol)
1359 : integer dopsf
1360 : complex weightgrid(nx, ny, npol, nchan)
1361 : integer doweightgrid
1362 :
1363 : double complex nvalue
1364 : double complex nweight
1365 : integer convsize, sampling
1366 : integer nconvplane, nconvchan, nconvpol
1367 : integer convplanemap(nrow)
1368 : integer convchanmap(nvischan)
1369 : integer convpolmap(nvispol)
1370 : complex convfunc(convsize, convsize, nconvpol, nconvchan,
1371 : $ nconvplane), cwt, crot
1372 : complex convweight(convsize, convsize, nconvpol, nconvchan,
1373 : $ nconvplane)
1374 : real sumsconv
1375 : real sumsconv2
1376 : real ratioofbeam
1377 :
1378 : real norm
1379 : real wt
1380 :
1381 : logical omos, doshift
1382 :
1383 : real pos(3)
1384 : integer loc(2), off(2), iloc(2)
1385 : integer iiloc(2)
1386 : integer rbeg, rend
1387 : integer ix, iy, iz, ipol, ichan, xind, yind
1388 : integer apol, achan, aconvplane, irow
1389 : integer aconvpol, aconvchan, xind2, yind2
1390 : double precision pi
1391 : data pi/3.14159265358979323846/
1392 : real maxsconv2, minsconv2
1393 0 : irow=rownum
1394 :
1395 0 : sumsconv=0
1396 0 : sumsconv2=0
1397 0 : ratioofbeam=1.0
1398 :
1399 0 : doshift=.FALSE.
1400 :
1401 0 : if(irow.ge.0) then
1402 0 : rbeg=irow+1
1403 0 : rend=irow+1
1404 : else
1405 0 : rbeg=1
1406 0 : rend=nrow
1407 : end if
1408 :
1409 : C write(*,*) 'max of sconvs ', maxsconv2, minsconv2, sampsupp,
1410 : C $ convsize
1411 :
1412 :
1413 :
1414 0 : xlast=0.0
1415 0 : ylast=0.0
1416 0 : do irow=rbeg, rend
1417 0 : aconvplane=convplanemap(irow)+1
1418 0 : if(rflag(irow).eq.0) then
1419 0 : do ichan=1, nvischan
1420 0 : achan=chanmap(ichan)+1
1421 0 : aconvchan=convchanmap(ichan)+1
1422 0 : if((achan.ge.1).and.(achan.le.nchan).and.
1423 0 : $ (weight(ichan,irow).ne.0.0)) then
1424 : call smos(uvw(1,irow), dphase(irow), freq(ichan), c,
1425 0 : $ scale, offset, sampling, pos, loc, off, phasor)
1426 0 : if (omos(nx, ny, loc, support)) then
1427 0 : do ipol=1, nvispol
1428 0 : apol=polmap(ipol)+1
1429 0 : aconvpol=convpolmap(ipol)+1
1430 : if((flag(ipol,ichan,irow).ne.1).and.
1431 0 : $ (apol.ge.1).and.(apol.le.npol)) then
1432 : C If we are making a PSF then we don't want to phase
1433 : C rotate but we do want to reproject uvw
1434 0 : if(dopsf.eq.1) then
1435 0 : nvalue=cmplx(weight(ichan,irow))
1436 : else
1437 : nvalue=weight(ichan,irow)*
1438 0 : $ (values(ipol,ichan,irow)*phasor)
1439 : end if
1440 0 : if(doweightgrid .gt. 0) then
1441 0 : nweight=cmplx(weight(ichan,irow))
1442 : end if
1443 :
1444 : C norm will be the value we would get for the peak
1445 : C at the phase center. We will want to normalize
1446 : C the final image by this term.
1447 0 : norm=0.0
1448 0 : if(sampling.eq.1) then
1449 0 : do iy=-support,support
1450 0 : yind=iy+(convsize)/2+1
1451 0 : do ix=-support,support
1452 0 : xind=ix+(convsize)/2+1
1453 : grid(loc(1)+ix,
1454 : $ loc(2)+iy,apol,achan)=
1455 : $ grid(loc(1)+ix,
1456 : $ loc(2)+iy,apol,achan)+
1457 : $ nvalue*convfunc(xind, yind,
1458 0 : $ aconvpol, aconvchan, aconvplane)
1459 :
1460 0 : if(doweightgrid .gt. 0) then
1461 0 : iloc(1)=nx/2+1+ix
1462 0 : iloc(2)=ny/2+1+iy
1463 : weightgrid(iloc(1),iloc(2),
1464 : $ apol,achan)= weightgrid(
1465 : $ iloc(1),iloc(2),apol,achan)
1466 : $ + nweight*convweight(xind, yind,
1467 0 : $ aconvpol, aconvchan, aconvplane)
1468 :
1469 : end if
1470 : end do
1471 : end do
1472 : else
1473 : C write(*,*)off
1474 0 : do iy=-support, support
1475 0 : iloc(2)=(sampling*iy)+off(2)
1476 0 : yind=iloc(2)+(convsize)/2+1
1477 0 : do ix=-support, support
1478 0 : iloc(1)=(sampling*ix)+off(1)
1479 0 : xind=iloc(1)+(convsize)/2+1
1480 : cwt=convfunc(xind, yind,
1481 0 : $ aconvpol, aconvchan, aconvplane)
1482 : C write(*,*) yind, xind, loc(1)+ix, loc(2)+iy, apol, achan,
1483 : C $ aconvpol, aconvchan, aconvplane
1484 : grid(loc(1)+ix,
1485 : $ loc(2)+iy,apol,achan)=
1486 : $ grid(loc(1)+ix,
1487 : $ loc(2)+iy,apol,achan)+
1488 0 : $ nvalue*cwt
1489 0 : if(doweightgrid .gt. 0) then
1490 0 : xind2=sampling*ix+(convsize)/2+1
1491 0 : yind2=sampling*iy+(convsize)/2+1
1492 : cwt=convweight(xind2,
1493 : $ yind2, aconvpol, aconvchan,
1494 0 : $ aconvplane)
1495 0 : iiloc(1)=nx/2+1+ix
1496 0 : iiloc(2)=ny/2+1+iy
1497 : weightgrid(iiloc(1),iiloc(2),
1498 : $ apol,achan)= weightgrid(
1499 : $ iiloc(1),iiloc(2),apol,achan)
1500 0 : $ + nweight*cwt
1501 :
1502 : end if
1503 : end do
1504 : end do
1505 : end if
1506 : sumwt(apol, achan)= sumwt(apol,achan)+
1507 0 : $ weight(ichan,irow)
1508 : end if
1509 : end do
1510 : end if
1511 : end if
1512 : end do
1513 : end if
1514 : end do
1515 0 : return
1516 : end
1517 : C
1518 :
1519 :
1520 :
1521 : C Degrid a number of visibility records
1522 : C
1523 0 : subroutine dmos (uvw, dphase, values, nvispol, nvischan,
1524 0 : $ flag, rflag,
1525 0 : $ nrow, rownum, scale, offset, grid, nx, ny, npol, nchan, freq,
1526 0 : $ c, support, convsize, sampling, convfunc,
1527 0 : $ chanmap, polmap, convplanemap, nconvplane)
1528 :
1529 : implicit none
1530 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
1531 : integer nconvplane
1532 : complex values(nvispol, nvischan, nrow)
1533 : complex grid(nx, ny, npol, nchan)
1534 : double precision uvw(3, nrow), freq(nvischan), c, scale(2),
1535 : $ offset(2)
1536 : double precision dphase(nrow), uvdist
1537 : double precision xlast, ylast
1538 : complex phasor
1539 : integer flag(nvispol, nvischan, nrow)
1540 : integer rflag(nrow)
1541 : integer rownum
1542 : integer support
1543 : integer chanmap(nchan), polmap(npol)
1544 : integer convplanemap(nrow)
1545 : complex nvalue
1546 :
1547 : integer convsize, sampling
1548 : complex convfunc(convsize, convsize, nconvplane), cwt, crot
1549 :
1550 : integer sampsupp
1551 :
1552 : complex sconv(-(support+1)*sampling:(support+1)*sampling,
1553 0 : $ -(support+1)*sampling:(support+1)*sampling, nconvplane)
1554 :
1555 : real norm, phase
1556 :
1557 : logical omos, doshift
1558 :
1559 : real pos(2)
1560 : integer loc(2), off(2), iloc(2)
1561 : integer rbeg, rend
1562 : integer ix, iy, iz, ipol, ichan
1563 : integer apol, achan, aconvplane, irow
1564 : real wt, wtx, wty
1565 : double precision pi
1566 : data pi/3.14159265358979323846/
1567 :
1568 0 : sampsupp=(support+1)*sampling
1569 0 : irow=rownum
1570 :
1571 0 : do iz=1, nconvplane
1572 0 : do iy=-sampsupp,sampsupp
1573 0 : iloc(2)=iy+convsize/2+1
1574 0 : do ix=-sampsupp,sampsupp
1575 0 : iloc(1)=ix+convsize/2+1
1576 0 : sconv(ix,iy,iz)=conjg(convfunc(iloc(1), iloc(2),iz))
1577 : end do
1578 : end do
1579 : end do
1580 0 : doshift=.FALSE.
1581 :
1582 0 : if(irow.ge.0) then
1583 0 : rbeg=irow+1
1584 0 : rend=irow+1
1585 : else
1586 0 : rbeg=1
1587 0 : rend=nrow
1588 : end if
1589 : C
1590 0 : xlast=0.0
1591 0 : ylast=0.0
1592 0 : do irow=rbeg, rend
1593 0 : aconvplane=convplanemap(irow)+1
1594 0 : if(rflag(irow).eq.0) then
1595 0 : do ichan=1, nvischan
1596 0 : achan=chanmap(ichan)+1
1597 0 : if((achan.ge.1).and.(achan.le.nchan)) then
1598 : call smos(uvw(1,irow), dphase(irow), freq(ichan), c,
1599 0 : $ scale, offset, sampling, pos, loc, off, phasor)
1600 0 : if (omos(nx, ny, loc, support)) then
1601 0 : do ipol=1, nvispol
1602 0 : apol=polmap(ipol)+1
1603 : if((flag(ipol,ichan,irow).ne.1).and.
1604 0 : $ (apol.ge.1).and.(apol.le.npol)) then
1605 0 : nvalue=0.0
1606 0 : norm=0.0
1607 0 : if(sampling.eq.1) then
1608 0 : do iy=-support,support
1609 0 : do ix=-support,support
1610 : nvalue=nvalue+(sconv(ix,iy,
1611 : $ aconvplane))*
1612 : $ grid(loc(1)+ix,loc(2)+iy,
1613 0 : $ apol,achan)
1614 : end do
1615 : end do
1616 : else
1617 0 : do iy=-support,support
1618 0 : iloc(2)=sampling*iy+off(2)
1619 0 : do ix=-support,support
1620 : iloc(1)=ix*sampling
1621 0 : $ +off(1)
1622 : cwt=sconv(iloc(1),
1623 0 : $ iloc(2),aconvplane)
1624 : nvalue=nvalue+cwt*
1625 : $ grid(loc(1)+ix,loc(2)+iy,
1626 0 : $ apol,achan)
1627 : end do
1628 : end do
1629 : end if
1630 : values(ipol,ichan,irow)=nvalue*conjg(
1631 0 : $ phasor)
1632 : end if
1633 : end do
1634 : end if
1635 : end if
1636 : end do
1637 : end if
1638 : end do
1639 0 : return
1640 : end
1641 : C
1642 :
1643 : C Degrid a number of visibility records
1644 : C
1645 0 : subroutine dmos2 (uvw, dphase, values, nvispol, nvischan,
1646 0 : $ flag, rflag,
1647 0 : $ nrow, rownum, scale, offset, grid, nx, ny, npol, nchan, freq,
1648 0 : $ c, support, convsize, sampling, convfunc,
1649 0 : $ chanmap, polmap, convplanemap, convchanmap, convpolmap,
1650 : $ nconvplane, nconvchan, nconvpol)
1651 :
1652 : implicit none
1653 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
1654 : integer nconvplane, nconvchan, nconvpol
1655 : complex values(nvispol, nvischan, nrow)
1656 : complex grid(nx, ny, npol, nchan)
1657 : double precision uvw(3, nrow), freq(nvischan), c, scale(2),
1658 : $ offset(2)
1659 : double precision dphase(nrow), uvdist
1660 : double precision xlast, ylast
1661 : complex phasor
1662 : integer flag(nvispol, nvischan, nrow)
1663 : integer rflag(nrow)
1664 : integer rownum
1665 : integer support
1666 : integer chanmap(nchan), polmap(npol)
1667 : integer convplanemap(nrow), convchanmap(nvischan)
1668 : integer convpolmap(nvispol)
1669 : complex nvalue
1670 :
1671 : integer convsize, sampling
1672 : complex convfunc(convsize, convsize, nconvpol, nconvchan,
1673 : $ nconvplane), cwt, crot
1674 :
1675 : integer sampsupp
1676 :
1677 : C complex sconv(-(support+1)*sampling:(support+1)*sampling,
1678 : C $ -(support+1)*sampling:(support+1)*sampling, nconvplane)
1679 :
1680 : real norm, phase
1681 :
1682 : logical omos, doshift
1683 :
1684 : real pos(2)
1685 : integer loc(2), off(2), iloc(2)
1686 : integer rbeg, rend
1687 : integer ix, iy, iz, ipol, ichan, xind, yind
1688 : integer apol, achan, aconvplane, irow
1689 : integer aconvchan, aconvpol
1690 : real wt, wtx, wty
1691 : double precision pi
1692 : data pi/3.14159265358979323846/
1693 :
1694 0 : sampsupp=(support+1)*sampling
1695 0 : irow=rownum
1696 :
1697 : C do iz=1, nconvplane
1698 : C do iy=-sampsupp,sampsupp
1699 : C iloc(2)=iy+convsize/2+1
1700 : c do ix=-sampsupp,sampsupp
1701 : C iloc(1)=ix+convsize/2+1
1702 : C sconv(ix,iy,iz)=conjg(convfunc(iloc(1), iloc(2),iz))
1703 : C end do
1704 : C end do
1705 : C end do
1706 0 : doshift=.FALSE.
1707 :
1708 0 : if(irow.ge.0) then
1709 0 : rbeg=irow+1
1710 0 : rend=irow+1
1711 : else
1712 0 : rbeg=1
1713 0 : rend=nrow
1714 : end if
1715 : C
1716 0 : xlast=0.0
1717 0 : ylast=0.0
1718 0 : do irow=rbeg, rend
1719 0 : aconvplane=convplanemap(irow)+1
1720 0 : if(rflag(irow).eq.0) then
1721 0 : do ichan=1, nvischan
1722 0 : achan=chanmap(ichan)+1
1723 0 : aconvchan=convchanmap(ichan)+1
1724 0 : if((achan.ge.1).and.(achan.le.nchan)) then
1725 : call smos(uvw(1,irow), dphase(irow), freq(ichan), c,
1726 0 : $ scale, offset, sampling, pos, loc, off, phasor)
1727 0 : if (omos(nx, ny, loc, support)) then
1728 0 : do ipol=1, nvispol
1729 0 : apol=polmap(ipol)+1
1730 0 : aconvpol=convpolmap(ipol)+1
1731 : if((flag(ipol,ichan,irow).ne.1).and.
1732 0 : $ (apol.ge.1).and.(apol.le.npol)) then
1733 : C write(*,*) 'aindices', aconvplane, aconvchan, aconvpol
1734 0 : nvalue=0.0
1735 0 : norm=0.0
1736 0 : if(sampling.eq.1) then
1737 0 : do iy=-support,support
1738 0 : yind=iy+(convsize)/2+1
1739 0 : do ix=-support,support
1740 0 : xind=ix+(convsize)/2+1
1741 : nvalue=nvalue+(convfunc(xind,yind,
1742 : $ aconvpol, aconvchan,
1743 : $ aconvplane))*
1744 : $ grid(loc(1)+ix,loc(2)+iy,
1745 0 : $ apol,achan)
1746 : end do
1747 : end do
1748 : else
1749 0 : do iy=-support,support
1750 0 : iloc(2)=sampling*iy+off(2)
1751 0 : yind=iloc(2)+(convsize)/2+1
1752 : C write(*,*) 'iloc(2)', iloc(2), off(2), yind
1753 0 : do ix=-support,support
1754 : iloc(1)=ix*sampling
1755 0 : $ +off(1)
1756 0 : xind=iloc(1)+(convsize)/2+1
1757 : C write(*,*) 'iloc(1)', iloc(1), off(1), xind
1758 : cwt=convfunc(xind, yind, aconvpol,
1759 0 : $ aconvchan,aconvplane)
1760 : nvalue=nvalue+cwt*
1761 : $ grid(loc(1)+ix,loc(2)+iy,
1762 0 : $ apol,achan)
1763 : end do
1764 : end do
1765 : end if
1766 : values(ipol,ichan,irow)=nvalue*conjg(
1767 0 : $ phasor)
1768 : end if
1769 : end do
1770 : end if
1771 : end if
1772 : end do
1773 : end if
1774 : end do
1775 0 : return
1776 : end
1777 : C
1778 :
1779 : C
1780 167172 : subroutine sectdmos2 (values, nvispol, nvischan,
1781 167172 : $ flag, rflag,
1782 167172 : $ nrow, grid, nx, ny, npol, nchan,
1783 167172 : $ support, convsize, sampling, convfunc,
1784 167172 : $ chanmap, polmap, convplanemap, convchanmap, convpolmap,
1785 167172 : $ nconvplane, nconvchan, nconvpol, rbeg,rend,loc,off,phasor)
1786 :
1787 : implicit none
1788 : integer, intent(in) :: nx, ny,npol,nchan,nvispol, nvischan, nrow
1789 : integer, intent(in) :: nconvplane, nconvchan, nconvpol
1790 : complex, intent(inout) :: values(nvispol, nvischan, nrow)
1791 : complex, intent(in) :: grid(nx, ny, npol, nchan)
1792 : complex, intent(in) :: phasor(nvischan, nrow)
1793 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
1794 : integer, intent(in) :: rflag(nrow)
1795 : integer, intent(in) :: support
1796 : integer, intent(in) :: chanmap(*), polmap(*)
1797 : integer, intent(in) :: convplanemap(nrow), convchanmap(nvischan)
1798 : integer, intent(in) :: convpolmap(nvispol)
1799 : complex :: nvalue
1800 :
1801 : integer, intent(in) :: convsize, sampling
1802 : complex, intent(in) :: convfunc(convsize, convsize, nconvpol,
1803 : $ nconvchan, nconvplane)
1804 : complex :: cwt, crot
1805 :
1806 : C complex sconv(-(support+1)*sampling:(support+1)*sampling,
1807 : C $ -(support+1)*sampling:(support+1)*sampling, nconvplane)
1808 :
1809 : real :: norm, phase
1810 :
1811 : logical :: omos
1812 : logical :: doconj
1813 :
1814 : integer, intent(in) :: loc(2, nvischan, nrow),
1815 : $ off(2,nvischan,nrow)
1816 : integer :: iloc(2)
1817 : integer, intent(in) :: rbeg, rend
1818 : integer :: ix, iy, iz, ipol, ichan, xind, yind
1819 : integer :: apol, achan, aconvplane, irow
1820 : integer :: aconvchan, aconvpol
1821 : real :: wt, wtx, wty
1822 :
1823 :
1824 : C write(*,*) 'polmap, ', polmap
1825 : C write(*,*) 'convpm,', convpolmap
1826 : C write(*,*) 'chanmp,', chanmap
1827 : C write(*,*) 'convcm,', convchanmap
1828 :
1829 3129054 : do irow=rbeg, rend
1830 2961882 : aconvplane=abs(convplanemap(irow))+1
1831 2961882 : doconj=(convplanemap(irow) < 0)
1832 3129054 : if(rflag(irow).eq.0) then
1833 9388050 : do ichan=1, nvischan
1834 6473874 : achan=chanmap(ichan)+1
1835 6473874 : aconvchan=convchanmap(ichan)+1
1836 9388050 : if((achan.ge.1).and.(achan.le.nchan)) then
1837 6473874 : if (omos(nx, ny, loc(1, ichan,irow),support)) then
1838 19410888 : do ipol=1, nvispol
1839 12940592 : apol=polmap(ipol)+1
1840 12940592 : aconvpol=convpolmap(ipol)+1
1841 : if((flag(ipol,ichan,irow).ne.1).and.
1842 19410888 : $ (apol.ge.1).and.(apol.le.npol)) then
1843 : C write(*,*) 'aindices', aconvplane, aconvchan, aconvpol
1844 12849548 : nvalue=0.0
1845 12849548 : norm=0.0
1846 :
1847 274323984 : do iy=-support,support
1848 261474436 : iloc(2)=sampling*iy+off(2, ichan, irow)
1849 261474436 : yind=iloc(2)+(convsize)/2+1
1850 : C write(*,*) 'iloc(2)', iloc(2), off(2), yind
1851 8127945004 : do ix=-support,support
1852 : iloc(1)=ix*sampling
1853 7853621020 : $ +off(1, ichan, irow)
1854 7853621020 : xind=iloc(1)+(convsize)/2+1
1855 : C write(*,*) 'iloc(1)', iloc(1), off(1), xind
1856 : cwt=convfunc(xind, yind, aconvpol,
1857 7853621020 : $ aconvchan,aconvplane)
1858 7853621020 : if(doconj) cwt=conjg(cwt)
1859 : nvalue=nvalue+cwt*
1860 : $ grid(loc(1, ichan, irow)+ix,
1861 : $ loc(2, ichan, irow)+iy,
1862 8115095456 : $ apol,achan)
1863 : end do
1864 : end do
1865 :
1866 : values(ipol,ichan,irow)=nvalue*conjg(
1867 12849548 : $ phasor(ichan, irow))
1868 : end if
1869 : end do
1870 : end if
1871 : end if
1872 : end do
1873 : end if
1874 : end do
1875 167172 : return
1876 : end
1877 : C
1878 : C same as sectdmos2 except with varying support
1879 0 : subroutine sectdmos3 (values, nvispol, nvischan,
1880 0 : $ flag, rflag,
1881 0 : $ nrow, grid, nx, ny, npol, nchan,
1882 0 : $ supports, convsize, sampling, convfunc,
1883 0 : $ chanmap, polmap, convplanemap, convchanmap, convpolmap,
1884 0 : $ nconvplane, nconvchan, nconvpol, rbeg,rend,loc,off,phasor)
1885 :
1886 : implicit none
1887 : integer, intent(in) :: nx, ny,npol,nchan,nvispol, nvischan, nrow
1888 : integer, intent(in) :: nconvplane, nconvchan, nconvpol
1889 : complex, intent(inout) :: values(nvispol, nvischan, nrow)
1890 : complex, intent(in) :: grid(nx, ny, npol, nchan)
1891 : complex, intent(in) :: phasor(nvischan, nrow)
1892 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
1893 : integer, intent(in) :: rflag(nrow)
1894 : integer, intent(in) :: chanmap(*), polmap(*)
1895 : integer, intent(in) :: convplanemap(nrow), convchanmap(nvischan)
1896 : integer, intent(in) :: convpolmap(nvispol)
1897 : complex :: nvalue
1898 :
1899 : integer, intent(in) :: convsize, sampling
1900 : complex, intent(in) :: convfunc(convsize, convsize, nconvpol,
1901 : $ nconvchan, nconvplane)
1902 : integer, intent(in) :: supports(nconvplane)
1903 : complex :: cwt, crot
1904 :
1905 : C complex sconv(-(support+1)*sampling:(support+1)*sampling,
1906 : C $ -(support+1)*sampling:(support+1)*sampling, nconvplane)
1907 0 : complex cfunc(convsize, convsize)
1908 : real :: norm, phase
1909 : integer :: support
1910 : logical :: omos
1911 : logical :: doconj
1912 :
1913 : integer, intent(in) :: loc(2, nvischan, nrow),
1914 : $ off(2,nvischan,nrow)
1915 : integer :: iloc(2)
1916 : integer, intent(in) :: rbeg, rend
1917 : integer :: ix, iy, iz, ipol, ichan, xind, yind
1918 : integer :: apol, achan, aconvplane, irow
1919 : integer :: aconvchan, aconvpol
1920 : real :: wt, wtx, wty
1921 :
1922 :
1923 : C write(*,*) 'polmap, ', polmap
1924 : C write(*,*) 'convpm,', convpolmap
1925 : C write(*,*) 'chanmp,', chanmap
1926 : C write(*,*) 'convcm,', convchanmap
1927 :
1928 0 : do irow=rbeg, rend
1929 0 : aconvplane=abs(convplanemap(irow))+1
1930 0 : support=supports(aconvplane)
1931 0 : doconj=(convplanemap(irow) < 0)
1932 0 : if(rflag(irow).eq.0) then
1933 0 : do ichan=1, nvischan
1934 0 : achan=chanmap(ichan)+1
1935 0 : aconvchan=convchanmap(ichan)+1
1936 0 : if((achan.ge.1).and.(achan.le.nchan)) then
1937 0 : if (omos(nx, ny, loc(1, ichan,irow),support)) then
1938 0 : do ipol=1, nvispol
1939 0 : apol=polmap(ipol)+1
1940 0 : aconvpol=convpolmap(ipol)+1
1941 : if((flag(ipol,ichan,irow).ne.1).and.
1942 0 : $ (apol.ge.1).and.(apol.le.npol)) then
1943 : C write(*,*) 'aindices', aconvplane, aconvchan, aconvpol
1944 0 : nvalue=0.0
1945 0 : norm=0.0
1946 0 : cfunc=convfunc(:,:,aconvpol, aconvchan, aconvplane)
1947 0 : do iy=-support,support
1948 0 : iloc(2)=sampling*iy+off(2, ichan, irow)
1949 0 : yind=iloc(2)+(convsize)/2+1
1950 : C write(*,*) 'iloc(2)', iloc(2), off(2), yind
1951 0 : do ix=-support,support
1952 : iloc(1)=ix*sampling
1953 0 : $ +off(1, ichan, irow)
1954 0 : xind=iloc(1)+(convsize)/2+1
1955 : C write(*,*) 'iloc(1)', iloc(1), off(1), xind
1956 0 : cwt=cfunc(xind, yind)
1957 0 : if(doconj) cwt=conjg(cwt)
1958 : nvalue=nvalue+cwt*
1959 : $ grid(loc(1, ichan, irow)+ix,
1960 : $ loc(2, ichan, irow)+iy,
1961 0 : $ apol,achan)
1962 : end do
1963 : end do
1964 :
1965 : values(ipol,ichan,irow)=nvalue*conjg(
1966 0 : $ phasor(ichan, irow))
1967 : end if
1968 : end do
1969 : end if
1970 : end if
1971 : end do
1972 : end if
1973 : end do
1974 0 : return
1975 : end
1976 : C
1977 :
1978 :
1979 : C Calculate gridded coordinates and the phasor needed for
1980 : C phase rotation.
1981 : C
1982 0 : subroutine smos (uvw, dphase, freq, c, scale, offset,
1983 : $ sampling, pos, loc, off, phasor)
1984 : implicit none
1985 : integer loc(2), off(2), sampling
1986 : double precision uvw(3), freq, c, scale(2), offset(2)
1987 : real pos(2)
1988 : double precision dphase, phase
1989 : complex phasor
1990 : integer idim
1991 : double precision pi
1992 : data pi/3.14159265358979323846/
1993 :
1994 0 : do idim=1,2
1995 : pos(idim)=scale(idim)*uvw(idim)*freq/c+
1996 0 : $ (offset(idim)+1)
1997 0 : loc(idim)=nint(pos(idim))
1998 0 : if(sampling.gt.1) then
1999 : C if((pos(idim)-loc(idim)) < 0.0)then
2000 : C loc(idim)=loc(idim)-1
2001 : C end if
2002 0 : off(idim)=nint((pos(idim)-real(loc(idim)))*real(-sampling))
2003 : C if(off(idim).eq.sampling) then
2004 : C off(idim)=0
2005 : C loc(idim)=loc(idim)+1
2006 : C end if
2007 : else
2008 0 : off(idim)=0
2009 : end if
2010 : end do
2011 0 : phase=-2.0D0*pi*dphase*freq/c
2012 0 : phasor=cmplx(cos(phase), sin(phase))
2013 0 : return
2014 : end
2015 :
2016 6473874 : logical function omos (nx, ny, loc, support)
2017 : implicit none
2018 : integer nx, ny, nw, loc(2), support
2019 : omos=(loc(1)-support.ge.1).and.(loc(1)+support.le.nx).and.
2020 6473874 : $ (loc(2)-support.ge.1).and.(loc(2)+support.le.ny)
2021 6473874 : return
2022 : end
2023 :
2024 :
2025 2179673859 : logical function onmosgrid(loc, nx, ny, nx0, ny0,
2026 : $ nxsub, nysub, support, msuppx, msuppy,
2027 : $ psuppx, psuppy, centerin)
2028 : implicit none
2029 : integer, intent(in) :: nx, ny, nx0, ny0, nxsub, nysub, loc(2),
2030 : $ support
2031 : logical, intent(out) :: centerin
2032 :
2033 : integer, intent(out) :: msuppx, msuppy, psuppx, psuppy
2034 : integer :: loc1sub, loc1plus, loc2sub, loc2plus
2035 2179673859 : msuppx=merge(-support, nx0-loc(1), loc(1)-support >= nx0)
2036 2179673859 : msuppy=merge(-support, ny0-loc(2), loc(2)-support >= ny0)
2037 : psuppx=merge(support, nx0+nxsub-loc(1)-1 , (loc(1)+support)
2038 2179673859 : $ < (nx0+nxsub))
2039 : psuppy=merge(support, ny0+nysub-loc(2)-1 , (loc(2)+support)
2040 2179673859 : $ < (ny0+nysub))
2041 : C write(*,*) 'ny0,nysub,loc(2), msuppy',ny0,nysub,loc(2), msuppy,
2042 : C $ support, ((loc(2)+support) < (ny0+nysub))
2043 2179673859 : loc1sub=loc(1)-support
2044 2179673859 : loc1plus=loc(1)+support
2045 2179673859 : loc2sub=loc(2)-support
2046 2179673859 : loc2plus=loc(2)+support
2047 :
2048 : centerin=(loc(1).ge.nx0) .and. (loc(1).lt. (nx0+nxsub)).and.
2049 2179673859 : $ (loc(2).ge.ny0).and.(loc(2) .lt. (ny0+nysub))
2050 : onmosgrid=(loc1plus .ge. nx0) .and. (loc1sub
2051 : $ .le. (nx0+nxsub)) .and.(loc2plus .ge. ny0) .and.
2052 2179673859 : $ (loc2sub .le. (ny0+nysub))
2053 :
2054 2179673859 : return
2055 : end
|
