C+
C NAME:
C       subroutine pixelcon
C PURPOSE:
C       To convert x y values from 
C
C CALLING SEQUENCE:
C       ./ashi_flat_field
C INPUTS:
C     it1               type of removal
C                       1 - a rectangular spot centered on x and y of half extent size
C                       2 - a circular spot centered on x and y of radius size
C                       3 - an angular extent with theta centered on x and y of angular extent size
C FUNCTIONS/SUBROUTINES:
C
C MODIFICATION HISTORY:
C     September 2023, Bernard Jackson (UCSD)
C-    

      subroutine pixelcon(it1,x,y,theax,theay,sizet,rdiam,dhmean,NXY,NOBB,NOBE,distcx,distcy)

        real dhmean   (NXY,NXY,NOBE)                     ! the mean value of the +-5 box that will be zeroed       

        if(it1.eq.1) then
C         Rectangle limit
          ix    = nint(x)
          iy    = nint(y)
          write(*,'(A,F10.1,A,F10.1,A,F5.1)') 'Blank a rectangle, x =',x,' y = ',y,' half size = ',sizet
          isize = nint(sizet)
          do K=NOBB,NOBE
            do J=1,NXY
              do I=1,NXY
                if(I.ge.(ix-isize).and.I.le.(ix+isize).and.J.ge.(iy-isize).and.J.le.(iy+isize)) then
                  dhmean(I,J,K) = 0.0
                end if                
              end do
            end do
          end do
        end if          
C       Circle limit
        if(it1.eq.2) then
          write(*,'(A,F10.1,A,F10.1,A,F5.1)') 'Blank a circle x =',x,' y = ',y,' radius = ',rdiam 
          do K=NOBB,NOBE
            do J=1,NXY
              do I=1,NXY
                xi = float(I)
                yj = float(J)
                rad = sqrt((xi-x)**2 + (yj-y)**2)
                if(rad.le.rdiam) then
                  dhmean(I,J,K) = 0.0
                end if                
              end do
            end do
          end do 
        end if                    

        if(it1.eq.3) then
C         Theta limit
          thea = atan2d((theay-distcy),(theax-distcx))

          write(*,'(A,F8.1,A,F6.1)')'Blank a sector. Theta =',thea,'   +-angle = ',sizet
          do K=NOBB,NOBE
            do J=1,NXY
              do I=1,NXY
                xi = float(I)
                yj = float(J)
                theta = atan2d((yj-distcy),(xi-distcx))
                if(theta.ge.(thea-sizet).and.theta.le.(thea+sizet)) then
                  dhmean(I,J,K) = 0.0                  
                end if
              end do
            end do
          end do
        end if
        
        if(it1.eq.4) then
C         Rectangle area check
          ix    = nint(x)
          iy    = nint(y)
          write(*,'(A,F10.1,A,F10.1,A,F5.1)') 'Check a rectangular patch, x =',x,' y = ',y,' Dist from llc = ',sizet
          isize = nint(sizet)
          do K=NOBB,NOBE
            do I=1,NXY
              do J=1,NXY
                if(J.lt.iy.or.J.gt.(iy+isize)) then
                  dhmean(I,J,K) = 0.0
                end if
                if(I.lt.ix.or.I.gt.(ix+isize)) then
                  dhmean(I,J,K) = 0.0
                end if
              end do
            end do
          end do
          value1 = 0.0
          avalue = 0.0
          do K=NOBB,NOBE
            do I=1,NXY
              do J=1,NXY
                if(dhmean(I,J,K).gt.200.0) then
                  value1 = value1 + dhmean(I,J,K)
                  avalue = avalue + 1.0
                end if
              end do
            end do
          end do
         print*,' '
C          write(*,'(3(A,F11.1))') ' Total removed = ',dmean0,' Total in all = ',2048.*2048.,' Diff = ',2048.*2048.-dmean0
          print*, 'Average of the total number kept', value1,avalue,value1/avalue
          print*,' '
        end if          

        return
      end