program resids
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  rflag,decs
	character*5  apa,asep,abess
	character*8  code1
	character*35 code2
	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='resids.in2',status='UNKNOWN')
	open(11,file='resids.out2',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
  901	format(a80,t1,f2.0,f3.1,a1,2f2.0)
	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(/1x,'Date.....',2x,'Codes (method, wt, n, ref, etc).',
     $     '............',4x,'Observed........',4x,'Calculated......',
     $     4x,'O-C (dtheta, drho, drho/rho, dx, dy).......')
c
  200	read(10,903,end=350) bess,pa,rflag,sep,code1,code2,
     $     abess,apa,asep
  903	format(t7,f10.5,t20,f7.3,t32,a1,f9.5,t50,a8,5x,a35,
     $     t7,a5,t20,a5,t35,a5)
	if (abess .eq. '     ') go to 350
	if (rflag  .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,rflag,sep,xneg,yneg,code1,code2)
	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)
	character*100 data
	character*8 refcode
        character*1 pcode,acode,tcode
	common pm,p,a,xi,xnode,t0,ecc,omega
	degrad=57.2957795
	read(10,901) data,p,pcode,a,acode,xi,xnode,t0,tcode,ecc,
     $     omega,refcode
  901	format(a100,t3,f12.6,a1,f10.5,a1,2f9.4,f13.6,a1,f9.6,f9.4,2x,a8)
     	write(11,902) data
  902	format(a100)
        if (pcode .eq. 'c') p=p*100.
        if (pcode .eq. 'd') p=p/365.2421987
	if (pcode .eq. 'h') p=p/365.2421987/24.
	if (pcode .eq. 'm') p=p/365.2421987/24./60.
        if (acode .eq. 'm') a=a/1000.
        if (acode .eq. 'M') a=a*60.
        if (tcode .eq. 'c') t0=t0*100.
        if (tcode .eq. 'd') t0=(t0-33282.423)/365.2421987+1950.
        if (tcode .eq. 'm') t0=(t0+0.5-33282.423)/365.2421987+1950.
	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,rflag,sep,xneg,yneg,code1,code2)
	common pm,p,a,xi,xnode,t0,ecc,omega
	character*8 code1
	character*35 code2
        character*1 rflag
c
	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
        sep2=sep
        if (rflag .eq. 'D') sep2=sep*3600.
        if (rflag .eq. 'M') sep2=sep*60.
        if (rflag .eq. 'm') sep2=sep/1000.
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)
        if (rflag .eq. 'D') rho=rho/3600.
        if (rflag .eq. 'M') rho=rho/60.
        if (rflag .eq. 'm') rho=rho*1000.
	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
c       remove 180deg differences between observed and calculated angles
c
        if ((theta-pa) .gt.  90.) pa=pa+180.
        if ((theta-pa) .gt.  90.) pa=pa+180.
        if ((theta-pa) .lt. -90.) pa=pa-180.
        if ((theta-pa) .lt. -90.) pa=pa-180.
        if (pa .gt. 360.) pa=pa-360.
        if (pa .lt.   0.) pa=pa+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,pa,rflag,sep,
     $     theta,rflag,rho,respa,rflag,ressep,perres,resx,resy
  901	format(1x,f9.4,2x,a8,1x,a35,3(f10.2,1x,a1,f8.4),3f9.4)
  400	return
  902	format(1x,f10.4,2x,a2,1x,a12,18x,f10.2,f8.3)
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc