; loop fitting a sphere to the data
;1. inpscan19 .. inputs data and clips using blmselectpnts.
;   51297801  before any clipping
;   17364752  after  using the clipping below
;   - keep points if:
;   - +/- 150 meters in x and y
;   -     150 meters in xy radius
;   -  -.7 to 47 meters in z
;
;2. when fitting sphere loop 10 times
;   - after each loop take difRadius=(measuredRadius - radius fit) for each point
;     measured radius is length of ([x0,y0,z0] - measuredpoint)j
;   - compute the rms
;   - and keep all points within 3 sigma .. 
;;  - then iterate
;3. do the above fitting for x0,y0,z0 and radius
;   then repeat fitting for x0,y0,z0 and keeping radius fixed at 265.176 meters 
;
; save files generated:
; fitiarall.sav .. fitiar[10,2] {}
; iikeep_scan19.sav.. iikeepsav,iikeepf3,iikeepf4,nsig,toloop
;  iikeepsav d[51e6]->xyz[17e6] via blmselectpnts
;  iikeepf4  xyz[*,iikeepf4] points after 10 iterations fit 4 params
;  iikeepf3  xyz[*,iikeepf3] points after 10 iterations fit 3 params, radius fixed at 265.176
; 
; 
useWeights=1
ao9todish=.7   ; the starting offset p50 to dish
;         x0 y0  z0                radius
a=size(xyz)
ntot=a[2]
; start with all
initCoef=[0d,0d,265.176- ao9ToDish,265.176d]
; j=0 fitall, 
; j=1 fixed radius
; on each loop keep points within 3 sigma then refit
; in case this is run from upper dir
savDir='/share/megs/phil/x101/p50/200311/fitsphere/'
; need the radius for computing the weights
rtp=blmxyztosph(xyz,/deg)
!p.multi=0
hor
ver
nsig=3.
toloop = 10
start=1
radxy=sqrt(xyz[0,*]^2 + xyz[1,*]^2)
for j=0,1 do begin
	fitA=(j eq 0)?[1,1,1,1]:[1,1,1,0]
	iikeep1=lindgen(ntot)
; 	fitradDif= FitRadius - pointRadius
	for iloop=0,toloop-1 do begin
		istat=fitsphere_cmpw(radxy[0,iikeep1],w)
		if useWeights eq 0 then begin
			w=w*0 + 1.
	    endif
	    print,"  computing fit " + systime(0)
		istat=blmfitsphere(xyz[*,iikeep1],fitA=fitA,fitIcf,initCoef=initCoef,w=w,fitRadDif=fitRadDif)
	    print,"  computing fit done " + systime(0)
		if start then begin
			fitIar=replicate(fiticf,toloop,2)
			start=0
		endif
		print,j,iloop,fiticf.npnts
        print,'   ',fiticf.fitcoef
		print,'   ',fiticf.sigmas
		print,'   ',fiticf.fitErr
		sig=fiticf.fitErr
		fitiar[iloop,j]=fiticf
		ii=where(abs(fitraddif) lt (nsig*sig),cntg)
		jj=where(abs(fitraddif) gt (nsig*sig),cntb)
		cnttot=n_elements(iikeep1)
		tit=string(format='("j=",i," iloop:",i," ntot:",i," badpnt:",i)',$
				j,iloop,cnttot,cntb)

		plot,xyz[0,iikeep1[jj]],xyz[2,iikeep1[jj]],psym=3,tit=tit,$
				xtit='xaxis',ytit='zaxis'
		empty
		iikeep1=iikeep1[ii]
	endfor
	if j eq 0 then iikeepf4=iikeep1
	if j eq 1 then iikeepf3=iikeep1
endfor
if useWeights eq 0 then begin
save,fitiar,toloop,nsig,file=savDir + 'fitiarAll_nw.sav'
save,iikeepsav,iikeepf4,iikeepf3,toloop,nsig,file=savDir+'iikeep_scan19_nw.sav'
endif else begin
save,fitiar,toloop,nsig,file=savDir + 'fitiarAll.sav'
save,iikeepsav,iikeepf4,iikeepf3,toloop,nsig,file=savDir+'iikeep_scan19.sav'
endelse
end