        implicit double precision (a-h,o-z)
	character*80 infile
	character*3 color(5)
	character*1 sign
	parameter(maxstars=250000)
	dimension ra(maxstars,2),dec(maxstars,2),phot(maxstars,5,2),
     &   phote(maxstars,5,2)
	dimension diff(maxstars,5),icol(5),data(maxstars),
     &  datae(maxstars),diffe(maxstars,5)
	data color/'V  ','B-V','U-B','V-R','R-I'/
*  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)
	write(6,1)
1	format('Enter name of catalog-1:',$)
	read(5,2) infile
	open(unit=1,file=infile,status='old')
	ix=0
2	format(a80)
101	read(1,102,end=190) ihr,imin,sec,sign,ideg,iamin,asec,
     &   v,ve,bmv,bmve,umb,umbe,vmr,vmre,rmi,rmie
102	format(i2,1x,i2,1x,f6.3,1x,a1,i2,1x,i2,1x,f5.2,
     &  10(1x,f6.3))
Cwrite(6,102) ihr,imin,sec,sign,ideg,iamin,asec,
C    &   v,ve,bmv,bmve,umb,umbe,vmr,vmre,rmi,rmie
	ix=ix+1
Cformat is fixed at this point...
	ra(ix,1)=stor(ihr,imin,sec)
	if(sign.eq.'-') then
	ideg=-ideg
	iamin=-iamin
	asec=-asec
	end if
	dec(ix,1)=astor(ideg,iamin,asec)
	phot(ix,1,1)=v
	phot(ix,2,1)=bmv
	phot(ix,3,1)=umb
	phot(ix,4,1)=vmr
	phot(ix,5,1)=rmi
	phote(ix,1,1)=ve
	phote(ix,2,1)=bmve
	phote(ix,3,1)=umbe
	phote(ix,4,1)=vmre
	phote(ix,5,1)=rmie
	go to 101
190	close(unit=1)
	inum1=ix
	write(6,191) inum1
191	format('Number of stars in 1st file=',i9)
	write(6,193)
193	format('Enter name of catalog-2:',$)
	read(5,2) infile
	open(unit=1,file=infile,status='old')
	ix=0
201	read(1,102,end=290) ihr,imin,sec,sign,ideg,iamin,asec,
     &   v,ve,bmv,bmve,umb,umbe,vmr,vmre,rmi,rmie
Cwrite(6,102) ihr,imin,sec,sign,ideg,iamin,asec,
C    &   v,ve,bmv,bmve,umb,umbe,vmr,vmre,rmi,rmie
	ix=ix+1
	ra(ix,2)=stor(ihr,imin,sec)
	if(sign.eq.'-') then
	ideg=-ideg
	iamin=-iamin
	asec=-asec
	end if
	dec(ix,2)=astor(ideg,iamin,asec)
	phot(ix,1,2)=v
	phot(ix,2,2)=bmv
	phot(ix,3,2)=umb
	phot(ix,4,2)=vmr
	phot(ix,5,2)=rmi
	phote(ix,1,2)=ve
	phote(ix,2,2)=bmve
	phote(ix,3,2)=umbe
	phote(ix,4,2)=vmre
	phote(ix,5,2)=rmie
	go to 201
290	close(unit=1)
	inum2=ix
	write(6,291) inum2
291	format('Number of stars in 2nd file:',i9)
	do 303 i=1,maxstars
	data(i)=0.
	datae(i)=0.
303	continue
C Now, find pairs
	do 302 k=1,5
	icol(k)=0
	do 302 i=1,maxstars
	diff(i,k)=0.
	diffe(i,k)=0.
302	continue
	do 300 i=1,inum1
	do 305 j=1,inum2
	call apa(ra(i,1),dec(i,1),ra(j,2),dec(j,2),dist,pa)
	d=dist/as2r
	if(d.lt.1.) then
C we have a winner!
C maybe...
	if(abs(phot(i,1,1)-phot(j,1,2)).gt.0.5) go to 305
C OK, now find the differences
	do 307 k=1,5
Cwrite(6,1001) k,phot(i,k,1),phot(j,k,2)
1001	format('Photometry for color:',i5,2f9.2)
	if((phot(i,k,1).gt.90.).or.(phot(j,k,2).gt.90.)) go to 307
	if((phote(i,k,1).gt.0.1).or.(phote(j,k,2).gt.0.1)) go to 307
C OK, count it!
	icol(k)=icol(k)+1
	diff(icol(k),k)=phot(i,k,1)-phot(j,k,2)
	diffe(icol(k),k)=sqrt(phote(i,k,1)**2+phote(j,k,2)**2)
307	continue
	end if
C
	
305	continue

300	continue
	do 310 k=1,5
	do 311 j=1,icol(k)
	data(j)=diff(j,k)
	datae(j)=diffe(j,k)
311	continue
	inum=icol(k)
Cwrite(6,*) inum,data(1)
        call xfit (data, datae, inum, 1, xmean, sigmam, sigma)
	call median(data,inum,amed)
	write(6,312) k,color(k),inum,amed,xmean
312	format('Median,mean difference for:',i3,1x,a3,1x,i6,2f8.3)
310	continue
	stop
	end	
	subroutine median(data,n,amed)
	implicit double precision(a-h,o-z)
	dimension data(*),indx(250000)
	if(n.eq.0) then
	amed=-99.
	return
	end if
	call sort(n,data,indx)
	if( ((n/2)*2).eq.n) then
	amed=(data(indx(n/2))+data(indx(n/2+1)))/2.
	else
	amed=data(indx(n/2+1))
	end if
	return
	end
        subroutine sort(n,arrin,indx)
	implicit double precision(a-h,o-z)
        dimension arrin(*),indx(*)
        do 11 j=1,n
        indx(j)=j
11      continue
        if(n.eq.1) return
        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
        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
c subroutine xfit.f
c
c source
c   Bevington, page 76.
c
c purpose
c   calculate the mean and estimated errors for a set of data points
c
c usage
c   call xfit (x, sigmax, npts, mode, xmean, sigmam, sigma)
c
c description of parameters
c   x      - array of data points
c   sigmax - array of standard deviations for data points
c   npts   - number of data points
c   mode   - determines method of weighting
c            +1 (instrumental) weight(i) = 1./sigmax(i)**2
c             0 (no weighting) weight(i) = 1.
c            -1 (statistical)  weight(i) = 1.
c   xmean  - weighted mean
c   sigmam - standard deviation of mean
c   sigma  - standard deviation of data
c
c subroutines and function subprograms required
c   none
c
	subroutine xfit (x,sigmax,npts,mode,xmean,sigmam,sigma)
        implicit double precision (a-h,o-z)
	dimension x(*),sigmax(*)
c
c accumulate weighted sums
c
11	sum=0.
	sumx=0.
	sigma=0.
	sigmam=0.
20	do 32 i=1,npts
21	if (mode) 22,22,24
22	weight=1.
	goto 31
24	weight=1./sigmax(i)**2
31	sum=sum+weight
32	sumx=sumx+weight*x(i)
c
c evaluate mean and standard deviations
c
41	xmean=sumx/sum
51	do 52 i=1,npts
52	sigma=sigma+(x(i)-xmean)**2
	free=npts-1
54	sigma=dsqrt(sigma/free)
61	if (mode) 62,64,66
62	sigmam=dsqrt(xmean/sum)
	goto 70
64	sigmam=sigma/dsqrt(sum)
	goto 70
66	sigmam=dsqrt(1./sum)
70	return
	end
        include '/tofu/home/massey/observing/astromsub.f'
