C	This program reads eq, np or sp*.grd files and
C	makes *.grd files of a combination for 
C	Viewing and further analysis
C
C
	program allskymap
	implicit none
	character*10 head(5)
	character*46 infile,infile1,infile0
	character*1 c,e,n,s,a(21),itype, three
	integer*4 i,j,k,imax,jmax,nskycombo(3600,1200),kmin,kmax,indx(19),karg,kkarg,imap
	integer*4 ixstars
	real*4 frame(3600,1200,19),skycombo(3600,1200),skymed(3600,1200),array(19)
	data infile0 /'c:\EngineeringXStars\cNXXZ_yyyy_doy_hhmmss.grd'/
	data infile1 /'c:\EngineeringOrbits\cNXXZ_yyyy_doy_hhmmss.grd'/
	data c,e,n,s,three /'c','e','n','s','3'/
	data skycombo /4320000*0./
	data nskycombo /4320000*0/
	data imax,jmax,kmax /3600,1200,19/
	data imap /1/	!imap=1 means "yes, output maps"
	data ixstars /0/ !ixstars=0 means: yes, do removed-star files in EngineeringXStars
c
	if(ixstars.eq.0)infile=infile0
	if(ixstars.eq.1)infile=infile1 	
c
	if(ixstars.eq.1)open (20,file='c:\EngineeringOrbits\filelist2eq3.txt',readonly)	
	if(ixstars.eq.0)open (20,file='c:\EngineeringXstars\filelist2sp3.txt',readonly)	
	do k=1,kmax
	read (20,20)(a(i),i=1,21)
20	format(39x,21a1)
c	c3eq3_2003_179_103910.grd		!bad camera 3 map
c	if(k.eq.3.and.a(2).eq.three)go to 96	!throw away cam 3 map: at wrong place in sky
	if(k.eq.1)then
	itype=a(3)
	if(itype.ne.e)then
		imax=800
		jmax=800
	endif
	endif
	if(a(1).ne.c)go to 98
	write(infile(22:42),'(21a1)')(a(i),i=1,21)
c
	print 15,k,a
15	format(i6,' Reading ',21a1)	
	open (10, file = infile, readonly)
	read (10,1)head
1	format(a10)	
	do j=1,jmax			!read in map
	read (10,10) (frame(i,j,k),i=1,imax)
10	format(20f8.1)	
	enddo
	close(10)
c
	do i=1,imax
	do j=1,jmax
c		if(frame(i,j,k).gt.60000.)frame(i,j,k)=0.	!remove the Sun marker
		if(frame(i,j,k).gt.0.)then
			skycombo(i,j)=skycombo(i,j)+frame(i,j,k)
			nskycombo(i,j)=nskycombo(i,j)+1
		endif
c		
	enddo
	enddo
c
96	enddo	!end of k loop, readin
c
c	Begin making maps
c
98	continue
	if(imap.ne.1)go to 95	
c
c	First calculate/output mean map
c
	do j=1,jmax
	do i=1,imax
		if(nskycombo(i,j).ne.0)skycombo(i,j)=skycombo(i,j)/float(nskycombo(i,j))
	enddo
	enddo
c			
	if(itype.eq.e)then
		if(ixstars.eq.1)open (11,file='c:\EngineeringOrbits\eq_allskymean.grd')
		if(ixstars.eq.0)open (11,file='c:\EngineeringXStars\eq_allskymean.grd')
		write(11,11)
11		format('DSAA',/,'3600 1200',/,'0 360',/,'-60 60',/,'0 1000')
	elseif(itype.eq.n)then
		if(ixstars.eq.1)open (11,file='c:\EngineeringOrbits\np_allskymean.grd')
		if(ixstars.eq.0)open (11,file='c:\EngineeringXStars\np_allskymean.grd')
		write(11,12)
12		format('DSAA',/, '800 800',/,'-40 40',/,'-40 40',/,'0 1000')
	elseif(itype.eq.s)then
		if(ixstars.eq.1)open (11,file='c:\EngineeringOrbits\sp_allskymean.grd')
		if(ixstars.eq.0)open (11,file='c:\EngineeringXStars\sp_allskymean.grd')
		write(11,12)
	endif
	do j=1,jmax
		write(11,10)(skycombo(i,j),i=1,imax)
	enddo
	close(11)
c
c	Now find/output minimum map
c	
	do j=1,jmax
	do i=1,imax
		kmin=0
		do k=1,kmax
			if(frame(i,j,k).gt.0.)then
			    if(frame(i,j,k).le.skycombo(i,j))then
				kmin=kmin+1
				skycombo(i,j)=frame(i,j,k)
			    endif
			endif
		enddo
		if(kmin.eq.0)then
c		print *,i,j,skycombo(i,j),((frame(i,j,k)),k=1,11)
		endif
	enddo
	enddo
c	
	if(itype.eq.e)then
		if(ixstars.eq.1)open (11,file='c:\EngineeringOrbits\eq_allskymin.grd')
		if(ixstars.eq.0)open (11,file='c:\EngineeringXStars\eq_allskymin.grd')
		write(11,11)
	elseif(itype.eq.n)then
		if(ixstars.eq.1)open (11,file='c:\EngineeringOrbits\np_allskymin.grd')
		if(ixstars.eq.0)open (11,file='c:\EngineeringXStars\np_allskymin.grd')
		write(11,12)
	elseif(itype.eq.s)then
		if(ixstars.eq.1)open (11,file='c:\EngineeringOrbits\sp_allskymin.grd')
		if(ixstars.eq.0)open (11,file='c:\EngineeringXStars\sp_allskymin.grd')
		write(11,12)
	endif
	do j=1,jmax
		write(11,10)(skycombo(i,j),i=1,imax)
	enddo
	close(11)
c
c	Now find/output median map
c	
95	continue	!jump to here if no maps are to be output
c	
	do j=1,jmax
	do i=1,imax
		kmin=0
		skycombo(i,j)=0.
		do k=1,kmax
			array(k)=0.
			if(frame(i,j,k).gt.0.)then
			    kmin=kmin+1
			    array(k)=frame(i,j,k)
			endif
		enddo
		if(kmin.gt.1)then
			karg=kmax-kmin/2
			call indexx(kmax,array,indx)
			skycombo(i,j)=array(indx(karg))
		endif
		if(kmin.le.1)then
c		print *,i,j,skycombo(i,j),kmin
		endif
		if(skycombo(i,j).gt.99999..or.skycombo(i,j).lt.-9999.)skycombo(i,j)=0.
		skymed(i,j)=skycombo(i,j)
	enddo
	enddo
	if(imap.ne.1)go to 94
c	
	if(itype.eq.e)then
		if(ixstars.eq.1)open (11,file='c:\EngineeringOrbits\eq_allskymed.grd')
		if(ixstars.eq.0)open (11,file='c:\EngineeringXStars\eq_allskymed.grd')
		write(11,11)
	elseif(itype.eq.n)then
		if(ixstars.eq.1)open (11,file='c:\EngineeringOrbits\np_allskymed.grd')
		if(ixstars.eq.0)open (11,file='c:\EngineeringXStars\np_allskymed.grd')
		write(11,12)
	elseif(itype.eq.s)then
		if(ixstars.eq.1)open (11,file='c:\EngineeringOrbits\sp_allskymed.grd')
		if(ixstars.eq.0)open (11,file='c:\EngineeringXStars\sp_allskymed.grd')
		write(11,12)
	endif
	do j=1,jmax
		write(11,10)(skycombo(i,j),i=1,imax)
	enddo
	close(11)
c
c	Now find/output sum-all-below-median map (note, some of this repeats above section, but tricky...)
c
94	continue	!jump to here if no maps are to be output
c	
	do j=1,jmax
	do i=1,imax
		kmin=0
		skycombo(i,j)=0.
		do k=1,kmax
			array(k)=0.
			if(frame(i,j,k).gt.0.)then
			    kmin=kmin+1
			    array(k)=frame(i,j,k)
			endif
		enddo
		if(kmin.gt.1)then
			karg=0
			kkarg=0
			call indexx(kmax,array,indx)
			do k=1,kmax
				if(array(indx(k)).gt.0.)then
					karg=karg+1
c
c					Here limit low contributers to sky map in faint sky, mitigates some overdone glare removal
c					
					if((array(indx(k)).gt.0.25*skymed(i,j)).or.skymed(i,j).gt.100.)then
						kkarg=kkarg+1
						skycombo(i,j)=skycombo(i,j)+array(indx(k))
					endif
					if(karg.ge.kmin/2)go to 97
				endif
			enddo
97			skycombo(i,j)=skycombo(i,j)/float(kkarg)
		endif
		if(kmin.le.1)then
c		print *,i,j,skycombo(i,j),kmin
		endif
		if(skycombo(i,j).gt.99999..or.skycombo(i,j).lt.-9999.)skycombo(i,j)=0.
	enddo
	enddo
	if(imap.ne.1)go to 99	!skip to end if no maps to be output
c	
	if(itype.eq.e)then
		if(ixstars.eq.1)open (11,file='c:\EngineeringOrbits\eq_allskymeed.grd')
		if(ixstars.eq.0)open (11,file='c:\EngineeringXStars\eq_allskymeed.grd')
		write(11,11)
	elseif(itype.eq.n)then
		if(ixstars.eq.1)open (11,file='c:\EngineeringOrbits\np_allskymeed.grd')
		if(ixstars.eq.0)open (11,file='c:\EngineeringXStars\np_allskymeed.grd')
		write(11,12)
	elseif(itype.eq.s)then
		if(ixstars.eq.1)open (11,file='c:\EngineeringOrbits\sp_allskymeed.grd')
		if(ixstars.eq.0)open (11,file='c:\EngineeringXStars\sp_allskymeed.grd')
		write(11,12)
	endif
	do j=1,jmax
		write(11,10)(skycombo(i,j),i=1,imax)
	enddo
	close(11)
99	continue
	stop 'done'
	end
c*********************************************************************************
      SUBROUTINE INDEXX(N,ARRIN,INDX)	!from Numerical Receipes...
      implicit none
      integer*4 N
      integer*4 INDX(N),I,J,L,IR,INDXT
      real*4 ARRIN(N),Q
      DO 11 J=1,N
        INDX(J)=J
11    CONTINUE
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
          ENDIF
          IF(Q.LT.ARRIN(INDX(J)))THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        INDX(I)=INDXT
      GO TO 10
      END