;
; look at the results from the fit
; restore,'fitiarAll.sav',/verb
; retstore,'iikeep_scan19.sav',/verb .. if you need the indices
;
; d[alldata] {}
; iikeepsav   d->xyz
; intensav[17e6]
; xyz[3,17e6] clipped in x,y,z, xyrad
; iikeepf4  apply to xyz fits for everything
; iikeepf3  apply to xyz, keep radius fixed.
;
; --> this file also generates the .png files 
; radialerr_pnts.png  radialerrs3.png  radialerrs4.png  xyraw_clip.png  xyusage_f4.png  xzraw_clip.png
; via screen shots
;
; plot how the fit progressed.
;
;goto,plotit1
goto,plothist
;
plotit1:
if hard then pscol,'fitsphere_plt.ps',/full
!p.multi=[0,1,2]
font=1
cs=1.4
csn=1.5
sym=-1
hor,0,11
nloop=n_elements(fitiar[*,0])
x=findgen(nloop) + 1
xtit='Fit iteration'
ytit='Number of points used'
tit='200311 p50 scan19. Iterate fitting sphere to data. Numpnts vs iteration.'
ver,14e6,16.5e6
plot,x,fitiar[*,0].npnts,psym=sym,$
	chars=cs,font=font,$
	xtit=xtit,ytit=ytit,$
	tit=tit
oplot,x,fitiar[*,1].npnts,psym=sym,col=colph[2] 
ln=2.5
xp=.3
scl=.7
note,ln      ,"Throw out points above 3 sigma after each iteration",chars=csn,font=font
xp=.04
ln=9
note,ln+0*scl,'fit x0,y0,z0,radius',chars=csn,font=font,xp=xp
note,ln+1*scl,'fit x0,y0,z0. Fixed radius=265.176 (m)',chars=csn,font=font,xp=xp,col=colph[2]
;
ver,.001,1
ytit='fit Sigma [meters]'
tit='fit sigma vs iteration'
plot,x,fitiar[*,0].fiterr,psym=sym,/ylog,$
	chars=cs,font=font,$
	xtit=xtit,ytit=ytit,$
	tit=tit
oplot,x,fitiar[*,1].fiterr,psym=sym,col=colph[2]
; ------------------------------------------------
; plot x0,y0,z0,radius
;
!p.multi=[0,1,3]
cs=2.2
sym1=-1
sym2=-2
ver,.01,.03
ytit='Meters'
tit='X0,Y0 fit coef vs fit iteration (offsets from scanner center)'
plot,x,fitiar[*,0].fitcoef[0],psym=sym1,$
	chars=cs,font=font,$
	xtit=xtit,ytit=ytit,$
	tit=tit
oplot,x,fitiar[*,1].fitcoef[0],psym=sym1,col=colph[2]
oplot,x,fitiar[*,0].fitcoef[1],psym=sym2
oplot,x,fitiar[*,1].fitcoef[1],psym=sym2,col=colph[2]
xp=.3 
ln=6.5
note,ln,"Coef from final iteration",chars=csn,font=font
for i=0,1 do begin &$
	lin=string(format='("Fit X0:",f6.4," Y0:",f6.4," Z0=",f8.4," Rad:",f8.4," (meters)")',$
		fitiar[9,i].fitcoef) &$
    note,ln+1+i*scl,lin,chars=csn,font=font,col=colph[i+1] &$
endfor
;
xp=.04
ln=5
note,ln      ,'+ X0',chars=csn,font=font,xp=xp
note,ln+1*scl,'* Y0',chars=csn,font=font,xp=xp
; z0
ver,264,265
tit='Z0 fit coef vs fit iteration'
plot,x,fitiar[*,0].fitcoef[2],psym=sym1,$
	chars=cs,font=font,$
	xtit=xtit,ytit=ytit,$
	tit=tit
oplot,x,fitiar[*,1].fitcoef[2],psym=sym1,col=colph[2]
ln=12
xp=.3
note,ln,"Scanner is about .7 meters above dish",chars=csn,font=font,xp=xp
; radius
ver,265.05,265.2
tit='Radius coef vs fit iteration'
plot,x,fitiar[*,0].fitcoef[3],psym=sym1,$
	chars=cs,font=font,$
	xtit=xtit,ytit=ytit,$
	tit=tit
oplot,x,fitiar[*,1].fitcoef[3],psym=sym1,col=colph[2]
;
; try make a histogram of the fit residuals
;
; cmp radius dif for final iteration
;
Radf3=fitiar[9,1].fitcoef[3]
cen3=reform(fitiar[9,1].fitcoef[0:2])
;
Radf4=fitiar[9,0].fitcoef[3]
cen4=reform(fitiar[9,0].fitcoef[0:2])
rdiff4=reform(blmcmpradius(xyz[*,iikeepf4],center=cen4)) - radf4
rdiff3=reform(blmcmpradius(xyz[*,iikeepf3],center=cen3)) - radf3

maxf4=max(rdiff4,min=minf4)
maxf3=max(rdiff3,min=minf3)
print,minf4,maxf4
print,minf3,maxf3
;
;make a histogram of the errors
;
minv=min(rdiff4)
maxv=max(rdiff4)
print,minv,maxv
maxv=1 
bin=.001
h4=histogram(rdiff4,bin=bin,loc=xlocf4)
h3=histogram(rdiff3,bin=bin,loc=xlocf3)
!p.multi=[0,1,2]
cs=1.4
ytit='Counts'
xtit='radial error [cm]'
tit='Histogram of radial errors (measured - fit)'
ver,0,1.4e6
hor,-1.7,1.7
sym=-1
plot,xlocf4*100.,h4,psym=sym,$
	chars=cs,font=font,$
	xtit=xtit,ytit=ytit,$
	tit=tit
oplot,xlocf3*100.,h3,col=colph[2],psym=sym
ls=2
flag,0,linestyle=ls,col=colph[3]
ln=2.5
xp=.04
scl=.7
note,ln      ,'4 param fit',chars=csn,font=font,xp=xp
note,ln+1*scl,'3 param fit',chars=csn,font=font,xp=xp,col=colph[2]
if hard then hardcopy
x
ldcolph
stop
;-------------------------------------------------------
; the fits are throwing away parts of the dish that are bad..
; keeping only points within 3 sigma give
;              4parmfit      3parmfit
; 3sigma was .0057*3=1.7cm   .0072*3= 2.2 cm
; we should probably keep points within say 5cm..
; recompute points to exclude using the 2 fits
;
; start with all data 
rad3=blmcmpradius(xyz,center=cen3)
rad4=blmcmpradius(xyz,center=cen4)
rdif3=rad3 - radF3
rdif4=rad4 - radF4
;
; now limit to 5 cm
rlim=.05
ii3=where(abs(rdif3) lt rlim,cnt3)
ii4=where(abs(rdif4) lt rlim,cnt4)
print,"cnt3:",cnt3," cnt4:",cnt4
;
plot,xyz[0L,ii3],xyz[1L,ii3],psym=3,/iso
plot,xyz[0L,ii4],xyz[1L,ii4],psym=3,/iso
;
;----------------------------------------------------
; 
stop
; --------------------------------------
; testing.. probably need to edit the  pallete to put a separate color
; for the no data location.
; or in iimage, set the min value to a little less than the min (without -999)
; histimg4 . look for missing value (0)
print,min(img4),max(img4)

bin=.0005
h=histogram(img4,bin=bin,loc=hloc,min=-rlim,max=rlim)
plot,hloc,h
; generate xy grid  so we can get locations for no data in inage
nx=n_elements(xout)
ny=n_elements(yout)
xx=fltarr(nx,nx)
yy=fltarr(nx,nx)
for i=0,nx-1 do xx[*,i]=xout
for i=0,ny-1 do yy[i,*]=yout
ii=where(img4 eq -999)
;
;  show the nodata on the image
;
plot,xx[ii],yy[ii],psym=3
; -----------------------------------------------
; histogram of points kept  vs radial xy distance
;
plothist:
rtp=blmxyztosph(xyz,/deg,/cwy)
;
bin=2
hr=histogram(reform(rtp[0,*]) ,binsize=bin,min=0,max=155,loc=hlocr)
hrf4=histogram(reform(rtp[0,iikeepf4]) ,binsize=bin,min=0,max=155,loc=hlocrf4)
hrf3=histogram(reform(rtp[0,iikeepf3]) ,binsize=bin,min=0,max=155,loc=hlocrf3)
;
if hard then pscol,'fitsphere_hist_rae.ps',/full
!p.multi=[0,1,3]
;-------------------------------------------------
; radius
cs=2.2
csn=1.5
sym=10
ylog=0
hor,0,150
ver,0,1.1e6
xtit='radial xy distance [meters]'
ytit='counts'
tit='scan19 histogram of points after fit exclusion  vs radial distance'
col3=2
col4=3

plot,hlocr,hr,psym=sym,ylog=ylog,$
	chars=cs,font=font,$
	xtit=xtit,ytit=ytit,$
	tit=tit
oplot,hlocrf3,hrf3,psym=sym,col=colph[col3]
oplot,hlocrf4,hrf4,psym=sym,col=colph[col4]
ln=4
xp=.5
scl=.7
note,ln,'Points before fit exclusion',chars=csn,font=font,col=colph[1],xp=xp
note,ln+1*scl,'Points after 3 param fit',chars=csn,font=font,col=colph[col3],xp=xp
note,ln+2*scl,'Points after 4 param fit',chars=csn,font=font,col=colph[col4],xp=xp
fl=[16.1,28.11]
ls=2
colfl8=4 
colfl8=10

flag,fl[0],col=colph[colfl8],linestyle=ls
flag,fl[1],col=colph[colfl8],linestyle=ls
xeps=.3
x1=fl[0]+xeps
y1=1.0e6
x2=fl[1]+xeps
y2=.90e6
lin=string(format='("radius hf 8 Mhz Dipoles: ",f6.3," m")',fl[0])
xyouts,x1,y1,lin,chars=csn,font=font,col=colph[colfl8]
lin=string(format='("radius hf 5 Mhz Dipoles: ",f6.3," m")',fl[1])
xyouts,x2,y2,lin,chars=csn,font=font,col=colph[colfl5]
;-------------------------------------------------
; hist vs az
;
minv=-90.
maxv=270.
bin=1
h=histogram(reform(rtp[1,*]) ,binsize=bin,min=minv,max=maxv,loc=hloc)
h4=histogram(reform(rtp[1,iikeepf4]) ,binsize=bin,min=minv,max=maxv,loc=hloc4)
h3=histogram(reform(rtp[1,iikeepf3]) ,binsize=bin,min=minv,max=maxv,loc=hloc3)
ii=where(hloc lt 0,cnt)
if cnt gt 0 then begin &$
	 hloc[ii]+=360. &$
	 hloc3[ii]+=360. &$
	 hloc4[ii]+=360. &$
endif
;
ii=sort(hloc)
tit='histogram of points after fit exclusion  vs scanner  azimuth'
xtit='az [deg scanner coord sys]'
hor,0,360
ver
plot,hloc[ii],h[ii],psym=sym,ylog=ylog,$
    chars=cs,font=font,$
    xtit=xtit,ytit=ytit,$
    tit=tit
oplot,hloc3[ii],h3[ii],psym=sym,col=colph[col3]
oplot,hloc4[ii],h4[ii],psym=sym,col=colph[col4]
;
ls=2
fl8=(360 - ([  -96.5136,30.509,147.001] - 90) - 1) mod 360.
flag,fl8,col=colph[colfl8],linestyle=ls
fl5=(360. -([  -146.982, -24.497, 96.995] -90) -1 ) mod 360.
flag,fl5,col=colph[colfl5],linestyle=ls
;
; flag large counts
;
fll=[168.38,348.42]-.5
flag,fll,linestyle=ls,col=colph[7]
;
xp=.04
ln=11.5
note,ln,'HF 8mhz dipoles',chars=csn,font=font,col=colph[colfl8],xp=xp
note,ln+1*scl,'HF 5mhz dipoles',chars=csn,font=font,col=colph[colfl5],xp=xp
note,ln+2*scl,'excessive counts',chars=csn,font=font,col=colph[7],xp=xp
; put hf dipole facing north az 0deg
;  
;-------------------------------------------------
; hist vs el
;
minv=-20.
maxv=18.
bin=.4
h=histogram(reform(rtp[2,*]) ,binsize=bin,min=minv,max=maxv,loc=hloc)
h4=histogram(reform(rtp[2,iikeepf4]) ,binsize=bin,min=minv,max=maxv,loc=hloc4)
h3=histogram(reform(rtp[2,iikeepf3]) ,binsize=bin,min=minv,max=maxv,loc=hloc3)
ii=sort(hloc)
xtit='el [deg scanner coord sys]'
tit='histogram of points after fit exclusion  vs scanner elevation'
hor,-15,18
ver,0,max(h)*1.05
plot,hloc[ii],h[ii],psym=sym,ylog=ylog,$
    chars=cs,font=font,$
    xtit=xtit,ytit=ytit,$
    tit=tit
oplot,hloc3[ii],h3[ii],psym=sym,col=colph[col3]
oplot,hloc4[ii],h4[ii],psym=sym,col=colph[col4]
if hard then hardcopy
x
ldcolph
stop
; ----------------------------------------------------------------------------------
;
; look at points in 3 param fit that were excluded from4 param fit
; used to create : radialerr_pnts.png
; import png:/share/megs/phil/x101/p50/200311/radialerr_pnts.png
;
n=n_elements(xyz[0,*])
ok=intarr(n)
; 3 in both
; 2 only in f4 4param
; 1 only in f3 3 param
ok[iikeepf4]=2
ok[iikeepf3]=ok[iikeepf3] + 1
;
ii=where(ok eq 2)
sym=3
!p.multi=[0,2,1] 
cs=1.4
csn=1.5
font=1
xtit='X [m] NorthWest (scannner coord)'
ytit='Y [m] NorthEast (scannner coord)'
tit='SouthEast'
hor,-150,150
ver,-150,150
plot,xyz[0,*],xyz[1,*],psym=sym,/iso,$
	chars=cs,/nodata,$
	xtit=xtit,ytit=ytit,$
	tit=tit
ln=.5
xp=.6
xpinc1=.12
xpinc2=.15
scl=.8
note,ln,"orig pnts",chars=csn,xp=xp
note,ln,'3paramFit pnts',chars=csn,xp=xp+xpinc1,col=colph[2]
note,ln,'4paramFit pnts',chars=csn,xp=xp+2*xpinc2,col=colph[3]
oplot,xyz[0,*],xyz[1,*],psym=sym
oplot,xyz[0,iikeepf3],xyz[1,iikeepf3],psym=sym,col=colph[2]
oplot,xyz[0,iikeepf4],xyz[1,iikeepf4],psym=sym,col=colph[3]
plot,xyz[0,ii],xyz[1,ii],psym=sym,/iso,$
	/nodata,$
	chars=cs,$
	xtit=xtit,ytit=ytit,$
	tit=tit
xp=.6
note,ln,'3paramFit pnts excluded by 4param fit',chars=csn,xp=xp,col=colph[2]
oplot,xyz[0,ii],xyz[1,ii],psym=sym,col=colph[2]  
stop
; -----------------------------------------------
; screen shots of 
; look at the x,y coverage after 4 param fit exclusion.

; -----------------------------------------------
; screen shot 1
; creates  xyusage_f4.png
; now using radialerr_pnts to show both
if 0 then begin
!p.multi=0
font=0
cs=1.4
hor,-150,150
ver,-150,150
xtit='x [meters]'
ytit='y [meters]'
tit='200311 4parm fit to sphere. resulting x,y pnts after fit clipping'
chth=2.
plot,xyz[0,iikeepf4],xyz[1,iikeepf4],psym=3,/iso,$
	chars=cs,charth=chth,$
	xtit=xtit,ytit=ytit,$
	tit=tit
endif
;plot,xyz[0,iikeepf3],xyz[1,iikeepf3],psym=3
; -----------------------------------------------
; screen shot of raw and clipped data
;
xlim=[-150.,150]
ylim=[-150.,150]
zlim=[-.7,47.]
xyrad=150.
xs=880
ys=829
window,0,xsize=xs,ysize=ys
xx=blmscale(hdr,d)
;
; display x,y to get xyraw_clip.png
; display y,z to get xzraw_clip.png
a=blmselectpnts(xx,iijunk,xlim=xlim,ylim=ylim,zlim=zlim,xyrad=xyrad)
stop
;-------------------
; compare scanner height above the dish
;
print,cen3- cen4

end