        program olresids
c
c       determine orbit and linear O-Cs for file of measures
c       Checks both orb6 and lin1 files for matches, so no prior
c       search of catalogs needed.
c
        dimension p(2500),a(2500),xi(2500),xn(2500),t(2500),e(2500),
     $     xo(2500),x0(2000),xa(2000),y0(2000),ya(2000),t0(2000)
        character*1 pcode,acode,tcode,junk,dsign,grade,informat,rflag
        character*5 compo(2500),comp1,compl(2500)
        character*7 ddo(2500),dd1,ddl(2000)
        character*8 oref(2500),lref(2000),ref1,ai,an,ae,ao
        character*9 aa
        character*10 wdso(2500),wds1,wdsl(2000)
        character*11 ap
        character*12 at
        character*104 data
c
        open(10,file='/home/wih/Orbits/Orb6/orb6orbits.txt',
     $       status='UNKNOWN')
        open(11,file='/home/wih/Orbits/Linear/LIN1/lelements.txt',
     $       status='UNKNOWN')
        open(12,file='olresids.in',status='UNKNOWN')
        open(13,file='olresids.out',status='UNKNOWN')
c
	degrad=0.0174533
	equ=2000.0
c
c       read orbits from orb6. First reads 7 lines from header, then
c       stars and elements. Rejects systems with incomplete elements,
c       as well as astrometric systems. Converts P,T, and a to years
c       and arcseconds if necessary
c
        do 100 n=1,7
  100   read(10,901) junk
  901   format(a1)
        norb=0
c
  200   read(10,902,end=300,err=900) wds1,dd1,comp1,p1,pcode,a1,acode,
     $     xi1,xn1,t1,tcode,e1,xo1,grade,ref1,ap,aa,ai,an,at,ae,ao
  902   format(t20,a10,1x,a7,a5,t81,f12.6,a1,t106,f9.5,a1,t126,f8.4,
     $     t144,f8.4,t163,f12.6,a1,t188,f8.6,t206,f8.4,t234,a1,3x,a8,
     $     t82,a11,t106,a9,t126,a8,t144,a8,t163,a12,t188,a8,t206,a8)
c
        if (grade .eq. '9') go to 200
        if (ap .eq. '    .      ') go to 200
        if (aa .eq. '   .     ') go to 200
        if (ai .eq. '   .    ') go to 200
        if (an .eq. '   .    ') go to 200
        if (at .eq. '     .      ') go to 200
        if (ae .eq. ' .     ') go to 200
        if (ao .eq. '   .    ') go to 200
        if (pcode .eq. 'c') p1=p1*100.
        if (pcode .eq. 'd') p1=p1/365.2422
        if (acode .eq. 'm') a1=a1/1000.
        if (tcode .eq. 'd') t1=1950.+((t1-33282.423)/365.2422)
        if (comp1 .eq. '     ') comp1='AB   '
        norb=norb+1
        wdso(norb)=wds1
        ddo(norb)=dd1
        compo(norb)=comp1
        oref(norb)=ref1
        p(norb)=p1
        a(norb)=a1
        xi(norb)=xi1
        xn(norb)=xn1
        t(norb)=t1
        e(norb)=e1
        xo(norb)=xo1
        go to 200
c
c       read elements from linear catalog (first read 6 lines of header)
c
  300   do 350 n=1,6
  350   read(11,901) junk
c
        nlin=0
c
  400   read(11,903,end=450,err=410) wds1,dd1,comp1,x01,xa1,y01,ya1,t01,
     $     ref1
  903   format(t21,a10,2x,a7,a5,t91,2f11.6,2x,2f11.6,f11.3,t189,a8//)
        if (comp1 .eq. '     ') comp1='AB   '
        nlin=nlin+1
        wdsl(nlin)=wds1
        ddl(nlin)=dd1
        compl(nlin)=comp1
        x0(nlin)=x01
        xa(nlin)=xa1
        y0(nlin)=y01
        ya(nlin)=ya1
        t0(nlin)=t01
	lref(nlin)=ref1
        go to 400
  410	write(6,994) wds1,dd1
  994	format(' lelements read error near ',a10,1x,a7)
	go to 999
c
c       read measures, match against each catalog in turn, derive
c       O-C values, and write out results
c
  450	write(6,995) norb,nlin
  995	format(' orbit and linear files read; norb, nlin = ',2i6//)
	informat='n'
  500   read(12,904,end=999) data,wds1,dd1,comp1,date,theta,rflag,rho,
     $     rah,ram,dsign,decd,decm
  904   format(a104,t14,a10,1x,a7,a5,f10.4,f10.2,1x,a1,f7.3,
     $     t14,f2.0,f3.1,a1,2f2.0)
c
        if (comp1 .eq. '     ') comp1='AB   '
	pm=0.
	ra=rah*15. + ram/4.
	dec=decd + decm/60.
	if(dsign .eq. '-') dec=-dec
	corr = (0.0056*sin(ra/degrad)/cos(dec/degrad))
     $     +0.00417*pm*sin(dec/degrad)
c
        iflag=0
        do 600 n=1,norb
        if (wds1 .ne. wdso(n)) go to 600
        if (dd1 .ne. ddo(n)) go to 600
        if (comp1 .ne. compo(n)) go to 600
        call oresid(date,theta,rho,corr,p(n),a(n),xi(n),xn(n),t(n),
     $     e(n),xo(n),dtheta,drho,thetac,rhoc)
	drr=100.*drho/rho
        if (iflag .eq. 0) write(13,905) data,dtheta,drho,drr,oref(n),
     $     thetac,rhoc
  905   format(a104,' ORB ',f9.2,2f10.4,3x,a8,2x,f9.2,f10.3)
        if (iflag .ne. 0) write(13,906) dtheta,drho,drr,oref(n),
     $     thetac,rhoc
  906   format(104x,' ORB ',f9.2,2f10.4,3x,a8,2x,f9.2,f10.3)
        iflag=iflag+1
  600   continue
c
        do 700 n=1,nlin
        if (wds1 .ne. wdsl(n)) go to 700
        if (dd1 .ne. ddl(n)) go to 700
        if (comp1 .ne. compl(n)) go to 700
        call lresid(date,theta,rho,corr,x0(n),xa(n),y0(n),ya(n),t0(n),
     $     dtheta,drho,thetac,rhoc)
	drr=100.*drho/rho
        if (iflag .eq. 0) write(13,907) data,dtheta,drho,drr,lref(n),
     $     thetac,rhoc
  907   format(a104,' LIN ',f9.2,2f10.4,3x,a8,2x,f9.2,f10.3)
        if (iflag .ne. 0) write(13,908) dtheta,drho,drr,lref(n),
     $     thetac,rhoc
  908   format(104x,' LIN ',f9.2,2f10.4,3x,a8,2x,f9.2,f10.3)
        iflag=iflag+1
  700   continue
c
c       if no match, write out data line
c
        if (iflag .eq. 0) write(13,909) data
  909   format(a104,29x,' ')
c
        go to 500
c
  900   write(6,991) norb,wds1
  991   format(i6,2x,a10)
c
  999   stop
        end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine oresid(date,pa,sep,corr,xxp,xxa,xxi,xxn,xxt,xxe,xxo,
     $     respa,ressep,theta,rho)
	degrad=57.2957795
        equ=2000.
	ee=sqrt((1.+xxe)/(1.-xxe))
	aee=xxa*(1.-xxe*xxe)
	xxo2=xxo/degrad
	xxn2=xxn/degrad
	cosi=cos(xxi/degrad)
	xmu=6.2831853/xxp
	em=xmu*(date-xxt)
	pa = pa + corr*(equ - date)
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=xxn2+atan(cosi*tan(xnu+xxo2))
	rho=r*cos(xnu+xxo2)/cos(theta-xxn2)
	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.
	if(respa .lt.  -90.) respa=respa+180.
	if(respa .gt.   90.) respa=respa-180.
	ressep=sep-rho
c
  	return
	end
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine lresid(date,theta,rho,corr,x0,xa,y0,ya,t0,
     $     dtheta,drho,tcalc,rcalc)
	degrad=57.2957795
        equ=2000.
c
	theta2 = theta + corr*(equ - date)
        xobs  =  sin(theta2/degrad) * rho
        yobs  = -cos(theta2/degrad) * rho
c
        xcalc = xa*(date-t0) + x0
        ycalc = ya*(date-t0) + y0
        rcalc = sqrt(xcalc*xcalc + ycalc*ycalc)
        tcalc = atan(-xcalc/ycalc) * degrad
        if (ycalc .gt. 0.) tcalc=tcalc-180.
        if (tcalc .lt.   0.) tcalc=tcalc+360.
        if (tcalc .gt. 360.) tcalc=tcalc-360.
c
        dx=xobs-xcalc
        dy=yobs-ycalc
        drho=rho-rcalc
        dtheta=theta-tcalc
        if (dtheta .lt. -180.) dtheta=dtheta+360.
        if (dtheta .gt.  180.) dtheta=dtheta-360.
        if (dtheta .lt.  -90.) dtheta=dtheta+180.
        if (dtheta .gt.   90.) dtheta=dtheta-180.
    	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
