	program ave2 
C rewritten 12 April 2002
	implicit double precision (a-h,o-z)
	parameter(maxstars=250000)
	parameter(maxfiles=20)
	character*80 inline
	character*50 filename
	integer irvec1(4),idvec1(4)
	character*1 sign
	dimension ra(maxstars,maxfiles),dec(maxstars,maxfiles),
     &  am(maxstars,maxfiles),ame(maxstars,maxfiles), 
     &  istat(maxstars,maxfiles),nstars(maxfiles) 
	dimension diff(maxstars,maxfiles),ndif(maxfiles)
	dimension arr(maxstars)
*  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/
	open(unit=1,file="tempzero",status='old')
C Get the name of the next file
	nfiles=0
	iz=0
1	read(1,2,end=11) inline
2	format(a80)
Cwrite(6,2) inline
	do 3 i=2,50
	if(inline(i:i).eq.' ') then
	read(inline(i+1:80),*) zeropt
	go to 4
	end if
3	continue
	stop "problem finding the end of the file name in ave"
4	nfiles=nfiles+1
Cwrite(6,144) nfiles
144	format('Reading in file',i5)
	if(nfiles.gt.maxfiles) stop "too many files"
	filename=''
	filename=inline(1:i-1)//'.rd'
	close(unit=2)
	istar=0
	open(unit=2,file=filename,status='old')
5	read(2,2,end=10) inline
Cwrite(6,2) inline
	if(inline(1:1).eq.'#') go to 5
Cif(inline(1:1).eq.'') go to 5
Cwrite(6,2) inline
	istar=istar+1
	if(istar.gt.maxstars) stop "too many stars"
C Now read in the ra, dec, mag, and merr:
C I THINK this will work....
	call bear(inline,ihr,imin,sec,ideg,iamin,asec)
C slightly risky!
	ra(istar,nfiles)=stor(ihr,imin,sec)
	dec(istar,nfiles)=astor(ideg,iamin,asec)
	read(inline(26:80),*) amag,amerr	
	amagp=amag-zeropt
	am(istar,nfiles)=amagp
	ame(istar,nfiles)=amerr
	istat(istar,nfiles)=0
	go to 5
10	nstars(nfiles)=istar
	write(6,146) nfiles,istar
146	format('In file',i5,' there were',i5,' stars')
	go to 1
C OK, we now have nfiles worth of nstars(i) stars each.  Let's go looking
C for matches!
11      continue
C Take a star in the first file and see if it is in any of the other files:
C for each star in each file
	do 150 iff=1,nfiles
	ndif(iff)=0
150	continue
	do 100 i=1,nfiles
	write(6,133) i,nstars(i)
133	format('Checking file',i5,' which contains',i5,' stars')
	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
	num=1
	sumerr=ame(j,i)
	if(ame(j,i).lt.0.005) then
	wt=1./.005**2
	else
	wt=1./ame(j,i)**2
	end if
	sumwt=wt
	sum=am(j,i)*wt
	istat(j,i)=1
	rapos=ra(j,i)
	decpos=dec(j,i)
C this is just a cheap way of having to deal with the final file separately.
	if(i.eq.nfiles) go to 410
	do 300 k=i+1,nfiles
Cwrite(6,301) k,nstars(k)
301	format('Against file',i5,' containing',i5,' stars')
	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 maybe a hit!
C See if the mags are in relatively OK agreement...
C	dm2=(am(j,i)-am(l,k))**2
C	err2=ame(j,i)**2+ame(l,k)**2
Cwrite(6,302) dist,hit
302	format('We have a hit!!!!',2f8.2)
	num=num+1
	sumerr=sumerr+ame(l,k)
	if(ame(l,k).lt.0.005) then
	wt=1./.005**2
	else
	wt=1./ame(l,k)**2	
	end if
	sum=sum+am(l,k)*wt
	sumwt=sumwt+wt
	istat(l,k)=1
	if(i.eq.1) then
	if((ame(j,i).lt.0.05).and.(ame(l,k).lt.0.05)) then
	ndif(k-1)=ndif(k-1)+1
	diff(ndif(k-1),k-1)=am(j,i)-am(l,k)
	write(10,123) am(j,i)-am(l,k)
123	format(f10.3)
	end if
	end if
	go to 300
	end if
400	continue
300	continue
410	ave=sum/sumwt
	aveerr=(sumerr/num)/sqrt(float(num))
	sigerr=sqrt(1./sumwt)
C Make the ra and dec look like numbers again:
	call dr2tf(3,rapos,sign,irvec1)
	call dr2af(2,decpos,sign,idvec1)
	write(9,501) irvec1,sign,idvec1,ave,sigerr,num
501	format(i2.2,':',i2.2,':',i2.2,'.',i3.3,1x,a1,
     &         i2.2,':',i2.2,':',i2.2,'.',i2.2,
     &        2f9.3,1x,i3)
200	continue
100	continue
	do 500 ji=1,nfiles-1
	do 490 ii=1,ndif(ji)
	arr(ii)=diff(ii,ji)
490	continue
	call median(arr,ndif(ji),amed)
	write(6,491) ji+1,amed
491	format('Median difference between files 1 and',i3,' =',f10.3)
500	continue
	stop
	end
	subroutine bear(inline,ihr,imin,sec,ideg,iamin,asec)
C The assumption is that the ra and dec have been printed as 12.3H 12.2h
	character*80 inline
	character*1 sign
	implicit double precision (a-h,o-z)
C step by step
	i=2
1	continue
	if(inline(i:i).eq.':') then
	read(inline(1:i-1),*) ihr
	else
	i=i+1
	go to 1
	end if
	read(inline(i+1:i+2),*) imin
	read(inline(i+4:i+9),*) sec
	k=i+10
2	if(inline(k:k).eq.':') then
C first see if there's a minus sign or not
	sign=''
	do 3 j=i+10,k
	if(inline(j:j).eq.'-') sign='-'
3	continue
	read(inline(i+10:k-1),*) ideg
	ideg=abs(ideg)
	read(inline(k+1:k+3),*) iamin
	read(inline(k+4:k+8),*) asec
	if(sign.eq.'-') then
	ideg=-ideg
	iamin=-iamin
	asec=-asec
	end if
	return
	else
	k=k+1
	go to 2
	end if
	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 'mediandbl.f'
	include '/tofu/home/massey/observing/astromsub.f'
