	program imtotal
C  April 14, 2002
C This program averages the photometry from chip to chip and produces
C a single output file from the 8 chips.
C It also compares the photometry from one chip to the next.
C Really?
	implicit double precision (a-h,o-z)
	parameter(maxstars=250000)
	parameter(nc=8)
	dimension ra(maxstars,nc),dec(maxstars,nc),
     &  phot(maxstars,5,nc),phote(maxstars,5,nc),
     &  nstars(nc),istat(maxstars,nc),
     &  num(5),sumerr(5),sumwt(5),sum(5),
     &  ave(5),aveerr(5),sigerr(5)
	dimension irvec1(4),idvec1(4)
	character*9 mosfil
	character*1 sign
*  SECONDS TO RADIANS
      PARAMETER (S2R=0.72722052166430398D-4)
*  ARC SECONDS TO RADIANS
      PARAMETER (AS2R=0.48481368110953599D-5)
C
*  SECONDS TO RADIANS
      STOR(RH,RM,RS)=S2R*(60D0*(60D0*RH+RM)+RS)
*  ARC SECONDS TO RADIANS
      ASTOR(DD,DM,DS)=AS2R*(60D0*(60D0*DD+DM)+DS)
C Distance to be considered the same star (in arcsec)
	data hit/0.75/

C OK, now let's read all the files
	do 10 ichip=1,nc
	write(mosfil,11) ichip
11	format('IM',i1,'.ubvri')
C see if this file exists
	nstars(ichip)=0
	istar=0
	open(unit=1,file=mosfil,status='old',err=12)
	go to 13
12	write(6,144) mosfil
144	format('This file does not exist: ',a9) 
	go to 10
13	istar=istar+1
C zero the most important counter!
	istat(istar,ichip)=0
	if(istar.gt.maxstars) stop 'too many stars in total'
	read(1,14,end=19) ihr,imin,sec,sign,ideg,iamin,asec,
     &  (phot(istar,j,ichip),phote(istar,j,ichip),j=1,5)
14	format(i2,1x,i2,1x,f6.3,1x,a1,i2,1x,i2,1x,f5.2,10f7.3)
	if(sign.eq.'-') then
	ideg=-ideg
	iamin=-iamin
	asec=-asec
	end if
	ra(istar,ichip)=stor(ihr,imin,sec)
	dec(istar,ichip)=astor(ideg,iamin,asec)
C Uh, great, I guess
	go to 13
19	nstars(ichip)=istar
	close(unit=1)
	write(6,18) mosfil,nstars(ichip)
18	format('In ',a9,' there were',i9,' stars')
10	continue
C OK, now look for matches, and all that...
	do 100 i=1,nc
	write(6,133) i,nstars(i)
133	format('Checking file',i5', which contains',i5,' stars')
	if(nstars(i).eq.0) go to 100
	do 200 j=1,nstars(i)
C Has this star been vetted?
	if(istat(j,i).ne.0) go to 200
C So now look at the subsequent files
	do 201 icol=1,5
	num(icol)=0
	if (phot(j,icol,i).gt.90.) go to 201
	num(icol)=1
	sumerr(icol)=phote(j,icol,i)
	if(phote(j,icol,i).lt..005) then
	wt=1./.005**2
	else
	wt=1./phote(j,icol,i)**2
	end if
	sumwt(icol)=wt
	sum(icol)=phot(j,icol,i)*wt
201	continue
	rapos=ra(j,i)
	decpos=dec(j,i)
C This is just a cheap way of not having to deal with the final file separately
	if(i.eq.j) go to 410
	do 300 k=i+1,8
	if(nstars(k).eq.0) go to 300
	do 400 l=1,nstars(k)
	if(istat(l,k).ne.0) go to 400
	if(abs(dec(l,k)-dec(j,i)).lt.1.0e-5) then
	call apa(ra(j,i),dec(j,i),ra(l,k),dec(l,k),sep,pa)
	else
	sep=1.
	end if
	dist=sep/as2r
	if(dist.lt.hit) then
C we have a winner!
	do 401 icol=1,5
	if(phot(l,icol,k).gt.90.) go to 401
	num(icol)=num(icol)+1
	sumerr(icol)=sumerr(icol)+phote(l,icol,k)
	if(phote(l,icol,k).lt.0.005) then
	st=1./.005**2
	else
	wt=1./phote(l,icol,k)**2
	end if
	sum(icol)=sum(icol)+phot(l,icol,k)*wt
	sumwt(icol)=sumwt(icol)+wt
401	continue
	istat(l,k)=1
	go to 300
	end if
400	continue
300	continue
410	continue
	do 411 icol=1,5
C allow for the posibility that a number was 90
	if (num(icol).gt.0) then
	ave(icol)=sum(icol)/sumwt(icol)
	aveerr(icol)=(sumerr(icol)/num(icol))/sqrt(float(num(icol)))
	sigerr(icol)=sqrt(1./sumwt(icol))
	else
	ave(icol)=99.999
	aveerr(icol)=99.999
	sigerr(icol)=99.999
	end if
411	continue
	call dr2tf(3,rapos,sign,irvec1)
	call dr2af(2,decpos,sign,idvec1)
        write(9,501) irvec1,sign,idvec1,
     &  (ave(icl),sigerr(icl),
     &  icl=1,5),(num(icol),icol=1,5)
501     format(i2.2,':',i2.2,':',i2.2,'.',i3.3,1x,a1,
     &         i2.2,':',i2.2,':',i2.2,'.',i2.2,
     &        10f7.3,5i2)
200	continue
100	continue
	stop
	end
        subroutine apa(rac,dcc,rao,dco,dist,pa)
c computes the distance (delt) and position (pa) along the
c sky FROM (rac,dec) to (rao,dco)
c written aug 14, 1987
C rewritten August 21, 1989, to do using "standard coords", rather
C than the old and peculiar (but possibly equivalent!) way
c plm
      implicit double precision (a-h,o-z)
        parameter (PI=3.141592653589793238462643)
C CONVERT TO PROPER COORDS
        CALL DS2TP(RAO,DCO,RAC,DCC,XI,ETA)
        pa=0.
        if(xi*eta.ne.0.)PA=DATAN2(XI,ETA)
        if(pa.le.0.) pa=pa+2.*pi
        DIST=DSQRT(XI**2+ETA**2)
        return
        end
        include '/tofu/home/massey/observing/astromsub.f'
        include 'mediandbl.f'
