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: 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: fpbwproj.f,v 1.13 2006/07/20 00:24:20 sbhatnag Exp $
29 : *-----------------------------------------------------------------------
30 :
31 0 : subroutine gfeij(uvw,area,ra1,dec1,ra2,dec2,dograd,pcwt,pdcwt1,
32 : $ pdcwt2,pa)
33 : implicit none
34 : double precision uvw(2), area, pa,pi,ra1,ra2,dec1,dec2
35 : integer dograd
36 : complex pcwt,pdcwt1,pdcwt2
37 : complex kij
38 : data pi/3.14159265358979323846/
39 :
40 0 : kij = cmplx(0,pi*(uvw(1)*(ra1+ra2) + uvw(2)*(dec1+dec2)))
41 0 : pcwt=exp(kij)/area
42 0 : pdcwt1 = cmplx(1,0)
43 0 : pdcwt2 = cmplx(1,0)
44 :
45 0 : if (dograd .gt. 0) then
46 0 : pdcwt1 = (cmplx(0,pi*uvw(1)))
47 0 : pdcwt2 = (cmplx(0,pi*uvw(2)))
48 : endif
49 0 : return
50 : end
51 :
52 0 : complex function accumulate(doavgpb,area, grid,gridorigin,func,
53 : $ nvalue, cnorm, rsupport,sampling, off,
54 : $ cfscale,
55 : $ scale,lambda, sDPA,cDPA,currentCFPA,convOrigin,convSize,
56 : $ wconvsize, polused, dopointingcorrection, dograd,nant,
57 0 : $ raoff,decoff, ConjPlane, PolnPlane, nx,ny,npol,nchan,
58 0 : $ apol,achan,ant1,ant2,uvw,nrow,irow)
59 :
60 :
61 : integer nrow, irow,doavgpb
62 : integer ant1(nrow),ant2(nrow)
63 : integer gridorigin(3)
64 : integer rsupport,sampling, off(3), convOrigin, convSize,wconvsize,
65 : $ polused
66 : integer dopointingcorrection, dograd,nant
67 : integer ix,iy,iloc(3),loc(3),origin(3),
68 : $ ConjPlane, PolnPlane,apol,achan
69 : double precision offset(3),cfscale,lambda,sDPA,cDPA,currentCFPA,
70 : $ area
71 : real raoff(2,1,nant), decoff(2,1,nant)
72 : double precision griduvw(2), scale(3),uvw(3,nrow),
73 : $ ra1,ra2,dec1,dec2
74 : complex cwt,pcwt,pdcwt1,pdcwt2,tcnorm
75 : complex func(convsize, convsize, wconvsize, polused)
76 : complex grid(nx, ny, npol, nchan),nvalue,cnorm
77 : integer ixr,iyr
78 : logical opbmos,mreindex
79 : external gcppeij
80 : double precision ts,tc,tmpnorm
81 : double precision scaledSampling
82 : integer scaledSupport
83 :
84 0 : scaledSampling=int(sampling*cfscale)
85 0 : scaledSupport=int(rsupport/cfscale+0.5)
86 :
87 0 : origin(1)=gridorigin(1)
88 0 : origin(2)=gridorigin(2)
89 0 : origin(3)=gridorigin(3)
90 0 : tcnorm=0
91 0 : iloc(3)=1
92 :
93 0 : do iy=-scaledSupport,scaledSupport
94 : c iloc(2)=(iy*sampling+off(2))* cfscale
95 0 : iloc(2)=iy * scaledSampling + off(2)
96 0 : iv = iloc(2)
97 0 : do ix=-scaledSupport,scaledSupport
98 : c iloc(1)=(ix*sampling+off(1))* cfscale
99 0 : iloc(1)=ix * scaledSampling + off(1)
100 0 : iu = iloc(1)
101 0 : ts=sDPA
102 0 : tc=cDPA
103 0 : if (mreindex(iu,iv,iloc(1),iloc(2),
104 0 : $ ts,tc,convOrigin,convSize)) then
105 0 : if (dopointingcorrection .eq. 1) then
106 0 : griduvw(1)=(iloc(1)-convOrigin)/(scale(1)*sampling)
107 0 : griduvw(2)=(iloc(2)-convOrigin)/(scale(2)*sampling)
108 0 : ra1 = raoff(apol,1,ant1(irow)+1)
109 0 : ra2 = raoff(apol,1,ant2(irow)+1)
110 0 : dec1= decoff(apol,1,ant1(irow)+1)
111 0 : dec2= decoff(apol,1,ant2(irow)+1)
112 : call gfeij(griduvw,area,ra1,dec1,ra2,dec2,
113 0 : $ dograd,pcwt,pdcwt1,pdcwt2,currentCFPA)
114 : else
115 0 : pcwt=cmplx(1.0,0.0)
116 : endif
117 0 : if(uvw(3,irow).gt.0.0) then
118 : cwt=conjg(func(iloc(1),iloc(2), iloc(3),
119 0 : $ ConjPlane))
120 : else
121 : cwt=(func(iloc(1),iloc(2), iloc(3),
122 0 : $ PolnPlane))
123 : end if
124 0 : cwt=(func(iloc(1),iloc(2), iloc(3), PolnPlane))
125 0 : cwt = (cwt/cnorm)
126 :
127 : c tcnorm=tcnorm+cwt*nvalue
128 : grid(origin(1)+ix,origin(2)+iy,apol,achan)=
129 : $ grid(origin(1)+ix,origin(2)+iy,apol,achan)+
130 0 : $ nvalue*cwt*(pcwt)
131 : c$$$ if ((doavgpb .ne. 0)) then
132 : c$$$ write(*,*) cnorm,tcnorm,nvalue,cwt
133 : c$$$ endif
134 : endif
135 : end do
136 : c$$$ write(*,*)
137 : end do
138 : c$$$ if ((doavgpb .eq. 1)) then
139 : c$$$ write(*,*) tcnorm,nvalue
140 : c$$$ endif
141 0 : accumulate = cnorm
142 : c$$$ stop
143 0 : end
144 : C
145 : C Compute area under the function func
146 : C
147 0 : complex function getarea(doavgpb,func,area,rsupport, sampling,
148 : $ off,
149 : $ cfscale,scale,lambda,sDPA, cDPA, currentCFPA, convOrigin,
150 : $ convSize, wconvsize, polused,dopointingcorrection, dograd,
151 0 : $ nant, raoff,decoff, ConjPlane, PolnPlane,ant1,ant2,uvw,
152 : $ nrow,irow)
153 :
154 : integer nrow,irow,doavgpb
155 : integer ant1(nrow),ant2(nrow)
156 : integer rsupport,sampling, off(3), convOrigin, convSize,wconvsize,
157 : $ polused
158 : integer dopointingcorrection, dograd,nant
159 : integer ix,iy,iloc(3),loc(3), ConjPlane, PolnPlane
160 : real raoff(2,1,nant), decoff(2,1,nant)
161 : double precision ra1,ra2,dec1,dec2
162 : double precision offset(3),cfscale,lambda,sDPA,cDPA,currentCFPA,
163 : $ area
164 : double precision griduvw(2), scale(3),uvw(3,nrow)
165 : complex cnorm,cwt,pcwt,pdcwt1,pdcwt2
166 : complex func(convsize, convsize, wconvsize, polused)
167 : integer ixr, iyr
168 : logical opbmos,mreindex
169 : external gcppeij
170 : integer OMP_GET_THREAD_NUM,TID
171 : integer scaledSupport
172 : double precision ts, tc, scaledSampling
173 :
174 0 : cnorm=0
175 0 : iloc(3)=1
176 0 : scaledSampling = int(sampling*cfscale);
177 0 : scaledSupport = int(rsupport/cfscale + 0.5)
178 :
179 0 : do iy=-scaledSupport,scaledSupport
180 : c iloc(2)=(iy*sampling+off(2))*cfscale
181 0 : iloc(2)=iy * scaledSampling + off(2)
182 0 : iv = iloc(2)
183 0 : do ix=-scaledSupport,scaledSupport
184 : c iloc(1)=(ix*sampling+off(1))*cfscale
185 0 : iloc(1)=ix * scaledSampling + off(1)
186 0 : iu=iloc(1)
187 0 : ts = sDPA
188 0 : tc = cDPA
189 0 : if (mreindex(iu,iv,iloc(1),iloc(2), ts,tc,
190 0 : $ convOrigin, convSize)) then
191 0 : if(uvw(3,irow).gt.0.0) then
192 : cwt=conjg(func(iloc(1),iloc(2), iloc(3),
193 0 : $ ConjPlane))
194 : else
195 0 : cwt=(func(iloc(1),iloc(2),iloc(3),PolnPlane))
196 : end if
197 0 : cwt=(func(iloc(1),iloc(2),iloc(3),PolnPlane))
198 0 : cnorm=cnorm+cwt
199 : endif
200 : enddo
201 : enddo
202 : c !$OMP ENDDO
203 : c !$OMP END PARALLEL
204 0 : getarea = cnorm
205 0 : end
206 : C
207 : C Grid a number of visibility records
208 : C
209 0 : subroutine gpbmos (uvw, dphase, values, nvispol, nvischan,
210 0 : $ dopsf, flag, rflag, weight, nrow, rownum,
211 0 : $ scale, offset, grid, nx, ny, npol, nchan, freq, c,
212 0 : $ support, convsize, convwtsize, sampling, wconvsize, convfunc,
213 0 : $ chanmap, polmap,polused,sumwt,
214 0 : $ ant1, ant2, nant, scanno, sigma,raoff, decoff,area,
215 0 : $ dograd,dopointingcorrection,npa,paindex,cfmap,conjcfmap,
216 0 : $ currentCFPA,actualPA,doavgpb,avgpb,cfRefFreq,convWts,
217 : $ wtsupport)
218 :
219 :
220 : implicit none
221 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow,polused
222 : integer npa,paindex
223 : integer convsize, sampling, wconvsize,dograd,dopointingcorrection
224 : complex values(nvispol, nvischan, nrow)
225 : complex grid(nx, ny, npol, nchan),avgpb(nx,ny,npol,nchan)
226 : double precision uvw(3, nrow), freq(nvischan), c, scale(3),
227 : $ offset(3),currentCFPA,actualPA,cfRefFreq
228 : double precision dphase(nrow), uvdist
229 : complex phasor
230 : integer flag(nvispol, nvischan, nrow)
231 : integer rflag(nrow)
232 : real weight(nvischan, nrow), phase
233 : double precision sumwt(npol, nchan)
234 : integer rownum
235 : integer support(wconvsize,polused,npa), rsupport
236 : integer chanmap(nchan), polmap(npol),cfmap(npol),conjcfmap(npol)
237 : integer dopsf, doavgpb,nodoavgpb
238 : integer wtsupport, convWtOrigin, convwtsize
239 : complex nvalue
240 :
241 : complex convfunc(convsize, convsize, wconvsize, polused),
242 : $ cwt,cwt2,convWts(convwtsize, convwtsize, wconvsize, polused)
243 :
244 : real norm
245 : complex cnorm,tcnorm,cnorm2
246 : real wt
247 :
248 : logical opbmos,mreindex,cfOK, cfwtOK
249 : external gcppeij
250 : complex getarea,accumulate
251 :
252 : real pos(3)
253 : integer loc(3), pboff(3), off(3), iloc(3),iu,iv,lloc(3)
254 : double precision pboffset(3)
255 : integer rbeg, rend
256 : integer ix, iy, ipol, ichan, tx,ty
257 : integer apol, achan, irow, PolnPlane, ConjPlane
258 : double precision pi
259 : data pi/3.14159265358979323846/
260 :
261 : double precision griduvw(2),dPA, cDPA,sDPA
262 : double precision mysigma, ra1,ra2,dec1,dec2
263 : double precision sigma,area,lambda,cfscale
264 : complex pcwt, pdcwt1, pdcwt2,tt
265 : integer nant, scanno, ant1(nrow), ant2(nrow)
266 : integer convOrigin
267 : real raoff(nant), decoff(nant)
268 : integer accumPB, chanPolPB
269 : complex tmpvalue
270 : integer OMP_GET_THREAD_NUM,TID
271 0 : logical chansDone(nvischan), polsDone(npol)
272 :
273 0 : chanPolPB=0
274 0 : irow=rownum
275 :
276 0 : if(irow.ge.0) then
277 0 : rbeg=irow+1
278 0 : rend=irow+1
279 : else
280 0 : rbeg=1
281 0 : rend=nrow
282 : end if
283 :
284 0 : dPA = -(currentCFPA - actualPA)
285 : c dPA = 0
286 0 : cDPA = cos(dPA)
287 0 : sDPA = sin(dPA)
288 :
289 0 : do ichan=1,nvischan
290 0 : chansDone(ichan) = .false.
291 : enddo
292 0 : do ichan=1,npol
293 0 : polsDone(ichan) = .false.
294 : enddo
295 :
296 0 : convOrigin = (convsize-1)/2
297 0 : convWtOrigin = (convwtsize-1)/2
298 0 : accumPB=1
299 0 : tcnorm=0.0
300 0 : do irow=rbeg, rend
301 0 : if(rflag(irow).eq.0) then
302 0 : do ichan=1, nvischan
303 0 : achan=chanmap(ichan)+1
304 :
305 0 : lambda = c/freq(ichan)
306 0 : cfscale = cfRefFreq/freq(ichan)
307 : c$$$ uvw(1,irow)=0
308 : c$$$ uvw(2,irow)=0
309 : c$$$ uvw(3,irow)=0
310 0 : if((achan.ge.1).and.(achan.le.nchan).and.
311 0 : $ (weight(ichan,irow).ne.0.0)) then
312 : call spbmos(uvw(1,irow), dphase(irow),
313 : $ freq(ichan), c, scale, offset, sampling,
314 0 : $ pos, loc, off, phasor)
315 0 : iloc(3)=max(1, min(wconvsize, loc(3)))
316 0 : rsupport=support(iloc(3),1,paindex)
317 :
318 : c wtsupport=rsupport
319 :
320 0 : cfOK=opbmos(nx, ny, wconvsize, loc, rsupport)
321 0 : cfwtOK=opbmos(nx, ny, wconvsize, loc, wtsupport)
322 0 : if (cfOK .or. cfwtOK) then
323 0 : PolnPlane=polused+1
324 0 : do ipol=1, nvispol
325 0 : apol=polmap(ipol)+1
326 :
327 : if((flag(ipol,ichan,irow).ne.1).and.
328 0 : $ (apol.ge.1).and.(apol.le.npol)) then
329 0 : ConjPlane = cfmap(ipol)+1
330 0 : PolnPlane = conjcfmap(ipol)+1
331 :
332 0 : if(dopsf.eq.1) then
333 0 : nvalue=cmplx(weight(ichan,irow))
334 : else
335 : nvalue=weight(ichan,irow)*
336 0 : $ (values(ipol,ichan,irow)*phasor)
337 : end if
338 :
339 0 : norm=0.0
340 0 : cnorm=0.0
341 0 : tcnorm=0.0
342 0 : cnorm2=0.0
343 : c accumPB=1
344 0 : if ((doavgpb .gt. 0) .and. (cfwtOK) .and.
345 : $ (accumPB .eq. 1)) then
346 : C
347 : C This is the center of the WT image is. Weight functions are
348 : C accumulated at the origin of the UV-grid.
349 : C
350 0 : lloc(1)=(nx)/2+1
351 0 : lloc(2)=(ny)/2+1
352 0 : lloc(3)=0
353 0 : pboff(1)=0
354 0 : pbOff(2)=0
355 0 : pboff(3)=off(3)
356 0 : convWtOrigin=convwtsize/2+1
357 : c tmpvalue=cmplx(1,0)
358 0 : tmpvalue=cmplx(weight(ichan,irow))
359 : cnorm2 = getarea(doavgpb,convWts,area,
360 : c cnorm2 = getarea(doavgpb,convfunc,area,
361 : $ wtsupport, sampling,pboff,
362 : $ cfscale,scale,lambda, sDPA, cDPA,
363 : $ currentCFPA, convWtOrigin,
364 : $ convwtsize,
365 : $ wconvsize, polused,
366 : $ dopointingcorrection,dograd,
367 : $ nant,raoff,decoff,ConjPlane,
368 0 : $ PolnPlane,ant1,ant2,uvw,nrow,irow)
369 : cnorm2=accumulate(doavgpb,area,avgpb,lloc,
370 : $ convWts,tmpvalue,cnorm2,wtsupport,
371 : c $ convfunc,tmpvalue,cnorm2,wtsupport,
372 : $ sampling,pboff,
373 : $ cfscale,scale,lambda,sDPA,cDPA,
374 : $ currentCFPA,convWtOrigin,convwtsize,
375 : $ wconvsize,polused,
376 : $ dopointingcorrection,dograd,nant,
377 : $ raoff,decoff,ConjPlane,PolnPlane,
378 : $ nx,ny,npol,nchan,apol,achan,
379 0 : $ ant1,ant2,uvw,nrow,irow)
380 0 : chansDone(ichan)=.true.
381 0 : polsDone(PolnPlane) = .true.
382 0 : tcnorm=tcnorm+weight(ichan,irow)
383 : c$$$ write(*,*)'PB: ',irow,tcnorm,ichan,apol
384 : c$$$ if (irow .eq. 1) then
385 : c$$$ write(*,*) 'PB:',ichan,PolnPlane,lloc,
386 : c$$$ $ pboff,cnorm2
387 : c$$$ endif
388 : endif
389 0 : if (cfOK) then
390 0 : convOrigin = convSize/2+1
391 : c$$$ nvalue=cmplx(1,0)
392 : c$$$ loc(1)=(nx)/2+1
393 : c$$$ loc(2)=(ny)/2+1
394 : c$$$ loc(3)=0
395 : c$$$ off(1)=0
396 : c$$$ off(2)=0
397 0 : nodoavgpb=0
398 : cnorm = getarea(nodoavgpb,convfunc,area,
399 : $ rsupport, sampling,
400 : $ off, cfscale,scale,lambda,
401 : $ sDPA, cDPA,currentCFPA,
402 : $ convOrigin, convSize,
403 : $ wconvsize, polused,
404 : $ dopointingcorrection,
405 : $ dograd, nant,raoff,decoff,
406 : $ ConjPlane,PolnPlane,ant1,ant2,
407 0 : $ uvw,nrow,irow)
408 : cnorm=accumulate(nodoavgpb,area,grid,loc,
409 : $ convfunc,nvalue,
410 : $ cnorm,rsupport,
411 : $ sampling,off,cfscale,
412 : $ scale,lambda,
413 : $ sDPA,cDPA,currentCFPA,convOrigin,
414 : $ convSize,wconvsize,polused,
415 : $ dopointingcorrection,dograd,nant,
416 : $ raoff,decoff,ConjPlane,PolnPlane,
417 : $ nx,ny,npol,nchan,apol,achan,
418 0 : $ ant1,ant2,uvw,nrow,irow)
419 : c$$$ write(*,*)'V: ',irow,weight(ichan,irow),
420 : c$$$ $ ichan,apol,nvalue
421 : c$$$ if (irow .eq. 1) then
422 : c$$$ write(*,*)'Vis:',ichan,PolnPlane,loc,
423 : c$$$ $ off,cnorm
424 : c$$$ endif
425 : sumwt(apol,achan)=sumwt(apol,achan)+
426 0 : $ weight(ichan,irow)
427 : c*real(cnorm)
428 : endif
429 : C
430 : end if
431 : C stop
432 : end do
433 : end if
434 : end if
435 : end do
436 0 : chanPolPB=0
437 0 : do ichan=1,npol
438 0 : if (polsDone(ichan) .eqv. .true.) chanPolPB=chanPolPB+1
439 : enddo
440 0 : if (chanPolPB .lt. npol) then
441 0 : accumPB=1
442 : else
443 0 : accumPB=0
444 : endif
445 : end if
446 : end do
447 0 : return
448 : end
449 : C
450 : C Degrid a number of visibility records
451 : C
452 0 : subroutine dpbmos0 (uvw, dphase, values, nvispol, nvischan,
453 0 : $ flag, rflag,
454 0 : $ nrow, rownum, scale, offset, grid, nx, ny, npol, nchan, freq,
455 0 : $ c, support, convsize, sampling, wconvsize, convfunc,
456 0 : $ chanmap, polmap,polused, ant1, ant2, nant,
457 0 : $ scanno, sigma, raoff, decoff,area,dograd,
458 0 : $ dopointingcorrection,npa,paindex,cfmap,conjcfmap,
459 : $ currentCFPA,actualPA,cfRefFreq)
460 :
461 : implicit none
462 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow,polused
463 : integer npa,paindex
464 : integer convsize, wconvsize, sampling
465 : complex values(nvispol, nvischan, nrow)
466 : complex grid(nx, ny, npol, nchan)
467 : double precision uvw(3, nrow), freq(nvischan), c, scale(3),
468 : $ offset(3),currentCFPA,actualPA, dPA, sDPA, cDPA,cfRefFreq
469 : double precision dphase(nrow), uvdist
470 : complex phasor
471 : integer flag(nvispol, nvischan, nrow)
472 : integer rflag(nrow),cfmap(npol),conjcfmap(npol)
473 : integer rownum
474 : integer support(wconvsize,polused,npa), rsupport
475 : integer chanmap(*), polmap(*)
476 :
477 : integer nant, scanno, ant1(nrow), ant2(nrow),dograd,
478 : $ dopointingcorrection
479 : real raoff(nant), decoff(nant)
480 : double precision sigma,area,lambda,cfscale
481 :
482 : complex nvalue
483 :
484 : complex convfunc(convsize, convsize, wconvsize, polused),
485 : $ cwt, pcwt,pdcwt1,pdcwt2
486 : double precision griduvw(2)
487 : double precision mysigma, ra1,ra2,dec1,dec2
488 :
489 : complex norm(4)
490 :
491 : logical opbmos,mreindex
492 : external gcppeij
493 :
494 : real pos(3)
495 : integer loc(3), off(3), iloc(3)
496 : integer rbeg, rend
497 : integer ix, iy, ipol, ichan, PolnPlane, ConjPlane
498 : integer convOrigin
499 : integer apol, achan, irow
500 : real wt, wtx, wty
501 : double precision pi
502 : data pi/3.14159265358979323846/
503 : integer ii,jj,iu,iv
504 :
505 : complex tmp
506 :
507 0 : irow=rownum
508 :
509 0 : if(irow.ge.0) then
510 0 : rbeg=irow+1
511 0 : rend=irow+1
512 : else
513 0 : rbeg=1
514 0 : rend=nrow
515 : end if
516 : C
517 :
518 0 : dPA = -(currentCFPA - actualPA)
519 : c dPA=0
520 0 : cDPA = cos(dPA)
521 0 : sDPA = sin(dPA)
522 0 : convOrigin = (convsize-1)/2
523 :
524 0 : do irow=rbeg, rend
525 0 : if(rflag(irow).eq.0) then
526 0 : do ichan=1, nvischan
527 0 : achan=chanmap(ichan)+1
528 :
529 0 : lambda = c/freq(ichan)
530 0 : cfscale = cfRefFreq/freq(ichan)
531 :
532 0 : if((achan.ge.1).and.(achan.le.nchan)) then
533 : call spbmos(uvw(1,irow), dphase(irow), freq(ichan), c,
534 0 : $ scale, offset, sampling, pos, loc, off, phasor)
535 0 : iloc(3)=max(1, min(wconvsize, loc(3)))
536 0 : rsupport=support(iloc(3),1,paindex)
537 0 : rsupport = nint( (rsupport / cfscale)+0.5 )
538 0 : if (opbmos(nx, ny, wconvsize, loc, rsupport)) then
539 0 : PolnPlane=0
540 :
541 0 : do ipol=1, nvispol
542 0 : apol=polmap(ipol)+1
543 : if((flag(ipol,ichan,irow).ne.1).and.
544 0 : $ (apol.ge.1).and.(apol.le.npol)) then
545 0 : ConjPlane = cfmap(ipol)+1
546 0 : PolnPlane = conjcfmap(ipol)+1
547 :
548 0 : nvalue=0.0
549 0 : norm(apol)=cmplx(0.0,0.0)
550 0 : pcwt=cmplx(1.0,0.0)
551 :
552 0 : do iy=-rsupport,rsupport
553 0 : iloc(2)=1+(iy*sampling+off(2))+convOrigin
554 0 : iloc(2) = iloc(2) * cfscale
555 0 : iv = (iy*sampling+off(2))
556 0 : do ix=-rsupport,rsupport
557 : iloc(1)=1+(ix*sampling+off(1))
558 0 : $ +convOrigin
559 0 : iloc(1) = iloc(1) * cfscale
560 0 : iu = (ix*sampling+off(1))
561 :
562 0 : if(mreindex(iu,iv,iloc(1),iloc(2),
563 : $ sDPA,cDPA, convOrigin, convSize))
564 0 : $ then
565 0 : if (dopointingcorrection .eq. 1)then
566 : griduvw(1)=(loc(1)-offset(1)
567 : $ +ix-1)/scale(1)-uvw(1,irow)
568 0 : $ /lambda
569 : griduvw(2)=(loc(2)-offset(2)
570 : $ +iy-1)/scale(2)-uvw(2,irow)
571 0 : $ /lambda
572 0 : ra1 = raoff(ant1(irow)+1)
573 0 : ra2 = raoff(ant2(irow)+1)
574 0 : dec1= decoff(ant1(irow)+1)
575 0 : dec2= decoff(ant2(irow)+1)
576 : call gfeij(griduvw,area,
577 : $ ra1,dec1,ra2,dec2,
578 : $ dograd,pcwt,pdcwt1,pdcwt2,
579 0 : $ currentCFPA)
580 : c$$$ call gcppeij(griduvw,area,
581 : c$$$ $ ra1,dec1,ra2,dec2,
582 : c$$$ $ dograd,pcwt,pdcwt1,pdcwt2,
583 : c$$$ $ currentCFPA)
584 : endif
585 :
586 0 : if(uvw(3,irow).gt.0.0) then
587 : cwt=conjg(convfunc(iloc(1),
588 0 : $ iloc(2), iloc(3),ConjPlane))
589 0 : pcwt = (pcwt)
590 : else
591 : cwt=(convfunc(iloc(1),
592 0 : $ iloc(2), iloc(3),PolnPlane))
593 0 : pcwt = (pcwt)
594 : endif
595 : nvalue=nvalue+(cwt)*(pcwt)*
596 : $ grid(loc(1)+ix,loc(2)+iy,apol,
597 0 : $ achan)
598 0 : norm(apol)=norm(apol)+(pcwt*cwt)
599 : endif
600 : end do
601 : end do
602 : values(ipol,ichan,irow)=
603 : $ nvalue*conjg(phasor)
604 0 : $ /norm(apol)
605 : end if
606 : end do
607 : end if
608 : end if
609 : end do
610 : end if
611 : end do
612 0 : return
613 : end
614 : C
615 : C Degrid a number of visibility records along with the grad. computations
616 : C
617 0 : subroutine dpbmosgrad (uvw, dphase, values, nvispol, nvischan,
618 0 : $ gazvalues, gelvalues, doconj,
619 0 : $ flag, rflag, nrow, rownum, scale, offset, grid,
620 0 : $ nx, ny, npol, nchan, freq, c, support, convsize, sampling,
621 0 : $ wconvsize, convfunc, chanmap, polmap,polused,ant1,ant2,nant,
622 0 : $ scanno, sigma, raoff, decoff,area,dograd,
623 0 : $ dopointingcorrection,npa,paindex,cfmap,conjcfmap,
624 : $ currentCFPA,actualPA,cfRefFreq)
625 :
626 : implicit none
627 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow,polused
628 : integer npa,paindex,doconj
629 : integer convsize, wconvsize, sampling
630 : complex values(nvispol, nvischan, nrow)
631 : complex gazvalues(nvispol, nvischan, nrow)
632 : complex gelvalues(nvispol, nvischan, nrow)
633 : complex grid(nx, ny, npol, nchan)
634 : double precision uvw(3, nrow), freq(nvischan), c, scale(3),
635 : $ offset(3),currentCFPA,actualPA, dPA, sDPA, cDPA,cfRefFreq
636 : double precision dphase(nrow), uvdist
637 : complex phasor
638 : integer flag(nvispol, nvischan, nrow)
639 : integer rflag(nrow),cfmap(npol),conjcfmap(npol)
640 : integer rownum
641 : integer support(wconvsize,polused,npa), rsupport
642 : integer chanmap(*), polmap(*)
643 :
644 : integer nant, scanno, ant1(nrow), ant2(nrow),dograd,
645 : $ dopointingcorrection
646 : real raoff(2,1,nant), decoff(2,1,nant)
647 : double precision sigma,area,lambda,cfscale
648 :
649 : complex nvalue,ngazvalue,ngelvalue
650 :
651 : complex convfunc(convsize, convsize, wconvsize, polused),
652 : $ cwt, pcwt,pdcwt1,pdcwt2
653 : double precision griduvw(2)
654 : double precision mysigma, ra1,ra2,dec1,dec2
655 :
656 : complex norm(4),gradaznorm(4),gradelnorm(4)
657 :
658 : logical opbmos,mreindex
659 : external gcppeij
660 :
661 : real pos(3)
662 : integer loc(3), off(3), iloc(3)
663 : integer rbeg, rend
664 : integer ix, iy, ipol, ichan, PolnPlane, ConjPlane
665 : integer convOrigin
666 : integer apol, achan, irow
667 : real wt, wtx, wty
668 : double precision pi, scaledSampling
669 : data pi/3.14159265358979323846/
670 : integer ii,jj,iu,iv
671 : integer scaledSupport
672 : complex tmp
673 :
674 0 : irow=rownum
675 :
676 0 : if(irow.ge.0) then
677 0 : rbeg=irow+1
678 0 : rend=irow+1
679 : else
680 0 : rbeg=1
681 0 : rend=nrow
682 : end if
683 : C
684 0 : dPA = -(currentCFPA - actualPA)
685 : c dPA = 0
686 0 : cDPA = cos(dPA)
687 0 : sDPA = sin(dPA)
688 : c convOrigin = (convsize-1)/2
689 0 : convOrigin = (convsize)/2
690 :
691 0 : do irow=rbeg, rend
692 0 : if(rflag(irow).eq.0) then
693 0 : do ichan=1, nvischan
694 0 : achan=chanmap(ichan)+1
695 :
696 0 : lambda = c/freq(ichan)
697 0 : cfscale = cfRefFreq/freq(ichan)
698 0 : scaledSampling = int(sampling*cfscale)
699 :
700 0 : if((achan.ge.1).and.(achan.le.nchan)) then
701 : call spbmos(uvw(1,irow), dphase(irow), freq(ichan), c,
702 0 : $ scale, offset, sampling, pos, loc, off, phasor)
703 0 : iloc(3)=max(1, min(wconvsize, loc(3)))
704 0 : rsupport=support(iloc(3),1,paindex)
705 0 : scaledSupport = int(rsupport/cfscale+0.5)
706 : c rsupport = nint( (rsupport / cfscale)+0.5 )
707 0 : if (opbmos(nx, ny, wconvsize, loc, rsupport)) then
708 0 : do ipol=1, nvispol
709 0 : apol=polmap(ipol)+1
710 : if((flag(ipol,ichan,irow).ne.1).and.
711 0 : $ (apol.ge.1).and.(apol.le.npol)) then
712 0 : ConjPlane = conjcfmap(ipol)+1
713 0 : PolnPlane = cfmap(ipol)+1
714 0 : ConjPlane = cfmap(ipol)+1
715 0 : PolnPlane = conjcfmap(ipol)+1
716 :
717 0 : nvalue=0.0
718 0 : ngazvalue = 0.0
719 0 : ngelvalue = 0.0
720 :
721 0 : norm(apol)=cmplx(0.0,0.0)
722 0 : gradaznorm(apol)=cmplx(0.0,0.0)
723 0 : gradelnorm(apol)=cmplx(0.0,0.0)
724 0 : pcwt=cmplx(1.0,0.0)
725 :
726 0 : do iy=-scaledSupport,scaledSupport
727 : c iloc(2) = (iy*sampling+off(2)) * cfscale
728 : iloc(2) = (1+(iy*sampling+off(2)))
729 0 : $ *cfscale + convOrigin
730 0 : iv=iy*scaledSampling+off(2)
731 0 : do ix=-scaledSupport,scaledSupport
732 : c iloc(1) = (ix*sampling+off(1))*cfscale
733 : iloc(1) = (1+(ix*sampling+off(1)))
734 0 : $ *cfscale + convOrigin
735 0 : iu=ix*scaledSampling+off(1)
736 :
737 0 : if(mreindex(iu,iv,iloc(1),iloc(2),
738 : $ sDPA,cDPA, convOrigin, convSize))
739 0 : $ then
740 0 : if (dopointingcorrection .eq. 1)then
741 : c$$$ griduvw(1)=(loc(1)-offset(1)
742 : c$$$ $ +ix-1)/scale(1)-uvw(1,irow)
743 : c$$$ $ /lambda
744 : c$$$ griduvw(2)=(loc(2)-offset(2)
745 : c$$$ $ +iy-1)/scale(2)-uvw(2,irow)
746 : c$$$ $ /lambda
747 0 : iu=iloc(1)-convOrigin
748 0 : iv=iloc(2)-convOrigin
749 0 : griduvw(1)=iu/(scale(1)*sampling)
750 0 : griduvw(2)=iv/(scale(2)*sampling)
751 0 : ii=apol
752 0 : ra1 = raoff(ii,1,ant1(irow)+1)
753 : c $ *cfscale
754 0 : ra2 = raoff(ii,1,ant2(irow)+1)
755 : c $ *cfscale
756 0 : dec1= decoff(ii,1,ant1(irow)+1)
757 : c $ *cfscale
758 0 : dec2= decoff(ii,1,ant2(irow)+1)
759 : c $ *cfscale
760 : call gfeij(griduvw,area,
761 : $ ra1,dec1,ra2,dec2,
762 : $ dograd,pcwt,pdcwt1,pdcwt2,
763 0 : $ currentCFPA)
764 : c$$$ call gcppeij(griduvw,area,
765 : c$$$ $ ra1,dec1,ra2,dec2,
766 : c$$$ $ dograd,pcwt,pdcwt1,pdcwt2,
767 : c$$$ $ currentCFPA)
768 : endif
769 : c$$$ if(uvw(3,irow).gt.0.0) then
770 : c$$$ cwt=conjg(convfunc(iloc(1),
771 : c$$$ $ iloc(2), iloc(3),ConjPlane))
772 : c$$$ pcwt = conjg(pcwt)
773 : c$$$ pdcwt1 = conjg(pdcwt1)
774 : c$$$ pdcwt2 = conjg(pdcwt2)
775 : c$$$ else
776 : c$$$ cwt=(convfunc(iloc(1),
777 : c$$$ $ iloc(2), iloc(3),PolnPlane))
778 : c$$$ endif
779 : cwt=(convfunc(iloc(1),
780 0 : $ iloc(2), iloc(3),PolnPlane))
781 0 : cwt = cwt*pcwt
782 : nvalue=nvalue+(cwt)*
783 : $ grid(loc(1)+ix,loc(2)+iy,apol,
784 0 : $ achan)
785 :
786 0 : if ((doconj .eq. 1) .and.
787 : $ (dograd .eq. 1)) then
788 0 : pdcwt1 = conjg(pdcwt1)
789 0 : pdcwt2 = conjg(pdcwt2)
790 : endif
791 0 : if (dograd .eq. 1) then
792 0 : pdcwt1 = pdcwt1*cwt
793 0 : pdcwt2 = pdcwt2*cwt
794 : ngazvalue=ngazvalue+(pdcwt1)*
795 : $ (grid(loc(1)+ix,loc(2)+iy,
796 0 : $ apol,achan))
797 : ngelvalue=ngelvalue+(pdcwt2)*
798 : $ (grid(loc(1)+ix,loc(2)+iy,
799 0 : $ apol,achan))
800 : endif
801 :
802 0 : norm(apol)=norm(apol)+(cwt)
803 : gradaznorm(apol)=gradaznorm(apol)+
804 0 : $ pdcwt1
805 : gradelnorm(apol)=gradelnorm(apol)+
806 0 : $ pdcwt1
807 : endif
808 : end do
809 : end do
810 : values(ipol,ichan,irow)=
811 : $ nvalue*conjg(phasor)
812 0 : $ /norm(apol)
813 0 : if (dograd .eq. 1) then
814 : gazvalues(ipol,ichan,irow)=ngazvalue*
815 0 : $ conjg(phasor)/norm(apol)
816 : gelvalues(ipol,ichan,irow)=ngelvalue*
817 0 : $ conjg(phasor)/norm(apol)
818 : endif
819 : end if
820 : end do
821 : end if
822 : end if
823 : end do
824 : end if
825 : end do
826 0 : return
827 : end
828 : C
829 : C Degrid a number of visibility records
830 : C
831 0 : subroutine dpbmos (uvw, dphase, values, nvispol, nvischan,
832 0 : $ flag, rflag,
833 0 : $ nrow, rownum, scale, offset, grid, nx, ny, npol, nchan, freq,
834 0 : $ c, support, convsize, sampling, wconvsize, convfunc,
835 0 : $ chanmap, polmap,polused, ant1, ant2, nant,
836 0 : $ scanno, sigma, raoff, decoff,area,dograd,
837 0 : $ dopointingcorrection,npa,paindex,cfmap,conjcfmap,
838 : $ currentCFPA,actualPA,cfRefFreq)
839 :
840 : implicit none
841 : integer nx, ny, npol, nchan, nvispol, nvischan, nrow,polused
842 : integer npa,paindex
843 : integer convsize, wconvsize, sampling
844 : complex values(nvispol, nvischan, nrow)
845 : complex grid(nx, ny, npol, nchan)
846 : double precision uvw(3, nrow), freq(nvischan), c, scale(3),
847 : $ offset(3),currentCFPA,actualPA, dPA, sDPA, cDPA,cfRefFreq
848 : double precision dphase(nrow), uvdist
849 : complex phasor
850 : integer flag(nvispol, nvischan, nrow)
851 : integer rflag(nrow),cfmap(npol),conjcfmap(npol)
852 : integer rownum
853 : integer support(wconvsize,polused,npa), rsupport
854 : integer chanmap(*), polmap(*)
855 :
856 : integer nant, scanno, ant1(nrow), ant2(nrow),dograd,
857 : $ dopointingcorrection
858 : real raoff(2,1,nant), decoff(2,1,nant)
859 : double precision sigma,area,lambda,cfscale
860 :
861 : complex nvalue, junk
862 :
863 : C complex convfunc(convsize/2-1, convsize/2-1, wconvsize, polused),
864 : complex convfunc(convsize, convsize, wconvsize, polused),
865 : $ cwt, pcwt,pdcwt1,pdcwt2
866 : double precision griduvw(2)
867 : double precision mysigma, ra1,ra2,dec1,dec2
868 :
869 : complex norm(4)
870 :
871 : logical opbmos,mreindex
872 : c external nwcppEij
873 : external gcppeij
874 :
875 : real pos(3)
876 : integer loc(3), off(3), iloc(3)
877 : integer rbeg, rend
878 : integer ix, iy, ipol, ichan, PolnPlane, ConjPlane
879 : integer convOrigin
880 : integer apol, achan, irow
881 : real wt, wtx, wty
882 : double precision pi, scaledSampling
883 : data pi/3.14159265358979323846/
884 : integer ii,jj,iu,iv
885 : integer scaledSupport
886 : complex tmp
887 :
888 0 : irow=rownum
889 :
890 0 : if(irow.ge.0) then
891 0 : rbeg=irow+1
892 0 : rend=irow+1
893 : else
894 0 : rbeg=1
895 0 : rend=nrow
896 : end if
897 : C
898 :
899 0 : dPA = -(currentCFPA - actualPA)
900 : c dPA = 0
901 0 : cDPA = cos(dPA)
902 0 : sDPA = sin(dPA)
903 0 : convOrigin = (convsize-1)/2
904 : c convOrigin = (convsize)/2
905 :
906 0 : do irow=rbeg, rend
907 0 : if(rflag(irow).eq.0) then
908 0 : do ichan=1, nvischan
909 0 : achan=chanmap(ichan)+1
910 :
911 0 : lambda = c/freq(ichan)
912 0 : cfscale = cfRefFreq/freq(ichan)
913 0 : scaledSampling = int(sampling*cfscale)
914 :
915 0 : if((achan.ge.1).and.(achan.le.nchan)) then
916 : call spbmos(uvw(1,irow), dphase(irow), freq(ichan), c,
917 0 : $ scale, offset, sampling, pos, loc, off, phasor)
918 0 : iloc(3)=max(1, min(wconvsize, loc(3)))
919 0 : rsupport=support(iloc(3),1,paindex)
920 0 : scaledSupport = int(rsupport/cfscale+0.5)
921 : c rsupport = nint( (rsupport / cfscale)+0.5 )
922 :
923 0 : if (opbmos(nx, ny, wconvsize, loc, rsupport)) then
924 0 : PolnPlane=0
925 0 : do ipol=1, nvispol
926 0 : apol=polmap(ipol)+1
927 : if((flag(ipol,ichan,irow).ne.1).and.
928 0 : $ (apol.ge.1).and.(apol.le.npol)) then
929 :
930 : c$$$ ConjPlane = conjcfmap(ipol)+1
931 : c$$$ PolnPlane = cfmap(ipol)+1
932 : C The following after feed_x -> -feed_x and PA -> PA + PI/2
933 0 : ConjPlane = cfmap(ipol)+1
934 0 : PolnPlane = conjcfmap(ipol)+1
935 :
936 0 : nvalue=0.0
937 0 : norm(apol)=cmplx(0.0,0.0)
938 0 : pcwt=cmplx(1.0,0.0)
939 :
940 0 : do iy=-scaledSupport,scaledSupport
941 : c iloc(2)=1+(iy*sampling+off(2))+convOrigin
942 0 : iv = iy*scaledSampling+off(2)
943 0 : do ix=-scaledSupport,scaledSupport
944 : c iloc(1)=1+(ix*sampling+off(1))+convOrigin
945 0 : iu = ix*scaledSampling+off(1)
946 :
947 0 : if(mreindex(iu,iv,iloc(1),iloc(2),
948 : $ sDPA,cDPA, convOrigin, convSize))
949 0 : $ then
950 0 : if (dopointingcorrection .eq. 1)
951 : $ then
952 : c$$$ griduvw(1) = (loc(1)-offset(1)
953 : c$$$ $ +ix-1)
954 : c$$$ $ /scale(1)-uvw(1,irow)/lambda
955 : c$$$ griduvw(2) = (loc(2)-offset(2)
956 : c$$$ $ +iy-1)
957 : c$$$ $ /scale(2)-uvw(2,irow)/lambda
958 : C
959 : C (iu,iv) are the non-PA-rotated co-ords. for the Conv. Func. Use
960 : C those to compute the phase gradient for pointing offset(squint).
961 : C Pointing errors must rotate with PA and for those, iloc must be
962 : C used to compute the appropriate phase gradient.
963 : C
964 : griduvw(1)=(iloc(1)-convOrigin)/
965 0 : $ (scale(1)*sampling)
966 : griduvw(2)=(iloc(2)-convOrigin)/
967 0 : $ (scale(2)*sampling)
968 0 : ii=apol
969 0 : ra1 = raoff(ii,1,ant1(irow)+1)
970 0 : ra2 = raoff(ii,1,ant2(irow)+1)
971 0 : dec1= decoff(ii,1,ant1(irow)+1)
972 0 : dec2= decoff(ii,1,ant2(irow)+1)
973 : call gfeij(griduvw,area,
974 : $ ra1,dec1,ra2,dec2,
975 : $ dograd,pcwt,pdcwt1,pdcwt2,
976 0 : $ currentCFPA)
977 : c$$$ call gcppeij(griduvw,area,
978 : c$$$ $ ra1,dec1,ra2,dec2,
979 : c$$$ $ dograd,pcwt,pdcwt1,pdcwt2,
980 : c$$$ $ currentCFPA)
981 : endif
982 :
983 0 : if(uvw(3,irow).gt.0.0) then
984 : cwt=(convfunc(iloc(1),
985 0 : $ iloc(2), iloc(3),PolnPlane))
986 : else
987 : cwt=conjg(convfunc(iloc(1),
988 0 : $ iloc(2), iloc(3),ConjPlane))
989 : endif
990 : C New code start
991 : cwt=(convfunc(iloc(1),
992 0 : $ iloc(2), iloc(3),PolnPlane))
993 0 : cwt=cwt*pcwt;
994 : nvalue=nvalue+(cwt)
995 : $ *grid(loc(1)+ix,loc(2)+iy,
996 0 : $ apol,achan)
997 : C New code end
998 :
999 : c norm(apol)=norm(apol)+(pcwt*cwt)
1000 0 : norm(apol)=norm(apol)+cwt
1001 :
1002 : c$$$ junk=grid(loc(1)+ix,loc(2)+iy,
1003 : c$$$ $ apol,achan)
1004 : c$$$ write(*,*)nvalue,
1005 : c$$$ $ abs(junk),
1006 : c$$$ $ atan2(imag(junk),
1007 : c$$$ $ real(junk))*57.29563,
1008 : c$$$ $ ix,iy,apol,loc(1),loc(2),
1009 : c$$$ $ iloc(1), iloc(2),
1010 : c$$$ $ (cwt),iu,iv,
1011 : c$$$ $ norm(apol)
1012 :
1013 : c$$$ write(*,*)abs(grid(loc(1)+ix,loc(2)+iy
1014 : c$$$ $ ,apol,achan)),ix,iy,apol,abs(cwt),
1015 : c$$$ $ (pcwt),dopointingcorrection
1016 : endif
1017 : end do
1018 : c write(*,*)
1019 : end do
1020 : c stop
1021 : c norm(apol) = norm(apol)/abs(norm(apol))
1022 : values(ipol,ichan,irow)=
1023 : $ nvalue*conjg(phasor)
1024 0 : $ /norm(apol)
1025 : c$$$ if ((ant1(irow)+1 .eq. 2) .and.
1026 : c$$$ $ (ant2(irow)+1 .eq. 4)) then
1027 : c$$$ write (*,*) ant1(irow)+1, ant2(irow)+1,
1028 : c$$$ $ ipol,ichan,irow,
1029 : c$$$ $ abs(values(ipol,ichan,irow)),
1030 : c$$$ $ atan2(imag(values(ipol,ichan,irow)),
1031 : c$$$ $ real(values(ipol,ichan,irow))),
1032 : c$$$ $ abs(nvalue),
1033 : c$$$ $ atan2(imag(nvalue),real(nvalue)),
1034 : c$$$ $ abs(norm(apol)),
1035 : c$$$ $ atan2(imag(norm(apol)),
1036 : c$$$ $ real(norm(apol)))
1037 : c$$$ endif
1038 : c stop
1039 : end if
1040 : end do
1041 : end if
1042 : end if
1043 : end do
1044 : end if
1045 : end do
1046 0 : return
1047 : end
1048 : C
1049 : C Calculate gridded coordinates and the phasor needed for
1050 : C phase rotation.
1051 : C
1052 0 : subroutine spbmos (uvw, dphase, freq, c, scale, offset,
1053 : $ sampling, pos, loc, off, phasor)
1054 : implicit none
1055 : integer loc(3), off(3), sampling
1056 : double precision uvw(3), freq, c, scale(3), offset(3)
1057 : real pos(3)
1058 : double precision dphase, phase
1059 : complex phasor
1060 : integer idim
1061 : double precision pi
1062 : data pi/3.14159265358979323846/
1063 :
1064 0 : pos(3)=sqrt(abs(scale(3)*uvw(3)*freq/c))+offset(3)+1.0
1065 0 : loc(3)=nint(pos(3))
1066 0 : off(3)=0
1067 :
1068 0 : do idim=1,2
1069 : pos(idim)=scale(idim)*uvw(idim)*freq/c+
1070 0 : $ (offset(idim)+1.0)
1071 0 : loc(idim)=nint(pos(idim))
1072 0 : off(idim)=nint((loc(idim)-pos(idim))*sampling)
1073 : end do
1074 :
1075 0 : phase=-2.0D0*pi*dphase*freq/c
1076 0 : phasor=cmplx(cos(phase), sin(phase))
1077 0 : return
1078 : end
1079 0 : logical function opbmos (nx, ny, nw, loc, support)
1080 : implicit none
1081 : integer nx, ny, nw, loc(3), support
1082 : opbmos=(support.gt.0).and.
1083 : $ (loc(1)-support.ge.1).and.(loc(1)+support.le.nx).and.
1084 : $ (loc(2)-support.ge.1).and.(loc(2)+support.le.ny).and.
1085 0 : $ (loc(3).ge.1).and.(loc(3).le.nw)
1086 0 : return
1087 : end
1088 :
1089 0 : logical function mreindex(inx,iny,outx,outy,sinDPA,cosDPA,
1090 : $ Origin, Size)
1091 : integer inx,iny,outx,outy, Origin, Size
1092 : double precision sinDPA, cosDPA
1093 : integer ix,iy
1094 :
1095 0 : if (sinDPA .eq. 0) then
1096 0 : ix = inx
1097 0 : iy = iny
1098 : else
1099 0 : ix = nint(cosDPA*inx + sinDPA*iny)
1100 0 : iy = nint(-sinDPA*inx + cosDPA*iny)
1101 : endif
1102 :
1103 0 : outx = ix+Origin
1104 0 : outy = iy+Origin
1105 :
1106 : mreindex=(outx .ge. 1 .and.
1107 : $ outx .le. Size .and.
1108 : $ outy .ge. 1 .and.
1109 0 : $ outy .le. Size)
1110 :
1111 0 : return
1112 : end
1113 :
|
