	program matchcheckfin
C  April 14, 2002
C compares the "total.ubvri" with 42-inch catalog
C "close the loop"
C Mildly rewritten version of "matchem.f"
C Assumes that the catalog name is still stored in "tempcat"
	implicit double precision (a-h,o-z)
	parameter(maxstars=250000)
	character*120 inline
	character*80 filein
	dimension ra(maxstars),dec(maxstars),v(maxstars),
     &  ve(maxstars),bmv(maxstars),bmve(maxstars),
     &  umb(maxstars),umbe(maxstars),
     &  vmr(maxstars),vmre(maxstars),
     &  rmi(maxstars),rmie(maxstars),
     &  ra42(maxstars),dec42(maxstars),v42(maxstars),
     &  ve42(maxstars),bmv42(maxstars),bmve42(maxstars),
     &  umb42(maxstars),umbe42(maxstars),
     &  vmr42(maxstars),vmre42(maxstars),
     &  rmi42(maxstars),rmie42(maxstars),
     &  difv(maxstars),difbmv(maxstars),difumb(maxstars),
     &  difvmr(maxstars),difrmi(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 first, get the name of the 42-inch catalog:
	open(unit=1,file='tempcat',status='old')
	read(1,22) filein
22	format(a80)
	close(unit=1)
	open(unit=1,file=filein,status='old')
C Read the survey catalog into memory
	write(6,1001) filein
1001	format('Reading the catalog into memory:',a80)
1	read(1,2,end=20) inline
2	format(a120)
	if(inline(1:1).eq.'#') go to 1
	call bear(inline,ihr,imin,sec,ideg,iamin,asec,isuc)
	if(isuc.eq.0) go to 1
	ik=ik+1
	if(ik.gt.maxstars) stop 'too many stars'
	ra42(ik)=stor(ihr,imin,sec)
	dec42(ik)=astor(ideg,iamin,asec)
Cwrite(6,221) ik,ra42(ik),dec42(ik)
221	format(i7,2g12.6)
C OK, find the first ":"
	do 3 j=1,50
	if(inline(j:j).eq.':') go to 4
3	continue
C give up
	write(6,33) inline
33	format('Minor problem here:',a120)
	go to 1
C Now get rid of the indef's, replacing them with 99.99
4	continue
	do 110 i=j+22,115
	if(inline(i:i+4).eq.'INDEF') write(inline(i:i+4),111)
111	format('99.99')
110	continue
	read(inline(j+22:120),*) v42(ik),ve42(ik),bmv42(ik),
     &  bmve42(ik),
     &  umb42(ik),umbe42(ik),vmr42(ik),vmre42(ik),
     &  rmi42(ik),rmie42(ik)
	go to 1
20	continue
	close(unit=1)
	num42=ik
	write(6,1003) num42
1003	format('Total number of stars in 42-inch catalog=',i5)
	open(unit=1,file='total.ubvri',status='old')
	open(unit=3,file='total.ubvridif',status='new')
	ik=0
1008    read(1,2,end=1100) inline
	call bear(inline,ihr,imin,sec,ideg,iamin,asec,isuc)
	if(isuc.eq.0) stop "surprising development"
	ik=ik+1
	if(ik.gt.maxstars) stop "too many stars in Mos cat"
	ra(ik)=stor(ihr,imin,sec)
	dec(ik)=astor(ideg,iamin,asec)
	read(inline(26:120),*) v(ik),ve(ik),bmv(ik),bmve(ik),
     &  umb(ik),umbe(ik),vmr(ik),vmre(ik),rmi(ik),rmie(ik)
	go to 1008
1100	num=ik
	write(6,1101) num
C zero the difference counters
	idifv=0
	idifbmv=0
	idifumb=0
	idifvmr=0
	idifrmi=0
1101	format('total.ubvri contains',i8,' stars.')
C Now, find an uncrowded star
C
	do 201 i=1,num
	ra1=ra(i)
	dec1=dec(i)
	do 202 j=1,num
	if(j.eq.i) go to 202
	ra2=ra(j)
	dec2=dec(j)
	call apa(ra1,dec1,ra2,dec2,sep,pa)
	dist=sep/as2r
	if(dist.gt.5.) go to 202
C Is the second star brigher than 2.5 mag fainter?
	dmag=v(i)-v(j)
C 	-5=    10.   - 15.
C       -2.5=  10.   - 12.5
C        0=    10.   - 10.0
C       2.5 =  10.   - 7.5
C so, don't worry about if dmag<-2.5
	if(dmag.lt.-2.5) go to 202
C otherwise, forget this star
	go to 201
202	continue
C OK, we're here, which means this star is well isolated.  Now try to find
C a match in the 42-inch catalog
	ibest=0
	closest=20000.
	do 205 k=1,num42
	ra2=ra42(k)
	dec2=dec42(k)
	call apa(ra1,dec1,ra2,dec2,sep,pa)
	dist=sep/as2r
	if(dist.lt.5.) then
C this star is a winner...we hope...
	closest=dist
	ibest=k
	end if
205	continue
	if(ibest.gt.0) then
	if((ve(i).lt.0.03).and.(ve42(ibest).lt.0.03)) then
	idifv=idifv+1
	difv(idifv)=v(i)-v42(ibest)
	end if
	if((bmve(i).lt.0.05).and.(bmve42(ibest).lt.0.05)) then
	idifbmv=idifbmv+1
	difbmv(idifbmv)=bmv(i)-bmv42(ibest)
	end if
	if((umbe(i).lt.0.05).and.(umbe42(ibest).lt.0.05)) then
	idifumb=idifumb+1
	difumb(idifumb)=umb(i)-umb42(ibest)
	end if
	if((vmre(i).lt.0.05).and.(vmre42(ibest).lt.0.05)) then
	idifvmr=idifvmr+1
	difvmr(idifvmr)=vmr(i)-vmr42(ibest)
	end if
	if((rmie(i).lt.0.05).and.(rmie42(ibest).lt.0.05)) then
	idifrmi=idifrmi+1
	difumb(idifrmi)=rmi(i)-rmi42(ibest)
	end if
	write(3,206) v(i),v(i)-v42(ibest),
     &              bmv(i),bmv(i)-bmv42(ibest),
     &              umb(i),umb(i)-umb42(ibest),
     &              vmr(i),vmr(i)-vmr42(ibest),
     &              rmi(i),rmi(i)-rmi42(ibest)
206	format(10f8.3)
	end if
201	continue
	call median(difv,idifv,amedv)
	call median(difbmv,idifbmv,amedbmv)
	call median(difumb,idifumb,amedumb)
	call median(difvmr,idifvmr,amedvmr)
	call median(difrmi,idifrmi,amedrmi)
	write(6,2011) idifv,amedv,
     &                idifbmv,amedbmv,
     &                idifumb,amedumb,
     &                idifvmr,amedvmr,
     &                idifrmi,amedrmi 
2011	format('For total ','V: num stars=',i7,' Delta =',f9.3/
     &         '   ',5x,'B-V: num stars=',i7,' Delta =',f9.3/
     &         '   ',5x,'U-B: num stars=',i7,' Delta =',f9.3/
     &         '   ',5x,'V-R: num stars=',i7,' Delta =',f9.3/
     &         '   ',5x 'R-I: num stars=',i7,' Delta =',f9.3/)
	close(unit=3)
	close(unit=1)
	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
		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)
Cwrite(6,12) inline(1:70)
12	format('Bear:',a70)
	isuc=0
C step by step
	i=2
1	continue
	if(inline(i:i).eq.':') then
	read(inline(i-2:i-1),*) ihr
	else
	i=i+1
	if(i.lt.30) go to 1
C return with failure dripping from our brow
	return
	end if
	isuc=1
	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
        include '/tofu/home/massey/observing/astromsub.f'
        include 'mediandbl.f'
