	program resid3
c
c	program reads "keeper" orbits in file orb6.master, then
c 	matches objects against a list of observations. Observed
c	systems with no matching orbits are written to file
c 	resid3.noorb, while systems with matching orbits are written
c	to file resid3.yesorb with O-C residuals appended.
c
	character*1 keep,junk,pcode,acode,tcode
	character*8 ref(2500),xref
        character*7 ystar1
        character*5 ystar2
	character*3 xgrade,grade(2500)
        character*10 xwds(2500),wds
	character*12 xstar(2500),xxstar,ystar
	character*100 data
	dimension xp(2500),xa(2500),xi(2500),xn(2500),xt(2500),
     $     xe(2500),xo(2500)
c
	open(10,file='orb6.master',  status='UNKNOWN')
	open(11,file='resid3.in',    status='UNKNOWN')
	open(12,file='resid3.noorb', status='UNKNOWN')
	open(13,file='resid3.yesorb',status='UNKNOWN')
c
c	read "header" info from orbit catalog
c
	do 100 n=1,4
  100	read(10,901) junk
  901	format(a1)
c
	norb=0
  200	read(10,902,end=250) xxwds,xxstar,xxp,pcode,xxa,acode,xxi,
     $     xxn,xxt,tcode,xxe,xxo,xgrade,keep,xref
  902	format(t20,a10,1x,a12,t82,f11.6,a1,t106,f9.5,a1,t126,f8.4,
     $     t144,f8.4,t163,f12.6,a1,t188,f8.6,t206,f8.4,t244,a3,
     $     t250,a1,t252,a8)
	if (keep .ne. 'y') go to 200
	norb=norb+1
	xstar(norb)=xxstar
        if (pcode .eq. 'c') xxp=xxp*100.
        if (pcode .eq. 'd') xxp=xxp/365.2421987
	if (pcode .eq. 'h') xxp=xxp/365.2421987/24.
	if (pcode .eq. 'm') xxp=xxp/365.2421987/24./60.
        if (acode .eq. 'm') xxa=xxa/1000.
        if (tcode .eq. 'd') xxt=(xxt-33282.423)/365.2421987+1950.
	xp(norb)=xxp
	xa(norb)=xxa
	xi(norb)=xxi
	xn(norb)=xxn
	xt(norb)=xxt
	xe(norb)=xxe
	xo(norb)=xxo
	grade(norb)=xgrade
	ref(norb)=xref
	go to 200
c
  250	write(6,991) norb
  991	format(' Finished reading orbit file; norb = ',i6)
c
c	start reading from observations list and check for matches
c	
  300	read(11,903,end=999) data,ystar1,ystar2,date,pa,sep
c  903	format(a100,t1,a12,f9.4,f7.2,f8.4)
  903	format(a100,t13,a7,1x,a5,f9.4,f7.2,f8.4)
        ystar=ystar1//ystar2
c
        if (date .lt. 2004.) pa=pa-1.84
        if (date .gt. 2004.) pa=pa-5.26
c
	iflag=0
  	do 400 n=1,norb
  	if (ystar .ne. xstar(n)) go to 400
c
c	Match - determine residual and write to file resid3.yesorb. 
c	Continue checking against orbits in case this star has more 
c	than one good orbit
c
        xxp=xp(n)
        xxa=xa(n)
        xxi=xi(n)
        xxn=xn(n)
        xxt=xt(n)
        xxe=xe(n)
        xxo=xo(n)
	call porbit(date,pa,sep,xxp,xxa,xxi,xxn,xxt,xxe,xxo,
     $     theta,rho,respa,ressep)
	iflag=iflag+1  	
  	write(13,904) data,theta,rho,respa,ressep,ref(n),grade(n),iflag
  904	format(a100,f8.2,f9.4,2x,f8.2,f9.4,2x,a8,2x,a3,i3)
  400	continue
c
c	If iflag=0 then no match - write to resid3.noorb, then cycle 
c	back for next observation
c
	if (iflag .eq. 0) write(12,905) data
  905	format(a100)
  	go to 300
c
  999	stop
  	end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine porbit(date,pa,sep,xxp,xxa,xxi,xxn,xxt,xxe,xxo,
     $     theta,rho,respa,ressep)
	degrad=57.2957795
	ee=sqrt((1.+xxe)/(1.-xxe))
	aee=xxa*(1.-xxe*xxe)
	xxo=xxo/degrad
	xxn=xxn/degrad
	cosi=cos(xxi/degrad)
	xmu=6.2831853/xxp
	em=xmu*(date-xxt)
c
	e0=em+xxe*sin(em)+0.5*xxe*xxe*sin(2.*em)
	do 200 l=1,10
	em0=e0-xxe*sin(e0)
	e0=e0+(em-em0)/(1.-xxe*cos(e0))
  200	continue
c
	xnu=2.*atan(ee*tan(0.5*e0))
	r=aee/(1.+xxe*cos(xnu))
	theta=xxn+atan(cosi*tan(xnu+xxo))
	rho=r*cos(xnu+xxo)/cos(theta-xxn)
	theta=theta*degrad
	if(rho .gt. 0.) go to 300
	rho=-rho
	theta=theta+180.
  300	if(theta .lt.   0.) theta=theta+360.
	if(theta .gt. 360.) theta=theta-360.
c
	respa=pa-theta
	if(respa .lt. -180.) respa=respa+360.
	if(respa .gt.  180.) respa=respa-360.
	ressep=sep-rho
c
  	return
	end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
