	program ntotal
C  April 15, 2002
C This program averages the photometry from the northern and southern pair
C of a 42-inch catalog
	implicit double precision (a-h,o-z)
	parameter(maxstars=250000)
	parameter(nc=2)
	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),idif(5),dif(5,maxstars),x(maxstars)
	dimension irvec1(4),idvec1(4)
	character*80 infile
	character*120 inline
	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,2
	write(6,11)
11	format('Input the name of a 42-inch catalog:',$)
	read(5,1111) infile
1111	format(a80)
	nstars(ichip)=0
	istar=0
	open(unit=1,file=infile,status='old')
15	read(1,14,end=19) inline
Cwrite(6,14) inline
14	format(a120)
	if(inline(1:1).eq.'#') go to 15
	call bear(inline,ihr,imin,sec,ideg,iamin,asec,isuc)
Cwrite(6,*) isuc,ihr,imin
	if(isuc.eq.0) go to 15
	istar=istar+1
	ra(istar,ichip)=stor(ihr,imin,sec)
	dec(istar,ichip)=astor(ideg,iamin,asec)
C Now find the photometry
C Now clean up "inline" and find the 3rd ":" while you're at it.
	ic=0
	ib=0
	do 110 i=1,115
	if(inline(i:i+4).eq.'INDEF') write(inline(i:i+4),111)
111	format('99.99')
	if(inline(i:i).eq.':') ic=ic+1
Cwrite(6,1123) i,inline(i:i),ic
1123	format('i,inline(i),ic=',i5,1x,a1,i5)
C this is obscure, but it is really finding the 3rd one.  think about it.
	if(ic.eq.2) ib=i
110	continue
	if(ib.eq.0) stop "problem here"
Cwrite(6,14) inline
Cwrite(6,*) ib
	read(inline(ib+10:120),*) (phot(istar,j,ichip),
     &                             phote(istar,j,ichip),
     &                             j=1,5)
Cwrite(14,1333) (phot(istar,j,ichip),phote(istar,j,ichip),j=1,5)
1333	format(10f8.3)
	go to 15
19	nstars(ichip)=istar
	close(unit=1)
	write(6,18) nstars(ichip)
18	format('There were',i9,' stars in this file')
10	continue
C OK, now look for matches, and all that...
	do 35 i=1,5
	idif(i)=0
35	continue
	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,2
	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
C Since there are only two files I can get away with this:
	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
	if(phote(l,icol,k).lt.0.10) then
	idif(icol)=idif(icol)+1
	dif(icol,idif(icol))=phot(j,icol,i)-phot(l,icol,k)
Cif(icol.eq.1) write(13,*) dif(icol,idif(icol))
	end if
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
	do 305 i=1,5
	if(idif(i).eq.0) go to 305
	do 306 ij=1,idif(i)
	x(ij)=dif(i,ij)
306	continue
	call median(x,idif(i),amed)
	write(6,307) i,idif(i),amed
307	format('For color:',i5,' there were',i5,' stars in common; med=',
     &  f8.3)
305	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'
	subroutine bear(inline,ihr,imin,sec,ideg,iamin,asec,isuc)
C The assumption is that the ra and dec have been printed as 12.3H 12.2h
	character*120 inline
	character*1 sign
	implicit double precision (a-h,o-z)
C step by step
	isuc=0
	i=2
1	continue
	if(inline(i:i).eq.':') then
	read(inline(i-2:i-1),*) ihr
	else
	i=i+1
	if(i.eq.120) return 
	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
	isuc=1
	return
	else
	k=k+1
	go to 2
	end if
	end
