C+
C NAME:
C	read_grd_image.f
C PURPOSE:
C	Read the ASHI image times and te centroid of the star from a file, and provide them to the main program as Julian dates
C     and x y pixel positions
C CATEGORY:
C	I/O
C CALLING SEQUENCE:
C	call read_grd_image(itype,cfilei,cfileo,NXY,dmax,cimghead,aimagexy,imagexy)
C INPUTS:
C     itype             integer     Type of file input and output
C                                      0 = read in and write out a grd file with <= 5 integer digits 
C                                      1 = read in a grd file with <= 5 integer digits; write out a grd file with F8.1 floating values
C                                      2 = read in a grd file with <= 5 integer digits; output a file with F8.1 floating values
C                                      3 = read in and write out a grd file with F8.1 floating values 
C                                      4 = input a grd file with <= 5 integer digits; output a file <= 5 integer digit values
C                                      5 = input a grd file with <= 5 integer digits; output a file with F8.1 floating values
C                                      6 = input a grd file with <= 5 integer digits; write out a grd file with <= 5 integer digit values
C                                      7 = input a grd file with <= 5 integer digits; write out a grd file with F8.1 floating values
C                                      8 = input a grd file with F8.1 floating values; write out a grd file with F8.1 floating values
C                                      9 = input a grd file with F8.1 floating values; write out a grd file with <= 5 integer digit values
C                                     10 = input a grd file with F8.1 floating values; output a file with F8.1 floating values
C                                     11 = input a grd file with F8.1 floating values; output a file with F8.1 floating values
C                                     12 = read in a grd file with F8.1 floating values and output the file - 1 value input short each line
C                                     13 = read in a grd file with F8.2 floating values and output the file
C                                     14 = read in and write out a grd file with F9.2 floating values 
C                                     15 = input and write out a grd file with F9.2 floating values 
C
C	cfilei            character	Name of the ASHI image grd input file
C	cfileo            character	Name of the ASHI image grd output file
C	NXY	            integer	Maximum number of image brightness values possible for both x and y in the file (generally 2048)
C     cimghead(5)       character   Grd file 5 line header
C     aimagexy(NXY,NXY)	real	      Floating grd image 
C     imagexy (NXY,NXY) integer	Fixed grd image     
C
C OUTPUTS:
C	aimagexy(NXY,NXY)	real	      Floating grd image
C     imagexy (NXY,NXY) integer	Fixed grd image  
C
C FUNCTIONS/SUBROUTINES:
C
C MODIFICATION HISTORY:
C	November-2022, Bernard Jackson (UCSD)
C-
	subroutine read_grd_image(itype,cfilei,cfileo,NXY,dmax,cimghead,aimagexy,imagexy)

        parameter       (ihead = 5)
        character 	cfilei*56
        character 	cfileo*57
        character       cimage1*9
        character       cimghead(ihead)*9

        real            aimagexy(NXY,NXY)
        real            aim     (NXY,NXY)
        real            aimt     (NXY*NXY)

        integer         imagexy(NXY,NXY)        
        
          distcx  = 1126.5                            ! center of the image in X as determined by Matthew 
          distcy  = 1000.5                            ! center of the image in Y as determined by Matthew

          print*,' '
          print*, 'Into read_grd_image.f itype =',itype
C          print*, cfilei

          if(itype.le.3.or.itype.eq.12.or.itype.eq.13.or.itype.eq.14) then
	      open (13, file=cfilei,status='old',recl=120,access='sequential',form='formatted',iostat=iReadashi)
	      if(iReadashi.ne.0) then
              close(13)
		  print *, 'There is no ',cfilei,' file in this subdirectory' 
		  stop
	      end if

            NXYM = NXY
            J2   = 1
            if(itype.eq.12) then
              NXYM = NXY-1
              J2   = 2
            end if
            do J=J2,NXYM
              if(J.eq.J2) then
                do K=1,5
                  read (13,'(A)',iostat=iReadashi) cimage1
                  if(K.eq.1) then
                    if(iReadashi.ne.0.or.cimage1.ne.'DSAA') then
                      close(13)
                      print*,'This is not a grd file - stop!'
                      stop
                    else
                      cimghead(K) = cimage1
                    end if
                  else
                    cimghead(K) = cimage1
                  end if
                end do
              end if              
              if(itype.le.2)  read (13,'(10F6.0)',iostat=iReadashi) (aimagexy(I,J),I=1,NXY)
              if(itype.eq.3)  read (13,'(10F8.2)',iostat=iReadashi) (aimagexy(I,J),I=1,NXY)
              if(itype.eq.12) read (13,'(10F8.2)',iostat=iReadashi) (aim(I,J),I=2,NXY)
              if(itype.eq.13) read (13,'(10F8.2)',iostat=iReadashi) (aimagexy(I,J),I=2,NXY)
              if(itype.eq.14) read (13,'(10F9.2)',iostat=iReadashi) (aimagexy(I,J),I=2,NXY)
	        if(iReadashi.ne.0) then
                close(13)
                print *, 'This file is corrupt - stop'
                stop
              end if
C              write (*,'(10F8.1)') (aimagexy(I,J),I=1,NXY)
	      end do
C            print*,' '
            print*, 'File just read: ',cfilei
          else
            if(itype.ge.4.and.itype.le.7) then
              do J=1,NXY
                do I=1,NXY
                  aimagexy(I,J) = float(imagexy(I,J))
                end do
              end do
              print*,' '
              print*, 'There was a floating file input and output'
            end if
          end if

          if(itype.eq.12) then
            IJ  = 1
            IJ1 = 0
            do J=1,NXYM
              do I=1,NXYM
                II = I+1
                JJ = J+1
                IIJJ = II+JJ
                IB = 0
                JB = 0
C                if(J.eq.NXYM) JB = 1
                if(I.eq.NXYM) IB = 2
                IJ   = IJ + 1 + IB + JB
                IJ1  = IJ1 + 1
                aimt(IJ)   = aim(II,JJ)
                aimagexy(I,J) = aimt(IJ1)
              end do
            end do 
          end if


          aimagel  = 0.0
          aimages  = 2000.0
          aimagess = 2000.0
          dinner   = 700.0
          nsmall   = 0
          do J=1,NXY
            do I=1,NXY
              distxy = sqrt((float(I)-distcx)**2 + (float(J)-distcy)**2)
              if(distxy.gt.dmax) aimagexy(I,J) = 0.0
              if(distxy.lt.dinner.and.aimagexy(I,J).ne.0.0.and.aimagexy(I,J).lt.aimages) aimages = aimagexy(I,J)
              if(distxy.lt.dinner.and.aimagexy(I,J).eq.0.0) nsmall = nsmall + 1
              if(aimagexy(I,J).gt.aimagel) aimagel = aimagexy(I,J)
              if(aimagexy(I,J).ne.0.0.and.aimagexy(I,J).lt.aimagess) aimagess = aimagexy(I,J)
            end do
          end do
          if(dmax.lt.2500.0) then
            print*,' '
            print*,'This is map rotated with the north pole at the top.'
            write(*,'(A,F10.2,A,F6.1)') ' The largest image value is', aimagel,'; The greatest radius mapped =',dmax
            write(*,'(A,F10.2,A,F6.1)') ' The smallest non-zero image value is', aimages,' inside the radius =',dinner
            write(*,'(A,I8,A,F6.1)')    ' There are ', nsmall,' zero values inside the radius =',dinner
          else
            print*,' '
            print*,'This is an image that fills to the edge of the 2048 pixel image'
            write(*,'(A,F10.2,A,F6.1)') ' The largest image value is', aimagel,'; The image is mapped to its edge.'
            write(*,'(A,F10.2,A,F6.1)') ' The smallest non-zero image value is',aimagess            
          end if
          if(itype.eq.0.or.itype.eq.1.or.itype.eq.3.or.itype.eq.6.or.itype.eq.7.or.itype.eq.8.or.itype.eq.9
     &       .or.itype.eq.15) then
            print*,' '
            print*, 'write out the grd file just read'
	      open (13, file=cfileo,status='new',recl=120,access='sequential',form='formatted',iostat=iReadashi)
            if(iReadashi.ne.0) then
              close(13)
              print*,'Could not open the ',cfileo,' file - stop!'
              stop
            end if
            do J=1,NXY         
              if(J.eq.1) then
                do K=1,5
                  write(13,'(A)') cimghead(K)
                end do
              end if
              if(itype.eq.0.or.itype.eq.6.or.itype.eq.9) then
                write(13,'(10I6)') (nint(aimagexy(I,J)),I=1,NXY)
C                write(*, '(10I6)') (nint(aimagexy(I,J)),I=1,NXY)
                if(J.eq.NXY) then
                  print*, cfileo,' written with <= 5 integer digit values'  
                end if
              end if
              if(itype.eq.1.or.itype.eq.3.or.itype.eq.7.or.itype.eq.8) then
                if(aimagexy(I,J).gt.99990.) aimagexy(I,J) = 99990.
                write(13,'(10F8.2)') (aimagexy(I,J),I=1,NXY)
C                write(*, '(10F8.1)') (aimagexy(I,J),I=1,NXY)
                if(J.eq.NXY) then
                  print*, cfileo,' was written with F8.2 floating values' 
                end if
              end if
              if(itype.eq.15) then
                if(aimagexy(I,J).gt.99990.) aimagexy(I,J) = 99990.
                write(13,'(10F9.2)') (aimagexy(I,J),I=1,NXY)
C                write(*, '(10F8.1)') (aimagexy(I,J),I=1,NXY)
                if(J.eq.NXY) then
                  print*, cfileo,' was written with F9.2 floating values' 
                end if
              end if
	      end do
          else  
C            print*,' '
            print*, 'This was an output rather than a write of whatever file was input'            
          end if
C          print*,' '
          print*, 'Out of read_grd_image.f'
	  return
	end