program bincombine

	parameter (maxfilt=2)
        character file*64, header*2880
	real mag(4), x(4), y(4), cnts(maxfilt), zz(maxfilt)
	integer id(4), nstar(4)
	parameter (imin = -200, imax = 200)
	integer*2 count(imin:imax,imin:imax)
	real rtot, rmatch, rcount(imin:imax,imin:imax)

C  Reads combine files for each chip where combine file
C    combines
C       1. fake star list
C       2. observed star list on fake frame
C       3. observed star list on original frame
C  Input is:
C       1. number observed stars on chip1
C       2. name of chip1 combine file
C       3. number observed stars on chip2
C       4. name of chip2 combine file
c       ...
C       8. name of chip4 combine file

C  Output is a historgram giving number of stars per
C    magnitude difference. Each chip is weighted by
C    the number of observed stars (currently magnitude
C    independent) on the chip. Also output the
C    total number of stars in the histogram and the
C    number of bins which are output

	read *, nfilt
	do i=1,nfilt
	   read *, cnts(i)
           cnts(i) = -2.5*log10(cnts(i)) + 25
	end do

	read '(a)', file
        l = index(file,' ') - 1
        istat = openc(ifile,file(1:l)//'.img'//char(0),1)

	open(7,file=file(1:l)//'.hdr',status='unknown')
        header = 'END '
        call lheadset('SIMPLE',.true.,header)
	call inheadset('BITPIX',16,header)
	call inheadset('BITPIX',16,header)
	call inheadset('NAXIS',2,header)
	call inheadset('NAXIS1',imax-imin+1,header)
	call inheadset('NAXIS2',imax-imin+1,header)
	call inheadset('CRVAL1',imin,header)
	call inheadset('CRVAL2',imin,header)
	call fheadset('CDELT1',1.d0,header)
	call fheadset('CDELT2',1.d0,header)
	call cheadset('OBJECT',file(1:l)//'  error distribution',header)
	call fheadset('MAG1',dble(cnts(1)),header)
	call fheadset('MAG2',dble(cnts(2)),header)

	jmin = imax
	jmax = imin
	rtot = 0
        rmatch = 0
        do i=imin,imax
          do j=imin,imax
            rcount(i,j) = 0.
          end do
	end do

	do ichip=1,4

C       Get number of observed stars on this chip for normalization
          read *, nstar(ichip), (zz(i),i=1,nfilt)

C       Get combine file name and open file
C        print '(1x,''Enter cmb file name: ''$)'
	  read '(a)', file
          l = index(file,' ') - 1
	  open(1,file=file(1:l)//'.cmb',status='old')
	  do i=1,3
	    read(1,*)
          end do

C        Loop over all the stars and bin the matched fake
C         stars, count all of the fake stars
1	  read(1,*,end=99) (id(i), mag(i), z, z, z, i=1,4)
	  read(1,*) (x(i),y(i),z,z,z,i=1,4)

C  Keep track of the total number of fake stars for total completeness
	  if (mag(3) .gt. 0) rtot = rtot + float(nstar(ichip))/nstar(1)

C  If we've matched a fake star, see if it's a good match.
C     Accept it if:
C       1. it doesn't match any star on the original list with
C           2 pixels
C       2. It does match, but the measured mag is closer to
C           the fake than to the original

	  if (mag(1) .gt. 0. .and. mag(3) .gt. 0) then

            if (mag(4) .eq. 0. .or. 
     &       abs(mag(1) - mag(3)) .lt. abs(mag(1)-mag(4))) then

               ibin = nint( (cnts(1)+zz(1)-mag(1))/0.01 )
               jbin = nint( (cnts(2)+zz(2)-mag(2))/0.01 )
               if (ibin .ge. imin .and. ibin .le. imax .and.
     &             jbin .ge. imin .and. jbin .le. imax) then
     	         rcount(ibin,jbin) = rcount(ibin,jbin) + 
     &                float(nstar(ichip))/nstar(1) 
	         jmin = min(jmin, ibin)
	         jmax = max(jmax, ibin)
                 rmatch = rmatch + float(nstar(ichip))/nstar(1)
C	        print *, ichip, id(1), (mag(i),i=1,4), (cnts(i),i=1,2),
C     &               ibin, jbin
	       else
C	        print *, id(1), (mag(i),i=1,4), (cnts(i),i=1,2)
               end if

            end if
          end if

	  goto 1

99	  continue
        end do

        ntot = nint(rtot)
        nmatch = nint(rmatch)
	call inheadset('NTOT',ntot,header)
	call inheadset('NMATCH',nmatch,header)


	ctot = 0
	ctot2 = 0	
	kernel = -12345
	do i=imin,imax
          do j=imin,imax
            count(j,i) = nint(rcount(j,i)+ran1(kernel)-0.5)
            ctot = ctot + count(j,i)
            ctot2 = ctot2 + rcount(j,i)
          end do
        end do
	call fheadset('CTOT',dble(ctot),header)
	call fheadset('CTOT2',dble(ctot2),header)

C  Write out the header
        LS = 1
        LE = LS + 79
        LH = LEN(HEADER)
51      CONTINUE
        IF (HEADER(LS:LE) .NE. ' ')
     &     WRITE(7,'(A80)',IOSTAT=IERR) HEADER(LS:LE)
        LS = LS + 80
        LE = LE + 80
        IF (LE .LE. LH .AND. HEADER(LS-80:LS+3-80) .NE. 'END ') GOTO 51

C  Close the header file 
        CLOSE(7)

C  Write out the image
	istat = writeint(ifile,count,2*(imax-imin+1)**2)
	istat = closec(ifile)

C  Output the total number of fake stars and the number of
C    fake stars found as a function of magnitude difference
	print *, ntot, jmax-jmin+1, nmatch
C	do ibin = jmin, jmax
C          diff = ibin*0.01
C	  print *, diff, count(ibin)
C        end do

	stop
	end