	program philmagselect
C      written by p. massey, Lowell Observatory, March 2002
C
C
C This goes through a final "calibrated" file coordinate file (blah.ubvrird)
C and selects out only the least crowded stars
C
	parameter(maxstar=200000)
        implicit double precision (a-h,o-z)
	dimension ra(maxstar),dec(maxstar),v(maxstar),ve(maxstar),
     &  b(maxstar),be(maxstar),u(maxstar),ue(maxstar),r(maxstar),
     &  re(maxstar),ai(maxstar),aie(maxstar)
	character*160 inline,inlinen(maxstar)
	character*80 filename
*  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 plate scale
	scale=1.13


	write(6,2001)
2001	format('Input name of ubvrird file:',$)
	read(5,2002) filename
2002	format(a80)
	last=lnblnk(filename)
	open(unit=1,file=filename,status='old')
        open(unit=2,file=filename(1:last)//'ed',status="new")
	idebug=1
	i=0
C Get another star and store all the facts
1	read(1,2,end=1001) inline
	if((inline(1:1).eq.'#').or.(inline(1:1).eq.'').or.
     &      (inline(1:1).eq.' ')) then
	write(2,2) inline
	go to 1
	end if
	if(i.gt.maxstar) stop 'TOO MANY STARS!'
2	format(a160)
C OK, store everything
	call getstuff(inline,ihr,imin,sec,ideg,iamin,asec,
     &  vret,veret,bret,beret,uret,ueret,rret,reret,airet,aieret,isuc)
	if(isuc.ne.1) go to 1
C In other words, jsut skip a star if anything is INDEF
	i=i+1
	inlinen(i)=inline
	if(i.gt.maxstar) stop 'TOO MANY STARS'
	ra(i)=STOR(ihr,imin,sec)
	dec(i)=ASTOR(ideg,iamin,asec)
	v(i)=vret
	ve(i)=veret
	u(i)=uret
	ue(i)=ueret
	b(i)=bret
	be(i)=beret
	re(i)=reret
	r(i)=rret
	ai(i)=airet
	aie(i)=aieret
	go to 1
1001	inum=i
	close(unit=1)
	write(6,1002) inum
1002	format('Read in',i5,' stars')
3000	continue
	iz=0
	do 1000 i=1,inum
	do 500 j=1,inum
	if(j.eq.i) go to 500
	call apa(ra(i),dec(i),ra(j),dec(j),dist,pa)
	d=dist/as2r
C d is now in arcseconds
C now make it pixels
	dp=d/scale
	if(d.gt.22.) go to 500
C OK, it's nearby.  Now what?
C How much fainter is amag2?
	dif=amag-amag2
C      -10=  15 -  25
C       -5=  15 -  20
C        0=  15 -  15
C        5=  15 -  10
C      
	difv=v(i)-v(j)
	difu=u(i)-u(j)
	difb=b(i)-b(j)
	difr=r(i)-r(j)
	difi=ai(i)-ai(j)
	dif=-1000
	if(difv.gt.dif) dif=difv
	if(difb.gt.dif) dif=difb
	if(difu.gt.dif) dif=difu
	if(difr.gt.dif) dif=difr
	if(difi.gt.dif) dif=difi
C
C
	if(dif.lt.-4) go to 500
	if((dp.gt.(11.+2.)).and.(dif.lt.-2.5)) go to 500
	if((dp.gt.(11.+5.)).and.(dif.lt.-1.5)) go to 500
	if((dp.gt.(11.+7.)).and.(dif.lt.-0.5)) go to 500
	if((dp.gt.(11.+9.)).and.(dif.lt.0.0)) go to 500
	if(debug.eq.1) write(6,37)
37	format('WE HAVE REJECTED THIS POOR STAR!')
	go to 1000
500	continue
88      write(2,2) inlinen(i)
	iz=iz+1
1000	continue
	write(6,22001) iz,inum
22001	format('We found',i5,' isolated stars out of ',i5)
	close(unit=2)
	stop
	end
	subroutine getstuff(inline,ihr,imin,sec,ideg,iamin,asec,
     &  vret,veret,bret,beret,uret,ueret,rret,reret,airet,aieret,isuc)
        implicit double precision (a-h,o-z)
	character*160 inline
C33CA-N-1 1:33:37.778 30:57:13.20  17.481  0.028   0.533   0.059   0.331   0.123
C  0.420   0.037   0.488   0.042
C        1 2  3      4  5  6     7
	idebug=0
	iflag1=0
	iflag2=0
	iflag3=0
	iflag4=0
	iflag5=0
	iflag6=0
	iflag7=0
	iiflag1=0
	iiflag2=0
	iiflag3=0
	iiflag4=0
	iiflag5=0
	iiflag6=0
	iiflag7=0
	do 10 i=1,160
	isuc=0
	if((inline(i:i).eq.' ').and.(iflag1.eq.0))then
	iflag1=1
	iiflag1=i
	go to 10
	end if
	if((inline(i:i).eq.':').and.(iflag2.eq.0)) then
	iflag2=1
	iiflag2=i
	go to 10
	end if
	if((inline(i:i).eq.':').and.(iflag3.eq.0)) then
	iflag3=1
	iiflag3=i
	go to 10
	end if
	if((inline(i:i).eq.' ').and.(iflag3.eq.1).and.(iflag4.eq.0)) then
	iflag4=1
	iiflag4=i
	go to 10
	end if
	if((inline(i:i).eq.':').and.(iflag5.eq.0)) then
	iflag5=1
	iiflag5=i
	go to 10
	end if
	if((inline(i:i).eq.':').and.(iflag6.eq.0)) then
	iflag6=1
	iiflag6=i
	go to 10
	end if
	if((inline(i:i).eq.' ').and.(iflag6.eq.1).and.(iflag7.eq.0)) then
	iflag7=1
	iiflag7=i
	end if
10	continue
	if(iflag1*iflag2*iflag3*iflag4*iflag5*iflag6*iflag7.ne.1) then
	write(6,11) inline
11	format('Problem with this line:',/a160)
	write(6,5) iflag1,iiflag1
	write(6,5) iflag2,iiflag2
	write(6,5) iflag3,iiflag3
	write(6,5) iflag4,iiflag4
	write(6,5) iflag5,iiflag5
	write(6,5) iflag6,iiflag6
	write(6,5) iflag7,iiflag7
	stop
	end if
C OK, now here goes!
	read(inline(iiflag7:160),*,err=999) v,ve,bmv,bmve,
     &  umb,umbe,vmr,vmre,rmi,rmie
C if we're here, then it worked.  Get ready to return
C V
	vret=v
	veret=ve
C B
	bret=bmv+v
	be2=bmve**2-ve**2
	if(be2.gt.0.) then
	beret=sqrt(be2)
	else
	beret=0.
	end if
	if(beret.lt.0.002) beret=0.002
C U
	uret=umb+bret
	ue2=umbe**2-beret**2
	if(ue2.gt.0.) then
	ueret=sqrt(ue2)
	else
	ueret=0.
	end if
	if(ueret.lt.0.002) ueret=0.002
C R
	rret=v-vmr
	re2=vmre**2-veret**2
	if(re2.gt.0.) then
	reret=sqrt(re2)
	else
	reret=0.
	end if
	if(reret.lt.0.002) reret=0.002

C I
	airet=rret-rmi
	aie2=rmie**2-reret**2
	if(aie2.gt.0.) then
	aieret=sqrt(aie2)
	else
	aieret=0.
	end if
	if(aieret.lt.0.002) aieret=0.002
C now get the RA:
	read(inline(iiflag1:iiflag2-1),*) ihr
	read(inline(iiflag2+1:iiflag3-1),*) imin
	read(inline(iiflag3+1:iiflag4),*) sec
	read(inline(iiflag4:iiflag5-1),*) ideg
	read(inline(iiflag5+1:iiflag6-1),*) iamin
	read(inline(iiflag6+1:iiflag7),*) asec
	do 3 i= iiflag4,iiflag5
	if(inline(i:i).eq.'-') then
	ideg=-abs(ideg)
	iamin=-iamin
	asec=-asec
	end if
3	continue
	isuc=1
	if(idebug.eq.1) then
	write(6,4) inline
	write(6,5) iflag1,iiflag1
	write(6,5) iflag2,iiflag2
	write(6,5) iflag3,iiflag3
	write(6,5) iflag4,iiflag4
	write(6,5) iflag5,iiflag5
	write(6,5) iflag6,iiflag6
	write(6,5) iflag7,iiflag7
4	format(a160)
5	format(2i5)
	write(6,6) ihr,imin,sec,ideg,iamin,asec
6	format(i2.2,1x,i2.2,1x,f6.3,1x,i3,1x,i3,1x,f6.2)
	write(6,7) vret,veret,bret,beret,uret,ueret,
     &             rret,reret,airet,aieret
7	format('Photometry: ',10f6.2)
	end if
999	return
	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 '/home/massey/observing/astromsub.f'
