;
; generate clipping region using scan 10
; this only needs to be run once.. 
; it will generate the parameters used to clip all of the
; other scans.
; you need to  run wedgeinp first
; then run things from this page using cut and paste.
;
zlim=[-.7,47.]
xyrad=150.

restore,'fitiarAll.sav',/verb
; 0 --> 4 param fit
; 1 --> 3 param fit r=265.176
; looks like 4parm fit distributes errors more evening vs range
ii=0
;ii=1
center=fitiar[9,ii].fitcoef[0:2]
R0=fitiar[9,ii].fitcoef[3]
;
; input all of the wedges..
;
.run wedgeinp
; try fitting the sphere.
;
i=0l
xyz=blmscale(dar[i].hdr,*dar[i].pd)
Rm=reform(blmcmpradius(xyz,center=center))
;
; call blmselectpnts to limit the points
hor
ver
a=blmselectpnts(xyz,iikeepT,zlim=zlim,xyrad=xyrad)
ok=intarr(dar[i].npnts)
ok[iikeepT]=1
Rdif=Rm - R0
;
plot,xyz[0,*],xyz[1,*],psym=3
plot,xyz[0,*],xyz[2,*],psym=3
plot,xyz[1,*],xyz[2,*],psym=3
;
;  look at spherical coord
;
rtp=blmxyztosph(xyz,/deg,/cwy)
plot,xyz[0,iikeepT],rdif[iikeepT],psym=3
plot,xyz[1,iikeepT],rdif[iikeepT],psym=3
plot,xyz[2,iikeepT],rdif[iikeepT],psym=3
;
hor
plot,rtp[1,iikeepT],rdif[iikeepT],psym=3
plot,rtp[1,iikeepT],rtp[2,iikeept],psym=3
;
; clip on elevation
el0=3
el1=17
ii=where((rtp[2,*] lt el0) or (rtp[2,*] gt el1),cnt)  
plot,xyz[0,ii],rdif[ii],psym=3
plot,xyz[1,ii],rdif[ii],psym=3
plot,xyz[2,ii],rdif[ii],psym=3
ok[ii]=0 
;
; limit y to  gt 35 m .. 1 km mode didn't work well close in
y0=36.
ii=where(xyz[1,*] lt y0,cnt)
ok[ii]=0


ii=where(ok eq 1,cnt)
;
plot,xyz[0,ii],rdif[ii],psym=3
plot,xyz[1,ii],rdif[ii],psym=3
plot,xyz[2,ii],rdif[ii],psym=3
plot,rtp[0,ii],rdif[ii],psym=3
oplot,[0,1000],[0,0],col=colph[2]
plot,rtp[1,ii],rdif[ii],psym=3
oplot,[-1000,1000],[0,0],col=colph[2]
plot,rtp[2,ii],rdif[ii],psym=3
oplot,[0,1000],[0,0],col=colph[2]
;
;
; look at 
; el 6-9
ok1=ok
ii=where((rtp[2,*] lt 6) or (rtp[2,*] gt 9),cnt)
ok1[ii]=0
;
;
; keep az -5.7 to 4
ii=where((rtp[1,*] lt -5.7) or (rtp[1,*] gt 4),cnt)
ok1[ii]=0
ii=where(ok1 eq 1,cnt)
;
hor
ver
plot,xyz[0,ii],xyz[1,ii],psym=3
plot,xyz[0,ii],xyz[2,ii],psym=3
;
; keep x=-6 to x=+4  so we have a rectangular region of data
; for smoothing
;
ii=where((xyz[0,*] lt -6) or (xyz[0,*] gt 4),cnt)
ok1[ii]=0
;
ii=where(ok1 eq 1,cnt)
;
xyzt=xyz[*,ii]
rtpT=rtp[*,ii]
rdift=rdif[ii]
plot,xyzt[0,*],xyzt[1,*],psym=3
maxx=max(xyzt[0,*],min=minx)
dx=maxx -minx
maxy=max(xyzt[1,*],min=miny)
dy=maxy -miny
;
stp=.3d
nx=long(dx/stp + .5)
ny=long(dy/stp + .5)
cntAr=lonarr(nx,ny)
rmsAr=dblarr(nx,ny)
avgAr=dblarr(nx,ny)
x0=minx*1d
y0=miny*1D
for ix=0,nx-1 do begin &$
	print,ix &$
	xs=x0 + ix*stp	 &$
    xe=xs + stp &$
	iix=where((xyzt[0,*] ge xs) and (xyzt[0,*] lt xe),cnt) &$
	for iy=0,ny-1 do begin &$
		ys=y0+iy*stp &$
		ye=ys + stp &$
	    ii=where((xyzt[1,iix] ge ys) and (xyzt[1,iix] lt ye),cnt) &$
		cntAr[ix,iy]=cnt &$
		if (cnt gt 2) then 	begin &$
			a=rms(rdift[iix[ii]],/quiet)	 &$
			avgAr[ix,iy]=a[0] &$
			rmsAr[ix,iy]=a[1] &$
		endif &$
	endfor &$
endfor
plot,cntar,psym=1
;
xx=dindgen(nx)*stp + x0
yy=dindgen(ny)*stp + y0
;
stripsxy,xx,cntar,0,0,/step
ver,0,.005
sym=1
stripsxy,xx,rmsar,0,0,/step ,psym =sym
stripsxy,yy,transpose(rmsar),0,0,/step
stripsxy,xx,avgar,0,0,/step
plot,avgar
end