LCOV - code coverage report
Current view: top level - synthesis/fortran - fgridft.f (source / functions) Hit Total Coverage
Test: casacpp_coverage.info Lines: 181 293 61.8 %
Date: 2024-11-06 17:42:47 Functions: 7 11 63.6 %

          Line data    Source code
       1             : *=======================================================================
       2             : *     Copyright (C) 1999-2012
       3             : *     Associated Universities, Inc. Washington DC, USA.
       4             : *
       5             : *     This library is free software; you can redistribute it and/or
       6             : *     modify it under the terms of the GNU Library General Public
       7             : *     License as published by the Free Software Foundation; either
       8             : *     version 2 of the License, or (at your option) any later version.
       9             : *
      10             : *     This library is distributed in the hope that it will be useful,
      11             : *     but WITHOUT ANY WARRANTY; without even the implied warranty of
      12             : *     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      13             : *     GNU Library General Public License for more details.
      14             : *
      15             : *     You should have received a copy of the GNU Library General Public
      16             : *     License along with this library; if not, write to the Free
      17             : *     Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,
      18             : *     MA 02139, USA.
      19             : *
      20             : *     Correspondence concerning AIPS++ should be addressed as follows:
      21             : *            Internet email: casa-feedback@nrao.edu.
      22             : *            Postal address: AIPS++ Project Office
      23             : *                            National Radio Astronomy Observatory
      24             : *                            520 Edgemont Road
      25             : *                            Charlottesville, VA 22903-2475 USA
      26             : *
      27             : *     $Id: fgridft.f 19685 2006-10-05 20:57:29Z rurvashi $
      28             : *-----------------------------------------------------------------------
      29             : C
      30             : C Grid a number of visibility records
      31             : C
      32           0 :       subroutine ggrid (uvw, dphase, values, nvispol, nvischan,
      33           0 :      $   dopsf, flag, rflag, weight, nrow, rownum,
      34           0 :      $   scale, offset, grid, nx, ny, npol, nchan, freq, c,
      35           0 :      $   support, sampling, convFunc, chanmap, polmap, sumwt)
      36             : 
      37             :       implicit none
      38             :       integer nx, ny, npol, nchan, nvispol, nvischan, nrow
      39             :       complex values(nvispol, nvischan, nrow)
      40             :       double complex grid(nx, ny, npol, nchan)
      41             :       double precision uvw(3, nrow), freq(nvischan), c, scale(2),
      42             :      $     offset(2)
      43             :       double precision dphase(nrow), uvdist
      44             :       complex phasor
      45             :       integer flag(nvispol, nvischan, nrow)
      46             :       integer rflag(nrow)
      47             :       real weight(nvischan, nrow)
      48             :       double precision sumwt(npol, nchan)
      49             :       integer rownum
      50             :       integer support, sampling
      51             :       integer chanmap(nchan), polmap(npol)
      52             :       integer dopsf
      53             : 
      54             :       double complex nvalue
      55             : 
      56             :       double precision convFunc(*)
      57             :       double precision norm
      58             :       double precision wt, wtx, wty
      59             : 
      60             :       logical ogrid
      61             : 
      62             :       double precision pos(2)
      63             :       integer loc(2), off(2), iloc(2)
      64             :       integer rbeg, rend
      65             :       integer ix, iy, ipol, ichan
      66             :       integer apol, achan, irow
      67             :       
      68           0 :       irow=rownum
      69             : 
      70           0 :       if(irow.ge.0) then
      71           0 :          rbeg=irow+1
      72           0 :          rend=irow+1
      73             :       else 
      74           0 :          rbeg=1
      75           0 :          rend=nrow
      76             :       end if
      77             : 
      78           0 :       do irow=rbeg, rend
      79           0 :          if(rflag(irow).eq.0) then 
      80           0 :          do ichan=1, nvischan
      81           0 :             achan=chanmap(ichan)+1
      82           0 :             if((achan.ge.1).and.(achan.le.nchan).and.
      83           0 :      $           (weight(ichan,irow).ne.0.0)) then
      84             :                call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c, 
      85           0 :      $              scale, offset, sampling, pos, loc, off, phasor)
      86           0 :                if (ogrid(nx, ny, loc, support)) then
      87           0 :                   do ipol=1, nvispol
      88           0 :                      apol=polmap(ipol)+1
      89             :                      if((flag(ipol,ichan,irow).ne.1).and.
      90           0 :      $                    (apol.ge.1).and.(apol.le.npol)) then
      91             : C If we are making a PSF then we don't want to phase
      92             : C rotate but we do want to reproject uvw
      93           0 :                         if(dopsf.eq.1) then
      94           0 :                            nvalue=cmplx(weight(ichan,irow))
      95             :                         else
      96             :                            nvalue=weight(ichan,irow)*
      97           0 :      $                        (values(ipol,ichan,irow)*phasor)
      98             :                         end if
      99           0 :                         norm=0.0
     100           0 :                         do iy=-support,support
     101           0 :                            iloc(2)=abs(sampling*iy+off(2))+1
     102           0 :                            wty=convFunc(iloc(2))
     103           0 :                            do ix=-support,support
     104           0 :                               iloc(1)=abs(sampling*ix+off(1))+1
     105           0 :                               wtx=convFunc(iloc(1))
     106           0 :                               wt=wtx*wty
     107             :                               grid(loc(1)+ix,loc(2)+iy,apol,achan)=
     108             :      $                             grid(loc(1)+ix,loc(2)+iy,apol,achan)+
     109           0 :      $                             nvalue*wt
     110           0 :                               norm=norm+wt
     111             :                            end do
     112             :                         end do
     113             :                         sumwt(apol,achan)=sumwt(apol,achan)+
     114           0 :      $                       weight(ichan,irow)*norm
     115             :                      end if
     116             :                   end do
     117             :                end if
     118             :             end if
     119             :          end do
     120             :          end if
     121             :       end do
     122           0 :       return
     123             :       end
     124             : C
     125             : C     Single precision gridder
     126           0 :       subroutine ggrids (uvw, dphase, values, nvispol, nvischan,
     127           0 :      $     dopsf, flag, rflag, weight, nrow, rownum,
     128           0 :      $     scale, offset, grid, nx, ny, npol, nchan, freq, c,
     129           0 :      $     support, sampling, convFunc, chanmap, polmap, sumwt)
     130             : 
     131             :       implicit none
     132             :       integer nx, ny, npol, nchan, nvispol, nvischan, nrow
     133             :       complex values(nvispol, nvischan, nrow)
     134             :       complex grid(nx, ny, npol, nchan)
     135             :       double precision uvw(3, nrow), freq(nvischan), c, scale(2),
     136             :      $     offset(2)
     137             :       double precision dphase(nrow), uvdist
     138             :       complex phasor
     139             :       integer flag(nvispol, nvischan, nrow)
     140             :       integer rflag(nrow)
     141             :       real weight(nvischan, nrow)
     142             :       double precision sumwt(npol, nchan)
     143             :       integer rownum
     144             :       integer support, sampling
     145             :       integer chanmap(nchan), polmap(npol)      
     146             :       integer dopsf
     147             : 
     148             :       double complex nvalue
     149             : 
     150             :       double precision convFunc(*)
     151             :       double precision norm
     152             :       double precision wt, wtx, wty
     153             : 
     154             :       logical ogrid
     155             : 
     156             :       double precision pos(2)
     157             :       integer loc(2), off(2), iloc(2)
     158             :       integer rbeg, rend
     159             :       integer ix, iy, ipol, ichan
     160             :       integer apol, achan, irow
     161             :       
     162           0 :       irow=rownum
     163             : 
     164           0 :       if(irow.ge.0) then
     165           0 :          rbeg=irow+1
     166           0 :          rend=irow+1
     167             :       else 
     168           0 :          rbeg=1
     169           0 :          rend=nrow
     170             :       end if
     171             : 
     172           0 :       do irow=rbeg, rend
     173           0 :          if(rflag(irow).eq.0) then 
     174           0 :          do ichan=1, nvischan
     175           0 :             achan=chanmap(ichan)+1
     176           0 :             if((achan.ge.1).and.(achan.le.nchan).and.
     177           0 :      $           (weight(ichan,irow).ne.0.0)) then
     178             :                call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c, 
     179           0 :      $              scale, offset, sampling, pos, loc, off, phasor)
     180           0 :                if (ogrid(nx, ny, loc, support)) then
     181           0 :                   do ipol=1, nvispol
     182           0 :                      apol=polmap(ipol)+1
     183             :                      if((flag(ipol,ichan,irow).ne.1).and.
     184           0 :      $                    (apol.ge.1).and.(apol.le.npol)) then
     185             : C If we are making a PSF then we don't want to phase
     186             : C rotate but we do want to reproject uvw
     187           0 :                         if(dopsf.eq.1) then
     188           0 :                            nvalue=cmplx(weight(ichan,irow))
     189             :                         else
     190             :                            nvalue=weight(ichan,irow)*
     191           0 :      $                        (values(ipol,ichan,irow)*phasor)
     192             :                         end if
     193           0 :                         norm=0.0
     194           0 :                         do iy=-support,support
     195           0 :                            iloc(2)=abs(sampling*iy+off(2))+1
     196           0 :                            wty=convFunc(iloc(2))
     197           0 :                            do ix=-support,support
     198           0 :                               iloc(1)=abs(sampling*ix+off(1))+1
     199           0 :                               wtx=convFunc(iloc(1))
     200           0 :                               wt=wtx*wty
     201             :                               grid(loc(1)+ix,loc(2)+iy,apol,achan)=
     202             :      $                             grid(loc(1)+ix,loc(2)+iy,apol,achan)+
     203           0 :      $                             nvalue*wt
     204           0 :                               norm=norm+wt
     205             :                            end do
     206             :                         end do
     207             :                         sumwt(apol,achan)=sumwt(apol,achan)+
     208           0 :      $                       weight(ichan,irow)*norm
     209             :                      end if
     210             :                   end do
     211             :                end if
     212             :             end if
     213             :          end do
     214             :          end if
     215             :       end do
     216           0 :       return
     217             :       end
     218             : 
     219             : 
     220    16617992 :       subroutine sectggridd(values, nvispol, nvischan,
     221    16617992 :      $   dopsf, flag, rflag, weight, nrow, 
     222    16617992 :      $   grid, nx, ny, npol, nchan, 
     223    16617992 :      $   support, sampling, convFunc, chanmap, polmap, sumwt, x0,
     224    16617992 :      $    y0, nxsub, nysub, rbeg, rend, loc, off, phasor)
     225             :       implicit none
     226             : 
     227             :       integer, intent(in) :: nx, ny, npol, nchan, nvispol, nvischan,nrow
     228             :       complex, intent(in) :: values(nvispol, nvischan, nrow)
     229             :       double complex, intent(inout) :: grid(nx, ny, npol, nchan)
     230             :       integer, intent(in) :: x0, y0, nxsub, nysub
     231             :       double precision, intent(in) :: convFunc(*)
     232             :       integer, intent(in) :: chanmap(nvischan), polmap(nvispol)
     233             :       integer, intent(in) ::  flag(nvispol, nvischan, nrow)
     234             :       integer, intent(in) ::  rflag(nrow)
     235             :       real, intent(in) :: weight(nvischan, nrow)
     236             :       double precision, intent(inout) :: sumwt(npol, nchan)
     237             :       integer, intent(in) :: support, sampling
     238             :       integer, intent(in) ::  dopsf, rbeg, rend
     239             : 
     240             : 
     241             :       integer, intent(in)  :: loc(2, nvischan, nrow)
     242             :       integer, intent(in) :: off(2, nvischan, nrow) 
     243             :       integer :: iloc(2)
     244             :       complex, intent(in) :: phasor(nvischan, nrow) 
     245             :       double complex :: nvalue
     246             : 
     247             :       double precision :: norm
     248             :       double precision :: wt, wtx, wty
     249             : 
     250             :       logical :: onmygrid
     251             : 
     252             :       double precision :: pos(2)
     253             :       integer :: ix, iy, ipol, ichan
     254             :       integer :: apol, achan, irow
     255             :       integer :: posx, posy
     256             :       integer :: msupporty, psupporty
     257             :       integer :: msupportx, psupportx
     258    16617992 :       double precision :: cfnow(-support:support, -support:support)
     259             : 
     260             : 
     261  5716392680 :       do irow=rbeg, rend
     262  5716392680 :          if(rflag(irow).eq.0) then 
     263 48196014660 :          do ichan=1, nvischan
     264 42651879956 :             achan=chanmap(ichan)+1
     265 42651879956 :             if((achan.ge.1).and.(achan.le.nchan).and.
     266  5544134704 :      $           (weight(ichan,irow).ne.0.0)) then
     267             : C               call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c, 
     268             : C     $              scale, offset, sampling, pos, loc, off, phasor)
     269             : C      write(*,*) loc(1, ichan, irow), loc(2, ichan, irow), irow,ichan
     270 40683206292 :                if (onmygrid(loc(1,ichan, irow), nx, ny, x0, y0, nxsub, 
     271             :      $          nysub, support)) then
     272 17249060517 :                   do ipol=1, nvispol
     273 11499998520 :                      apol=polmap(ipol)+1
     274             :                      if((flag(ipol,ichan,irow).ne.1).and.
     275 17249060517 :      $                    (apol.ge.1).and.(apol.le.npol)) then
     276             : C If we are making a PSF then we don't want to phase
     277             : C rotate but we do want to reproject uvw
     278 10991852776 :                         if(dopsf.eq.1) then
     279  4584949468 :                            nvalue=cmplx(weight(ichan,irow))
     280             :                         else
     281             :                            nvalue=weight(ichan,irow)*
     282  6406903308 :      $                    (values(ipol,ichan,irow)*phasor(ichan, irow))
     283             :                         end if
     284 10991852776 :                         norm=0.0
     285 10991852776 :                         msupporty=-support
     286 10991852776 :                         psupporty=support
     287 10991852776 :                         psupportx=support
     288 10991852776 :                         msupportx=-support
     289 87934822208 :                         do iy=-support, support
     290 76942969432 :                            iloc(2)=abs(sampling*iy+off(2,ichan,irow))+1
     291 76942969432 :                            wty=convFunc(iloc(2))
     292 >62653*10^7 :                            do ix=-support, support
     293             :                               iloc(1)=abs(sampling*ix+off(1,ichan, 
     294 >53860*10^7 :      $                             irow))+1
     295 >53860*10^7 :                               wtx=convFunc(iloc(1))
     296 >53860*10^7 :                               cfnow(ix, iy)=wtx*wty
     297 >61554*10^7 :                               norm=norm+cfnow(ix,iy)
     298             :                            end do
     299             :                         end do
     300 87934822208 :                         do iy=-support, support
     301 >62653*10^7 :                            do ix=-support, support
     302 >61554*10^7 :                               cfnow(ix, iy)=cfnow(ix, iy)/norm
     303             :                            end do
     304             :                         end do 
     305 10991852776 :                         if((loc(1, ichan, irow)-support) .lt. x0) 
     306  5309909806 :      $                       msupportx=  -(loc(1, ichan, irow)-x0)
     307 10991852776 :                         if((loc(1, ichan, irow)+support).ge.(x0+nxsub)) 
     308  6195845206 :      $                       psupportx=  x0+nxsub-loc(1, ichan, irow)-1   
     309 10991852776 :                         if((loc(2, ichan, irow)-support) .lt. y0) 
     310  5230516868 :      $                       msupporty=  -(loc(2, ichan, irow)-y0)
     311 10991852776 :                         if((loc(2, ichan, irow)+support).ge.(y0+nysub)) 
     312  6162964310 :      $                       psupporty=  y0+nysub-loc(2, ichan, irow)-1   
     313             : C     write(*,*) msupportx, psupportx, msupporty, psupporty
     314 10991852776 :                         norm=0.0
     315 44863344750 :                         do iy=msupporty,psupporty
     316 33871491974 :                            posy=loc(2, ichan, irow)+iy
     317             : C           if( (posy .lt. (y0+nysub)) .and. 
     318             : C     $        (posy.ge. y0)) then
     319             : C                            iloc(2)=abs(sampling*iy+off(2,ichan,irow))+1
     320             : C                            wty=convFunc(iloc(2))
     321 >17072*10^7 :                             do ix=msupportx,psupportx
     322 >12586*10^7 :                                posx=loc(1, ichan, irow)+ix
     323             : C           if( (posx .lt. (x0+nxsub)) .and. 
     324             : C     $       (posx .ge. x0)) then
     325             : C            write(*,*) posx, posy, loc(1), loc(2), x0, y0, nxsub, nysub
     326             : C                                iloc(1)=abs(sampling*ix+off(1,ichan, 
     327             : C     $                                 irow))+1
     328             : C                                  wtx=convFunc(iloc(1))
     329             : C                                  wt=wtx*wty
     330             :                                   grid(posx,posy,apol,achan)=
     331             :      $                             grid(posx, posy,apol,achan)+
     332 >12586*10^7 :      $                                   nvalue*cfnow(ix,iy)
     333 >15973*10^7 :                                     norm=norm+cfnow(ix,iy)
     334             : C              write(*,*) iloc(1), iloc(2), posx, posy
     335             : C               end if
     336             :                               end do
     337             : C               end if
     338             :                         end do
     339             :                         sumwt(apol,achan)=sumwt(apol,achan)+
     340 10991852776 :      $                       weight(ichan,irow)*norm
     341             :                      end if
     342             :                   end do
     343             :                end if
     344             : C          if onmygrid
     345             :             end if
     346             :          end do
     347             :          end if
     348             :       end do
     349    16617992 :       return 
     350             :       end
     351             : 
     352             : C      subroutine sectggrids(uvw, dphase, values, nvispol, nvischan,
     353             : C     $     dopsf, flag, rflag, weight, nrow, 
     354             : C     $     scale, offset, grid, nx, ny, npol, nchan, freq, c,
     355             : C     $     support, sampling, convFunc, chanmap, polmap, sumwt, x0,
     356             : C     $     y0, nxsub, nysub, rbeg, rend)
     357             : 
     358       77336 :        subroutine sectggrids(values, nvispol, nvischan,
     359       77336 :      $   dopsf, flag, rflag, weight, nrow, 
     360       77336 :      $   grid, nx, ny, npol, nchan, 
     361       77336 :      $   support, sampling, convFunc, chanmap, polmap, sumwt, x0,
     362       77336 :      $    y0, nxsub, nysub, rbeg, rend, loc, off, phasor)
     363             : 
     364             : 
     365             : 
     366             :       implicit none
     367             :       
     368             :       integer, intent(in) ::  nx,ny,npol,nchan, nvispol, nvischan, nrow
     369             :       complex, intent(in) :: values(nvispol, nvischan, nrow)
     370             :       complex, intent(inout) :: grid(nx, ny, npol, nchan)
     371             :       integer, intent(in) ::  x0, y0, nxsub, nysub
     372             :       double precision, intent(in) ::  convFunc(*)
     373             :       integer, intent(in) :: chanmap(nvischan), polmap(nvispol)
     374             :       integer, intent(in) ::  flag(nvispol, nvischan, nrow)
     375             :       integer, intent(in) ::  rflag(nrow)
     376             :       real, intent(in) ::  weight(nvischan, nrow)
     377             :       double precision, intent(inout) ::  sumwt(npol, nchan)
     378             :       integer, intent(in) :: support, sampling
     379             :       integer, intent(in) ::  dopsf, rbeg, rend
     380             : 
     381             :        integer, intent(in)  :: loc(2, nvischan, nrow)
     382             :       integer, intent(in) :: off(2, nvischan, nrow) 
     383             :       integer :: iloc(2)
     384             :       complex, intent(in) :: phasor(nvischan, nrow) 
     385             :       double complex :: nvalue
     386             : 
     387             :       double precision :: norm
     388             :       double precision :: wt, wtx, wty
     389             : 
     390             :       logical :: onmygrid
     391             : 
     392             :       double precision :: pos(2)
     393             :       integer :: ix, iy, ipol, ichan
     394             :       integer :: apol, achan, irow
     395             :       integer :: posx, posy
     396             :       integer :: msupporty, psupporty
     397             :       integer :: msupportx, psupportx
     398       77336 :       double precision :: cfnow(-support:support, -support:support)
     399             : 
     400             : 
     401    25091312 :       do irow=rbeg, rend
     402    25091312 :          if(rflag(irow).eq.0) then 
     403  1070526672 :          do ichan=1, nvischan
     404  1046038456 :             achan=chanmap(ichan)+1
     405  1046038456 :             if((achan.ge.1).and.(achan.le.nchan).and.
     406    24488216 :      $           (weight(ichan,irow).ne.0.0)) then
     407             : C               call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c, 
     408             : C     $              scale, offset, sampling, pos, loc, off, phasor)
     409   123928760 :                if (onmygrid(loc(1,ichan, irow), nx, ny, x0, y0, nxsub, 
     410             :      $          nysub, support)) then
     411    97195632 :                   do ipol=1, nvispol
     412    65318562 :                      apol=polmap(ipol)+1
     413             :                      if((flag(ipol,ichan,irow).ne.1).and.
     414    97195632 :      $                    (apol.ge.1).and.(apol.le.npol)) then
     415             : C If we are making a PSF then we don't want to phase
     416             : C rotate but we do want to reproject uvw
     417    51146388 :                         if(dopsf.eq.1) then
     418    49581966 :                            nvalue=cmplx(weight(ichan,irow))
     419             :                         else
     420             :                            nvalue=weight(ichan,irow)*
     421     1564422 :      $                    (values(ipol,ichan,irow)*phasor(ichan, irow))
     422             :                         end if
     423    51146388 :                         norm=0.0
     424   111679308 :                         do iy=-support, support
     425    60532920 :                            iloc(2)=abs(sampling*iy+off(2,ichan,irow))+1
     426    60532920 :                            wty=convFunc(iloc(2))
     427   237917952 :                            do ix=-support, support
     428             :                               iloc(1)=abs(sampling*ix+off(1,ichan, 
     429   126238644 :      $                             irow))+1
     430   126238644 :                               wtx=convFunc(iloc(1))
     431   126238644 :                               cfnow(ix, iy)=wtx*wty
     432   186771564 :                               norm=norm+cfnow(ix,iy)
     433             :                            end do
     434             :                         end do
     435   111679308 :                         do iy=-support, support
     436   237917952 :                            do ix=-support, support
     437   186771564 :                               cfnow(ix, iy)=cfnow(ix, iy)/norm
     438             :                            end do
     439             :                         end do
     440    51146388 :                         norm=0.0
     441    51146388 :                         msupporty=-support
     442    51146388 :                         psupporty=support
     443    51146388 :                         psupportx=support
     444    51146388 :                         msupportx=-support
     445    51146388 :                         if((loc(1, ichan, irow)-support) .lt. x0) 
     446      658284 :      $                       msupportx=  -(loc(1, ichan, irow)-x0)
     447    51146388 :                         if((loc(1, ichan, irow)+support).ge.(x0+nxsub)) 
     448     4879990 :      $                       psupportx=  x0+nxsub-loc(1, ichan, irow)-1   
     449    51146388 :                         if((loc(2, ichan, irow)-support) .lt. y0) 
     450      659190 :      $                       msupporty=  -(loc(2, ichan, irow)-y0)
     451    51146388 :                         if((loc(2, ichan, irow)+support).ge.(y0+nysub)) 
     452     4693580 :      $                       psupporty=  y0+nysub-loc(2, ichan, irow)-1   
     453   102266944 :                         do iy=msupporty,psupporty
     454    51120556 :                            posy=loc(2, ichan, irow)+iy
     455   164891000 :                            do ix=msupportx,psupportx
     456    62624056 :                               posx=loc(1, ichan, irow)+ix
     457             :                               grid(posx,posy,apol,achan)=
     458             :      $                             grid(posx, posy,apol,achan)+
     459    62624056 :      $                             nvalue*cfnow(ix, iy)
     460   113744612 :                               norm=norm+cfnow(ix,iy)
     461             :                            end do
     462             :                         end do
     463             :                         sumwt(apol,achan)=sumwt(apol,achan)+
     464    51146388 :      $                       weight(ichan,irow)*norm
     465             :                      end if
     466             :                   end do
     467             :                end if
     468             :             end if
     469             :          end do
     470             :          end if
     471             :       end do
     472       77336 :       return 
     473             :       end
     474             : 
     475             : 
     476             : 
     477             : C
     478             : C Degrid a number of visibility records
     479             : C
     480           0 :       subroutine dgrid (uvw, dphase, values, nvispol, nvischan,
     481           0 :      $     flag, rflag,
     482           0 :      $     nrow, rownum, scale, offset, grid, nx, ny, npol, nchan, freq,
     483             :      $     c, support, sampling, convFunc, chanmap, polmap)
     484             : 
     485             :       implicit none
     486             :       integer nx, ny, npol, nchan, nvispol, nvischan, nrow
     487             :       complex values(nvispol, nvischan, nrow)
     488             :       complex grid(nx, ny, npol, nchan)
     489             :       double precision uvw(3, nrow), freq(nvischan), c, scale(2),
     490             :      $     offset(2)
     491             :       double precision dphase(nrow), uvdist
     492             :       complex phasor
     493             :       integer flag(nvispol, nvischan, nrow)
     494             :       integer rflag(nrow)
     495             :       integer rownum
     496             :       integer support, sampling
     497             :       integer chanmap(*), polmap(*)
     498             : 
     499             :       complex nvalue
     500             : 
     501             :       double precision convFunc(*)
     502             :       real norm
     503             : 
     504             :       logical ogrid
     505             : 
     506             :       double precision pos(2)
     507             :       integer loc(2), off(2), iloc(2)
     508             :       integer rbeg, rend
     509             :       integer ix, iy, ipol, ichan
     510             :       integer apol, achan, irow
     511             :       real wt, wtx, wty
     512             : 
     513           0 :       irow=rownum
     514             : 
     515           0 :       if(irow.ge.0) then
     516           0 :          rbeg=irow+1
     517           0 :          rend=irow+1
     518             :       else 
     519           0 :          rbeg=1
     520           0 :          rend=nrow
     521             :       end if
     522             : 
     523           0 :       do irow=rbeg, rend
     524           0 :          if(rflag(irow).eq.0) then
     525           0 :          do ichan=1, nvischan
     526           0 :             achan=chanmap(ichan)+1
     527           0 :             if((achan.ge.1).and.(achan.le.nchan)) then
     528             :                call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c,
     529           0 :      $              scale, offset, sampling, pos, loc, off, phasor)
     530           0 :                if (ogrid(nx, ny, loc, support)) then
     531           0 :                   do ipol=1, nvispol
     532           0 :                      apol=polmap(ipol)+1
     533             :                      if((flag(ipol,ichan,irow).ne.1).and.
     534           0 :      $                    (apol.ge.1).and.(apol.le.npol)) then
     535           0 :                         nvalue=0.0
     536           0 :                         norm=0.0
     537           0 :                         do iy=-support,support
     538           0 :                            iloc(2)=abs(sampling*iy+off(2))+1
     539           0 :                            wty=convFunc(iloc(2))
     540           0 :                            do ix=-support,support
     541           0 :                               iloc(1)=abs(sampling*ix+off(1))+1
     542           0 :                               wtx=convFunc(iloc(1))
     543           0 :                               wt=wtx*wty
     544           0 :                               norm=norm+wt
     545             :                               nvalue=nvalue+wt*
     546           0 :      $                             grid(loc(1)+ix,loc(2)+iy,apol,achan)
     547             :                            end do
     548             :                         end do
     549             :                         values(ipol,ichan,irow)=(nvalue*conjg(phasor))
     550           0 :      $                       /norm
     551             :                      end if
     552             :                   end do
     553             :                end if
     554             :             end if
     555             :          end do
     556             :          end if
     557             :       end do
     558           0 :       return
     559             :       end
     560             : 
     561             : 
     562      965194 :       subroutine sectdgrid (values, nvispol, nvischan,
     563      965194 :      $     flag, rflag,
     564      965194 :      $     nrow, grid, nx, ny, npol, nchan, 
     565             :      $     support, sampling, convFunc, chanmap, polmap, rbeg, rend, 
     566      965194 :      $     loc,off,phasor)
     567             : 
     568             :       implicit none
     569             :       integer, intent(in) :: nx,ny,npol,nchan,nvispol,nvischan,nrow
     570             :       complex, intent(inout) :: values(nvispol, nvischan, nrow)
     571             :       complex, intent(in) :: grid(nx, ny, npol, nchan)
     572             :       integer, intent(in) :: flag(nvispol, nvischan, nrow)
     573             :       integer, intent(in) :: rflag(nrow)
     574             :       integer, intent(in) :: support, sampling
     575             :       integer, intent(in) ::  chanmap(*), polmap(*)
     576             : 
     577             :       complex :: nvalue
     578             : 
     579             :       double precision, intent(in) :: convFunc(*)
     580             :       real norm
     581             : 
     582             :       logical ogrid
     583             :       integer, intent(in)  :: loc(2, nvischan, nrow)
     584             :       integer, intent(in) :: off(2, nvischan, nrow) 
     585             :       complex, intent(in) :: phasor(nvischan, nrow) 
     586             :       integer :: iloc(2)
     587             :       integer, intent(in) :: rbeg, rend
     588             :       integer ix, iy, ipol, ichan
     589             :       integer apol, achan, irow
     590             :       real wt, wtx, wty
     591             : 
     592             : 
     593    56068336 :       do irow=rbeg, rend
     594    56068336 :          if(rflag(irow).eq.0) then
     595   521365933 :          do ichan=1, nvischan
     596   467600784 :             achan=chanmap(ichan)+1
     597   521365933 :             if((achan.ge.1).and.(achan.le.nchan)) then
     598             : C               call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c,
     599             : C     $              scale, offset, sampling, pos, loc, off, phasor)
     600   460314024 :                if (ogrid(nx, ny, loc(1,ichan, irow) , support)) then
     601  1380731073 :                   do ipol=1, nvispol
     602   920615582 :                      apol=polmap(ipol)+1
     603             :                      if((flag(ipol,ichan,irow).ne.1).and.
     604  1380731073 :      $                    (apol.ge.1).and.(apol.le.npol)) then
     605   895248406 :                         nvalue=0.0
     606   895248406 :                         norm=0.0
     607  7161987248 :                         do iy=-support,support
     608             :                            iloc(2)=abs(sampling*iy+off(2,ichan, 
     609  6266738842 :      $                                 irow))+1
     610  6266738842 :                            wty=convFunc(iloc(2))
     611 51029159142 :                            do ix=-support,support
     612             :                               iloc(1)=abs(sampling*ix+off(1,ichan,
     613 43867171894 :      $                             irow))+1
     614 43867171894 :                               wtx=convFunc(iloc(1))
     615 43867171894 :                               wt=wtx*wty
     616 43867171894 :                               norm=norm+wt
     617             :                               nvalue=nvalue+wt*
     618             :      $                             grid(loc(1, ichan, irow)+ix,
     619 50133910736 :      $                             loc(2, ichan, irow)+iy,apol,achan)
     620             :                            end do
     621             :                         end do
     622             :                         values(ipol,ichan,irow)=(nvalue*conjg(
     623             :      $                       phasor(ichan, irow)))
     624   895248406 :      $                       /norm
     625             :                      end if
     626             :                   end do
     627             :                end if
     628             :             end if
     629             :          end do
     630             :          end if
     631             :       end do
     632      965194 :       return
     633             :       end
     634             : 
     635        1088 :       subroutine sectdgridjy (values, nvispol, nvischan,
     636        1088 :      $     flag, rflag,
     637        1088 :      $     nrow, grid, nx, ny, npol, nchan, 
     638             :      $     support, sampling, convFunc, chanmap, polmap, rbeg, rend, 
     639        1088 :      $     loc,off,phasor, scalechan)
     640             : 
     641             :       implicit none
     642             :       integer, intent(in) :: nx,ny,npol,nchan,nvispol,nvischan,nrow
     643             :       complex, intent(inout) :: values(nvispol, nvischan, nrow)
     644             :       complex, intent(in) :: grid(nx, ny, npol, nchan)
     645             :       integer, intent(in) :: flag(nvispol, nvischan, nrow)
     646             :       integer, intent(in) :: rflag(nrow)
     647             :       integer, intent(in) :: support, sampling
     648             :       integer, intent(in) ::  chanmap(*), polmap(*)
     649             :       complex :: nvalue
     650             : 
     651             :       double precision, intent(in) :: convFunc(*), scalechan(*)
     652             :       real norm
     653             : 
     654             :       logical ogrid
     655             :       integer, intent(in)  :: loc(2, nvischan, nrow)
     656             :       integer, intent(in) :: off(2, nvischan, nrow) 
     657             :       complex, intent(in) :: phasor(nvischan, nrow) 
     658             :       integer :: iloc(2)
     659             :       integer, intent(in) :: rbeg, rend
     660             :       integer ix, iy, ipol, ichan
     661             :       integer apol, achan, irow
     662             :       real wt, wtx, wty
     663             : 
     664             : 
     665       96560 :       do irow=rbeg, rend
     666       96560 :          if(rflag(irow).eq.0) then
     667     7368192 :          do ichan=1, nvischan
     668     7272720 :             achan=chanmap(ichan)+1
     669     7368192 :             if((achan.ge.1).and.(achan.le.nchan)) then
     670             : C               call sgrid(uvw(1,irow), dphase(irow), freq(ichan), c,
     671             : C     $              scale, offset, sampling, pos, loc, off, phasor)
     672      322920 :                if (ogrid(nx, ny, loc(1,ichan, irow) , support)) then
     673     1109160 :                   do ipol=1, nvispol
     674      786240 :                      apol=polmap(ipol)+1
     675             :                      if((flag(ipol,ichan,irow).ne.1).and.
     676     1109160 :      $                    (apol.ge.1).and.(apol.le.npol)) then
     677      645840 :                         nvalue=0.0
     678      645840 :                         norm=0.0
     679     5166720 :                         do iy=-support,support
     680             :                            iloc(2)=abs(sampling*iy+off(2,ichan, 
     681     4520880 :      $                                 irow))+1
     682     4520880 :                            wty=convFunc(iloc(2))
     683    36812880 :                            do ix=-support,support
     684             :                               iloc(1)=abs(sampling*ix+off(1,ichan,
     685    31646160 :      $                             irow))+1
     686    31646160 :                               wtx=convFunc(iloc(1))
     687    31646160 :                               wt=wtx*wty
     688    31646160 :                               norm=norm+wt
     689             :                               nvalue=nvalue+wt*scalechan(ichan)*
     690             :      $                             grid(loc(1, ichan, irow)+ix,
     691    36167040 :      $                             loc(2, ichan, irow)+iy,apol,achan)
     692             :                            end do
     693             :                         end do
     694             :                         values(ipol,ichan,irow)=(nvalue*conjg(
     695             :      $                       phasor(ichan, irow)))
     696      645840 :      $                       /norm
     697             :                      end if
     698             :                   end do
     699             :                end if
     700             :             end if
     701             :          end do
     702             :          end if
     703             :       end do
     704        1088 :       return
     705             :       end
     706             : 
     707             : 
     708             : 
     709             : C
     710             : C Calculate gridded coordinates and the phasor needed for
     711             : C phase rotation.
     712             : C
     713           0 :       subroutine sgrid (uvw, dphase, freq, c, scale, offset, sampling, 
     714             :      $     pos, loc, off, phasor)
     715             :       implicit none
     716             :       integer sampling
     717             :       integer loc(2), off(2)
     718             :       double precision uvw(3), freq, c, scale(2), offset(2)
     719             :       double precision pos(2)
     720             :       double precision dphase, phase
     721             :       complex phasor
     722             :       integer idim
     723             :       double precision pi
     724             :       data pi/3.14159265358979323846/
     725             : 
     726           0 :       do idim=1,2
     727           0 :          pos(idim)=scale(idim)*uvw(idim)*freq/c+(offset(idim)+1.0)
     728           0 :          loc(idim)=nint(pos(idim))
     729           0 :          off(idim)=nint((loc(idim)-pos(idim))*sampling)
     730             :       end do
     731           0 :       phase=-2.0D0*pi*dphase*freq/c
     732           0 :       phasor=cmplx(cos(phase), sin(phase))
     733           0 :       return 
     734             :       end
     735             : C
     736             : C Is this on the grid?
     737             : C
     738   460636944 :       logical function ogrid (nx, ny, loc, support)
     739             :       implicit none
     740             :       integer nx, ny, loc(2), support
     741             :       ogrid=(loc(1)-support.ge.1).and.(loc(1)+support.le.nx).and.
     742   460636944 :      $        (loc(2)-support.ge.1).and.(loc(2)+support.le.ny)
     743   460636944 :       return
     744             :       end
     745             : 
     746 40807135052 :       logical function onmygrid(loc, nx, ny, nx0, ny0, nxsub, nysub, 
     747             :      $     support)
     748             :       implicit none
     749             :       integer nx, ny, nx0, ny0, nxsub, nysub, loc(2), support
     750             :       logical ogrid
     751             : C      logical onmygrid
     752             : C      onmygrid=ogrid(nx, ny, loc, support)
     753             : C      if(.not. onmygrid) then
     754             : C         return 
     755             : C      end if
     756             : C      onmygrid=(loc(1)-nx0 .ge. (-support)) .and. ((loc(1)-nx0-nxsub) 
     757             : C     $     .le. (support)) .and.((loc(2)-ny0) .ge. (-support)) .and. 
     758             : C     $ ((loc(2)-ny0-nysub) .le. (support))  
     759             : C----------------------------
     760             :       integer loc1sub, loc1plus, loc2sub, loc2plus
     761 40807135052 :       loc1sub=loc(1)-support
     762 40807135052 :       loc1plus=loc(1)+support
     763 40807135052 :       loc2sub=loc(2)-support
     764 40807135052 :       loc2plus=loc(2)+support
     765             : C----------------
     766             : C      onmygrid=(loc1plus .ge. nx0) .and. (loc1sub 
     767             : C     $     .le. (nx0+nxsub)) .and.(loc2plus .ge. ny0) .and. 
     768             : C     $     (loc2sub .le. (ny0+nysub))
     769             : C--------------
     770             :       onmygrid=(loc1plus .ge. nx0) .and. (loc1sub 
     771             :      $     .le. (nx0+nxsub)) .and.(loc2plus .ge. ny0) .and. 
     772             :      $     (loc2sub .le. (ny0+nysub)) .and. (loc1sub.ge.1)
     773             :      $     .and.(loc1plus.le.nx).and.
     774 40807135052 :      $     (loc2sub.ge.1).and.(loc2plus.le.ny)
     775             :      
     776 40807135052 :       return 
     777             :       end
     778   249422610 :       subroutine locuvw(uvw, dphase, freq, nvchan, scale, offset,
     779   249422610 :      $     sampling, loc, off, phasor, irow, doW, Cinv)
     780             :       implicit none
     781             :       double precision,  intent(in) :: uvw(3, *), dphase(*), freq(*)
     782             :       integer, intent(in) :: nvchan, sampling, irow, doW
     783             :       double precision, intent(in) :: scale(3), offset(3), Cinv
     784             :       integer, intent(inout) :: loc(2+doW, nvchan, *), off(2+doW,
     785             :      $     nvchan, *)
     786             :       complex, intent(inout) :: phasor(nvchan, *)
     787             :       integer :: f, k, row
     788             :       double precision :: phase, pos
     789             :       double precision :: pi
     790             :       data pi/3.14159265358979323846/
     791             : 
     792   249422610 :       row=irow+1
     793             : 
     794  2278728907 :       do f=1, nvchan
     795  6087918891 :          do k=1,2
     796  4058612594 :             pos=scale(k)*uvw(k,row)*freq(f)*Cinv+(offset(k)+1.0)
     797  4058612594 :             loc(k, f, row)=nint(pos)
     798  6087918891 :             off(k, f, row)=nint((loc(k, f, row)-pos)*sampling)
     799             :          end do
     800  2029306297 :          phase=-2.0D0*pi*dphase(row)*freq(f)*Cinv
     801  2029306297 :          phasor(f,row)=cmplx(cos(phase), sin(phase))
     802  2278728907 :          if(doW .eq. 1) then
     803             :             pos=sqrt(abs(scale(3)*uvw(3, row)*freq(f)*Cinv))+offset(3)
     804    18268146 :      $           +1.0
     805    18268146 :             loc(3,f, row)=nint(pos)
     806    18268146 :             off(3, f, row)=0
     807             :          end if
     808             :       end do
     809   249422610 :       return
     810             :       end

Generated by: LCOV version 1.16