PROGRAM RESID
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  ORB,ELT,decs
	CHARACTER*3  CODE1
	CHARACTER*5  APA,ASEP
	CHARACTER*12 CODE2
	CHARACTER*20 file01,file02
	character*30 star
C
	DEGRAD=0.0174533
	PM=0.0
	EQU=2000.0
C	Enter file names for i/o
C
	write(6,931)
  931	format(/' Enter input data file name: ',$)
	read(5,932) file01
  932	format(a20)
	open(10,file=file01)
	write(6,933)
  933	format(' Enter output file name: ',$)
	read(5,932) file02
	open(11,file=file02)
C
C	READ STAR NAME AND ORB. IF ORB=YES THEN READ ORBITAL ELEMENTS
C
	READ(10,901) STAR,RAH,RAM,DECS,DECD,DECM
  901	FORMAT(a30,T56,F2.0,F3.1,A1,2F2.0)
	READ(10,902) ORB
  902	FORMAT(A1)
	IF(ORB .NE. 'Y') GO TO 150
	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,906) STAR,RAH,RAM,DECS,DECD,DECM,PM
  906	FORMAT(1X,a30,10X,F3.0,F5.1,1X,A1,F3.0,F4.0,F9.4)
	WRITE(11,904)
  904	FORMAT(/4X,'Date',6X,'Codes',13X,'Observed',8X,
     $     'Calculated',10X,'O-C (theta, rho, rho%, x, y)')
C
C
  200	READ(10,903,END=350) BESS,PA,ELT,SEP,CODE1,CODE2,APA,ASEP
  903	FORMAT(T7,F9.4,T20,F6.2,T34,A1,F6.4,T50,A3,1X,A12,T20,A5,T35,A5)
	IF(ELT .EQ. '<') GO TO 200
	IF(APA .EQ. ' ' .OR. ASEP .EQ. ' ') GO TO 200
	NPT=NPT+1
	IF(ORB .EQ. 'Y') PA=PA + CORR*(EQU-BESS)
	IF(ORB .EQ. 'Y') CALL PORBIT(STAR,BESS,PA,SEP,XNEG,YNEG,
     $     CODE1,CODE2)
	GO TO 200
C
  350	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,PM
	WRITE(11,901) P,A,XI,XNODE,T0,ECC,OMEGA,PM
  901	FORMAT(/' Orbital elements are as follows:'/
     $     '    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(STAR,DATE,PA,SEP,XNEG,YNEG,CODE1,CODE2)
	COMMON PM,P,A,XI,XNODE,T0,ECC,OMEGA
	CHARACTER*30 star
	CHARACTER*3  CODE1
	CHARACTER*12 CODE2
	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,PA,SEP,
     $     THETA,RHO,RESPA,RESSEP,PERRES,RESX,RESY
  901	FORMAT(1X,F9.4,2X,A3,1X,A12,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