C program to read in a .mag file and a specified magnitude and bin width,
C   and write out to fort.7 the number of stars found in this bin and
C   the median mags in the other colors

	real data(22), sum(50000,6), fmed(6)
	integer nsum(6)
	character file*64

	print '(1x,''Enter file name to get fake mags on: ''$)'
	read '(a)', file
	open(1,file=file,status='old')
        print '(1x,''Enter desired magnitude, bin width: ''$)'
	read *, amag, abin

	read(1,*)
	read(1,*)
	read(1,*) nfilt
        nfilt = abs(nfilt)

        a1 = amag - abin/2.
        a2 = amag + abin/2.

	n = 0
	do i=1,nfilt
          nsum(i) = 0
        end do
1	read(1,*,end=99) (data(i),i=1,3+2*nfilt)
        if (data(4) .ge. a1 .and. data(4) .lt. a2) then
          n = n + 1
	if (n .gt. 50000) then
          print *, 'n too large: ', n
          pause
        end if
          do i=2,nfilt
            if (data(2+i*2) .lt. 90) then
              nsum(i) = nsum(i) + 1
              sum(nsum(i),i) = data(4)-data(2+i*2)
            end if
          end do
        end if
	goto 1

99	continue

  	do i=2,nfilt
          if (nsum(i) .gt. 10) then
	    call median(sum(1,i),nsum(i),fmed(i))
          else
            fmed(i) = amag
          end if
  	end do

	write(7,*) n
	write(7,*) (amag,amag-fmed(i),i=2,nfilt)

	stop
	end

	
        SUBROUTINE MEDIAN(X,NPTS,FMEDIA)

C       This routine finds the median in an array of numbers.

C       It uses a quicksort algorithm from the text :
C       'Algorithms + Data Structures = Programs' by Wirth (1976) page 84.
C       In the worst case it completely sorts the input array.
C       (Time is proportional to n * log n.)
C       Sorting is done on the indices of the array, so as not to
C       disturb its original ordering.

C       Input:  X       The array of points to be examined
C               NPTS    The length of the array
C               KEY     The keys to be partially sorted

C       Output: FMEDIA  The median value

C       Authors:Mike Fich, Berkeley, April 1981.
C               B.F.

        REAL*4 X(NPTS)
        INTEGER KEY(50000)

C       Initialize the keys.

        DO 8701 I = 1, NPTS
                KEY(I) = I
8701    CONTINUE

        ILEFT   = 1
        IRITE   = NPTS
        MIDDLE  = (NPTS+1)/2

C       Find the middle element, sort everything by exchanging everything
C       bigger on the left with everything smaller on the right end of the
C       array.

8702         IF (ILEFT .LT. IRITE) THEN
                XMID    = X(KEY(MIDDLE))

C       INDEX1 moves from left to right, INDEX2 moves from right to left
C       INDEX1 points to something bigger than the key being sorted on.
C       INDEX2 points to something smaller.

                INDEX1  = ILEFT
                INDEX2  = IRITE
8703                 IF (INDEX1 .LE. INDEX2) THEN
                        IF (X(KEY(INDEX1)) .LT. XMID) THEN
8704                                 IF (X(KEY(INDEX1)) .LT. XMID) THEN
                                        INDEX1  = INDEX1+1
                                GO TO 8704
                                END IF
                        END IF

                        IF (X(KEY(INDEX2)) .GT. XMID) THEN
8705                                 IF (X(KEY(INDEX2)) .GT. XMID) THEN
                                        INDEX2  = INDEX2-1
                                GO TO 8705
                                END IF
                        END IF

C       If INDEX1 is still to the left of INDEX2 then exchange the
C       keys found out of place and increment pointers.

                        IF (INDEX1 .LE. INDEX2) THEN
                                KTEMP   = KEY(INDEX1)
                                KEY(INDEX1)     = KEY(INDEX2)
                                KEY(INDEX2)     = KTEMP
                                INDEX1  = INDEX1+1
                                INDEX2  = INDEX2-1
                        END IF

C       If not sorted for initial key go back.

                GO TO 8703
                END IF

C       Choose new limits for sorting remaining block.

                IF (INDEX2 .LT. MIDDLE) ILEFT = INDEX1
                IF (MIDDLE .LT. INDEX1) IRITE = INDEX2
        GO TO 8702
        END IF

        FMEDIA  = X(KEY(MIDDLE))

        RETURN
        END