program grader1
c
c	The program reads a list of binary star data (separation   
c	and position angle) and a set of orbital elements (from 
c	neworb4.new), then determines O-C residuals and other 
c	statistics. Purpose is to determine orbit grade criterion.
c	             [ orbgrades.5* --> orbgrades.6* ]
c
	common       pm,p,a,xi,xnode,t0,ecc,omega,ee,aee,cosi,xmu
	character*1  decs
	character*13 ref
	character*20 infile, outfile
	character*22 star
	dimension date(2000),theta(2000),rho(2000),wt(2000),
     $     diff(2000),theta2(2000),phase(2000)
c
	degrad=0.0174533
	equ=2000.0
c
	write(6,931)
  931	format(' Enter input file: ',$)
  	read(5,932) infile
  932	format(a20)
  	write(6,933)
  933	format(' Enter output file: ',$)
  	read(5,932) outfile
c
	open(10,file=infile, status='UNKNOWN')
	open(11,file=outfile,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(a22,t1,f2.0,f3.1,a1,2f2.0,t63,f4.3)
c	iflag=iflag+1
	read(10,904) p,a,xi,xnode2,t0,ecc,omega2,ref,dlast
  904	format(f15.6,f10.7,f8.5,f9.5,f11.6,f8.6,f11.7,a13,t81,f4.0)
c
	ee=sqrt((1.+ecc)/(1.-ecc))
	aee=a*(1.-ecc*ecc)
	omega=omega2*degrad
	cosi=cos(xi*degrad)
	xnode=xnode2*degrad
	xmu=6.2831853/p
c
	ra=15.*rah + ram/4.
	dec=decd + decm/60.
	if(decs .eq. '-') dec=-dec
 	corr = (0.0056*sin(ra*degrad)/cos(dec*degrad))
     $     + 0.00417*pm*sin(dec*degrad)
c
c	read data points, determine residuals
c
	npt=0
  200	read(10,903,end=300) bess,tobs,robs,weight
  903	format(t7,f9.4,t20,f6.2,t34,f9.6,t54,f4.1)
	if (bess .eq. 0.0) go to 300
	npt=npt+1
	tobs=tobs + corr*(equ-bess)
	theta(npt)=tobs
	rho(npt)=robs
	wt(npt)=weight
	date(npt)=bess
	call porbit(bess,tobs,robs,tres,rres,rresr,xyres)
	diff(npt)=xyres
	go to 200
c
c	determine weighted rms residuals, throw out outliers and recompute
c
  300	call rms(npt,wt,diff,sumwt,rmsdiff)
  	do 320 n=1,npt
  320	if (diff(n) .gt. 3.0*rmsdiff) wt(n)=0.
  	call rms(npt,wt,diff,sumwt2,rmsdiff2)
c
c	simple one - determine number of revolutions covered
c	up to date of orbit publication
c
	if (date(npt) .lt. dlast) dlast=date(npt)
	enrev=(dlast-date(1))/p
c
c	theta and phase coverage: extract theta and phase
c	values for observations up to publication date
c
	npt2=0
	do 350 n=1,npt
	if (wt(n) .eq. 0.) go to 350
	if (date(n) .gt. dlast) go to 350
	npt2=npt2+1
	theta2(npt2)=theta(n)
	phase(npt2)=(date(n)-t0)/p + 1000.
	phase(npt2)=phase(npt2)-float(int(phase(npt2)))
  350	continue
c
c	sort on theta, determine rms difference in
c	theta between successive observations
c	also determine maximum gap in theta
c
	do 400 n=1,npt2-1
	do 400 l=n+1,npt2
	if (theta2(l) .ge. theta2(n)) go to 400
	temp=theta2(l)
	theta2(l)=theta2(n)
	theta2(n)=temp
  400	continue
c		
	dth=0.
	dthmax=0.
	do 420 n=1,npt2-1
	tdiff=abs(theta2(n+1)-theta2(n))
  	dth=dth+tdiff*tdiff
	if (tdiff .gt. dthmax) dthmax=tdiff
  420	continue
	tdiff=abs(theta2(1)+360.-theta2(npt2))
	dth=dth+tdiff*tdiff
	if (tdiff .gt. dthmax) dthmax=tdiff
	dth=sqrt(dth/float(npt2))
c
c	sort on phase, determine rms difference in
c	phase between successive observations
c	also determine maximum gap in phase
c
	do 500 n=1,npt2-1
	do 500 l=n+1,npt2
	if (phase(l) .ge. phase(n)) go to 500
	temp=phase(l)
	phase(l)=phase(n)
	phase(n)=temp
  500	continue
c		
	dph=0.
	dphmax=0.
	do 520 n=1,npt2-1
	pdiff=abs(phase(n+1)-phase(n))
  	dph=dph+pdiff*pdiff
	if (pdiff .gt. dphmax) dphmax=pdiff
  520	continue
	pdiff=abs(phase(1)+1.-phase(npt2))
	if (pdiff .gt. dphmax) dphmax=pdiff
	dph=dph+pdiff*pdiff
	dph=sqrt(dph/float(npt2))
	dph=dph*360.
	dphmax=dphmax*360.
c
	write(11,902) star,p,a,xi,xnode2,t0,ecc,omega2,ref,
     $     npt2,rmsdiff2,rmsdiff2/a,dth,dph,dthmax,dphmax,enrev 
  902	format(a22,f11.4,f8.4,2f7.2,f10.4,f6.3,f7.2,a13,i5,
     $     2f8.4,4f7.2,f9.3)
c
  	go to 100
  999	stop
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine porbit(bess,tobs,robs,tres,rres,rresr,xyres)
	common pm,p,a,xi,xnode,t0,ecc,omega,ee,aee,cosi,xmu
	degrad=57.2957795
	em=xmu*(bess-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))
	tcalc=xnode+atan(cosi*tan(xnu+omega))
	rcalc=r*cos(xnu+omega)/cos(tcalc-xnode)
	tcalc=tcalc*degrad
	if(rcalc .gt. 0.) go to 300
	rcalc=-rcalc
	tcalc=tcalc+180.
  300	if(tcalc .lt.   0.) tcalc=tcalc+360.
	if(tcalc .gt. 360.) tcalc=tcalc-360.
c
	xcalc=-sin((tcalc+180.)/degrad)*rcalc
	ycalc= cos((tcalc+180.)/degrad)*rcalc
	xobs =-sin((tobs+180.)/degrad)*robs
	yobs = cos((tobs+180.)/degrad)*robs
	xres = xobs-xcalc
	yres = yobs-ycalc
	xyres= sqrt(xres*xres+yres*yres)
c
	tres=tobs-tcalc
	if(tres .lt. -180.) tres=tres+360.
	if(tres .gt.  180.) tres=tres-360.
	rres=robs-rcalc
	rresr=rres/rcalc
c
  400	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
	subroutine rms(npt,wt,diff,sumwt,rmsdiff)
	dimension wt(2000),diff(2000)
	sumwt=0.
	sumdiff=0.
	do 100 n=1,npt
	sumwt=sumwt+wt(n)
	sumdiff=sumdiff+wt(n)*diff(n)*diff(n)
  100	continue
  	rmsdiff=sqrt(sumdiff/sumwt)
  	return
  	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc