Line data Source code
1 : *=======================================================================
2 : * -*- FORTRAN -*-
3 : * Copyright (C) 1999,2001,2002
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: fpbwproj.f,v 1.13 2006/07/20 00:24:20 sbhatnag Exp $
29 : *-----------------------------------------------------------------------
30 : C
31 : C Grid a number of visibility records
32 : C
33 0 : subroutine gpbwproj (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, wconvsize, convfunc,
37 0 : $ chanmap, polmap,polused,sumwt,
38 0 : $ ant1, ant2, nant, scanno, sigma,raoff, decoff,area,
39 0 : $ dograd,dopointingcorrection,npa,paindex,cfmap,conjcfmap,
40 : $ currentCFPA,actualPA,cfRefFreq)
41 :
42 :
43 : implicit none
44 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow,polused
45 : integer npa,paindex
46 : integer convsize, sampling, wconvsize,dograd,dopointingcorrection
47 : complex values(nvispol, nvischan, nrow)
48 : complex grid(nx, ny, npol, nchan)
49 : double precision uvw(3, nrow), freq(nvischan), c, scale(3),
50 : $ offset(3),currentCFPA,actualPA,cfRefFreq
51 : double precision dphase(nrow), uvdist
52 : complex phasor
53 : integer flag(nvispol, nvischan, nrow)
54 : integer rflag(nrow)
55 : real weight(nvischan, nrow), phase
56 : double precision sumwt(npol, nchan)
57 : integer rownum
58 : integer support(wconvsize,polused,npa), rsupport
59 : integer chanmap(nchan), polmap(npol),cfmap(npol),conjcfmap(npol)
60 : integer dopsf
61 :
62 : complex nvalue
63 :
64 : C complex convfunc(convsize/2-1, convsize/2-1, wconvsize, polused),
65 : complex convfunc(convsize, convsize, wconvsize, polused),
66 : $ cwt
67 :
68 : real norm
69 : complex cnorm,tcnorm
70 : real wt
71 :
72 : logical opbwproj,reindex
73 : external nwcppeij
74 :
75 : real pos(3)
76 : integer loc(3), off(3), iloc(3),iu,iv
77 : integer rbeg, rend
78 : integer ix, iy, ipol, ichan
79 : integer apol, achan, irow, PolnPlane, ConjPlane
80 : double precision pi,tc,ts,cfscale,ixr,iyr
81 : data pi/3.14159265358979323846/
82 :
83 : double precision griduvw(2),dPA, cDPA,sDPA
84 : double precision mysigma, ra1,ra2,dec1,dec2
85 : double precision sigma,area,lambda
86 : complex pcwt, pdcwt1, pdcwt2
87 : integer nant, scanno, ant1(nrow), ant2(nrow)
88 : integer convOrigin
89 : real raoff(nant), decoff(nant)
90 :
91 :
92 : integer ii,jj,kk,imax,jmax
93 : real tmp
94 :
95 0 : irow=rownum
96 :
97 0 : if(irow.ge.0) then
98 0 : rbeg=irow+1
99 0 : rend=irow+1
100 : else
101 0 : rbeg=1
102 0 : rend=nrow
103 : end if
104 0 : cfscale=1.0
105 0 : dPA = -(currentCFPA - actualPA)
106 : c dPA = 0
107 0 : cDPA = cos(dPA)
108 0 : sDPA = sin(dPA)
109 :
110 0 : convOrigin = (convsize-1)/2
111 : c convOrigin = (convsize+1)/2
112 : c$$$ write(*,*) convOrigin, convsize, (convsize/2),
113 : c$$$ $ convfunc(convOrigin,convOrigin,1,1)
114 0 : do irow=rbeg, rend
115 0 : if(rflag(irow).eq.0) then
116 0 : do ichan=1, nvischan
117 0 : achan=chanmap(ichan)+1
118 :
119 0 : lambda = c/freq(ichan)
120 :
121 0 : if((achan.ge.1).and.(achan.le.nchan).and.
122 0 : $ (weight(ichan,irow).ne.0.0)) then
123 : call spbwproj(uvw(1,irow), dphase(irow),
124 : $ freq(ichan), c, scale, offset, sampling,
125 0 : $ pos, loc, off, phasor)
126 0 : iloc(3)=max(1, min(wconvsize, loc(3)))
127 0 : rsupport=support(iloc(3),1,paindex)
128 0 : if (opbwproj(nx, ny, wconvsize, loc, rsupport)) then
129 0 : PolnPlane=polused+1
130 :
131 0 : do ipol=1, nvispol
132 0 : apol=polmap(ipol)+1
133 :
134 : if((flag(ipol,ichan,irow).ne.1).and.
135 0 : $ (apol.ge.1).and.(apol.le.npol)) then
136 :
137 : c$$$ ConjPlane = conjcfmap(ipol)+1
138 : c$$$ PolnPlane = cfmap(ipol)+1
139 : C
140 : C The following after feed_x -> -feed_x and PA -> PA + PI/2
141 0 : ConjPlane = cfmap(ipol)+1
142 0 : PolnPlane = conjcfmap(ipol)+1
143 : C
144 : C If we are making a PSF then we don't want to phase rotate but we do want
145 : C to reproject uvw
146 : C
147 0 : if(dopsf.eq.1) then
148 0 : nvalue=cmplx(weight(ichan,irow))
149 : else
150 : nvalue=weight(ichan,irow)*
151 0 : $ (values(ipol,ichan,irow)*phasor)
152 : end if
153 : C
154 : C norm will be the value we would get for the peak at the phase center.
155 : C We will want to normalize the final image by this term.
156 : C
157 0 : norm=0.0
158 0 : cnorm=cmplx(1.0,0.0)
159 0 : cnorm=0.0
160 0 : tcnorm=0.0
161 0 : do iy=-rsupport,rsupport
162 : c$$$ iloc(2)=1+(iy*sampling+off(2))
163 : c$$$ $ +convOrigin
164 : c$$$ ii=iloc(2)
165 0 : iloc(2)=(iy*sampling+off(2))
166 0 : iv = iloc(2)
167 0 : do ix=-rsupport,rsupport
168 : c$$$ iloc(1)=1+(ix*sampling+off(1))
169 : c$$$ $ +convOrigin
170 : c$$$ jj=iloc(1)
171 0 : iloc(1)=(ix*sampling+off(1))
172 :
173 0 : iu = iloc(1)
174 :
175 : c$$$ ixr = ix*sampling*cfscale
176 : c$$$ iyr = iy*sampling*cfscale
177 : c$$$ iu = nint(cDPA*ixr + sDPA*iyr)
178 : c$$$ iv = nint(-sDPA*ixr + cDPA*iyr)
179 : c$$$ iu=iu+off(1)*cfscale
180 : c$$$ iv=iv+off(2)*cfscale
181 :
182 0 : ixr = (ix*sampling+off(1))*cfscale
183 0 : iyr = (iy*sampling+off(2))*cfscale
184 0 : iu = nint(cDPA*ixr + sDPA*iyr)
185 0 : iv = nint(-sDPA*ixr + cDPA*iyr)
186 :
187 0 : ts=0.0
188 0 : tc=1.0
189 0 : if (reindex(iu,iv,iloc(1),iloc(2),
190 : $ ts,tc,
191 0 : $ convOrigin, convSize)) then
192 0 : if (dopointingcorrection .eq. 1)
193 : $ then
194 : c$$$ griduvw(1) = (loc(1)-offset(1)+
195 : c$$$ $ ix-1)
196 : c$$$ $ /scale(1)-uvw(1,irow)/lambda
197 : c$$$ griduvw(2) = (loc(2)-offset(2)+
198 : c$$$ $ iy-1)
199 : c$$$ $ /scale(2)-uvw(2,irow)/lambda
200 0 : iu=iu-off(1)*cfscale
201 0 : iv=iv-off(2)*cfscale
202 0 : griduvw(1)=(iu)/(scale(1)*sampling)
203 0 : griduvw(2)=(iv)/(scale(2)*sampling)
204 0 : ra1 = raoff(ant1(irow)+1)
205 0 : ra2 = raoff(ant2(irow)+1)
206 0 : dec1= decoff(ant1(irow)+1)
207 0 : dec2= decoff(ant2(irow)+1)
208 : call nwcppeij(griduvw,area,
209 : $ ra1,dec1,ra2,dec2,
210 : $ dograd,pcwt,pdcwt1,pdcwt2,
211 0 : $ currentCFPA)
212 : else
213 0 : pcwt=cmplx(1.0,0.0)
214 : endif
215 : c$$$ if(uvw(3,irow).gt.0.0) then
216 : c$$$ cwt=conjg(convfunc(iloc(1),
217 : c$$$ $ iloc(2), iloc(3),
218 : c$$$ $ ConjPlane))
219 : c$$$ pcwt=conjg(pcwt)
220 : c$$$ else
221 : c$$$ cwt=(convfunc(iloc(1),
222 : c$$$ $ iloc(2), iloc(3),
223 : c$$$ $ PolnPlane))
224 : c$$$c pcwt=conjg(pcwt)
225 : c$$$ end if
226 :
227 : cwt=(convfunc(iloc(1),
228 : $ iloc(2), iloc(3),
229 0 : $ PolnPlane))
230 0 : cwt = cwt * pcwt
231 0 : if (dopsf .eq. 1) then
232 : C cnorm=cnorm+real(cwt)
233 0 : cnorm=cnorm+cwt
234 : else
235 0 : cnorm=cnorm+cwt
236 : endif
237 : endif
238 : enddo
239 : enddo
240 : C stop
241 : c$$$ else
242 : c$$$ cnorm=cmplx(1.0,0.0)
243 : c$$$ endif
244 :
245 :
246 0 : do iy=-rsupport,rsupport
247 0 : iloc(2)=(iy*sampling+off(2))
248 0 : iv = iloc(2)
249 0 : do ix=-rsupport,rsupport
250 0 : iloc(1)=(ix*sampling+off(1))
251 0 : iu = iloc(1)
252 :
253 : c$$$ ixr = ix*sampling*cfscale
254 : c$$$ iyr = iy*sampling*cfscale
255 : c$$$ iu = nint(cDPA*ixr + sDPA*iyr)
256 : c$$$ iv = nint(-sDPA*ixr + cDPA*iyr)
257 : c$$$ iu=iu+off(1)*cfscale
258 : c$$$ iv=iv+off(2)*cfscale
259 :
260 0 : ixr = (ix*sampling+off(1))*cfscale
261 0 : iyr = (iy*sampling+off(2))*cfscale
262 0 : iu = nint(cDPA*ixr + sDPA*iyr)
263 0 : iv = nint(-sDPA*ixr + cDPA*iyr)
264 :
265 0 : ts=0.0
266 0 : tc=1.0
267 0 : if (reindex(iu,iv,iloc(1),iloc(2),
268 : $ ts,tc,
269 0 : $ convOrigin,convSize)) then
270 : C
271 : C Compute the pointing offset term
272 : C
273 0 : if (dopointingcorrection .eq. 1)
274 : $ then
275 : c$$$ iu=iu-off(1)*cfscale
276 : c$$$ iv=iv-off(2)*cfscale
277 : C
278 : C Use the rotated CF co-oridinates to compute the phase grad. This
279 : C effectively rotates the pointing vector with the PA
280 : C
281 0 : iu = iloc(1)-convOrigin
282 0 : iv = iloc(2)-convOrigin
283 : griduvw(1)=(iu)/
284 0 : $ (scale(1)*sampling)
285 : griduvw(2)=(iv)/
286 0 : $ (scale(2)*sampling)
287 0 : ra1 = raoff(ant1(irow)+1)
288 0 : ra2 = raoff(ant2(irow)+1)
289 0 : dec1= decoff(ant1(irow)+1)
290 0 : dec2= decoff(ant2(irow)+1)
291 : call nwcppeij(griduvw,area,
292 : $ ra1,dec1,ra2,dec2,
293 : $ dograd,pcwt,pdcwt1,pdcwt2,
294 0 : $ currentCFPA)
295 : else
296 0 : pcwt=cmplx(1.0,0.0)
297 : endif
298 :
299 : c$$$ if(uvw(3,irow).gt.0.0) then
300 : c$$$ cwt=conjg(convfunc(iloc(1),
301 : c$$$ $ iloc(2), iloc(3),
302 : c$$$ $ ConjPlane))
303 : c$$$ pcwt=conjg(pcwt)
304 : c$$$ else
305 : c$$$ cwt=(convfunc(iloc(1),
306 : c$$$ $ iloc(2), iloc(3),
307 : c$$$ $ PolnPlane))
308 : c$$$ pcwt=(pcwt)
309 : c$$$ end if
310 : cwt=(convfunc(iloc(1),
311 : $ iloc(2), iloc(3),
312 0 : $ PolnPlane))
313 0 : cwt = cwt * pcwt
314 : c$$$ if (dopsf .eq. 1) then
315 : c$$$ cwt = real(cwt)
316 : c$$$ endif
317 : grid(loc(1)+ix,
318 : $ loc(2)+iy,apol,achan)=
319 : $ grid(loc(1)+ix,
320 : $ loc(2)+iy,apol,achan)+
321 : $ nvalue*cwt
322 0 : $ /cnorm
323 : c norm=norm+real(cwt/cnorm)
324 : c$$$ write(*,*) 'G ',abs(cwt),
325 : c$$$ $ ix,iy,iloc(1),iloc(2)
326 0 : tcnorm=tcnorm+(cwt/cnorm)
327 : c cnorm=cnorm+cwt
328 : endif
329 : end do
330 : c$$$ write(*,*)
331 : end do
332 : c$$$ stop
333 0 : norm=real(tcnorm)
334 : sumwt(apol,achan)=sumwt(apol,achan)+
335 0 : $ weight(ichan,irow)*norm
336 : end if
337 : end do
338 : c else
339 : end if
340 : end if
341 : end do
342 : end if
343 : end do
344 : c write(*,*)'Gridded ',irow,' visibilities'
345 : c$$$ close(11)
346 : c$$$ close(12)
347 0 : return
348 : end
349 : C
350 : C Degrid a number of visibility records
351 : C
352 0 : subroutine dpbwproj (uvw, dphase, values, nvispol, nvischan,
353 0 : $ flag, rflag,
354 0 : $ nrow, rownum, scale, offset, grid, nx, ny, npol, nchan, freq,
355 0 : $ c, support, convsize, sampling, wconvsize, convfunc,
356 0 : $ chanmap, polmap,polused, ant1, ant2, nant,
357 0 : $ scanno, sigma, raoff, decoff,area,dograd,
358 0 : $ dopointingcorrection,npa,paindex,cfmap,conjcfmap,
359 : $ currentCFPA,actualPA,cfRefFreq)
360 :
361 : implicit none
362 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow,polused
363 : integer npa,paindex
364 : integer convsize, wconvsize, sampling
365 : complex values(nvispol, nvischan, nrow)
366 : complex grid(nx, ny, npol, nchan)
367 : double precision uvw(3, nrow), freq(nvischan), c, scale(3),
368 : $ offset(3),currentCFPA,actualPA, dPA, sDPA, cDPA,cfRefFreq
369 : double precision dphase(nrow), uvdist
370 : complex phasor
371 : integer flag(nvispol, nvischan, nrow)
372 : integer rflag(nrow),cfmap(npol),conjcfmap(npol)
373 : integer rownum
374 : integer support(wconvsize,polused,npa), rsupport
375 : integer chanmap(*), polmap(*)
376 :
377 : integer nant, scanno, ant1(nrow), ant2(nrow),dograd,
378 : $ dopointingcorrection
379 : real raoff(nant), decoff(nant)
380 : double precision sigma,area,lambda
381 :
382 : complex nvalue
383 :
384 : C complex convfunc(convsize/2-1, convsize/2-1, wconvsize, polused),
385 : complex convfunc(convsize, convsize, wconvsize, polused),
386 : $ cwt, pcwt,pdcwt1,pdcwt2
387 : double precision griduvw(2)
388 : double precision mysigma, ra1,ra2,dec1,dec2
389 :
390 : complex norm(4)
391 :
392 : logical opbwproj,reindex
393 : external nwcppEij
394 :
395 : real pos(3)
396 : integer loc(3), off(3), iloc(3)
397 : integer rbeg, rend
398 : integer ix, iy, ipol, ichan, PolnPlane, ConjPlane
399 : integer convOrigin
400 : integer apol, achan, irow
401 : real wt, wtx, wty
402 : double precision pi,ts,tc,cfscale,ixr,iyr
403 : data pi/3.14159265358979323846/
404 : integer ii,jj,iu,iv
405 :
406 : complex tmp
407 :
408 0 : irow=rownum
409 :
410 0 : if(irow.ge.0) then
411 0 : rbeg=irow+1
412 0 : rend=irow+1
413 : else
414 0 : rbeg=1
415 0 : rend=nrow
416 : end if
417 : C
418 :
419 0 : cfscale=1.0
420 0 : dPA = -(currentCFPA - actualPA)
421 : c dPA=0
422 0 : cDPA = cos(dPA)
423 0 : sDPA = sin(dPA)
424 0 : convOrigin = (convsize+1)/2
425 0 : convOrigin = (convsize-1)/2
426 0 : convOrigin = (convsize)/2
427 :
428 0 : do irow=rbeg, rend
429 0 : if(rflag(irow).eq.0) then
430 0 : do ichan=1, nvischan
431 0 : achan=chanmap(ichan)+1
432 :
433 0 : lambda = c/freq(ichan)
434 :
435 0 : if((achan.ge.1).and.(achan.le.nchan)) then
436 : call spbwproj(uvw(1,irow), dphase(irow), freq(ichan), c,
437 0 : $ scale, offset, sampling, pos, loc, off, phasor)
438 0 : iloc(3)=max(1, min(wconvsize, loc(3)))
439 0 : rsupport=support(iloc(3),1,paindex)
440 : c$$$ off(1)=0
441 : c$$$ off(2)=0
442 0 : if (opbwproj(nx, ny, wconvsize, loc, rsupport)) then
443 0 : PolnPlane=0
444 : if (ant1(irow).eq. 1 .and.
445 : $ ant2(irow).eq. 2) then
446 : c$$$ write(*,*) ichan,off(1), off(2),loc(1),loc(2),
447 : c$$$ $ sampling,scale(1),uvw(1,irow),uvw(2,irow)
448 : endif
449 :
450 0 : do ipol=1, nvispol
451 0 : apol=polmap(ipol)+1
452 : if((flag(ipol,ichan,irow).ne.1).and.
453 0 : $ (apol.ge.1).and.(apol.le.npol)) then
454 :
455 0 : ConjPlane = cfmap(ipol)+1
456 0 : PolnPlane = conjcfmap(ipol)+1
457 :
458 0 : nvalue=0.0
459 0 : norm(apol)=cmplx(0.0,0.0)
460 0 : pcwt=cmplx(1.0,0.0)
461 0 : ii=0
462 0 : do iy=-rsupport,rsupport
463 0 : iloc(2)=1+(iy*sampling+off(2))+convOrigin
464 0 : iv = (iy*sampling+off(2))
465 0 : do ix=-rsupport,rsupport
466 0 : iloc(1)=1+(ix*sampling+off(1))+convOrigin
467 0 : iu = (ix*sampling+off(1))
468 :
469 : c$$$ ixr = ix*sampling*cfscale
470 : c$$$ iyr = iy*sampling*cfscale
471 : c$$$ iu = nint(cDPA*ixr + sDPA*iyr)
472 : c$$$ iv = nint(-sDPA*ixr + cDPA*iyr)
473 : c$$$ iu=iu+off(1)*cfscale
474 : c$$$ iv=iv+off(2)*cfscale
475 :
476 0 : ixr = (ix*sampling+off(1))*cfscale
477 0 : iyr = (iy*sampling+off(2))*cfscale
478 : c$$$ iu = nint(cDPA*ixr + sDPA*iyr)
479 : c$$$ iv = nint(-sDPA*ixr + cDPA*iyr)
480 0 : iu=ixr
481 0 : iv=iyr
482 0 : ts=0.0
483 0 : tc=1.0
484 0 : ts=sDPA
485 0 : tc=cDPA
486 0 : if(reindex(iu,iv,iloc(1),iloc(2),
487 : $ ts,tc, convOrigin, convSize))
488 0 : $ then
489 :
490 0 : if (dopointingcorrection .eq. 1) then
491 : c$$$ iu = iu - off(1)*cfscale
492 : c$$$ iv = iv - off(2)*cfscale
493 : C
494 : C Use the rotated CF co-oridinates to compute the phase grad. This
495 : C effectively rotates the pointing vector with the PA
496 : C
497 0 : iu = iloc(1)-convOrigin
498 0 : iv = iloc(2)-convOrigin
499 0 : griduvw(1)=(iu)/(scale(1)*sampling)
500 0 : griduvw(2)=(iv)/(scale(2)*sampling)
501 0 : ra1 = raoff(ant1(irow)+1)
502 0 : ra2 = raoff(ant2(irow)+1)
503 0 : dec1= decoff(ant1(irow)+1)
504 0 : dec2= decoff(ant2(irow)+1)
505 : call nwcppEij(griduvw,area,
506 : $ ra1,dec1,ra2,dec2,
507 : $ dograd,pcwt,pdcwt1,pdcwt2,
508 0 : $ currentCFPA)
509 : endif
510 :
511 : c$$$ if(uvw(3,irow).gt.0.0) then
512 : c$$$ cwt=conjg(convfunc(iloc(1),
513 : c$$$ $ iloc(2), iloc(3),ConjPlane))
514 : c$$$ pcwt = conjg(pcwt)
515 :
516 :
517 : c$$$ else
518 : c$$$ cwt=(convfunc(iloc(1),
519 : c$$$ $ iloc(2), iloc(3),PolnPlane))
520 : c$$$c pcwt = conjg(pcwt)
521 : c$$$ endif
522 : cwt=(convfunc(iloc(1),
523 0 : $ iloc(2), iloc(3),PolnPlane))
524 :
525 : nvalue=nvalue+(cwt)*(pcwt)*
526 0 : $ grid(loc(1)+ix,loc(2)+iy,apol,achan)
527 0 : norm(apol)=norm(apol)+(pcwt*cwt)
528 : c ii=ii+1
529 : c if (irow .eq. 2) then
530 : c$$$ tmp=grid(loc(1)+ix,loc(2)+iy,
531 : c$$$ $ apol,achan)
532 : c$$$ if (ant1(irow) .eq. 1 .and.
533 : c$$$ $ ant2(irow) .eq. 2) then
534 : c$$$ write(*,*)abs(nvalue),
535 : c$$$ $ atan2(imag(nvalue),
536 : c$$$ $ real(nvalue)),
537 : c$$$ $ abs(tmp),
538 : c$$$ $ atan2(imag(tmp),real(tmp)),
539 : c$$$ $ ix,iy,
540 : c$$$ $ iloc(1), iloc(2),
541 : c$$$ $ abs(cwt),
542 : c$$$ $ norm(apol),
543 : c$$$ $ ant1(irow),ant2(irow),
544 : c$$$ $ irow,achan,apol,
545 : c$$$ $ off(1), off(2)
546 : c$$$ endif
547 : endif
548 : end do
549 : c$$$ write(*,*)
550 : end do
551 : c norm(apol) = norm(apol)/abs(norm(apol))
552 : c$$$ write (*,*) nvalue, norm(apol)
553 : values(ipol,ichan,irow)=
554 : $ nvalue*conjg(phasor)
555 0 : $ /norm(apol)
556 : c$$$ if ((ant1(irow).eq.1).and.
557 : c$$$ $ (ant2(irow).eq.2)) then
558 : c$$$ write(*,*) irow,apol,ipol,
559 : c$$$ $ abs(values(ipol,ichan,irow)),
560 : c$$$ $ atan2(imag(values(ipol,ichan,irow)),
561 : c$$$ $ real(values(ipol,ichan,irow)))*57.2956,
562 : c$$$ $ real(norm(apol)),imag(norm(apol)),ii
563 : c$$$ endif
564 : end if
565 : end do
566 : c$$$ write(*,*) "====================================="
567 : c$$$ stop
568 : end if
569 : end if
570 : end do
571 : end if
572 : end do
573 0 : return
574 : end
575 : C
576 : C Degrid a number of visibility records along with the grad. computations
577 : C
578 0 : subroutine dpbwgrad (uvw, dphase, values, nvispol, nvischan,
579 0 : $ gazvalues, gelvalues, doconj,
580 0 : $ flag, rflag, nrow, rownum, scale, offset, grid,
581 0 : $ nx, ny, npol, nchan, freq, c, support, convsize, sampling,
582 0 : $ wconvsize, convfunc, chanmap, polmap,polused,ant1,ant2,nant,
583 0 : $ scanno, sigma, raoff, decoff,area,dograd,
584 0 : $ dopointingcorrection,npa,paindex,cfmap,conjcfmap,
585 : $ currentCFPA,actualPA,cfRefFreq)
586 :
587 : implicit none
588 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow,polused
589 : integer npa,paindex,doconj
590 : integer convsize, wconvsize, sampling
591 : complex values(nvispol, nvischan, nrow)
592 : complex gazvalues(nvispol, nvischan, nrow)
593 : complex gelvalues(nvispol, nvischan, nrow)
594 : complex grid(nx, ny, npol, nchan)
595 : double precision uvw(3, nrow), freq(nvischan), c, scale(3),
596 : $ offset(3),currentCFPA,actualPA, dPA, sDPA, cDPA,cfRefFreq
597 : double precision dphase(nrow), uvdist
598 : complex phasor
599 : integer flag(nvispol, nvischan, nrow)
600 : integer rflag(nrow),cfmap(npol),conjcfmap(npol)
601 : integer rownum
602 : integer support(wconvsize,polused,npa), rsupport
603 : integer chanmap(*), polmap(*)
604 :
605 : integer nant, scanno, ant1(nrow), ant2(nrow),dograd,
606 : $ dopointingcorrection
607 : real raoff(2,1,nant), decoff(2,1,nant)
608 : double precision sigma,area,lambda
609 :
610 : complex nvalue,ngazvalue,ngelvalue
611 :
612 : C complex convfunc(convsize/2-1, convsize/2-1, wconvsize, polused),
613 : complex convfunc(convsize, convsize, wconvsize, polused),
614 : $ cwt, pcwt,pdcwt1,pdcwt2
615 : double precision griduvw(2)
616 : double precision mysigma, ra1,ra2,dec1,dec2
617 :
618 : complex norm(4),gradaznorm(4),gradelnorm(4)
619 :
620 : logical opbwproj,reindex
621 : external nwcppEij
622 :
623 : real pos(3)
624 : integer loc(3), off(3), iloc(3)
625 : integer rbeg, rend
626 : integer ix, iy, ipol, ichan, PolnPlane, ConjPlane
627 : integer convOrigin
628 : integer apol, achan, irow
629 : real wt, wtx, wty
630 : double precision pi,ts,tc,cfscale,ixr,iyr
631 : data pi/3.14159265358979323846/
632 : integer ii,iu,iv,tt
633 :
634 : complex tmp
635 0 : tt=0
636 0 : irow=rownum
637 :
638 0 : if(irow.ge.0) then
639 0 : rbeg=irow+1
640 0 : rend=irow+1
641 : else
642 0 : rbeg=1
643 0 : rend=nrow
644 : end if
645 : C
646 :
647 0 : cfscale=1.0
648 0 : dPA = -(currentCFPA - actualPA)
649 : c dPA = 0
650 0 : cDPA = cos(dPA)
651 0 : sDPA = sin(dPA)
652 : c convOrigin = (convsize+1)/2
653 0 : convOrigin = (convsize-1)/2
654 0 : convOrigin = (convsize)/2
655 :
656 0 : do irow=rbeg, rend
657 0 : if(rflag(irow).eq.0) then
658 0 : do ichan=1, nvischan
659 0 : achan=chanmap(ichan)+1
660 :
661 0 : lambda = c/freq(ichan)
662 :
663 0 : if((achan.ge.1).and.(achan.le.nchan)) then
664 : call spbwproj(uvw(1,irow), dphase(irow), freq(ichan), c,
665 0 : $ scale, offset, sampling, pos, loc, off, phasor)
666 0 : iloc(3)=max(1, min(wconvsize, loc(3)))
667 0 : rsupport=support(iloc(3),1,paindex)
668 0 : if (opbwproj(nx, ny, wconvsize, loc, rsupport)) then
669 :
670 0 : do ipol=1, nvispol
671 0 : apol=polmap(ipol)+1
672 : if((flag(ipol,ichan,irow).ne.1).and.
673 0 : $ (apol.ge.1).and.(apol.le.npol)) then
674 0 : ConjPlane = conjcfmap(ipol)+1
675 0 : PolnPlane = cfmap(ipol)+1
676 0 : ConjPlane = cfmap(ipol)+1
677 0 : PolnPlane = conjcfmap(ipol)+1
678 0 : nvalue=0.0
679 0 : ngazvalue = 0.0
680 0 : ngelvalue = 0.0
681 :
682 0 : norm(apol)=cmplx(0.0,0.0)
683 0 : gradaznorm(apol)=cmplx(0.0,0.0)
684 0 : gradelnorm(apol)=cmplx(0.0,0.0)
685 0 : pcwt=cmplx(1.0,0.0)
686 :
687 0 : do iy=-rsupport,rsupport
688 0 : iloc(2)=1+(iy*sampling+off(2))+convOrigin
689 0 : iv=(iy*sampling+off(2))
690 0 : do ix=-rsupport,rsupport
691 0 : iloc(1)=1+(ix*sampling+off(1))+convOrigin
692 0 : iu=(ix*sampling+off(1))
693 :
694 0 : ixr = (ix*sampling+off(1))*cfscale
695 0 : iyr = (iy*sampling+off(2))*cfscale
696 0 : ts=sDPA
697 0 : tc=cDPA
698 0 : iu=ixr
699 0 : iv=iyr
700 0 : if(reindex(iu,iv,iloc(1),iloc(2),
701 : $ ts,tc, convOrigin, convSize))
702 : c$$$ if(reindex(ixr,iyr,iloc(1),iloc(2),
703 : c$$$ $ ts,tc, convOrigin, convSize))
704 0 : $ then
705 0 : if (dopointingcorrection .eq. 1) then
706 : c$$$ iu=iu-off(1)*cfscale
707 : c$$$ iv=iv-off(2)*cfscale
708 : C
709 : C Use the rotated CF co-oridinates to compute the phase grad. This
710 : C effectively rotates the pointing vector with the PA
711 : C
712 0 : iu=iloc(1)-convOrigin
713 0 : iv=iloc(2)-convOrigin
714 0 : griduvw(1)=(iu)/(scale(1)*sampling)
715 0 : griduvw(2)=(iv)/(scale(2)*sampling)
716 0 : ii = PolnPlane
717 : c ii=apol
718 0 : ra1 = raoff(ii,1,ant1(irow)+1)
719 0 : ra2 = raoff(ii,1,ant2(irow)+1)
720 0 : dec1= decoff(ii,1,ant1(irow)+1)
721 0 : dec2= decoff(ii,1,ant2(irow)+1)
722 : call nwcppEij(griduvw,area,
723 : $ ra1,dec1,ra2,dec2,
724 : $ dograd,pcwt,pdcwt1,pdcwt2,
725 0 : $ currentCFPA)
726 : endif
727 : c$$$ if(uvw(3,irow).gt.0.0) then
728 : c$$$ cwt=conjg(convfunc(iloc(1),
729 : c$$$ $ iloc(2), iloc(3),ConjPlane))
730 : c$$$ pcwt = conjg(pcwt)
731 : c$$$ pdcwt1 = conjg(pdcwt1)
732 : c$$$ pdcwt2 = conjg(pdcwt2)
733 : c$$$ else
734 : c$$$ cwt=(convfunc(iloc(1),
735 : c$$$ $ iloc(2), iloc(3),PolnPlane))
736 : c$$$ endif
737 : cwt=(convfunc(iloc(1),
738 0 : $ iloc(2), iloc(3),PolnPlane))
739 :
740 0 : cwt = cwt*pcwt
741 : nvalue=nvalue+(cwt)*
742 : $ grid(loc(1)+ix,loc(2)+iy,apol,
743 0 : $ achan)
744 :
745 0 : if ((doconj .eq. 1) .and.
746 : $ (dograd .eq. 1)) then
747 0 : pdcwt1 = conjg(pdcwt1)
748 0 : pdcwt2 = conjg(pdcwt2)
749 : endif
750 0 : if (dograd .eq. 1) then
751 0 : pdcwt1 = pdcwt1*cwt
752 0 : pdcwt2 = pdcwt2*cwt
753 : ngazvalue=ngazvalue+(pdcwt1)*
754 : $ (grid(loc(1)+ix,loc(2)+iy,
755 0 : $ apol,achan))
756 : ngelvalue=ngelvalue+(pdcwt2)*
757 : $ (grid(loc(1)+ix,loc(2)+iy,
758 0 : $ apol,achan))
759 : endif
760 :
761 0 : norm(apol)=norm(apol)+(cwt)
762 : gradaznorm(apol)=gradaznorm(apol)+
763 0 : $ pdcwt1
764 : gradelnorm(apol)=gradelnorm(apol)+
765 0 : $ pdcwt1
766 : c$$$ if (apol .eq. 1) then
767 : c$$$ write(*,*)ix,iy,
768 : c$$$ $ abs(cwt),atan2(imag(cwt),real(cwt)),
769 : c$$$ $ abs(pdcwt1),
770 : c$$$ $ atan2(imag(pdcwt1),real(pdcwt1)),
771 : c$$$ $ abs(pdcwt2),
772 : c$$$ $ atan2(imag(pdcwt2),real(pdcwt2))
773 : c$$$ endif
774 : endif
775 : end do
776 : c$$$ write(*,*)
777 : end do
778 : c norm(apol) = norm(apol)/abs(norm(apol))
779 : values(ipol,ichan,irow)=
780 : $ nvalue*conjg(phasor)
781 0 : $ /norm(apol)
782 : c$$$ if (ant1(irow) .eq. 3 .and.
783 : c$$$ $ ant2(irow) .eq. 7) then
784 : c$$$ write(*,*)irow,ant1(irow),
785 : c$$$ $ ant2(irow),ipol,ichan,
786 : c$$$ $ values(ipol,ichan,irow)
787 : c$$$ stop
788 : c$$$ endif
789 0 : if (dograd .eq. 1) then
790 : gazvalues(ipol,ichan,irow)=ngazvalue*
791 0 : $ conjg(phasor)/norm(apol)
792 : C $ conjg(phasor)/gradaznorm(apol)
793 : gelvalues(ipol,ichan,irow)=ngelvalue*
794 0 : $ conjg(phasor)/norm(apol)
795 : C $ conjg(phasor)/gradelnorm(apol)
796 : endif
797 : end if
798 : end do
799 0 : tt=0
800 : c$$$ stop
801 : end if
802 : end if
803 : end do
804 : end if
805 : end do
806 0 : return
807 : end
808 : C
809 : C Calculate gridded coordinates and the phasor needed for
810 : C phase rotation.
811 : C
812 0 : subroutine spbwproj (uvw, dphase, freq, c, scale, offset,
813 : $ sampling, pos, loc, off, phasor)
814 : implicit none
815 : integer loc(3), off(3), sampling
816 : double precision uvw(3), freq, c, scale(3), offset(3)
817 : real pos(3)
818 : double precision dphase, phase
819 : complex phasor
820 : integer idim
821 : double precision pi
822 : data pi/3.14159265358979323846/
823 :
824 : C pos(3)=(scale(3)*uvw(3)*freq/c)*(scale(3)*uvw(3)*freq/c)
825 : C $ +offset(3)+1.0
826 : C pos(3)=(scale(3)*uvw(3)*freq/c)+offset(3)+1.0
827 0 : pos(3)=sqrt(abs(scale(3)*uvw(3)*freq/c))+offset(3)+1.0
828 0 : loc(3)=nint(pos(3))
829 0 : off(3)=0
830 :
831 0 : do idim=1,2
832 : pos(idim)=scale(idim)*uvw(idim)*freq/c+
833 0 : $ (offset(idim)+1.0)
834 0 : loc(idim)=nint(pos(idim))
835 0 : off(idim)=nint((loc(idim)-pos(idim))*sampling)
836 : end do
837 :
838 0 : phase=-2.0D0*pi*dphase*freq/c
839 0 : phasor=cmplx(cos(phase), sin(phase))
840 0 : return
841 : end
842 0 : logical function opbwproj (nx, ny, nw, loc, support)
843 : implicit none
844 : integer nx, ny, nw, loc(3), support
845 : opbwproj=(support.gt.0).and.
846 : $ (loc(1)-support.ge.1).and.(loc(1)+support.le.nx).and.
847 : $ (loc(2)-support.ge.1).and.(loc(2)+support.le.ny).and.
848 0 : $ (loc(3).ge.1).and.(loc(3).le.nw)
849 0 : return
850 : end
851 :
852 0 : logical function reindex(inx,iny,outx,outy,sinDPA,cosDPA,
853 : $ Origin, Size)
854 : integer inx,iny,outx,outy, Origin, Size
855 : double precision sinDPA, cosDPA
856 : integer ix,iy
857 :
858 0 : ix = nint(cosDPA*inx + sinDPA*iny)+1
859 0 : iy = nint(-sinDPA*inx + cosDPA*iny)+1
860 : c$$$ ix = nint(cosDPA*inx + sinDPA*iny)
861 : c$$$ iy = nint(-sinDPA*inx + cosDPA*iny)
862 0 : outx = ix+Origin
863 0 : outy = iy+Origin
864 :
865 : c$$$ outx=ix
866 : c$$$ outy=iy
867 : reindex=(outx .ge. 1 .and.
868 : $ outx .le. Size .and.
869 : $ outy .ge. 1 .and.
870 0 : $ outy .le. Size)
871 :
872 0 : return
873 : end
|
