;
; analyzing the errors.
; make sure we have the errors computed
; th will rotate so y points north
thRotD=-129.25996D
mtocm=100.
ftom=.3048d
mtof=1d/ftom
cent=center3
radius=designRadiusM
goto,doit
;
; throw out points beyond 10cm
;
rtp=blmxyztosph(xyz,cent=center3,/deg)
radErrCm=(reform(rtp[0,*]) -radius)*100.
ii=where(abs(radErrCm) gt 10,cnt)
plot,xyz[0,ii],xyz[2,ii],psym=3
ii=where(abs(radErrCm) lt 10,cnt)
xyzr=xyz[*,ii]
; rotate to north.. xyz, and center
xyzR=rotvec(xyzR,thRotD,axis=3,/usedouble)
centR=rotvec(center3,thRotD,axis=3,/usedouble)
radErrCm=cmperrradius(xyzr,radius,cenDef=centR)
; now plot points > 2cm
; red is too long, blue is too low
;
hor,-500,500
ver,-500,500
lim=2.
ii=where(abs(radErrcm) gt lim,cnt)
;plot,[0,1],[0,1],/nodata,
ii=where(abs(radErrCm) le 10,cnt)
xyzR=xyzr[*,ii]
radERrCm=radErrCm[ii]
;
; now compute the errors for a number of radii
;
doit:
radToUseAr=[0.,-.03,-.07] + designRadiusM
th=2. 
charth=2
cs=1.5
colb='ffffff'xul ; white
col ='000000'xul ; black
errLimCm=2.
;
for i=0,n_elements(radToUseAr) - 1 do begin
	radToUse=radToUseAr[i]
	radErrCm=cmperrradius(xyzr,radToUse,cenDef=centR)
hor,-500,500
ver,-500,500
ii=where(abs(radErrCm) gt errLimCm,cnt)
xtit='west - X - east  [ft]'
ytit='Y ft'
lin=string(format='("Radial errors [cm] > ",f4.1," cm (using Radius = ",f8.3," m)")',errLimCm,radToUse)
tit=lin

plot,[0,1],[0,1],/nodata,background=colb,color=col,/iso,$
	chars=cs,xthick=th,ythick=th,charth=charth,$
	xtit=xtit,ytit=ytit,$
	tit=tit
ii1=where(radErrCm[ii] gt 0)
oplot,xyzr[0,ii[ii1]]*mtof,xyzr[1,ii[ii1]]*mtof,psym=3,col=colph[2]
ii2=where(radErrCm[ii] lt 0)
oplot,xyzr[0,ii[ii2]]*mtof,xyzr[1,ii[ii2]]*mtof,psym=3,col=colph[4]
;
; import png:radErrsStd.png
	print,"import png:radErrs_dr_xcm.png"
	key=checkkey(/wait)
endfor
stop


; look at the rotated offset
;rr=dist3([cent[0],cent[1],0d])*mtof*12
;drawcircle,0,0.,rr,color=2
;oplot,[0,cent[0]]*mtof*12.,[0,cent[1]]*mtof*12.
;oplot,[0,centr[0]]*mtof*12,[0,centr[1]]*mtof*12,col=colph[3]
;
rtprf=blmxyztosph(xyzrf,cent=centrf,/deg)
; now rotate az to align with north
plot,xyzrf[0,*],xyzrf[1,*],psym=3
; meausured - ref
radErrRCm=reform((rtprf[0,*]*ftom-(radius))*mtocm)
radErrRCm=reform(radErrRCm)
;
nn=n_elements(raderrRcm)
cliplev=10.
ii=where(abs(raderrRcm) gt cliplev,cnt)
print,'10cm:',cnt,cnt/(nn*1.)
!p.multi=[0,1,2]
hor,-500,500
ver,-500,500
plot,xyzrf[0,ii],xyzrf[1,ii],psym=3,/iso
ver,-1,47
plot,xyzrf[0,ii],xyzrf[2,ii],psym=3
;
; look at the points we are excluding
;
ii=where(abs(raderrRcm) lt cliplev,cnt)
clipLevcm=1.8 
print,cnt
!p.multi=0
hor,-500,500
ver,-500,500
jj=where(abs(radErrRCm[ii]) gt cliplevcm,cnt)
kk=where(radErrRcm[ii[jj]] gt 0,cnth)
colb='ffffff'xul
col ='000000'xul
plot,xyzrf[0,ii[jj[kk]]],xyzrf[1,ii[jj[kk]]],psym=0,/nodata,background=colb,color=col &$
oplot,xyzrf[0,ii[jj[kk]]],xyzrf[1,ii[jj[kk]]],psym=3,col=colph[2]
kk=where(radErrRcm[ii[jj]] lt 0,cntl)
oplot,xyzrf[0,ii[jj[kk]]],xyzrf[1,ii[jj[kk]]],psym=3,col=colph[4]
print,cnt,cnth,cntl


plot,[0,cent[0]],[0,cent[1]]
end