	program matchem
C  April 4, 2002
C
C Assumes that the catalog name is stored in "tempcat", and that
C the survey file is in "temprd"
	implicit double precision (a-h,o-z)
	parameter(maxstars=250000)
	character*120 inline
	character*80 filein
	character*1 sign,filt
	dimension ra(maxstars),dec(maxstars),amag(maxstars),
     &  ame(maxstars),
     &  ra42(maxstars),dec42(maxstars),v(maxstars),ve(maxstars),
     &  bmv(maxstars),bmve(maxstars),umb(maxstars),umbe(maxstars),
     &  vmr(maxstars),vmre(maxstars),rmi(maxstars),rmie(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)

	open(unit=3,file='tempmatchem.dat',status='new')
	open(unit=1,file='temprd',status='old')
C Read the survey catalog into memory
	write(6,1001)
1001	format('Reading the survey into memory')
1	read(1,2,end=20) inline
2	format(a120)
	if(inline(1:1).eq.'#') go to 1
C Now go looking for the first ":'
	do 3 i=1,50
	if(inline(i:i).eq.':') go to 5
3	continue
C give up
	go to 1
5	read(inline(i-2:i+9),6) ihr,imin,sec
Cwrite(6,2) inline
Cwrite(6,6) ihr,imin,sec
Cwrite(6,2001) ihr,imin,sec
6	format(i2,1x,i2,1x,f6.3)
C Now find the declination:
	do 7 j=i+10,i+20
	if(inline(j:j).eq.':') go to 8
7	continue
	write (6,9) inline(1:50)
	stop 'PROBLEM IN MATCHEM!!! HELP!!'
9	format('Problem with formatting:',a50)
8	read(inline(j-3:j-3),10) sign
10	format(a1)
	read(inline(j-2:j+8),11) ideg,iamin,asec
Cwrite(6,11) ideg,iamin,asec
Cwrite(6,2001) ideg,iamin,asec
11	format(i2,1x,i2,1x,f5.2)
	ik=ik+1
	if(ik.gt.maxstars) stop 'too many stars'
	ra(ik)=stor(ihr,imin,sec)
	if(sign.eq.'-') then
	ideg=-ideg
	iamin=-iamin
	asec=-asec
	end if
	dec(ik)=astor(ideg,iamin,asec)
	read(inline(j+10:120),*) am,ae
	amag(ik)=am
	ame(ik)=ae
C determine filter for the first file
	if(ik.eq.1) then
	filt=''
	do 15 ij=j+10,120
	if(inline(ij:ij).eq.'B') then 
	filt='B'
	go to 1
	end if
	if(inline(ij:ij).eq.'V') then
	filt='V'
	go to 1
	end if
	if(inline(ij:ij).eq.'R') then
	filt='R'
	go to 1
	end if
	if(inline(ij:ij).eq.'I') then
	if(inline(ij+1:ij+2).eq.'ND') stop "PROBLEM IN MATCHEM"
	filt='I'
	go to 1
	end if
	if(inline(ij:ij).eq.'U') then
	filt='U'
	go to 1
	end if
15	continue
	end if
	go to 1
20	num=ik
	write(6,1002) num,filt
1002	format('Total number of stars=',i5,' Filter=',a1)
C OK, now read the 42inch catalog into memory:
	close(unit=1)
	ik=0
	open(unit=2,file='tempcat',status='old')
	read(2,22) filein
22	format(a80)
	close(unit=2)
	open(unit=1,file=filein,status='old')
	write(6,10011)
10011	format('Reading the 42-inch catalog into memory')
100	read(1,2,end=200) inline
	if(inline(1:1).eq.'#') go to 100
C Now go looking for the first ":'
	do 30 i=1,50
	if(inline(i:i).eq.':') go to 50
30	continue
C give up
	go to 100
50	read(inline(i-2:i+9),6) ihr,imin,sec
Cwrite(6,2001) ihr,imin,sec
C Now find the declination:
	do 70 j=i+10,i+20
	if(inline(j:j).eq.':') go to 80
70	continue
	write (6,9) inline(1:50)
	stop
80	read(inline(j-3:j-3),10) sign
	read(inline(j-2:j+8),11) ideg,iamin,asec
Cwrite(6,2001) ideg,iamin,asec
2001	format(2i3,f8.3)
	ik=ik+1
	if(ik.gt.maxstars) stop 'too many stars'
	ra42(ik)=stor(ihr,imin,sec)
	if(sign.eq.'-') then
	ideg=-ideg
	iamin=-iamin
	asec=-asec
	end if
	dec42(ik)=astor(ideg,iamin,asec)
C Now get rid of the indef's, replacing them with 99.99
	do 110 i=j+10,115
	if(inline(i:i+4).eq.'INDEF') write(inline(i:i+4),111)
111	format('99.99')
110	continue
	read(inline(j+10:120),*) v(ik),ve(ik),bmv(ik),bmve(ik),
     &  umb(ik),umbe(ik),vmr(ik),vmre(ik),rmi(ik),rmie(ik)
	go to 100
200	continue
	close(unit=1)
	num42=ik
	write(6,1003) num42
1003	format('Total number of stars in 42-inch catalog=',i5)
C Now, find an uncrowded star
	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=amag(i)-amag(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
	write(3,206) filt,closest,amag(i),ame(i),v(ibest),ve(ibest),
     &   bmv(ibest),bmve(ibest),
     &   umb(ibest),umbe(ibest),
     &   vmr(ibest),vmre(ibest),
     &   rmi(ibest),rmie(ibest)
206	format(a1,1x,f6.1,12f7.3)
	end if
201	continue
	close(unit=3)
	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'
