C+
C NAME:
C	ashi_img_rmv.f
C PURPOSE:
C	Given the ASHI image location cfileo the image pixel parameters, its DSAA values, and the floating image (a itype 14 image)
C     or more normally the file provided (a type 15 file) with the icirx and iciry parameters, write out the image for use in Surfer
C     as an F9.2 file.
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
C	call ashi_img_rmv(NX,NY2,icirx,iciry,aimage,aimages)
C INPUTS:
C	NX	              integer	Number of image brightness values possible for x in the file (generally 2048)
C	NY2	              integer	Number of image brightness values possible for y in the file (usually 1024)
C	icirx	              integer	Maximum number of image brightness values possible for x in the file (generally 1024)
C	icirx	              integer	Maximum number of image brightness values possible for y in the file (generally 512)
C     aimage(NX,NY2)	  real	Floating grd image
C OUTPUTS:
C     aimage(icirx,iciry) real	Floating grd image
C
C FUNCTIONS/SUBROUTINES:
C
C MODIFICATION HISTORY:
C	October-2023, Bernard Jackson (UCSD)
C-
	subroutine ashi_img_rmv(NX,NY2,icirx,iciry,aimage,aimages)

        real            aimage  (NX,NY2)              ! image brightness values input generally x = 2048, y = 1024)  
        real            aimagem (NX,NY2)
        real            aimaged (NX,NY2)
        real            aimagedd(NX,NY2)
        real            aimagess(icirx,iciry)  
        real            aimages (icirx,iciry)         !  image brightness values output generally x = 1024, y = 512)   
        real            anumm   (icirx,iciry) 
        
          print*, 'Into ashi_img_rmv.f',NX,NY2,icirx,iciry
C         1024 j values = 180 deg input is about 5.7 pixels / degree
C         2048 i values = 360 deg input is about 5.7 pixels / degree

C          NMJ = 30                                                ! J step size to determine mean values (don't throw out the baby)
C          NMI = 60                                                ! X step size to determine mean values (don't throw out the baby)         
          NMJ = 20                                                ! J step size to determine mean values (don't throw out the baby)
          NMI = 40                                                ! X step size to determine mean values (don't throw out the baby)         

          NBJ = 5                                                 ! J step size to determine values removed
          NBI = 11                                                ! X step size to determine values removed

          aimgmax =  10000.0
          aimgmin = -10000.0
          CONRV = 5.0

          do I=1,NX
            do J=1,NY2
              aimagedd(I,J) = 0.0
              aimaged(I,J)  = 0.0
              aimagem(I,J)  = 0.0
            end do
          end do
          do K=1,3
            write(*,'(A,F12.2,A)') 'Maximum +- value from mean of image kept =',aimgmax,' ADU' 
            do I=1,NX
              do J=1,NY2
                dnump  =  0.0
                anumpv =  0.0
                anummm = 0.0
                anummv = 0.0
                dhum   = 0.0
                anum   = 0.0
                do JJ=J-NMJ,J+NMJ                  
                  do II=I-NMI,I+NMI
                    if(JJ.ge.1.and.JJ.le.NY2.and.II.ge.1.and.II.le.NX) then
                      if(aimage(II,JJ).gt.-aimgmax.and.aimage(II,JJ).lt.aimgmax) then
                        dhum = dhum + aimage(II,JJ)
                        anum = anum + 1.0
                      else
                        if(aimage(II,JJ).ge.aimgmax) then
                          if(aimage(II,JJ).gt.anumpv) anumpv = aimage(II,JJ)                  
                          anump = anump + 1.0
                        end if
                        if(aimage(II,JJ).le.-aimgmax) then
                          if(aimage(II,JJ).le.anummv) anummv = aimage(II,JJ)                 
                          anummn = anummn + 1.0
                        end if
                      end if
                    end if
                  end do
                end do
                if(anum.ne.0) then
                  aimagem(I,J) = dhum/anum                              ! mean value of the +-11 and +-5 box +-10 ADU
                end if
              end do
            end do                                                      ! This provides a mean of values throughout the input image
            aimgmax = aimgmax/10.0
            write(*,'(F12.1,A,F10.2)'),anump, ' values are above the maximum with the greatest =',anumpv
            write(*,'(F12.1,A,F10.2)'),anummm,' values are below the minimum with the smallest =',anummv
          end do

          do I=1,NX
            do J=1,NY2
              aimaged(I,J) = aimage(I,J) - aimagem(I,J)                 ! remove the mean from the original input image
            end do
          end do
          write(*,'(A,F7.2,A)') 'Maximum +- value from mean of image now kept at end =',aimgmax,' ADU' 

C         Now, provide mean values from this new image (aimaged) - anything over aimgmax ADU away from the norm is removed in the final image
          do I=1,NX
            do J=1,NY2
              dnump  =  0.0
              anumpv =  0.0
              anummm  = 0.0
              anummv = 0.0
              dhum   = 0.0
              anum   = 0.0
              do JJ=J-NBJ,J+NBJ                  
                do II=I-NBI,I+NBI
                  if(JJ.ge.1.and.JJ.le.NY2.and.II.ge.1.and.II.le.NX) then
                    if(aimaged(II,JJ).gt.-aimgmax.and.aimaged(II,JJ).lt.aimgmax) then
                      dhum = dhum + aimaged(II,JJ)
                      anum = anum + 1.0                        
                    else

                      if(aimaged(II,JJ).ge.aimgmax) then
                        if(aimaged(II,JJ).gt.anumpv) anumpv = aimaged(II,JJ)                  
                        anump = anump + 1.0
                        
                      end if
                      if(aimaged(II,JJ).le.-aimgmax) then
                        if(aimaged(II,JJ).le.anummv) anummv = aimaged(II,JJ)                 
                        anummm = anummm + 1.0 
                      end if
                    end if
                  end if
                end do
              end do
              if(anum.ne.0) then
                aimagedd(I,J) = (dhum/anum)                             ! new smoothed image of the +-11 and +-5 box +-10 ADU
              else
                aimagedd(I,J) = -99999.9                !  something wrong here
              end if
            end do
          end do
          write(*,'(F12.1,A,F10.2)'),anump, ' values are above the maximum with the greatest =',anumpv
          write(*,'(F12.1,A,F10.2)'),anummm,' values are below the minimum with the smallest =',anummv
C
C         This is the output image (Gridsphere runs vary slowly with large images and at least this works on Leela) 
C
          do II=1,icirx                                                 ! Now compress the image by a factor of 2 (area of 4)
            do JJ=1,iciry
              aimagess(II,JJ) = 0.0
              anumm   (II,JJ) = 0.0
            end do
          end do
            do II=1,icirx
              do JJ=1,iciry
                do I=II*2-1,II*2
                  do J=JJ*2-1,JJ*2
                    if(aimagedd(I,J).gt.-aimgmax.and.aimagedd(I,J).lt.aimgmax) then ! if any of the 4 pixels > aimgmax its mean, remove. Else ave.
                      aimagess(II,JJ) = aimagess(II,JJ) + aimagedd(I,J)
                      anumm(II,JJ) = anumm(II,JJ) + 1.0
                    end if
                  end do
                end do
              end do
            end do
            do II=1,icirx
              do JJ=1,iciry
                if(anumm(II,JJ).ne.0.0) then                
                  aimages(II,JJ) = aimagess(II,JJ)/anumm(II,JJ)
                end if
              end do
            end do
C
C           Now run gridsphere2D to smooth over the poles (Gridsphere is designed to only run on cylindrical images)
C           Gridsphere will not run on the original input images because they are too big
C           CONRV is the gaussian image smoother that at a value of 5.0 smooths images at 1/e to about 2.0 degrees circular everywhere
C             (This means that no stars, even polar ones, input very much structue over the poles)
C                 
            igrid = 1
            if(igrid.eq.1) then
              print *,'Gridsphere is being used 1/e +- CONRV =',CONRV,' degrees'
              call GridSphere2D(1.0,icirx,iciry,1,aimages,CONRV,0,0.0,0.0)               ! Smooth the final image by Gridsphere2D 
            end if

          print*, 'Out of ashi_img_rmv.f'
	  return
	end