program orbstats1
c
c	checks a few possible statistics for grading
c	orbits: number of orbital revolutions, density 
c	of the coverage
c
	character*90 title
	dimension date(1999),the(1999),rho(1999),wt(1999)
c
	open(10,file='orbtest.3',status='UNKNOWN')
	open(11,file='orbstats1.out',status='UNKNOWN')
c
  100	read(10,901,end=999) title
  901	format(a90)
c  	write(11,901) title
	read(10,902) title,p,a,xi,xnode,t0,ecc,omega,
     $     igrade,refcode
  902	format(a90,t1,f15.6,f10.7,f9.5,f7.3,2x,f11.6,f8.6,
     $     f9.5,3x,i1,2x,a8)
c  	write(11,901) title
c
	nobs=1
  200	read(10,903) date(nobs),the(nobs),rho(nobs),wt(nobs)
  903	format(t7,f9.4,t20,f6.2,t34,f7.4,t54,f4.1)
 	if (date(nobs) .eq. 0.0) go to 300
 	nobs=nobs+1
	if (nobs .gt. 1999) write(6,991) nobs
  991	format(' nobs =',i4)
 	go to 200
c
  300	nobs=nobs-1
  	xorbits=(date(nobs)-date(1))/p
c
c	sort observations in order of theta
c
	do 400 n=1,nobs-1
	do 400 l=n+1,nobs
	if (the(l) .ge. the(n)) go to 400
	ttemp=the(l)
	the(l)=the(n)
	the(n)=ttemp
	dtemp=date(l)
	date(l)=date(n)
	date(n)=dtemp
	rtemp=rho(l)
	rho(l)=rho(n)
	rho(n)=rtemp
	wtemp=wt(l)
	wt(l)=wt(n)
	wt(n)=wtemp
  400	continue
c
	trms=0.
	rrms=0.
c	write(6,991) nobs
	do 500 n=2,nobs
	trms=trms+(the(n)-the(n-1))*(the(n)-the(n-1))
	rrms=rrms+(rho(n)-rho(n-1))*(rho(n)-rho(n-1))
  500	continue
  	trms=trms+(the(1)-the(nobs)+360.)
     $     *(the(1)-the(nobs)+360.)
  	rrms=rrms+(rho(1)-rho(nobs))*(rho(1)-rho(nobs))
	trms=sqrt(trms/real(nobs))
	rrms=sqrt(rrms/real(nobs))/a
c
	write(11,904) p,a,xi,xnode,t0,ecc,omega,igrade,
     $     nobs,xorbits,trms,rrms
  904	format(f10.4,f10.4,f9.3,f9.3,f11.3,f9.5,
     $     f9.3,3x,i1,i5,f7.2,f12.2,f12.6)
c
	go to 100
  999	stop
  	end