program mean
c
c	determine mean drho and theta values versus date and theta
c
	dimension date(5045),theta(5045),rho(5045),dtheta(5045),
     $     drho(5045),nth(36),tdt(36),tdr(36),
     $     ddt(5045),ddr(5045),ddate(5045)
	real mdate,mdrh,mdth
c
	open(10,file='resids.wsi',status='UNKNOWN')
	open(11,file='meanresids.wsi',status='UNKNOWN')
c
c	read in residuals, discard outliers
c
	nmeas=0
	dtall=0.
	drall=0.
  100	read(10,901,end=150) date1,theta1,rho1,dtheta1,drho1
  901	format(f10.4,t76,f6.2,f9.4,f10.2,t110,f9.4)
	if (abs(dtheta1) .gt. 5.) go to 100
	if (abs(drho1) .gt. 0.2) go to 100
	nmeas=nmeas+1
	date(nmeas)=date1
	theta(nmeas)=theta1
	rho(nmeas)=rho1
	dtheta(nmeas)=dtheta1
	drho(nmeas)=drho1
	dtall=dtall+dtheta1
	drall=drall+drho1
	go to 100
c
  150	dtall=dtall/float(nmeas)
	drall=drall/float(nmeas)
	write(6,991) nmeas
  991	format(' nmeas (after removing outliers) = ',i6)	
	write(6,992) dtall,drall
  992	format(' overall mean dtheta and drho = ',f10.2,f10.4)	
c
c	determine mean theta and rho residuals versus theta. Bin in 10deg bins
c
	do 200 n=1,36
	nth(n)=0
	tdt(n)=0.
	tdr(n)=0.
  200	continue
c
	do 300 n=1,nmeas
	nn=int(theta(n)/10.)+1
	nth(nn)=nth(nn)+1
	tdt(nn)=tdt(nn)+dtheta(n)
	tdr(nn)=tdr(nn)+drho(n)
  300	continue
c
	do 400 n=1,36
	theta1=float(n*10)-5.
	tdt(n)=tdt(n)/float(nth(n))
	tdr(n)=tdr(n)/float(nth(n))
	write(11,902) theta1,tdt(n),tdr(n),nth(n)
  902	format(f5.0,f10.2,f9.4,i6)
  400	continue
c
c	determine running theta and rho residuals versus date, Bin in 50-measure bins
c
	dsum=0.
	tsum=0.
	rdum=0.
	do 500 n=1,50
	dsum=dsum+date(n)
	tsum=tsum+dtheta(n)
	rsum=rsum+drho(n)
  500	continue
	date1=dsum/50.
	dtheta1=tsum/50.
	drho1=rsum/50.
	write(11,903) date1,dtheta1,drho1
  903	format(f10.4,f10.2,f9.4)
c
	do 600 n=51,nmeas
	dsum=dsum+date(n)-date(n-50)
	tsum=tsum+dtheta(n)-dtheta(n-50)
	rsum=rsum+drho(n)-drho(n-50)
	date1=dsum/50.
	dtheta1=tsum/50.
	drho1=rsum/50.
	write(11,903) date1,dtheta1,drho1
  600	continue
c
  999	stop
	end