; 
; try to compute the gain loss
; do this with all errors < 5cm
; and restricted to just errors < 2cm (get rid of td errors)
;
; use the gridded img with 5cm error max
; generated in doit.
;
; xout[2000].,yout[2000]
; find points withing 150 meter radius
;
; --> 2nd try .. exclude all points > 2cm.. to look at after td fix
; added phase histograms
; --> run twice . 
;     once wit doclip=0 , then doclip=1
;
; 	goto,plotit2
   	goto,plotit
	doclip=1
	clipH=1.2
	rmsRefCm=.22	; use 2.2 mm for error in reference 
	cliprmsVal=.22 ; replace with rms 2.2 mm
	clipToVal=0. ;  replace with 1cm
	radDf=500.
	verb=1
	nodataVal=min(img)
	ftom=.3048d
; use 225 meter apperture
	radft=225/2./ftom
	ngrid=2000
	xx=fltarr(ngrid,ngrid)
	yy=fltarr(ngrid,ngrid)
	for i=0,ngrid-1 do begin &$
 		xx[*,i]=xout &$
		yy[i,*]=yout &$
	endfor
	if doclip then begin
		imgclip=img
		ii=where(imgclip gt clipH,nclipped)
 		imgclip[ii]=randomN(0,nclipped)*rmsRefCm + clipToVal
	endif
;
;  step thru y and x computing the gain change for a set of wavelengths
;
	ystart= -400. 
	xstart= -400.
	step=200.
	nsteps=abs(ystart)*2/step +1
	lambdaAr=[21.,12.6,6.,3.] ; all in cm 
	collambdaAr=[1,2,3,4]
	nlambda=n_elements(lambdaAr)
	ntot=nsteps*nsteps*nlambda
;
	icur=0l
	for iy=0,nsteps - 1 do begin
		y0=ystart + iy*step
 		for ix=0,nsteps -1 do begin
			x0=xstart + ix*step
;           max sure x0 is on the dish
			xedge=sqrt(radDf^2 - y0^2)
			if (x0 lt 0) and (x0 lt (-xedge)) then x0=-xedge
			if (x0 gt 0) and (x0 gt (xedge)) then x0=xedge
			if doclip then begin
			istat=gainlosscmp(imgclip,xx,yy,x0,y0,lambdaAr,gI,radF=radFt,$
					    rmsRefCm=rmsRefCm,nodataVal=nodataVal)
			endif else begin
			istat=gainlosscmp(img,xx,yy,x0,y0,lambdaAr,gI,radF=radFt,$
					    rmsRefCm=rmsRefCm,nodataVal=nodataVal)
			endelse
			if icur eq 0 then gIAr=replicate(gI[0],ntot)
			i0=icur
			i1=i0 + nlambda -1
			giar[i0:i1]=gI
			if verb then begin
				for i=0,nlambda -1 do begin
					j=i0 + i
				    ln=string(format='("x0,y0:",f6.1,1x,f6.1," lambda:",f4.1," npnts:",i7," gainLoss:",f5.3)',$
                        x0,y0,giar[j].lambdaCm,giar[j].npnts,giar[j].gainloss)

					print,ln
		        endfor
			endif
			icur+=nlambda
		endfor
	endfor
	if doclip then begin
	giarC=giar ; clipped
	endif else begin
	giarN=giar ; normal
	endelse
	stop
;
; make clipped img
;
; force image to go -5 to 5 cm
;
asave=max(imgclip,ind)
imgclip[ind]=5.
;.run ../displayimg
imgclip[ind]=asave
;
a=rms(imgclip)
print,a
; -1.3802026       2.1375222
;
;   nclipped/ntot=186158l/(2000d*2000) = 4.6%
;
; now plot out the results
;
plotit:
	 if hard then pscol,'gainloss.ps',/full
	for ic=0,1 do begin
		doclip=ic eq 1
		hor,-500,500
	xtit='X ft [west,east]'
	ytit='gain loss [db]'
	tit0=(doclip)?'11mar20 gain loss from dish errors(<1.2cm) vs x,y  and lambda.':$
	              '11mar20 gain loss from dish errors(<5cm) vs x,y  and lambda.'
	giar=(doclip)?giarC:giarN
	ver,-30,1
	font=1
	cs=2.2 
	csn=1.5
	!p.multi=[0,1,5]
	sym=-2
	for iy=nsteps-1 ,0,-1 do begin &$
		y0=ystart + iy*step &$
		ii=where(gIar.y0f eq y0,cnt) &$
		gtmp=gIar[ii] &$
		tit=tit0 + string(format='(" Y= ",i4," ft")',long(y0))  &$
		plot,[0,1],[0,1],/nodata,$
			chars=cs,font=font,$
			xtit=xtit,ytit=ytit,$
			title=tit  &$
		for il=0,nlambda-1 do begin &$
			l=lambdaar[il] &$
			icol=collambdaar[il] &$
			ii=where(gtmp.lambdaCm eq l) &$
			y=alog10(gtmp[ii].gainloss)*10. &$
			oplot,gtmp[ii].x0F,y,col=colph[icol],psym=sym &$
		endfor &$
		tit0=''
	endfor
	xp=.04
	csn=1.5
	stp=.8
	ln=1.5
	for i=0,nlambda-1 do begin &$
		l=lambdaAr[i] &$
	    note,ln + i*stp,string(format='(f4.1," cm")',l),chars=csn,font=font,col=colph[i+1],xp=xp &$
	endfor
	ln=2.3
	xp=.3
	scl=.7
	note,ln,"refRms: .22cm",chars=csn,font=font,xp=xp
	if ic eq 1 then begin
		lin=string(format='("errs>",f3.1,"cmm set to ",f4.2,"cm rms noise")',$
				clipH,clipRmsVal)
	    note,ln+1*scl,lin,chars=csn,font=font,xp=xp
	endif
	xp=.55
	note,ln,'225m beam',chars=csn,font=font,xp=xp
		
	ln=1.5
	xp=.85
	note,ln,"North",chars=csn,font=font,xp=xp
	ln=26  
	note,ln,"South",chars=csn,font=font,xp=xp
	endfor
;
; now plot the ratio of gains
;
plotit2:
	ver,0,5. 
	hor,-500,500
    font=1
    cs=2.2
    csn=1.5
    !p.multi=[0,1,5]
    sym=-2
	tit0='(Gain<1.2cmErr)/(Gain<5.cmErr).  '
	ytit='GainIncrease [db]'
    for iy=nsteps-1 ,0,-1 do begin &$
        y0=ystart + iy*step &$
        ii=where(gIarN.y0f eq y0,cnt) &$
        gtmpN=gIarN[ii] &$
        gtmpC=gIarC[ii] &$
        tit=tit0 + string(format='(" Y= ",i4," ft")',long(y0))  &$
        plot,[0,1],[0,1],/nodata,$
            chars=cs,font=font,$
            xtit=xtit,ytit=ytit,$
            title=tit  &$
        for il=0,nlambda-1 do begin &$
            l=lambdaar[il] &$
            icol=collambdaar[il] &$
            ii=where(gtmpN.lambdaCm eq l) &$
            y=alog10(gtmpC[ii].gainloss/gtmpN[ii].gainloss)*10. &$
            oplot,gtmpN[ii].x0F,y,col=colph[icol],psym=sym &$
        endfor &$
        tit0='' &$
    endfor
    xp=.04
    csn=1.5
    stp=.8
    ln=1.5
    for i=0,nlambda-1 do begin &$
        l=lambdaAr[i] &$
        note,ln + i*stp,string(format='(f4.1," cm")',l),chars=csn,font=font,col=colph[i+1],xp=xp &$
    endfor
	ln=1.5
	lin=string(format='("Errs >",f3.1,"cm set to ",f4.2,"cm rms noise")',$
  		clipH,cliprmsVal)
	note,ln,lin,chars=csn,font=font,xp=.2
	if hard then hardcopy
	x
	ldcolph
;
end