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: aips2-request@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 42195 : subroutine gmoswgtd (nvispol, nvischan,
450 126585 : $ flag, rflag, weight, nrow,
451 : $ nx, ny, npol, nchan,
452 : $ support, convsize, sampling,
453 84390 : $ chanmap, polmap,
454 84390 : $ weightgrid, sumwt, convweight, convplanemap,
455 84390 : $ convchanmap, convpolmap,
456 : $ nconvplane, nconvchan, nconvpol, rbeg,
457 42195 : $ 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 :
489 : integer :: iloc(2)
490 : integer :: iiloc(2)
491 : integer, intent(in) :: rbeg, rend
492 : integer :: ix, iy, iz, ipol, ichan, xind, yind
493 : integer :: apol, achan, aconvplane, irow
494 : integer :: aconvpol, aconvchan, xind2, yind2
495 : integer :: posx, posy, msupportx, msupporty, psupportx, psupporty
496 : logical :: centin
497 :
498 :
499 :
500 :
501 11273493 : do irow=rbeg, rend
502 11231298 : aconvplane=convplanemap(irow)+1
503 11273493 : if(rflag(irow).eq.0) then
504 38277780 : do ichan=1, nvischan
505 27159114 : achan=chanmap(ichan)+1
506 27159114 : aconvchan=convchanmap(ichan)+1
507 27159114 : if((achan.ge.1).and.(achan.le.nchan).and.
508 11118666 : $ (weight(ichan,irow).ne.0.0)) then
509 24394320 : if (onmosgrid(loc(1, ichan, irow), nx, ny, 1, 1,
510 : $ nx, ny, support, msupportx, msupporty,
511 : $ psupportx, psupporty, centin)) then
512 73157148 : do ipol=1, nvispol
513 48771432 : apol=polmap(ipol)+1
514 48771432 : aconvpol=convpolmap(ipol)+1
515 : if((flag(ipol,ichan,irow).ne.1).and.
516 73157148 : $ (apol.ge.1).and.(apol.le.npol)) then
517 : C If we are making a PSF then we don't want to phase
518 : C rotate but we do want to reproject uvw
519 :
520 :
521 44916624 : nweight=cmplx(weight(ichan,irow))
522 : sumwt(apol, achan)=sumwt(apol, achan)+
523 44916624 : $ weight(ichan, irow)
524 :
525 : C norm will be the value we would get for the peak
526 : C at the phase center. We will want to normalize
527 : C the final image by this term.
528 44916624 : norm=0.0
529 : C write(*,*)off
530 1377664644 : do iy=msupporty, psupporty
531 44608931490 : do ix=msupportx, psupportx
532 :
533 43231266846 : xind2=sampling*ix+(convsize)/2+1
534 43231266846 : yind2=sampling*iy+(convsize)/2+1
535 : cwt=convweight(xind2,
536 43231266846 : $ yind2, aconvpol, aconvchan, aconvplane)
537 43231266846 : iiloc(1)=nx/2+1+ix
538 43231266846 : iiloc(2)=ny/2+1+iy
539 : weightgrid(iiloc(1),iiloc(2),
540 : $ apol,achan)= weightgrid(
541 : $ iiloc(1),iiloc(2),apol,achan)
542 44564014866 : $ + nweight*cwt
543 :
544 :
545 : end do
546 : end do
547 :
548 : end if
549 : end do
550 : end if
551 : end if
552 : end do
553 : end if
554 : end do
555 42195 : return
556 : end
557 :
558 : C Single precision weight grid image...Damn you fortran...no templates
559 0 : subroutine gmoswgts (nvispol, nvischan,
560 0 : $ flag, rflag, weight, nrow,
561 : $ nx, ny, npol, nchan,
562 : $ support, convsize, sampling,
563 0 : $ chanmap, polmap,
564 0 : $ weightgrid, sumwt, convweight, convplanemap,
565 0 : $ convchanmap, convpolmap,
566 : $ nconvplane, nconvchan, nconvpol, rbeg,
567 0 : $ rend, loc, off, phasor)
568 :
569 : implicit none
570 : integer, intent(in) :: nx,ny,npol,nchan, nvispol, nvischan, nrow
571 :
572 :
573 : integer, intent(in) :: loc(2, nvischan, nrow)
574 : integer, intent(in) :: off(2, nvischan, nrow)
575 : complex, intent(in) :: phasor(nvischan, nrow)
576 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
577 : integer, intent(in) :: rflag(nrow)
578 : real, intent(in) :: weight(nvischan, nrow)
579 : integer, intent(in) :: support
580 : integer, intent(in) :: chanmap(nchan), polmap(npol)
581 : complex, intent(inout) :: weightgrid(nx, ny, npol, nchan)
582 : double precision, intent(inout) :: sumwt(npol, nchan)
583 : complex :: nweight
584 : integer, intent(in) :: convsize, sampling
585 : integer, intent(in) :: nconvplane, nconvchan, nconvpol
586 : integer, intent(in) :: convplanemap(nrow)
587 : integer, intent(in) :: convchanmap(nvischan)
588 : integer, intent(in) :: convpolmap(nvispol)
589 : complex :: cwt
590 : complex, intent(in) :: convweight(convsize, convsize, nconvpol,
591 : $ nconvchan, nconvplane)
592 :
593 :
594 : real :: norm
595 : real :: wt
596 :
597 : logical :: onmosgrid
598 :
599 : integer :: iloc(2)
600 : integer :: iiloc(2)
601 : integer, intent(in) :: rbeg, rend
602 : integer :: ix, iy, iz, ipol, ichan, xind, yind
603 : integer :: apol, achan, aconvplane, irow
604 : integer :: aconvpol, aconvchan, xind2, yind2
605 : integer :: posx, posy, msupportx, msupporty, psupportx, psupporty
606 : logical :: centin
607 :
608 :
609 :
610 :
611 0 : do irow=rbeg, rend
612 0 : aconvplane=convplanemap(irow)+1
613 0 : if(rflag(irow).eq.0) then
614 0 : do ichan=1, nvischan
615 0 : achan=chanmap(ichan)+1
616 0 : aconvchan=convchanmap(ichan)+1
617 0 : if((achan.ge.1).and.(achan.le.nchan).and.
618 0 : $ (weight(ichan,irow).ne.0.0)) then
619 0 : if (onmosgrid(loc(1, ichan, irow), nx, ny, 1, 1,
620 : $ nx, ny, support, msupportx, msupporty,
621 : $ psupportx, psupporty, centin)) then
622 0 : do ipol=1, nvispol
623 0 : apol=polmap(ipol)+1
624 0 : aconvpol=convpolmap(ipol)+1
625 : if((flag(ipol,ichan,irow).ne.1).and.
626 0 : $ (apol.ge.1).and.(apol.le.npol)) then
627 : C If we are making a PSF then we don't want to phase
628 : C rotate but we do want to reproject uvw
629 :
630 :
631 0 : nweight=cmplx(weight(ichan,irow))
632 : sumwt(apol, achan)=sumwt(apol, achan)+
633 0 : $ weight(ichan, irow)
634 :
635 : C norm will be the value we would get for the peak
636 : C at the phase center. We will want to normalize
637 : C the final image by this term.
638 0 : norm=0.0
639 : C write(*,*)off
640 0 : do iy=msupporty, psupporty
641 0 : do ix=msupportx, psupportx
642 :
643 0 : xind2=sampling*ix+(convsize)/2+1
644 0 : yind2=sampling*iy+(convsize)/2+1
645 : cwt=convweight(xind2,
646 0 : $ yind2, aconvpol, aconvchan, aconvplane)
647 0 : iiloc(1)=nx/2+1+ix
648 0 : iiloc(2)=ny/2+1+iy
649 : weightgrid(iiloc(1),iiloc(2),
650 : $ apol,achan)= weightgrid(
651 : $ iiloc(1),iiloc(2),apol,achan)
652 0 : $ + nweight*cwt
653 :
654 :
655 : end do
656 : end do
657 :
658 : end if
659 : end do
660 : end if
661 : end if
662 : end do
663 : end if
664 : end do
665 0 : return
666 : end
667 :
668 :
669 125242 : subroutine sectgmosd2 (values, nvispol, nvischan,
670 187863 : $ dopsf, flag, rflag, weight, nrow,
671 62621 : $ grid, nx, ny, npol, nchan,
672 62621 : $ support, convsize, sampling, convfunc,
673 125242 : $ chanmap, polmap,
674 125242 : $ sumwt, convplanemap,
675 125242 : $ convchanmap, convpolmap,
676 : $ nconvplane, nconvchan, nconvpol, x0, y0, nxsub, nysub, rbeg,
677 125242 : $ rend, loc, off, phasor)
678 :
679 : implicit none
680 : integer, intent(in) :: nx,ny,npol,nchan, nvispol, nvischan, nrow
681 : complex, intent(in) :: values(nvispol, nvischan, nrow)
682 : double complex, intent(inout) :: grid(nx, ny, npol, nchan)
683 :
684 : integer, intent(in) :: x0, y0, nxsub, nysub
685 : integer, intent(in) :: loc(2, nvischan, nrow)
686 : integer, intent(in) :: off(2, nvischan, nrow)
687 : complex, intent(in) :: phasor(nvischan, nrow)
688 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
689 : integer, intent(in) :: rflag(nrow)
690 : real, intent(in) :: weight(nvischan, nrow)
691 : double precision, intent(inout) :: sumwt(npol, nchan)
692 : integer, intent(in) :: support
693 : integer, intent(in) :: chanmap(nchan), polmap(npol)
694 : integer, intent(in) :: dopsf
695 :
696 : double complex :: nvalue
697 : integer, intent(in) :: convsize, sampling
698 : integer, intent(in) :: nconvplane, nconvchan, nconvpol
699 : integer, intent(in) :: convplanemap(nrow)
700 : integer, intent(in) :: convchanmap(nvischan)
701 : integer, intent(in) :: convpolmap(nvispol)
702 : complex, intent(in) :: convfunc(convsize, convsize, nconvpol,
703 : $ nconvchan, nconvplane)
704 : complex :: cwt
705 :
706 :
707 : real :: norm
708 : real :: wt
709 :
710 : logical :: onmosgrid
711 :
712 : integer :: iloc(2)
713 : integer :: iiloc(2)
714 : integer, intent(in) :: rbeg, rend
715 : integer :: ix, iy, iz, ipol, ichan, xind, yind
716 : integer :: apol, achan, aconvplane, irow
717 : integer :: aconvpol, aconvchan, xind2, yind2
718 : integer :: posx, posy, msupportx, msupporty, psupportx, psupporty
719 : logical :: centin
720 :
721 :
722 :
723 :
724 13725185 : do irow=rbeg, rend
725 13662564 : aconvplane=convplanemap(irow)+1
726 13725185 : if(rflag(irow).eq.0) then
727 45332740 : do ichan=1, nvischan
728 31792244 : achan=chanmap(ichan)+1
729 31792244 : aconvchan=convchanmap(ichan)+1
730 31792244 : if((achan.ge.1).and.(achan.le.nchan).and.
731 13540496 : $ (weight(ichan,irow).ne.0.0)) then
732 29020181 : if (onmosgrid(loc(1, ichan, irow), nx, ny, x0, y0,
733 : $ nxsub, nysub, support, msupportx, msupporty,
734 : $ psupportx, psupporty, centin)) then
735 87020355 : do ipol=1, nvispol
736 58013570 : apol=polmap(ipol)+1
737 58013570 : aconvpol=convpolmap(ipol)+1
738 : if((flag(ipol,ichan,irow).ne.1).and.
739 87020355 : $ (apol.ge.1).and.(apol.le.npol)) then
740 : C If we are making a PSF then we don't want to phase
741 : C rotate but we do want to reproject uvw
742 54158762 : if(dopsf.eq.1) then
743 24368476 : nvalue=cmplx(weight(ichan,irow))
744 : else
745 : nvalue=weight(ichan,irow)*
746 29790286 : $ (values(ipol,ichan,irow)*phasor(ichan, irow))
747 : end if
748 :
749 :
750 : C norm will be the value we would get for the peak
751 : C at the phase center. We will want to normalize
752 : C the final image by this term.
753 54158762 : norm=0.0
754 1667673680 : do iy=msupporty, psupporty
755 : iloc(2)=(sampling*iy)+
756 1613514918 : $ off(2, ichan, irow)
757 1613514918 : yind=iloc(2)+(convsize)/2+1
758 53984658888 : do ix=msupportx, psupportx
759 : iloc(1)=(sampling*ix)+
760 52316985208 : $ off(1, ichan, irow)
761 52316985208 : xind=iloc(1)+(convsize)/2+1
762 : cwt=convfunc(xind, yind,
763 52316985208 : $ aconvpol, aconvchan, aconvplane)
764 : C write(*,*) support, iloc
765 : C write(*,*) loc(1, ichan, irow)+ix,loc(2, ichan, irow)+iy,xind,yind
766 : grid(loc(1, ichan, irow)+ix,
767 : $ loc(2, ichan, irow)+iy,apol,achan)=
768 : $ grid(loc(1, ichan, irow)+ix,
769 : $ loc(2, ichan, irow)+iy,apol,achan)+
770 53930500126 : $ nvalue*cwt
771 : end do
772 : end do
773 54158762 : if(centin) then
774 : sumwt(apol, achan)= sumwt(apol,achan)+
775 54154754 : $ weight(ichan,irow)
776 : endif
777 : end if
778 : end do
779 : C if onmos
780 : end if
781 : end if
782 : end do
783 : end if
784 : end do
785 62621 : return
786 : end
787 : C Single Precision version
788 0 : subroutine sectgmoss2 (values, nvispol, nvischan,
789 0 : $ dopsf, flag, rflag, weight, nrow,
790 0 : $ grid, nx, ny, npol, nchan,
791 0 : $ support, convsize, sampling, convfunc,
792 0 : $ chanmap, polmap,
793 0 : $ sumwt, convplanemap,
794 0 : $ convchanmap, convpolmap,
795 : $ nconvplane, nconvchan, nconvpol, x0, y0, nxsub, nysub, rbeg,
796 0 : $ rend, loc, off, phasor)
797 :
798 : implicit none
799 : integer, intent(in) :: nx,ny,npol,nchan, nvispol, nvischan, nrow
800 : complex, intent(in) :: values(nvispol, nvischan, nrow)
801 : complex, intent(inout) :: grid(nx, ny, npol, nchan)
802 :
803 : integer, intent(in) :: x0, y0, nxsub, nysub
804 : integer, intent(in) :: loc(2, nvischan, nrow)
805 : integer, intent(in) :: off(2, nvischan, nrow)
806 : complex, intent(in) :: phasor(nvischan, nrow)
807 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
808 : integer, intent(in) :: rflag(nrow)
809 : real, intent(in) :: weight(nvischan, nrow)
810 : double precision, intent(inout) :: sumwt(npol, nchan)
811 : integer, intent(in) :: support
812 : integer, intent(in) :: chanmap(nchan), polmap(npol)
813 : integer, intent(in) :: dopsf
814 :
815 : complex :: nvalue
816 : integer, intent(in) :: convsize, sampling
817 : integer, intent(in) :: nconvplane, nconvchan, nconvpol
818 : integer, intent(in) :: convplanemap(nrow)
819 : integer, intent(in) :: convchanmap(nvischan)
820 : integer, intent(in) :: convpolmap(nvispol)
821 : complex, intent(in) :: convfunc(convsize, convsize, nconvpol,
822 : $ nconvchan, nconvplane)
823 : complex :: cwt
824 :
825 :
826 : real :: norm
827 : real :: wt
828 :
829 : logical :: onmosgrid
830 :
831 : integer :: iloc(2)
832 : integer :: iiloc(2)
833 : integer, intent(in) :: rbeg, rend
834 : integer :: ix, iy, iz, ipol, ichan, xind, yind
835 : integer :: apol, achan, aconvplane, irow
836 : integer :: aconvpol, aconvchan, xind2, yind2
837 : integer :: posx, posy, msupportx, msupporty, psupportx, psupporty
838 : logical :: centin
839 :
840 :
841 :
842 :
843 0 : do irow=rbeg, rend
844 0 : aconvplane=convplanemap(irow)+1
845 0 : if(rflag(irow).eq.0) then
846 0 : do ichan=1, nvischan
847 0 : achan=chanmap(ichan)+1
848 0 : aconvchan=convchanmap(ichan)+1
849 0 : if((achan.ge.1).and.(achan.le.nchan).and.
850 0 : $ (weight(ichan,irow).ne.0.0)) then
851 0 : if (onmosgrid(loc(1, ichan, irow), nx, ny, x0, y0,
852 : $ nxsub, nysub, support, msupportx, msupporty,
853 : $ psupportx, psupporty, centin)) then
854 0 : do ipol=1, nvispol
855 0 : apol=polmap(ipol)+1
856 0 : aconvpol=convpolmap(ipol)+1
857 : if((flag(ipol,ichan,irow).ne.1).and.
858 0 : $ (apol.ge.1).and.(apol.le.npol)) then
859 : C If we are making a PSF then we don't want to phase
860 : C rotate but we do want to reproject uvw
861 0 : if(dopsf.eq.1) then
862 0 : nvalue=cmplx(weight(ichan,irow))
863 : else
864 : nvalue=weight(ichan,irow)*
865 0 : $ (values(ipol,ichan,irow)*phasor(ichan, irow))
866 : end if
867 :
868 :
869 : C norm will be the value we would get for the peak
870 : C at the phase center. We will want to normalize
871 : C the final image by this term.
872 0 : norm=0.0
873 0 : do iy=msupporty, psupporty
874 : iloc(2)=(sampling*iy)+
875 0 : $ off(2, ichan, irow)
876 0 : yind=iloc(2)+(convsize)/2+1
877 0 : do ix=msupportx, psupportx
878 : iloc(1)=(sampling*ix)+
879 0 : $ off(1, ichan, irow)
880 0 : xind=iloc(1)+(convsize)/2+1
881 : cwt=convfunc(xind, yind,
882 0 : $ aconvpol, aconvchan, aconvplane)
883 : C write(*,*) support, iloc
884 : grid(loc(1, ichan, irow)+ix,
885 : $ loc(2, ichan, irow)+iy,apol,achan)=
886 : $ grid(loc(1, ichan, irow)+ix,
887 : $ loc(2, ichan, irow)+iy,apol,achan)+
888 0 : $ nvalue*cwt
889 : end do
890 : end do
891 0 : if(centin) then
892 : sumwt(apol, achan)= sumwt(apol,achan)+
893 0 : $ weight(ichan,irow)
894 : end if
895 : end if
896 : end do
897 : end if
898 : end if
899 : end do
900 : end if
901 : end do
902 0 : return
903 : end
904 : C
905 :
906 :
907 0 : subroutine gmoss (uvw, dphase, values, nvispol, nvischan,
908 0 : $ dopsf, flag, rflag, weight, nrow, rownum,
909 0 : $ scale, offset, grid, nx, ny, npol, nchan, freq, c,
910 0 : $ support, convsize, sampling, convfunc,
911 0 : $ chanmap, polmap,
912 0 : $ sumwt, weightgrid, convweight, doweightgrid, convplanemap,
913 : $ nconvplane)
914 :
915 : implicit none
916 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
917 : complex values(nvispol, nvischan, nrow)
918 : complex grid(nx, ny, npol, nchan)
919 :
920 : double precision uvw(3, nrow), freq(nvischan), c, scale(2),
921 : $ offset(2)
922 : double precision dphase(nrow), uvdist
923 : double precision xlast, ylast
924 : complex phasor
925 : integer flag(nvispol, nvischan, nrow)
926 : integer rflag(nrow)
927 : real weight(nvischan, nrow), phase
928 : double precision sumwt(npol, nchan)
929 : integer rownum
930 : integer support
931 : integer chanmap(nchan), polmap(npol)
932 : integer dopsf
933 : complex weightgrid(nx, ny, npol, nchan)
934 : integer doweightgrid
935 :
936 : complex nvalue
937 : complex nweight
938 : integer convsize, sampling
939 : integer nconvplane
940 : integer convplanemap(nrow)
941 : complex convfunc(convsize, convsize, nconvplane), cwt, crot
942 : complex convweight(convsize, convsize, nconvplane)
943 :
944 :
945 : complex sconv(-(support+1)*sampling:(support+1)*sampling,
946 0 : $ -(support+1)*sampling:(support+1)*sampling, nconvplane)
947 : complex sconv2(-(support+1)*sampling:(support+1)*sampling,
948 0 : $ -(support+1)*sampling:(support+1)*sampling, nconvplane)
949 : real sumsconv
950 : real sumsconv2
951 : real ratioofbeam
952 :
953 : real norm
954 : real wt
955 :
956 : logical omos, doshift
957 :
958 : real pos(3)
959 : integer loc(2), off(2), iloc(2)
960 : integer iiloc(2)
961 : integer rbeg, rend
962 : integer ix, iy, iz, ipol, ichan
963 : integer apol, achan, aconvplane, irow
964 : double precision pi
965 : data pi/3.14159265358979323846/
966 : integer sampsupp
967 0 : sampsupp=(support+1)*sampling
968 0 : irow=rownum
969 :
970 0 : sumsconv=0
971 0 : sumsconv2=0
972 0 : ratioofbeam=1.0
973 0 : do iz=1, nconvplane
974 0 : do iy=-sampsupp,sampsupp
975 0 : iloc(2)=iy+convsize/2+1
976 0 : do ix=-sampsupp,sampsupp
977 0 : iloc(1)=ix+convsize/2+1
978 0 : sconv(ix,iy,iz)=(convfunc(iloc(1), iloc(2),iz))
979 0 : sconv2(ix,iy,iz)=convweight(iloc(1), iloc(2),iz)
980 : end do
981 : end do
982 : end do
983 0 : doshift=.FALSE.
984 :
985 0 : if(irow.ge.0) then
986 0 : rbeg=irow+1
987 0 : rend=irow+1
988 : else
989 0 : rbeg=1
990 0 : rend=nrow
991 : end if
992 :
993 0 : xlast=0.0
994 0 : ylast=0.0
995 0 : do irow=rbeg, rend
996 0 : aconvplane=convplanemap(irow)+1
997 0 : if(rflag(irow).eq.0) then
998 0 : do ichan=1, nvischan
999 0 : achan=chanmap(ichan)+1
1000 0 : if((achan.ge.1).and.(achan.le.nchan).and.
1001 0 : $ (weight(ichan,irow).ne.0.0)) then
1002 : call smos(uvw(1,irow), dphase(irow), freq(ichan), c,
1003 0 : $ scale, offset, sampling, pos, loc, off, phasor)
1004 0 : if (omos(nx, ny, loc, support)) then
1005 0 : do ipol=1, nvispol
1006 0 : apol=polmap(ipol)+1
1007 : if((flag(ipol,ichan,irow).ne.1).and.
1008 0 : $ (apol.ge.1).and.(apol.le.npol)) then
1009 : C If we are making a PSF then we don't want to phase
1010 : C rotate but we do want to reproject uvw
1011 0 : if(dopsf.eq.1) then
1012 0 : nvalue=cmplx(weight(ichan,irow))
1013 : else
1014 : nvalue=weight(ichan,irow)*
1015 0 : $ (values(ipol,ichan,irow)*phasor)
1016 : end if
1017 0 : if(doweightgrid .gt. 0) then
1018 0 : nweight=cmplx(weight(ichan,irow))
1019 : end if
1020 :
1021 : C norm will be the value we would get for the peak
1022 : C at the phase center. We will want to normalize
1023 : C the final image by this term.
1024 0 : norm=0.0
1025 0 : if(sampling.eq.1) then
1026 0 : do iy=-support,support
1027 0 : do ix=-support,support
1028 : grid(loc(1)+ix,
1029 : $ loc(2)+iy,apol,achan)=
1030 : $ grid(loc(1)+ix,
1031 : $ loc(2)+iy,apol,achan)+
1032 0 : $ nvalue*sconv(ix,iy, aconvplane)
1033 :
1034 0 : if(doweightgrid .gt. 0) then
1035 0 : iloc(1)=nx/2+1+ix
1036 0 : iloc(2)=ny/2+1+iy
1037 : weightgrid(iloc(1),iloc(2),
1038 : $ apol,achan)= weightgrid(
1039 : $ iloc(1),iloc(2),apol,achan)
1040 0 : $ + nweight*sconv2(ix,iy,aconvplane)
1041 :
1042 : end if
1043 : end do
1044 : end do
1045 : else
1046 0 : do iy=-support,support
1047 0 : iloc(2)=(sampling*iy+off(2))+1
1048 0 : do ix=-support, support
1049 0 : iloc(1)=(sampling*ix+off(1))+1
1050 : cwt=sconv(iloc(1),
1051 0 : $ iloc(2),aconvplane)
1052 : grid(loc(1)+ix,
1053 : $ loc(2)+iy,apol,achan)=
1054 : $ grid(loc(1)+ix,
1055 : $ loc(2)+iy,apol,achan)+
1056 0 : $ nvalue*cwt
1057 0 : if(doweightgrid .gt. 0) then
1058 : cwt=sconv2(iloc(1), iloc(2),
1059 0 : $ aconvplane)
1060 0 : iiloc(1)=nx/2+1+ix
1061 0 : iiloc(2)=ny/2+1+iy
1062 : weightgrid(iiloc(1),iiloc(2),
1063 : $ apol,achan)= weightgrid(
1064 : $ iiloc(1),iiloc(2),apol,achan)
1065 0 : $ + nweight*cwt
1066 :
1067 : end if
1068 : end do
1069 : end do
1070 : end if
1071 : sumwt(apol, achan)= sumwt(apol,achan)+
1072 0 : $ weight(ichan,irow)
1073 : end if
1074 : end do
1075 : end if
1076 : end if
1077 : end do
1078 : end if
1079 : end do
1080 0 : return
1081 : end
1082 : C
1083 0 : subroutine gmoss2 (uvw, dphase, values, nvispol, nvischan,
1084 0 : $ dopsf, flag, rflag, weight, nrow, rownum,
1085 0 : $ scale, offset, grid, nx, ny, npol, nchan, freq, c,
1086 0 : $ support, convsize, sampling, convfunc,
1087 0 : $ chanmap, polmap,
1088 0 : $ sumwt, weightgrid, convweight, doweightgrid, convplanemap,
1089 0 : $ convchanmap, convpolmap,
1090 : $ nconvplane, nconvchan, nconvpol)
1091 :
1092 : implicit none
1093 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
1094 : complex values(nvispol, nvischan, nrow)
1095 : complex grid(nx, ny, npol, nchan)
1096 :
1097 : double precision uvw(3, nrow), freq(nvischan), c, scale(2),
1098 : $ offset(2)
1099 : double precision dphase(nrow), uvdist
1100 : double precision xlast, ylast
1101 : complex phasor
1102 : integer flag(nvispol, nvischan, nrow)
1103 : integer rflag(nrow)
1104 : real weight(nvischan, nrow), phase
1105 : double precision sumwt(npol, nchan)
1106 : integer rownum
1107 : integer support
1108 : integer chanmap(nchan), polmap(npol)
1109 : integer dopsf
1110 : complex weightgrid(nx, ny, npol, nchan)
1111 : integer doweightgrid
1112 :
1113 : double complex nvalue
1114 : double complex nweight
1115 : integer convsize, sampling
1116 : integer nconvplane, nconvchan, nconvpol
1117 : integer convplanemap(nrow)
1118 : integer convchanmap(nvischan)
1119 : integer convpolmap(nvispol)
1120 : complex convfunc(convsize, convsize, nconvpol, nconvchan,
1121 : $ nconvplane), cwt, crot
1122 : complex convweight(convsize, convsize, nconvpol, nconvchan,
1123 : $ nconvplane)
1124 : real sumsconv
1125 : real sumsconv2
1126 : real ratioofbeam
1127 :
1128 : real norm
1129 : real wt
1130 :
1131 : logical omos, doshift
1132 :
1133 : real pos(3)
1134 : integer loc(2), off(2), iloc(2)
1135 : integer iiloc(2)
1136 : integer rbeg, rend
1137 : integer ix, iy, iz, ipol, ichan, xind, yind
1138 : integer apol, achan, aconvplane, irow
1139 : integer aconvpol, aconvchan, xind2, yind2
1140 : double precision pi
1141 : data pi/3.14159265358979323846/
1142 : real maxsconv2, minsconv2
1143 0 : irow=rownum
1144 :
1145 0 : sumsconv=0
1146 0 : sumsconv2=0
1147 0 : ratioofbeam=1.0
1148 :
1149 0 : doshift=.FALSE.
1150 :
1151 0 : if(irow.ge.0) then
1152 0 : rbeg=irow+1
1153 0 : rend=irow+1
1154 : else
1155 0 : rbeg=1
1156 0 : rend=nrow
1157 : end if
1158 :
1159 : C write(*,*) 'max of sconvs ', maxsconv2, minsconv2, sampsupp,
1160 : C $ convsize
1161 :
1162 :
1163 :
1164 0 : xlast=0.0
1165 0 : ylast=0.0
1166 0 : do irow=rbeg, rend
1167 0 : aconvplane=convplanemap(irow)+1
1168 0 : if(rflag(irow).eq.0) then
1169 0 : do ichan=1, nvischan
1170 0 : achan=chanmap(ichan)+1
1171 0 : aconvchan=convchanmap(ichan)+1
1172 0 : if((achan.ge.1).and.(achan.le.nchan).and.
1173 0 : $ (weight(ichan,irow).ne.0.0)) then
1174 : call smos(uvw(1,irow), dphase(irow), freq(ichan), c,
1175 0 : $ scale, offset, sampling, pos, loc, off, phasor)
1176 0 : if (omos(nx, ny, loc, support)) then
1177 0 : do ipol=1, nvispol
1178 0 : apol=polmap(ipol)+1
1179 0 : aconvpol=convpolmap(ipol)+1
1180 : if((flag(ipol,ichan,irow).ne.1).and.
1181 0 : $ (apol.ge.1).and.(apol.le.npol)) then
1182 : C If we are making a PSF then we don't want to phase
1183 : C rotate but we do want to reproject uvw
1184 0 : if(dopsf.eq.1) then
1185 0 : nvalue=cmplx(weight(ichan,irow))
1186 : else
1187 : nvalue=weight(ichan,irow)*
1188 0 : $ (values(ipol,ichan,irow)*phasor)
1189 : end if
1190 0 : if(doweightgrid .gt. 0) then
1191 0 : nweight=cmplx(weight(ichan,irow))
1192 : end if
1193 :
1194 : C norm will be the value we would get for the peak
1195 : C at the phase center. We will want to normalize
1196 : C the final image by this term.
1197 0 : norm=0.0
1198 0 : if(sampling.eq.1) then
1199 0 : do iy=-support,support
1200 0 : yind=iy+(convsize)/2+1
1201 0 : do ix=-support,support
1202 0 : xind=ix+(convsize)/2+1
1203 : grid(loc(1)+ix,
1204 : $ loc(2)+iy,apol,achan)=
1205 : $ grid(loc(1)+ix,
1206 : $ loc(2)+iy,apol,achan)+
1207 : $ nvalue*convfunc(xind, yind,
1208 0 : $ aconvpol, aconvchan, aconvplane)
1209 :
1210 0 : if(doweightgrid .gt. 0) then
1211 0 : iloc(1)=nx/2+1+ix
1212 0 : iloc(2)=ny/2+1+iy
1213 : weightgrid(iloc(1),iloc(2),
1214 : $ apol,achan)= weightgrid(
1215 : $ iloc(1),iloc(2),apol,achan)
1216 : $ + nweight*convweight(xind, yind,
1217 0 : $ aconvpol, aconvchan, aconvplane)
1218 :
1219 : end if
1220 : end do
1221 : end do
1222 : else
1223 : C write(*,*)off
1224 0 : do iy=-support, support
1225 0 : iloc(2)=(sampling*iy)+off(2)
1226 0 : yind=iloc(2)+(convsize)/2+1
1227 0 : do ix=-support, support
1228 0 : iloc(1)=(sampling*ix)+off(1)
1229 0 : xind=iloc(1)+(convsize)/2+1
1230 : cwt=convfunc(xind, yind,
1231 0 : $ aconvpol, aconvchan, aconvplane)
1232 : C write(*,*) yind, xind, loc(1)+ix, loc(2)+iy, apol, achan,
1233 : C $ aconvpol, aconvchan, aconvplane
1234 : grid(loc(1)+ix,
1235 : $ loc(2)+iy,apol,achan)=
1236 : $ grid(loc(1)+ix,
1237 : $ loc(2)+iy,apol,achan)+
1238 0 : $ nvalue*cwt
1239 0 : if(doweightgrid .gt. 0) then
1240 0 : xind2=sampling*ix+(convsize)/2+1
1241 0 : yind2=sampling*iy+(convsize)/2+1
1242 : cwt=convweight(xind2,
1243 : $ yind2, aconvpol, aconvchan,
1244 0 : $ aconvplane)
1245 0 : iiloc(1)=nx/2+1+ix
1246 0 : iiloc(2)=ny/2+1+iy
1247 : weightgrid(iiloc(1),iiloc(2),
1248 : $ apol,achan)= weightgrid(
1249 : $ iiloc(1),iiloc(2),apol,achan)
1250 0 : $ + nweight*cwt
1251 :
1252 : end if
1253 : end do
1254 : end do
1255 : end if
1256 : sumwt(apol, achan)= sumwt(apol,achan)+
1257 0 : $ weight(ichan,irow)
1258 : end if
1259 : end do
1260 : end if
1261 : end if
1262 : end do
1263 : end if
1264 : end do
1265 0 : return
1266 : end
1267 : C
1268 :
1269 :
1270 :
1271 : C Degrid a number of visibility records
1272 : C
1273 0 : subroutine dmos (uvw, dphase, values, nvispol, nvischan,
1274 0 : $ flag, rflag,
1275 0 : $ nrow, rownum, scale, offset, grid, nx, ny, npol, nchan, freq,
1276 0 : $ c, support, convsize, sampling, convfunc,
1277 0 : $ chanmap, polmap, convplanemap, nconvplane)
1278 :
1279 : implicit none
1280 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
1281 : integer nconvplane
1282 : complex values(nvispol, nvischan, nrow)
1283 : complex grid(nx, ny, npol, nchan)
1284 : double precision uvw(3, nrow), freq(nvischan), c, scale(2),
1285 : $ offset(2)
1286 : double precision dphase(nrow), uvdist
1287 : double precision xlast, ylast
1288 : complex phasor
1289 : integer flag(nvispol, nvischan, nrow)
1290 : integer rflag(nrow)
1291 : integer rownum
1292 : integer support
1293 : integer chanmap(nchan), polmap(npol)
1294 : integer convplanemap(nrow)
1295 : complex nvalue
1296 :
1297 : integer convsize, sampling
1298 : complex convfunc(convsize, convsize, nconvplane), cwt, crot
1299 :
1300 : integer sampsupp
1301 :
1302 : complex sconv(-(support+1)*sampling:(support+1)*sampling,
1303 0 : $ -(support+1)*sampling:(support+1)*sampling, nconvplane)
1304 :
1305 : real norm, phase
1306 :
1307 : logical omos, doshift
1308 :
1309 : real pos(2)
1310 : integer loc(2), off(2), iloc(2)
1311 : integer rbeg, rend
1312 : integer ix, iy, iz, ipol, ichan
1313 : integer apol, achan, aconvplane, irow
1314 : real wt, wtx, wty
1315 : double precision pi
1316 : data pi/3.14159265358979323846/
1317 :
1318 0 : sampsupp=(support+1)*sampling
1319 0 : irow=rownum
1320 :
1321 0 : do iz=1, nconvplane
1322 0 : do iy=-sampsupp,sampsupp
1323 0 : iloc(2)=iy+convsize/2+1
1324 0 : do ix=-sampsupp,sampsupp
1325 0 : iloc(1)=ix+convsize/2+1
1326 0 : sconv(ix,iy,iz)=conjg(convfunc(iloc(1), iloc(2),iz))
1327 : end do
1328 : end do
1329 : end do
1330 0 : doshift=.FALSE.
1331 :
1332 0 : if(irow.ge.0) then
1333 0 : rbeg=irow+1
1334 0 : rend=irow+1
1335 : else
1336 0 : rbeg=1
1337 0 : rend=nrow
1338 : end if
1339 : C
1340 0 : xlast=0.0
1341 0 : ylast=0.0
1342 0 : do irow=rbeg, rend
1343 0 : aconvplane=convplanemap(irow)+1
1344 0 : if(rflag(irow).eq.0) then
1345 0 : do ichan=1, nvischan
1346 0 : achan=chanmap(ichan)+1
1347 0 : if((achan.ge.1).and.(achan.le.nchan)) then
1348 : call smos(uvw(1,irow), dphase(irow), freq(ichan), c,
1349 0 : $ scale, offset, sampling, pos, loc, off, phasor)
1350 0 : if (omos(nx, ny, loc, support)) then
1351 0 : do ipol=1, nvispol
1352 0 : apol=polmap(ipol)+1
1353 : if((flag(ipol,ichan,irow).ne.1).and.
1354 0 : $ (apol.ge.1).and.(apol.le.npol)) then
1355 0 : nvalue=0.0
1356 0 : norm=0.0
1357 0 : if(sampling.eq.1) then
1358 0 : do iy=-support,support
1359 0 : do ix=-support,support
1360 : nvalue=nvalue+(sconv(ix,iy,
1361 : $ aconvplane))*
1362 : $ grid(loc(1)+ix,loc(2)+iy,
1363 0 : $ apol,achan)
1364 : end do
1365 : end do
1366 : else
1367 0 : do iy=-support,support
1368 0 : iloc(2)=sampling*iy+off(2)
1369 0 : do ix=-support,support
1370 : iloc(1)=ix*sampling
1371 0 : $ +off(1)
1372 : cwt=sconv(iloc(1),
1373 0 : $ iloc(2),aconvplane)
1374 : nvalue=nvalue+cwt*
1375 : $ grid(loc(1)+ix,loc(2)+iy,
1376 0 : $ apol,achan)
1377 : end do
1378 : end do
1379 : end if
1380 : values(ipol,ichan,irow)=nvalue*conjg(
1381 0 : $ phasor)
1382 : end if
1383 : end do
1384 : end if
1385 : end if
1386 : end do
1387 : end if
1388 : end do
1389 0 : return
1390 : end
1391 : C
1392 :
1393 : C Degrid a number of visibility records
1394 : C
1395 0 : subroutine dmos2 (uvw, dphase, values, nvispol, nvischan,
1396 0 : $ flag, rflag,
1397 0 : $ nrow, rownum, scale, offset, grid, nx, ny, npol, nchan, freq,
1398 0 : $ c, support, convsize, sampling, convfunc,
1399 0 : $ chanmap, polmap, convplanemap, convchanmap, convpolmap,
1400 : $ nconvplane, nconvchan, nconvpol)
1401 :
1402 : implicit none
1403 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow
1404 : integer nconvplane, nconvchan, nconvpol
1405 : complex values(nvispol, nvischan, nrow)
1406 : complex grid(nx, ny, npol, nchan)
1407 : double precision uvw(3, nrow), freq(nvischan), c, scale(2),
1408 : $ offset(2)
1409 : double precision dphase(nrow), uvdist
1410 : double precision xlast, ylast
1411 : complex phasor
1412 : integer flag(nvispol, nvischan, nrow)
1413 : integer rflag(nrow)
1414 : integer rownum
1415 : integer support
1416 : integer chanmap(nchan), polmap(npol)
1417 : integer convplanemap(nrow), convchanmap(nvischan)
1418 : integer convpolmap(nvispol)
1419 : complex nvalue
1420 :
1421 : integer convsize, sampling
1422 : complex convfunc(convsize, convsize, nconvpol, nconvchan,
1423 : $ nconvplane), cwt, crot
1424 :
1425 : integer sampsupp
1426 :
1427 : C complex sconv(-(support+1)*sampling:(support+1)*sampling,
1428 : C $ -(support+1)*sampling:(support+1)*sampling, nconvplane)
1429 :
1430 : real norm, phase
1431 :
1432 : logical omos, doshift
1433 :
1434 : real pos(2)
1435 : integer loc(2), off(2), iloc(2)
1436 : integer rbeg, rend
1437 : integer ix, iy, iz, ipol, ichan, xind, yind
1438 : integer apol, achan, aconvplane, irow
1439 : integer aconvchan, aconvpol
1440 : real wt, wtx, wty
1441 : double precision pi
1442 : data pi/3.14159265358979323846/
1443 :
1444 0 : sampsupp=(support+1)*sampling
1445 0 : irow=rownum
1446 :
1447 : C do iz=1, nconvplane
1448 : C do iy=-sampsupp,sampsupp
1449 : C iloc(2)=iy+convsize/2+1
1450 : c do ix=-sampsupp,sampsupp
1451 : C iloc(1)=ix+convsize/2+1
1452 : C sconv(ix,iy,iz)=conjg(convfunc(iloc(1), iloc(2),iz))
1453 : C end do
1454 : C end do
1455 : C end do
1456 0 : doshift=.FALSE.
1457 :
1458 0 : if(irow.ge.0) then
1459 0 : rbeg=irow+1
1460 0 : rend=irow+1
1461 : else
1462 0 : rbeg=1
1463 0 : rend=nrow
1464 : end if
1465 : C
1466 0 : xlast=0.0
1467 0 : ylast=0.0
1468 0 : do irow=rbeg, rend
1469 0 : aconvplane=convplanemap(irow)+1
1470 0 : if(rflag(irow).eq.0) then
1471 0 : do ichan=1, nvischan
1472 0 : achan=chanmap(ichan)+1
1473 0 : aconvchan=convchanmap(ichan)+1
1474 0 : if((achan.ge.1).and.(achan.le.nchan)) then
1475 : call smos(uvw(1,irow), dphase(irow), freq(ichan), c,
1476 0 : $ scale, offset, sampling, pos, loc, off, phasor)
1477 0 : if (omos(nx, ny, loc, support)) then
1478 0 : do ipol=1, nvispol
1479 0 : apol=polmap(ipol)+1
1480 0 : aconvpol=convpolmap(ipol)+1
1481 : if((flag(ipol,ichan,irow).ne.1).and.
1482 0 : $ (apol.ge.1).and.(apol.le.npol)) then
1483 : C write(*,*) 'aindices', aconvplane, aconvchan, aconvpol
1484 0 : nvalue=0.0
1485 0 : norm=0.0
1486 0 : if(sampling.eq.1) then
1487 0 : do iy=-support,support
1488 0 : yind=iy+(convsize)/2+1
1489 0 : do ix=-support,support
1490 0 : xind=ix+(convsize)/2+1
1491 : nvalue=nvalue+(convfunc(xind,yind,
1492 : $ aconvpol, aconvchan,
1493 : $ aconvplane))*
1494 : $ grid(loc(1)+ix,loc(2)+iy,
1495 0 : $ apol,achan)
1496 : end do
1497 : end do
1498 : else
1499 0 : do iy=-support,support
1500 0 : iloc(2)=sampling*iy+off(2)
1501 0 : yind=iloc(2)+(convsize)/2+1
1502 : C write(*,*) 'iloc(2)', iloc(2), off(2), yind
1503 0 : do ix=-support,support
1504 : iloc(1)=ix*sampling
1505 0 : $ +off(1)
1506 0 : xind=iloc(1)+(convsize)/2+1
1507 : C write(*,*) 'iloc(1)', iloc(1), off(1), xind
1508 : cwt=convfunc(xind, yind, aconvpol,
1509 0 : $ aconvchan,aconvplane)
1510 : nvalue=nvalue+cwt*
1511 : $ grid(loc(1)+ix,loc(2)+iy,
1512 0 : $ apol,achan)
1513 : end do
1514 : end do
1515 : end if
1516 : values(ipol,ichan,irow)=nvalue*conjg(
1517 0 : $ phasor)
1518 : end if
1519 : end do
1520 : end if
1521 : end if
1522 : end do
1523 : end if
1524 : end do
1525 0 : return
1526 : end
1527 : C
1528 :
1529 : C
1530 39162 : subroutine sectdmos2 (values, nvispol, nvischan,
1531 39162 : $ flag, rflag,
1532 19581 : $ nrow, grid, nx, ny, npol, nchan,
1533 19581 : $ support, convsize, sampling, convfunc,
1534 58743 : $ chanmap, polmap, convplanemap, convchanmap, convpolmap,
1535 58743 : $ nconvplane, nconvchan, nconvpol, rbeg,rend,loc,off,phasor)
1536 :
1537 : implicit none
1538 : integer, intent(in) :: nx, ny,npol,nchan,nvispol, nvischan, nrow
1539 : integer, intent(in) :: nconvplane, nconvchan, nconvpol
1540 : complex, intent(inout) :: values(nvispol, nvischan, nrow)
1541 : complex, intent(in) :: grid(nx, ny, npol, nchan)
1542 : complex, intent(in) :: phasor(nvischan, nrow)
1543 : integer, intent(in) :: flag(nvispol, nvischan, nrow)
1544 : integer, intent(in) :: rflag(nrow)
1545 : integer, intent(in) :: support
1546 : integer, intent(in) :: chanmap(nchan), polmap(npol)
1547 : integer, intent(in) :: convplanemap(nrow), convchanmap(nvischan)
1548 : integer, intent(in) :: convpolmap(nvispol)
1549 : complex :: nvalue
1550 :
1551 : integer, intent(in) :: convsize, sampling
1552 : complex, intent(in) :: convfunc(convsize, convsize, nconvpol,
1553 : $ nconvchan, nconvplane)
1554 : complex :: cwt, crot
1555 :
1556 : C complex sconv(-(support+1)*sampling:(support+1)*sampling,
1557 : C $ -(support+1)*sampling:(support+1)*sampling, nconvplane)
1558 :
1559 : real :: norm, phase
1560 :
1561 : logical :: omos
1562 :
1563 :
1564 : integer, intent(in) :: loc(2, nvischan, nrow),
1565 : $ off(2,nvischan,nrow)
1566 : integer :: iloc(2)
1567 : integer, intent(in) :: rbeg, rend
1568 : integer :: ix, iy, iz, ipol, ichan, xind, yind
1569 : integer :: apol, achan, aconvplane, irow
1570 : integer :: aconvchan, aconvpol
1571 : real :: wt, wtx, wty
1572 :
1573 :
1574 2964291 : do irow=rbeg, rend
1575 2944710 : aconvplane=convplanemap(irow)+1
1576 2964291 : if(rflag(irow).eq.0) then
1577 9353790 : do ichan=1, nvischan
1578 6456744 : achan=chanmap(ichan)+1
1579 6456744 : aconvchan=convchanmap(ichan)+1
1580 9353790 : if((achan.ge.1).and.(achan.le.nchan)) then
1581 6456744 : if (omos(nx, ny, loc(1, ichan,irow),support)) then
1582 19359837 : do ipol=1, nvispol
1583 12906558 : apol=polmap(ipol)+1
1584 12906558 : aconvpol=convpolmap(ipol)+1
1585 : if((flag(ipol,ichan,irow).ne.1).and.
1586 19359837 : $ (apol.ge.1).and.(apol.le.npol)) then
1587 : C write(*,*) 'aindices', aconvplane, aconvchan, aconvpol
1588 12815514 : nvalue=0.0
1589 12815514 : norm=0.0
1590 :
1591 273831044 : do iy=-support,support
1592 261015530 : iloc(2)=sampling*iy+off(2, ichan, irow)
1593 261015530 : yind=iloc(2)+(convsize)/2+1
1594 : C write(*,*) 'iloc(2)', iloc(2), off(2), yind
1595 8121239326 : do ix=-support,support
1596 : iloc(1)=ix*sampling
1597 7847408282 : $ +off(1, ichan, irow)
1598 7847408282 : xind=iloc(1)+(convsize)/2+1
1599 : C write(*,*) 'iloc(1)', iloc(1), off(1), xind
1600 : cwt=convfunc(xind, yind, aconvpol,
1601 7847408282 : $ aconvchan,aconvplane)
1602 : nvalue=nvalue+cwt*
1603 : $ grid(loc(1, ichan, irow)+ix,
1604 : $ loc(2, ichan, irow)+iy,
1605 8108423812 : $ apol,achan)
1606 : end do
1607 : end do
1608 :
1609 : values(ipol,ichan,irow)=nvalue*conjg(
1610 12815514 : $ phasor(ichan, irow))
1611 : end if
1612 : end do
1613 : end if
1614 : end if
1615 : end do
1616 : end if
1617 : end do
1618 19581 : return
1619 : end
1620 : C
1621 :
1622 :
1623 : C Calculate gridded coordinates and the phasor needed for
1624 : C phase rotation.
1625 : C
1626 0 : subroutine smos (uvw, dphase, freq, c, scale, offset,
1627 : $ sampling, pos, loc, off, phasor)
1628 : implicit none
1629 : integer loc(2), off(2), sampling
1630 : double precision uvw(3), freq, c, scale(2), offset(2)
1631 : real pos(2)
1632 : double precision dphase, phase
1633 : complex phasor
1634 : integer idim
1635 : double precision pi
1636 : data pi/3.14159265358979323846/
1637 :
1638 0 : do idim=1,2
1639 : pos(idim)=scale(idim)*uvw(idim)*freq/c+
1640 0 : $ (offset(idim)+1)
1641 0 : loc(idim)=nint(pos(idim))
1642 0 : if(sampling.gt.1) then
1643 : C if((pos(idim)-loc(idim)) < 0.0)then
1644 : C loc(idim)=loc(idim)-1
1645 : C end if
1646 0 : off(idim)=nint((pos(idim)-real(loc(idim)))*real(-sampling))
1647 : C if(off(idim).eq.sampling) then
1648 : C off(idim)=0
1649 : C loc(idim)=loc(idim)+1
1650 : C end if
1651 : else
1652 0 : off(idim)=0
1653 : end if
1654 : end do
1655 0 : phase=-2.0D0*pi*dphase*freq/c
1656 0 : phasor=cmplx(cos(phase), sin(phase))
1657 0 : return
1658 : end
1659 :
1660 6456744 : logical function omos (nx, ny, loc, support)
1661 : implicit none
1662 : integer nx, ny, nw, loc(2), support
1663 : omos=(loc(1)-support.ge.1).and.(loc(1)+support.le.nx).and.
1664 6456744 : $ (loc(2)-support.ge.1).and.(loc(2)+support.le.ny)
1665 6456744 : return
1666 : end
1667 :
1668 :
1669 53414501 : logical function onmosgrid(loc, nx, ny, nx0, ny0,
1670 : $ nxsub, nysub, support, msuppx, msuppy,
1671 : $ psuppx, psuppy, centerin)
1672 : implicit none
1673 : integer, intent(in) :: nx, ny, nx0, ny0, nxsub, nysub, loc(2),
1674 : $ support
1675 : logical, intent(out) :: centerin
1676 :
1677 : integer, intent(out) :: msuppx, msuppy, psuppx, psuppy
1678 : integer :: loc1sub, loc1plus, loc2sub, loc2plus
1679 53414501 : msuppx=merge(-support, nx0-loc(1), loc(1)-support >= nx0)
1680 53414501 : msuppy=merge(-support, ny0-loc(2), loc(2)-support >= ny0)
1681 : psuppx=merge(support, nx0+nxsub-loc(1)-1 , (loc(1)+support)
1682 53414501 : $ < (nx0+nxsub))
1683 : psuppy=merge(support, ny0+nysub-loc(2)-1 , (loc(2)+support)
1684 53414501 : $ < (ny0+nysub))
1685 : C write(*,*) 'ny0,nysub,loc(2), msuppy',ny0,nysub,loc(2), msuppy,
1686 : C $ support, ((loc(2)+support) < (ny0+nysub))
1687 53414501 : loc1sub=loc(1)-support
1688 53414501 : loc1plus=loc(1)+support
1689 53414501 : loc2sub=loc(2)-support
1690 53414501 : loc2plus=loc(2)+support
1691 :
1692 : centerin=(loc(1).ge.nx0) .and. (loc(1).lt. (nx0+nxsub)).and.
1693 53414501 : $ (loc(2).ge.ny0).and.(loc(2) .lt. (ny0+nysub))
1694 : onmosgrid=(loc1plus .ge. nx0) .and. (loc1sub
1695 : $ .le. (nx0+nxsub)) .and.(loc2plus .ge. ny0) .and.
1696 53414501 : $ (loc2sub .le. (ny0+nysub))
1697 :
1698 53414501 : return
1699 : end
|
