C+
C NAME:
C       subroutine ashi_img
C PURPOSE:
C       To read the ASHI image grd files, rotate, stretch and combine them together
C CALLING SEQUENCE:
C       call ashi_img(iproc,IMGIN,heliolat,cpfilei,crfileo,cdfileo,NX,NY,dmax,cimghead,aimage,aimageout,aimagedif)
C INPUTS:
C       iproc           integer           Type operation to perform as a process
C                                               1 sequential image differences ecliptic coordinates
C                                               2 differences from the overall image base ecliptic coordinates
C                                               3 sequential image differences ecliptic coordinates
C                                               4 differences from the overall image base heliographic coordinates
C                                               5 
C       IMGIN                    integer           Number of composite images to input
C       cpfilei(IMGIN)           character         name of the image file to be input                
C       crfileo(IMGIN)           character         name of the reduced image file to be output 
C       cdfileo(IMGIN)           character         name of the difference image file to be output 
C       NX                       integer           size of the image file input (presumed 2048 x 2048)
C       NY                       integer           size of the somoothed output files 
C       dmax                     real              the cutoff size of the input image circle
C OUTPUTS
C       cimghead                 character         the image header of the output files
C       aimage   (NX,NY,IMGIN)     real            the images read in 
C       aimageout(NX,NY,IMGIN)     real            the images output      
C       aimagedif(NX/2,NY/4,IMGIN) real            the different images output
CC      aimagetst(NX,NY/2,IMGIN)   real            the test images output
C
C FUNCTIONS/SUBROUTINES:
C       read_grd_image
C
C MODIFICATION HISTORY:
C       November 2022, Bernard Jackson (UCSD)
C-    
        subroutine ashi_img(iproc,IMGIN,cpfilei,crfileo,cdfileo,NX,NY,dmax,cimghead,aimagei,aimageout,aimagedif)

          parameter     (ihead = 5)                                     ! # of lines in the GRD file header
          parameter     (NRA   = 3600)                                  ! size in RA of images
          parameter     (NDEC  = 1800)                                  ! size in Dec of images

          character     cpfileii (IMGIN)*66                             ! Dummy file to fill to input to write_image_grd_p
          character     cpfilei  (IMGIN)*47                             ! Location and name of an input file if there is one     
          character     cfileoo*66                                      ! Dummy file to fill to input to write_image_grd_p
          character     crfileo  (IMGIN)*47                             ! Location and name of a reduced output file if there is one 
          character     cdfileo  (IMGIN)*47                             ! Location and name of a differenced output file if there is one
          character     chfileo  (IMGIN)*47                             ! Location and name of a reduced output file if there is one 
          character     crfilet*46
          character     cimghead*9
          character     cpfilebs*55 /'./ASHI_grd_files/ASHI_img_0641-2790_2901-4540_StarT.grd'/
          character     cpfilebf*55 /'./ASHI_grd_files/ASHI_img_0641-2790_2901-4540_basee.grd'/
          
          real          aimage    (NX,NY/2,IMGIN),                      ! input images to use
     &                  aimageout (NX,NY/2,IMGIN),                      ! output images
     &                  aimagetst (NX,NY/2,IMGIN),                      ! dummy image used for tests
     &                  aimagedif (NX/2,NY/4,IMGIN),                    ! output difference images
     &                  anumout   (NX,NY/2,IMGIN),                      ! number of output images     
     &                  aimagebase(NX,NY/2)                             ! image base output
C     
          real          aimagei   (NX,NY),                              ! original image input       
     &                  aimageb   (NX,NY),                              ! image from the whole image data set used as a base
     &                  aimages   (NX,NY/2),                            ! dummy input image used for gridsphere2D
     &                  aimagehas (NX/2,NY/4)                           ! dummy image used for hammer-atioff input images differences and smoothed

          real          aimagess  (NX/2,NY/4,IMGIN)                     ! image that is smoothed, differenced, and gridsphered and changed to K=IMGIN
          real          aimagesf  (NX/2,NY/2,IMGIN)                     ! image format for a fisheye view
          real          anumf     (NX/2,NY/2,IMGIN)                     ! number of image fisheye view format values
          real          aimagehass(NX/2,NY/4,IMGIN)                     ! image used for hammer-atioff differenced and smoothed output
          real          anumhass  (NX/2,NY/4,IMGIN)                     ! image number used for hammer-atioff differenced and smoothed output

          real          aimagesss (NX/2,NY/4)                           ! image that is smoothed, differenced, and gridsphered from ashi_img_rmv.f
          real          aimagesb  (NX/2,NY/4)                           ! image that is smoothed, differenced, and gridsphered, and used as a base
          real          aimagesff (NX/2,NY/2)                           ! image format for a fisheye view
          real          anumff    (NX/2,NY/2)                           ! number of image fisheye view format values
          real          aimagehasss(NX/2,NY/4)
          real          anumhasss (NX/2,NY/4)
          real          aimval    (41,IMGIN)
          real          S10v      (41)

C         Dummy array
          integer       image    (NX,NY)

          print*, ' '
          print*, 'Into subroutine ashi_img.f',NX,NY

          dmax = 3000.0                                                 ! For input into read_grd_image
          RCONST = float(NX)/360.0                                      ! pixels per degree
          CONRV = 2.0*RCONST                                            ! width of the smoothing in degrees*pixels/degrees (pixels)
          CONRV = 2.0
          icirx = NX
          iciry = NY/2
          print*, icirx,iciry

          do K=1,IMGIN
            cpfileii(K) = '                                                                     '          
            cpfileii(K) = cpfilei(K)
            print*, K,'cpfileii(K) ' ,cpfileii(K)
            if(K.eq.1) chfileo(K) =  './ASHI_grd_files/ASHI_differenced_img_0641h.grd' 
            if(K.eq.2) chfileo(K) =  './ASHI_grd_files/ASHI_differenced_img_1941h.grd'
            if(K.eq.3) chfileo(K) =  './ASHI_grd_files/ASHI_differenced_img_3241h.grd'
          end do

          do K=1,IMGIN
            do I=1,NX
              do J=1,NX
                aimagei(I,J) = 0.0        
                image  (I,J) = 0
              end do
            end do
            cpfileii(K) = '                                                                     '  
            cpfileii(K) = cpfilei(K)
            call read_grd_image(14,cpfileii(K),cfileo,NX,dmax,cimghead,aimagei,image)     ! These are half files in cylindrical coordinates X,Y=2048 Y=1024 active
            do I=1,icirx
              do J=1,iciry
                aimage(I,J,K) = aimagei(I,J)                                              ! This loads aimage with the H-A image and zeros
              end do
            end do
C            cfileoo = '                                                                     '  
C            cfileoo = crfileo(K)
C            call write_image_grd_p(cfileoo,icirx,iciry,aimage(1,1,K))                 ! These are files in cylindrical coordinates X=2048 Y=1024
          end do

C         Now change the hammer-aitoff aimage image into cylindrical coordinates if needed
          ihammer = 0
          if(ihammer.eq.1) then
            sq2 = sqrt(2.0)
            do K=1,IMGIN
              do I=1,NX
                do J=1,NY
                  X    = float(I - iricx)/float(iricx)
                  X2   = X*X
                  Y    = float(J - iricy)/float(iricx)
                  Y2   = Y*Y
                  Z2   = 1.0 - X2/2.0 - Y2/2.0
                  Z    = sqrt(Z2)
                  ALON = 2.0*atan((sq2*X*Z)/(2.0*Z2 - 1.0))
                  ALAT = asin(sq2*Y*Z)
                  JC   = (NXYD2+(NXYD2/2))-1
                  if(I.eq.1.and.J.eq.JC)        print*,I,J,ALON,ALAT
                  if(I.eq.NXYD2.and.J.eq.JC)    print*,I,J,ALON,ALAT
                  if(I.eq.NXYD2/2.and.J.eq.1)   print*,I,J,ALON,ALAT
                  if(I.eq.NXYD2/2.and.J.eq.NXY) print*,I,J,ALON,ALAT
                end do
              end do
              call write_image_grd_p(crfileo,icirx,iciry,aimage(1,1,K)) 
            end do

          end if                          ! End Hammer-Aitoff image to cylindrical
C
C         There are two ways to difference images. One is by making consecutive differences, and the other way 
C         is by making a base from all the images and subtracting the previous from the current to learn the differences
C         To attempt the first way set ibase = 0. 
C         To attempt the second way set ibase = 1
C          ibase = 0
          ibase = 1
C          ibase = 2                                                           ! This provides NO reduced images and aimage is just passed
          if(ibase.eq.0) then
            idif = 1
C           Difference each cylindrical image from the earlier
            do I=1,icirx
              do J=1,iciry
                aimages(I,J) = aimage(i,j,1)
              end do
            end do
C           Smooth and reduce the first original image in size and gridsphere2D it
            call ashi_img_rmv(icirx,iciry,icirx/2,iciry/2,aimages,aimagesss)
            do I=1,icirx/2
              do J=1,iciry/2
                aimagess(I,J,1) = aimagesss(I,J)
              end do
            end do
            do K=2,IMGIN
             do I=1,icirx
               do J=1,iciry
                 aimages(I,J) = aimage(i,j,K)
               end do
             end do
             call ashi_img_rmv(icirx,iciry,icirx/2,iciry/2,aimages,aimagesss)
             do I=1,icirx/2
               do J=1,iciry/2
                 aimagess(I,J,K) = aimagesss(I,J)
               end do
             end do
             do I=1,icirx/2
               do J=1,iciry/2
                 aimagedif(I,J,K-1) = 0.0
                 if(aimagess(I,J,K-1).ne.0.0.and.aimagess(I,J,K).ne.0.0)              
     &             aimagedif(I,J,K-idif) = aimagess(I,J,K) - aimagess(I,J,K-idif)
               end do
             end do
           end do          
C           do K=2,IMGIN
C             cfileoo = '                                                                     '  
C             cfileoo = cdfileo(K)
C             call write_image_grd_p(cfileoo,icirx/2,iciry/2,aimagedif(1,1,K-1)) ! fisheye image difference from the last
C           end do

          end if

          if(ibase.eq.1) then
            idif = 0
C           Either read in a composite image (1) or put the cylindrical images together as a base image (0)
C            imagebase = 0
            imagebase = 1
            if(imagebase.EQ.1) then                                                   ! An image base is read in
              do I=1,NX
                do J=1,NX
                  aimagei(I,J) = 0.0        
                end do
              end do
              call read_grd_image(14,cpfilebs,cfileo,NX,dmax,cimghead,aimageb,image)  ! These are half files in cylindrical coordinates X,Y=2048 Y=1024 active
              do K=1,IMGIN                                                         
                do I=1,icirx
                  do J=1,iciry
                    aimagebase(I,J) = aimageb(I,J)
                  end do
                end do
              end do
            else                                                                      ! An image base is averaged from the three individual values
              do K=1,IMGIN                                                         
                do I=1,icirx
                  do J=1,iciry
                    aimagebase(I,J) = aimagebase(I,J) + aimageout(I,J,K)*(1.0/float(IMGIN)) 
                  end do
                end do
              end do
            end if
            print*,' '
            print*, 'The base has been produced now smooth,reduce, and gridsphere it'
C           Smooth and reduce the image base in size and gridsphere2D it
            call ashi_img_rmv(icirx,iciry,icirx/2,iciry/2,aimagebase,aimagesss)
            do I=1,icirx/2
              do J=1,iciry/2
                aimagesb(I,J) = aimagesss(I,J)
              end do
            end do
C           print*, 'After base determination, now write out the base',icirx,iciry
C           call write_image_grd_p(cpfilebf,icirx/2,iciry/2,aimagesb)
C           Now subtract the base image from each individual 
            do K=1,IMGIN                                                         
              do I=1,icirx
                do J=1,iciry
                  aimages(I,J) = aimage(I,J,K)
                end do
              end do
              print*,' '
              write(*,'(A,I2,A)') 'Now smooth, reduce, and gridsphere image', K,' and subtract the base'
              call ashi_img_rmv(icirx,iciry,icirx/2,iciry/2,aimages,aimagesss)
              do I=1,icirx/2
                do J=1,iciry/2
                  aimagedif(I,J,K) = 0.0
                  if(aimagesss(I,J).ne.0.0.and.aimagesb(I,J).ne.0.0) aimagedif(I,J,K) = aimagesss(I,J) - aimagesb(I,J)
                end do
              end do
            end do
            do K=1,IMGIN
              cfileoo = '                                                                     '  
              cfileoo = cdfileo(K)
C              call write_image_grd_p(cfileoo,icirx/2,iciry/2,aimagedif(1,1,K)) ! fisheye image difference from base
            end do
C      STOP
          end if
      print*, 'End of 3 image base removal'
          ifish = 0
          if(ifish.eq.1) then
C          Now change the aimagedif into a fisheye view

           do K=1,IMGIN 
             do I=1,icirx                                                       ! This provides a cylindrical image for a fisheye image test
               do J=1,iciry
                 aimagetst(I,J,K) = aimage(I,J,K)
               end do
             end do
C            print*,'to here'
C            cfileoo = '                                                                     '  
C            cfileoo = crfileo(K)
C            call write_image_grd_p(cfileoo,icirx,iciry,aimagetst(1,1,K))                 ! These are files in cylindrical coordinates X=2048 Y=1024 
           end do

           print*, ' '
           print*, 'Number of images IMGIN =', IMGIN
           print*, ' '

           icirc = NX/2  
C          icirc = aperture*RCONST/2.0
           ipdif  = icirc/2                                                     ! difference in circle pixels from the center     
           print *,' '
           print *, 'before fisheye image maker',icirc,icirx,iciry,ipdif
           print *,' '
           pi = 3.14159265

           do K=1,IMGIN
            do I=1,icirc                                                       ! This provides a fisheye image
              do J=1,icirc            
                aimagesf (I,J,K) = 0.0
                anumf    (I,J,K) = 0.0
                aimagesff(I,J)   = 0.0
                anumff   (I,J)   = 0.0
              end do
            end do
            do I=1,icirx                                                        ! This provides a fisheye image (must go over the size of the test image)
              do J=1,iciry                                                      ! Burke in degrees

                iBurke = 1
                if(iBurke.eq.1) then                                            ! Original version that should be OK

C                  aperture = 1.0                
                  aperture = 0.75                                               ! 1.0 should put aperture in degrees from +-180 deg by +-90 deg
                  alonx = (float(I) - 1024.0)/1024.0
                  alaty = (float(J) -  512.0)/512.0
                  alonp = alonx*pi
                  alatp = alaty*pi/2.0

                  sinlat = sin(alatp)
                  coslat = cos(alatp)
                  sinlon = sin(alonp)
                  coslon = cos(alonp)
            
                  PXX  = coslat*coslon
                  PYY  = coslat*sinlon
                  PZZ  = sinlat
                  RR     =  2.0*atan2(sqrt(PXX*PXX + PZZ*PZZ),PYY)/aperture
                  BARERR = (2.0*atan2(sqrt(PXX*PXX + PZZ*PZZ),PYY))
                  THEDD  = atan2(PZZ,PXX)
C                  if(I.gt.1024.and.J.gt.512) THEDD = THEDD                                ! quadrant I
C                  if(I.le.1024.and.J.gt.512) THEDD = THEDD                                ! quadrant II
C                  if(I.le.1024.and.J.le.512)  THEDD = 360.0 + THEDD                       ! Quadrant III OLD
C                  if(I.gt.1024.and.J.le.512)  THEDD = 360.0 + THEDD                       ! Quadrant IV  OLD
C                  if(I.le.1024.and.J.le.512)  THEDD = THEDD + 180.0                       ! Quadrant III
C                  if(I.gt.1024.and.J.le.512)  THEDD = THEDD - 180.0                       ! Quadrant IV
                  AZI    = THEDD*180.0/pi - 90.0
                  if(AZI.lt.0.0) AZI = AZI + 360.0
                  XX     = RR*cos(THEDD)*180.0/pi
                  YY     = RR*sin(THEDD)*180.0/pi
                
C                  alond = alonp*180.0/pi
C                  alatd = alatp*180.0/pi

                  alond = alonx
                  alatd = alaty

                end if                  ! Used version of Burke to test

                IX   = nint(XX) + ipdif                           ! IX AND JY +- from Burke in pixels from 512 (in pixel degrees)
                JY   = nint(YY) + ipdif

                if(IX.lt.1) then
                  IX = 1
                  print*,'Below IX inner bounds', IX,JY
                end if
                if(JY.lt.1) then
                  JY = 1
                  print*,'Below JY inner bounds', IX,JY
                end if
                if(IX.gt.icirc) then
                  IX = icirc
                  print*,'Above IX upper bounds', IX,JY
                end if

                if(JY.gt.icirc) then
                  JY = icirc
                  print*,'Above JY upper bounds', IX,JY
                end if
C                aimagesf(IX,JY,K) = aimagesf(IX,JY,K) + aimagetst (I,J,K)
C                anumf   (IX,JY,K) = anumf(IX,JY,K) + 1.0
                aimagesff(IX,JY)   = aimagesff(IX,JY) + aimagetst (I,J,K)
                anumff   (IX,JY)   = anumff(IX,JY) + 1.0

                if(THEDD.lt.0.0) THEDD = THEDD + 2.0*pi

C                if(IX.eq.256.and.JY.eq.256) write(*,'(4I5,3F10.2)'), I,J,IX,JY,THEDD*180./pi,aimagesff(IX,JY),anumff(IX,JY)
C                if(IX.eq.256.and.JY.eq.768) write(*,'(4I5,3F10.2)'), I,J,IX,JY,THEDD*180./pi,aimagesff(IX,JY),anumff(IX,JY)
C                if(IX.eq.768.and.JY.eq.256) write(*,'(4I5,3F10.2)'), I,J,IX,JY,THEDD*180./pi,aimagesff(IX,JY),anumff(IX,JY)
C                if(IX.eq.768.and.JY.eq.768) write(*,'(4I5,3F10.2)'), I,J,IX,JY,THEDD*180./pi,aimagesff(IX,JY),anumff(IX,JY)

                if(IX.eq.384.and.JY.eq.384) write(*,'(4I5,3F10.2)'), I,J,IX,JY,THEDD*180./pi,aimagesff(IX,JY),anumff(IX,JY)
                if(IX.eq.384.and.JY.eq.640) write(*,'(4I5,3F10.2)'), I,J,IX,JY,THEDD*180./pi,aimagesff(IX,JY),anumff(IX,JY)
                if(IX.eq.640.and.JY.eq.384) write(*,'(4I5,3F10.2)'), I,J,IX,JY,THEDD*180./pi,aimagesff(IX,JY),anumff(IX,JY)
                if(IX.eq.640.and.JY.eq.640) write(*,'(4I5,3F10.2)'), I,J,IX,JY,THEDD*180./pi,aimagesff(IX,JY),anumff(IX,JY)

C                if(I.eq.icirx-512-iXi.and.J.eq.iciry-256-iYi) print*,'Q4',I,J,IX,JY,aimagesff(IX,JY),anumff(IX,JY),aimagedif (I,J,K)
C                if(IX.eq.518.and.JY.eq.599) print*,'Q4',I,J,IX,JY,aimagesff(IX,JY),anumff(IX,JY),aimagetst (I,J,K)

C                if(anumff(IX,JY).eq.2) print*, 'val',I,J,alon,alat,PX,PY,PZ,R,THE,IX,JY,aimagesff(IX,JY),anumff(IX,JY)
C                if(BARER.lt.2.0*3.14160.and.BARER.gt.-2.0.and.J.eq.512)
C                if(BARER.lt.3.15.and.BARER.gt.-2.0.and.I.eq.1024) 
C     &       print*,'val barer',I,J,alon,alat,PX,PY,PZ,R,THE,IX,JY,aimagesff(IX,JY),anumff(IX,JY),BARER
              end do
            end do
            anu = 0.0
            do I=1,icirc
              do J=1,icirc
                if(anumff(I,J).ne.0.0) then
                  aimagesf(I,J,K)  = aimagesff(I,J)/anumff(I,J)                        ! This loads aimage with the circular image and zeros
C                  if(aimagesf(I,J,K).gt.0.0.and.anumff(I,J).eq.2) print*, 'all 2',I,J,anumff(I,J),aimagesff(I,J),aimagesf(I,J,K)
                  anu = anu + 1.0
                end if

                iet = 0
                if(iet.eq.1) then                                                   ! This can be used to restrict the radius to less than 360 deg
                  XII = float(I - ipdif)
                  YJJ = float(J - ipdif)
                  IRR = nint(SQRT(XII*XII + YJJ*YJJ))
                  if(IRR.lt.ipdif) then
                    aimagesf(I,J,K) = 800.0
                  else
                    aimagesf(I,J,K) = 0.0
                  end if

                end if
              end do
            end do
            print*,'The number of locations with data are',anu
            cfileoo = '                                                                     '  
            cfileoo = cdfileo(K)
C            if(K.eq.1) then
              print*, ' Before image write icirc = ',K, icirc
              call write_image_grd_p(cfileoo,icirc,icirc,aimagesf(1,1,K))        ! output fisheye difference image
C            end if
           end do
          end if
C
C         Provide a smoothed hammer-aitoff differenced image  
C
C          ihammer=0
          ihammer=1
          print *, 'Before ihammer make', ihammer
          if(ihammer.eq.1) then
            print *,' '
            if(idif.eq.0) write(*,'(A,3I5)'), 'provide a smoothed h-a K image with a base removed',K,icirx/2,iciry/2
            if(idif.eq.1) write(*,'(A,3I5)'), 'provide a smoothed h-a K differenced image maker',K,icirx/2,iciry/2
            print *,' '           

            S10pADU = 0.035
            do K=1,IMGIN-idif
              do I=1,icirx/2                                               ! This zeros the potential Kth smoothed composite h-a image
                do J=1,iciry/2            
                  aimagehass(I,J,K) = 0.0
                  anumhass  (I,J,K) = 0.0
                end do
              end do                                                       ! Load up the smoothed Kth heliographic or ecliptic image
              do I=1,icirx/2
                do J=1,iciry/2
                  aimagehas(I,J) = aimagedif(I,J,K)*S10pADU
                end do
              end do

              cfileoo = '                                                                     '  
              cfileoo = chfileo(K)
              print*, ' Loaded smoothed image to make a s10 h-a write K,icirx/2,iciry/2 = ',K,icirx/2,iciry/2
              call write_image_grd_hp(cfileoo,icirx/2,iciry/2,aimagehas) ! write smoothed hammer-aitoff image difference
C
C             Provide a hammer-aitoff image from the smoothed Kth differenced image
C
              anu = 0.0
              DEGPP = 360.0/float(icirx/2)                                        ! Degrees per pixel
              dDSun = 0.0
              NVAL  = 41

              print*, 'DEGPP,dDSun = ', DEGPP,dDSun,icirx/2,iciry/2
              do I=1,icirx/2
                do J=1,iciry/2
                  iheliographic = 1
                  if(iheliographic.eq.1) DECD  = float(J)*DEGPP - 90.0
                  RADPM = float(I)*DEGPP
                  RADPM = RADPM - 180.0
                  sqtopcc = SQRT(1.0 + cosd(DECD)*cosd(RADPM/2.0))
                  YJ = ((((sind(DECD))/sqtopcc)*90.0) + 90.0 - dDSun)/DEGPP
                  if(YJ.le.0.5)    YJ = 1.0
                  if(YJ.ge.float(iciry/2)+0.5) YJ = float(iciry/2)
                  XI = ((((cosd(DECD)*sind(RADPM/2.0))/sqtopcc)*180.0) + 180.0)/DEGPP
                  if(XI.le.0.05)   XI = 1.0
                  if(XI.ge.float(icirx/2)+0.5) XI = float(icirx/2)
                  IX = nint(XI)
                  if(IX.lt.1) IX = 1
                  if(IX.gt.icirx/2) IX = icirx/2
                  JY = nint(YJ)
                  if(JY.lt.1) JY = 1
                  if(JY.gt.iciry/2) JY = iciry/2
C                  if(J.eq.1.and.I.lt.20)     write(*,'(3(I5,F10.2),F10.2,F12.4)'),I,RADPM,J,DECD,K,XI,YJ,aimagess(IX,JY,K)
C                  if(J.eq.256.and.I.lt.20)   write(*,'(3(I5,F10.2),F10.2,F12.4)'),I,RADPM,J,DECD,K,XI,YJ,aimagess(IX,JY,K)
C                  if(J.eq.512.and.I.lt.20)   write(*,'(3(I5,F10.2),F10.2,F12.4)'),I,RADPM,J,DECD,K,XI,YJ,aimagess(IX,JY,K)
C                  if(J.eq.1.and.I.gt.1004)   write(*,'(3(I5,F10.2),F10.2,F12.4)'),I,RADPM,J,DECD,K,XI,YJ,aimagess(IX,JY,K)
C                  if(J.eq.256.and.I.gt.1004) write(*,'(3(I5,F10.2),F10.2,F12.4)'),I,RADPM,J,DECD,K,XI,YJ,aimagess(IX,JY,K)
C                  if(J.eq.512.and.I.gt.1004) write(*,'(3(I5,F10.2),F10.2,F12.4)'),I,RADPM,J,DECD,K,XI,YJ,aimagess(IX,JY,K)
C                  if(aimagehas(IX,JY).lt.-0.001.or.aimagehas(IX,JY).gt.0.001) then
                    aimagehass(IX,JY,K) = aimagess(IX,JY,K) + aimagehas(I,J)
                    anumhass  (IX,JY,K) = anumhass (IX,JY,K) + 1.0
C                  end if
                end do
              end do
              do I=1,icirx/2
                do J=1,iciry/2
                  if(anumhass(IX,JY,K).gt.-99990.0) then
                    aimagehass(IX,JY,K) = aimagehass(IX,JY,K)/anumhass(IX,JY,K)
                    anu = anu +1.0
                  end if
                end do
              end do
C
C
C             attempt to smooth out the more pattern
              imore = 0
              if(imore.eq.1) then
                print*,'Attempt to smooth the more pattern'
                do I=1,icirx/2
                  do J=1,iciry/2
                    aimagehasss(I,J) = 0.0
                    anumhasss(I,J)   = 0.0
                    do III=I-3,I+3
                      do JJJ=J-3,J+3
                        if(aimagehass(III,JJJ,K).ne.0.0) then
                          aimagehasss(I,J) = aimagehasss(I,J) + aimagehass(III,JJJ,K)
                          anumhasss(I,J)   = anumhasss(I,J) + 1.0
                        end if
                      end do
                    end do
                  end do
                end do
                do I=1,icirx/2
                  do J=1,iciry/2
                    if(anumhasss(I,J).ne.0.0.and.aimagehass(I,J,K).ne.0.0) aimagehasss(I,J) = aimagehasss(I,J)/anumhasss(I,J)
                  end do
                end do
              end if
C             end attempt to smooth out the more pattern
C
              print*,'The number of locations with h-a data are',anu
              cfileoo = '                                                                     '  
              cfileoo = cdfileo(K)
              print*, ' Before h-a image write K,icirx/2,iciry/2 = ',K,icirx/2,iciry/2
              if(imore.ne.1) call write_image_grd_hp(cfileoo,icirx/2,iciry/2,aimagehass(1,1,K)) ! write smoothed hammer-aitoff image difference
              if(imore.eq.1) call write_image_grd_hp(cfileoo,icirx/2,iciry/2,aimagehasss) ! write an even better smoothed hammer-aitoff image difference
C
C             Now make histograms
              do I=1,icirx/2
                do J=1,iciry/2
                  S10Mx = 0.10
                  dval  = S10Mx/(float((NVAL-1)/2))
                  dval2 = dval/2.0
                  do II=1,NVAL 
                    avalm = -S10Mx - dval2
                    avalp = avalm + dval

                    ihammero = 0
                    if(ihammero.eq.1) then                                                            ! make histograms of the output
                      if(I.eq.1.and.j.eq.1) print*,'Histogram of H-A output file' 
                      if(aimagehass(I,j,K).gt.(S10Mx+dval2))  aimval(NVAL,K) = aimval(NVAL,K) + 1.0
                      if(aimagehass(I,j,K).lt.(-S10Mx-dval2)) aimval(1,K) = aimval(1,K) + 1.0
                      if(aimagehass(I,j,K).gt.avalm.and.aimagehass(I,J,K).lt.avalp) aimval(II,K) = aimval(II,K) + 1.0
                    else                                                                              ! make histograms of the input
                      if(I.eq.1.and.j.eq.1) print*,'Histogram of H-A input file' 
                      if(aimagehas(I,j).gt.(S10Mx+dval2))  aimval(NVAL,K) = aimval(NVAL,K) + 1.0
                      if(aimagehas(I,j).lt.(-S10Mx-dval2)) aimval(1,K) = aimval(1,K) + 1.0
                      if(aimagehas(I,j).gt.avalm.and.aimagehas(I,J).lt.avalp) aimval(II,K) = aimval(II,K) + 1.0
                    end if

                    S10v(II) = -S10Mx
                    S10Mx = S10Mx - dval
                  end do
                end do
              end do

            end do
            do II=1,NVAL
              if(II.eq.1) print*,                      '  Histogram analysis         S10 value  DIf 1     Dif 2     Dif 3'
              if(idif.eq.0) write(*,'(A,F7.3,3F10.1)') 'Hammer-Atioff S10 histogram   ',S10v(II),aimval(II,1),aimval(II,2),
     &                                                 aimval(II,3)
              if(idif.eq.1) write(*,'(A,F7.3,3F10.1)') 'Hammer-Atioff S10 histogram   ',S10v(II),aimval(II,1),aimval(II,2)
            end do
          end if

      STOP

        return
      end