C This is FORTRAN source code for the Sigma-CLEAN deconvolution
c (described in an earlier version by Keel 1991 PASP 103, 723)
c The following is a sample parameter file, used to make background
c execution a little friendlier (program expects to find it in
c the file clean.par)
40.1                                Readout noise, ADU (real WFC)
4.1                                 Gain, electrons/ADU (real WFC)
n6240vstar.txt                       PSF file
21                                  X dimension of PSF file
21                                  Y dimension of PSF file
n6240vx.txt                         Data file to be deconvolved
30                                  X dimension of image
30                                  Y dimension of image
0.01                                Loop gain
10000                               Number of iterations
Y                                   Use image as noise model?
n6240vx.txt                         Noise model image file (same size as image)
3.5                                 Cutoff in image sigma to stop cleaning


	PROGRAM HSTCLN
C	  APPLIES CLEAN ALGORITHM TO OPTICAL DATA
C	  SMOOTHS AFTER EACH ITERATION TO SUPPRESS NOISE EFFECTS
C	  VERSION TO PRINT ONLY SELECTED CLEAN COMPONENTS
c	modifications 2/89 to read/write IRAF-compatible text format
c	further modifications started 9/10/90 for CRAY use on WFC data
c	selects most significant (not largest) peak for clean component
c	More improvements: input from a parameter file HSTCLN.PAR
c	option for different noise-model image than input data
c	address image arrays as one-dimensional instead of 2D
c	insert a cutoff in sigma (like 3.0) as sigmacut
	CHARACTER*35 FILEIN,FDATA,EXT,FOUT,FCC
	character*1 letter,bigy,littley
	data bigy,littley/'Y','y'/
c	  Dimensions for image arrays: data to 256x256 or same area
c	  PSF to 80x80 or same area
c	DIMENSION DATA(65536),SMTH(65536),CLN(65536)
c	dimension sigma0(65536),snr(65536)
c	DIMENSION PSF(6400)
	DIMENSION DATA(256,256),SMTH(256,256),CLN(256,256)
	dimension sigma0(256,256),snr(256,256)
	DIMENSION PSF(80,80)
	logical n100,stopflag

C	SET UP ARRAY DIMENSIONS TO BE READ BY PEAK
	NPSF=80
	MPSF=80
	NDATA=256
	MDATA=256
	open (unit=7,file='clean.par',form='formatted',status='old',
     1      access='sequential',carriagecontrol='list')
	open (unit=4,file='clean.log',form='formatted',status='new',
     1      access='sequential',carriagecontrol='list')
c	   set up parameters of noise model
c	   ron = effective readout noise IN ELECTRONS
c	   gain = electrons/ADU so Poisson noise drops by gain**1/2
	read (7,*) ron
	read (7,*) ccdgain
	sqrtg=sqrt(ccdgain)
	ccmin=1.e8
	WRITE (4,900)
  900   FORMAT (2X,' HSTCLN RUN LOG '/
     1   '  PSF FILE: ')
	read (7,901) filein
  901   format (a35)
	read (7,*) nx
	read (7,*) ny
	OPEN (UNIT=8,FILE=FILEIN,FORM='FORMATTED',ACCESS='SEQUENTIAL',
     1    STATUS='OLD')
	IPSF=NX
	JPSF=NY
	 READ (8,*) ((psf(i,j),i=1,ipsf),j=1,jpsf)
	CLOSE(8)
	WRITE (4,903) filein,ipsf,jpsf
  903   FORMAT (2x, a35,' FOUND AND READ IN '/2x,'NX=',i3,'  NY=',i3)
	read (7,901) FDATA
	read (7,*) nx
	read (7,*) ny
	OPEN (UNIT=9,FILE=FDATA,form='formatted',ACCESS='sequential',
     1   STATUS='OLD')
	 READ (9,*) ((data(i,j),i=1,nx),j=1,ny)
	CLOSE (9)

	CALL PEAK (DATA,NDATA,MDATA,NX,NY,IPK,JPK,DMAX)
	WRITE (4,904) fdata,NX,NY,DMAX,IPK,JPK
  904   FORMAT ('  DATA FILE:',a35/'  NX=',I3,' NY=',I3/
     1    2X,'MAXIMUM=',E12.6,' AT I=',I3,' J=',I3)
c	scale input data to electrons (rather than ADU) for consistency
c	with exact noise level
	do i=1,nx
	do j=1,ny
		data(i,j)=data(i,j)*ccdgain
	end do
	end do
c	   set up initial sigma image; sigma=sqrt(ron**2+image/gain)
c	   use the same image for cleaning and noise model?
	READ (7,*) gain
	read (7,*) NLOOPS
	read (7,905) letter
  905   format (a1)
	if (letter.ne.bigy.and.letter.ne.littley) go to 55
c	   case for sigma image the same as input image
c	   skip name of noise-model image first
	read (7,905) letter
	do 50 i=1,nx
	do 45 j=1,ny
	sigma0(i,j)=sqrt((ron)**2+data(i,j))
   45   continue
   50   continue
	go to 58
   55   read (7,901) filein
	open (unit=9,file=filein,form='formatted',status='old',
     1    access='sequential',carriagecontrol='list')
	read (9,*) ((sigma0(i,j),i=1,nx),j=1,ny)
	close(9)
	write (4,907) filein
  907   format (2x,'Noise model from ',a35)
	do i=1,nx
	do j=1,ny
		sigma0(i,j)=sqrt((ron)**2+ccdgain*sigma0(i,j))
	end do
	end do
C	FIND PEAK OF PSF AND ITS POSITION; NORMALIZE
   58   CALL PEAK (PSF,NPSF,MPSF,IPSF,JPSF,IPSPK,JPSPK,PSFMAX)
	PSUM=0.0
	DO 70 J=1,JPSF
	DO 60 I=1,IPSF
	  PSF(I,J)=PSF(I,J)/PSFMAX
	  PSUM=PSUM+PSF(I,J)
   60   CONTINUE
   70   CONTINUE
	read (7,*) sigmacut
	WRITE (4,908) PSFMAX,IPSPK,JPSPK,IPSF,JPSF,PSUM
  908   FORMAT ('  PSF PEAK = ',F10.4,' IPK=',I3,' JPK=',I3,' NX=',I4,
     1   ' NY=',I4,/' TOTAL FLUX=',E12.6)
	FTOTAL=0.0
	MI1=1
	MI2=NX
	MJ1=1
	MJ2=NY
	NLINE=1
	n100=.false.

C	MAIN CLEANING LOOP
	DO 500 NITR=1,NLOOPS
		do jsig=1,ny
		do isig=1,nx
		snr(isig,jsig)=data(isig,jsig)/sigma0(isig,jsig)
 		end do
		end do
	CALL PEAK (snr,NDATA,MDATA,NX,NY,IPK,JPK,VMAX)
	if ((abs(vmax)).lt.sigmacut) go to 550
	vmax=data(ipk,jpk)
	CCFLUX=GAIN*VMAX*PSUM
	CLN(IPK,JPK)=CLN(IPK,JPK)+CCFLUX
	SFLUX=GAIN*VMAX
	FTOTAL=FTOTAL+CCFLUX
	IF (NLINE.EQ.100.OR.NITR.EQ.1)
     1  WRITE (4,925) NITR,IPK,JPK,CCFLUX,FTOTAL
c	   if abs(ccflux) is going up, terminate clean
	if (nline.eq.100) n100=.true.
	if (.not.n100) go to 395
        if (abs(ccflux).lt.1.e-4) stopflag=.true.
c	if (abs(ccflux).gt.3.0*ccmin) stopflag=.true.
	if (stopflag) go to 550
	if (abs(ccflux).lt.ccmin) ccmin=abs(ccflux)
  395	NLINE=NLINE+1
	IF (NLINE.EQ.101) NLINE=1
C	  IPK,JPK,CCFLUX IS NITR'TH CLEAN COMPONENT
C	  SUBTRACT SCALED PSF
	  DO 450 J=1,JPSF
	  DO 440 I=1,IPSF
	  INDEX=IPK-IPSPK+I
	  JNDEX=JPK-JPSPK+J
	  IF (JNDEX.LT.1.OR.JNDEX.GT.NY) GO TO 450
	  IF (INDEX.LT.1.OR.INDEX.GT.NX) GO TO 440
	  DATA(INDEX,JNDEX)=DATA(INDEX,JNDEX)-SFLUX*PSF(I,J)
  440     CONTINUE
  450     CONTINUE
  500     CONTINUE

C	  WRITE RESIDUAL MAP TO FILE
c	convert back to ADU from electrons
  550	do i=1,nx
	do j=1,ny
		data(i,j)=data(i,j)/ccdgain
	end do
	end do
        OPEN (UNIT=8,FILE='clean.res',ACCESS='sequential',
     1    STATUS='NEW',form='formatted',carriagecontrol='list')
	WRITE (8,997) ((data(i,j),i=1,nx),j=1,ny)
  997   format (1x,f11.2)
	CLOSE(8)
C	  WRITE CLEAN MODEL TO FILE, in ADU again
	do i=1,nx
	do j=1,ny
		cln(i,j)=cln(i,j)/ccdgain
	end do
	end do
  690   OPEN (UNIT=3,FILE='clean.dat',form='formatted',
     1   ACCESS='sequential',STATUS='NEW',carriagecontrol='list')
	WRITE (3,997) ((cln(i,j),i=1,nx),j=1,ny)
	CLOSE (3)
  925   FORMAT (1x,'Iteration ',I7,' I=',I3,' J=',I3,' Flux=',
     1     E13.6,' Flux sum=',E13.6)
	STOP
	END

	SUBROUTINE PEAK (ARRAY,IDIM,JDIM,NX,NY,IPK,JPK,VMAX)
C	  FINDS PEAK (ABSOLUTE VALUE) IN ARRAY
C	  ARRAY IS IDIM*JDIM; NX*NY ACTUALLY FILLED
	DIMENSION ARRAY(IDIM,JDIM)
	VMAX=0.0
	AVMAX=0.0
	DO 100 J=1,NY
	  DO 50 I=1,NX
	  IF (ABS(ARRAY(I,J)).LT.AVMAX) GO TO 50
	  VMAX=ARRAY(I,J)
	  AVMAX=ABS(VMAX)
	  IPK=I
	  JPK=J
   50     CONTINUE
  100   CONTINUE
	RETURN
	END