program combine_eq_pl
	implicit none
c
c	This program combines eq, np, and sp maps into one
c	Note also, that the output format is changed, for the 2x2 binned output, to 20f10.2
c
	character*10 head(5)
	integer*4 i,j,imax,jmax,iflip,ix,iy,isum,ii,jj,iarg,jarg
	real*4 skycombo(3600,1800),skyeq(3600,1200),skynp(800,800),skysp(800,800),skycombo2(1800,900)
	real*4 x,y,r,theta,sum
	data iflip /1/	!iflip=1 means the input files (as is usual) are in Norton's coords.
c
	open(10,file='c:\EngineeringXStars\eq_allskymeed.grd',readonly)
	read (10,1)head
1	format(a10)
	imax=3600
	jmax=1200	
	do j=1,jmax			!read in map
		if(iflip.eq.0)read (10,10) (skyeq(i,j),i=1,imax)
		if(iflip.eq.1)read (10,10) (skyeq(3601-i,j),i=1,imax)
10		format(20f8.1)	
	enddo
	close(10)
c	
	open(10,file='c:\EngineeringXStars\np_allskymeed.grd',readonly)
	read (10,1)head
	imax=800
	jmax=800	
	do j=1,jmax			!read in map
		if(iflip.eq.0)read (10,10) (skynp(i,j),i=1,imax)
		if(iflip.eq.1)read (10,10) (skynp(i,801-j),i=1,imax)
	enddo
	close(10)
c	
	open(10,file='c:\EngineeringXStars\sp_allskymeed.grd',readonly)
	read (10,1)head
	do j=1,jmax			!read in map
		if(iflip.eq.0)read (10,10) (skysp(i,j),i=1,imax)
		if(iflip.eq.1)read (10,10) (skysp(801-i,801-j),i=1,imax)
	enddo
	close(10)
c
c	transfer equatorial plot
c
	imax=3600
	jmax=1200
	do j=1,jmax
	do i=1,imax
		skycombo(i,j+300)=skyeq(i,j)
	enddo
	enddo
c
c	transfer north & south polar plots, go to nearest map value
c
	do j=1,300
	    r=-0.05+0.1*float(j)
   	    do i=1,imax
		theta=-0.05+0.1*float(i)
		y=10.*(r*sind(theta))+400.
		x=10.*(r*cosd(theta))+400.
c		if(mod(i,10).eq.0)print *,i,x,y
		ix=nint(x)
		iy=nint(y)
		skycombo(i,j)=skysp(ix,iy)
		skycombo(i,1801-j)=skynp(ix,iy)
   	    enddo
	enddo
c
c	Output combined map
c
	open(10,file='c:\EngineeringXStars\comb_allskymeed.grd')
	write(10,110)
110	format('DSAA',/,'3600 1800',/,'0 360',/,'-90 90',/,'0 1000')
	do j=1,1800
		write(10,10)(skycombo(i,j),i=1,3600)
	enddo
	close(10)
c
c	Output combined map 2x2 binned for Bernie's program
c
	do i=1,1800
	do j=1,900
		sum=0.
		isum=0
		do ii=1,2
		do jj=1,2
			iarg=2*(i-1)+ii
			jarg=2*(j-1)+jj
			if(skycombo(iarg,jarg).gt.0.)then
				sum=sum+skycombo(iarg,jarg)
				isum=isum+1
			else
				print *,iarg,jarg,skycombo(iarg,jarg)
			endif
		enddo
		enddo
		skycombo2(i,j)=0.
		if(isum.gt.0)skycombo2(i,j)=sum/float(isum)
	enddo
	enddo				
	open(10,file='c:\EngineeringXStars\comb2_allskymeed.grd')
	write(10,111)
111	format('DSAA',/,'1800 900',/,'0 360',/,'-90 90',/,'0 1000')
	do j=1,900
		write(10,11)(skycombo2(i,j),i=1,1800)
11		format(20f10.2)	
	enddo
	close(10)
	stop	'done'
	end