program resid2
c
c	The program reads a list of binary star data (separation and 
c	position angle) and a set of orbital elements, then determines
c	O-C residuals.
c
	common       pm,p,a,xi,xnode,t0,ecc,omega
	character*1  elt,decs
	character*3  code1
	character*5  apa,asep
	character*4  code2
	character*5  code3
	character*80 star
c
	degrad=0.0174533
	pm=0.0
	equ=2000.0
c
c	enter file names for I/O
c
	open(10,file='calib.plt',status='UNKNOWN')
	open(11,file='calib.resids',status='UNKNOWN')
c
c	read star name and orbital elements
c
  100	read(10,901,end=999,err=999) star,rah,ram,decs,decd,decm,pm
  901	format(a80,t1,f2.0,f3.1,a1,2f2.0,t63,f4.3)
	write(11,908) star
  908	format(a80)
	call orbit(rx,ry,norb)
	ra=15.*rah + ram/4.
	dec=decd + decm/60.
	if(decs .eq. '-') dec=-dec
	corr=0.
	if(pm .ne. 0.0) corr = (0.0056*sin(ra*degrad)/cos(dec*degrad))
     $                       + 0.00417*pm*sin(dec*degrad)
c
c	read data points, determine residuals
c
  150	write(11,904)
  904	format(/4x,'Date',10x,'Codes',19x,'Observed',9x,
     $     'Calculated',12x,'O-C (theta, rho, rho%, x, y)')
c
c
  200	read(10,903,end=350) bess,pa,elt,sep,code1,code2,code3,
     $     apa,asep
  903	format(t7,f9.4,t20,f6.2,t34,a1,f6.4,t50,a3,1x,a4,6x,a5,
     $     t20,a5,t35,a5)
	if (bess .eq. 0.0) go to 350
	if(elt .eq. '<') go to 200
	if(apa .eq. ' ' .or. asep .eq. ' ') go to 200
	npt=npt+1
	pa=pa + corr*(equ-bess)
	call porbit(bess,pa,sep,xneg,yneg,code1,code2,code3)
	go to 200
c
  350	write(11,907)
  907	format(80x)
  	go to 100
  999	stop
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine orbit(x,y,norb)
	common pm,p,a,xi,xnode,t0,ecc,omega
	degrad=57.2957795
	read(10,*) p,a,xi,xnode,t0,ecc,omega
c	write(11,901) p,a,xi,xnode,t0,ecc,omega,pm
  901	format('    P = ',f12.4,'   a = ',f12.4,
     $         '    i = ',f12.4,'   ln = ',f12.4/
     $         '    T = ',f12.4,'   e = ',f12.4,
     $         '   lp = ',f12.4,'   pm = ',f12.4)
	ee=sqrt((1.+ecc)/(1.-ecc))
	aee=a*(1.-ecc*ecc)
	omega=omega/degrad
	cosi=cos(xi/degrad)
	xnode=xnode/degrad
	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine porbit(date,pa,sep,xneg,yneg,code1,code2,code3)
	common pm,p,a,xi,xnode,t0,ecc,omega
	character*3 code1
	character*4 code2
	character*5 code3
	degrad=57.2957795
	ee=sqrt((1.+ecc)/(1.-ecc))
	aee=a*(1.-ecc*ecc)
	cosi=cos(xi/degrad)
	xmu=6.2831853/p
	em=xmu*(date-t0)
c
	e0=em+ecc*sin(em)+0.5*ecc*ecc*sin(2.*em)
	do 200 l=1,10
	em0=e0-ecc*sin(e0)
	e0=e0+(em-em0)/(1.-ecc*cos(e0))
  200	continue
c
	xnu=2.*atan(ee*tan(0.5*e0))
	r=aee/(1.+ecc*cos(xnu))
	theta=xnode+atan(cosi*tan(xnu+omega))
	rho=r*cos(xnu+omega)/cos(theta-xnode)
	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
	xneg=-sin((theta+180.)/degrad)*rho
	yneg= cos((theta+180.)/degrad)*rho
	xobs=-sin((pa+180.)/degrad)*sep
	yobs= cos((pa+180.)/degrad)*sep
	resx=xobs-xneg
	resy=yobs-yneg
c
	respa=pa-theta
	if(respa .lt. -180.) respa=respa+360.
	if(respa .gt.  180.) respa=respa-360.
	ressep=sep-rho
	perres=ressep/rho
	if(sep .eq. 0.0) go to 400
c
	write(11,901) date,code1,code2,code3,pa,sep,
     $     theta,rho,respa,ressep,perres,resx,resy
  901	format(1x,f9.4,2x,a3,2x,a4,2x,a5,3(f10.2,f8.4),3f8.4)
  400	return
  902	format(1x,f10.4,2x,a2,1x,a12,18x,f10.2,f8.3)
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc