;
; compute the gain loss over an area of the dish using the measurement errors
;SYNTAX:istat=gainlosscmp(imgCm,xxF,yyF,x0F,y0F,lambdaArCm,gIAr,$
;                         radf=radf,rmsRefCM=rmsRefcm,nodataVal=nodataVal)
;ARGS:
;imgCm[n,n] : radial errors from griddata . units cm.
;xxf[n,n]  : fltarr . x position (feet) for each point of img
;yyf[n,n]  : fltarr . y position (feet) for each point of img
;x0F         : flt    . x center [ft] (in x,y plane) for the computation
;y0F         : flt    . y center [ft] (in x,y plane) for the computation
;lambdaArcm[m]:fltarr . wavelengths to compute the gain loss. (in cm).
;KEYWORDS:
; radF:    float       radius (feet) to use for beam. default is 225m/2/ftom=370 feet.
; rmsRefcm: float      rms phase error to use for the reference. def: .22cm
; nodataVal: double    value for nodata in img. this is the min value in the array.
;                      if you pass it in, the program will not bother to search for it.
;RETURNS:
;istat:    :  0 ok, -1 error
;gIar[m]   :  {}    struct contianing loss info
;DESCRIPTION:
; 	compute the gain loss for a set of wavelengths give the measured radial errors on the dish
;The user specifies the x,y center (in feet) where the value will be computed.
;The steps are:
; 1. Find the minimum value in the imgCm array. This will be the no data value (caused by 
;    things like the shadows from the hf dipoles. These points will not be included in the
;    reference or the computed loss.
; 2. find all points within radF of the x0,y0 location that have value measurements:
;    - when there is spillover, we are currently using a smaller beam area (that lies on the dish)
;      the reference is computed over the smaller area.. so the losses are measuring good dish vs bad
;      dish and not accounting for the decrease in gain caused by spillover.
; 3. for each wavelength
;    -compute the phase error for the reference (using the rmsDefmm)
;    -compute the phase error using the measured errors.
;    - integrate the E field for the ref and measurement
;    - compute the intensity and take the ratio.
;
; The gIAr[m] struct contains:
;** Structure <16bab458>, 7 tags, length=56, data length=52, refs=1:
;   LAMBDAcm  DOUBLE  0.0000000 ;  lambda, cm
;   RMSREFcm  DOUBLE  0.0000000 ; rms cm used for reference
;   RMSMEAScm DOUBLE  0.0000000 ; rms measured over the area.
;   X0F       DOUBLE  0.0000000 ; x center ft
;   Y0F       DOUBLE  0.0000000 ; y center ft
;   GAINLOSS  DOUBLE  0.0000000 ; linear
;   NPNTS     LONG    0         ; in the area
;
; 11may20.. changed phase errro err/lambda*2*pi to 2*err/lambda *2*pi
;-
function gainlosscmp,imgCm,xxF,yyF,x0F,y0F,lambdaArCm,gIAr,radf=radf,$
				       rmsRefCM=rmsRefcm,nodataVal=nodataVal
;
	if n_elements(nodataVal) eq 0 then nodataVal=min(imgCm)
	a=size(imgcm)
	gridLen=a[1]
	ftom=.3048d
	nlambda=n_elements(lambdaArcm)
	rmsRefCm=(n_elements(rmsRefcm) eq 0)?.22:rmsRefCm*1d
	radF    =(n_elements(radF) eq 0)? 225d/2d/ftom: radF*1d
	a={lambdaCm: 0D,$   ; lambda, cm
       rmsRefCm:0d,$  ; rms cm used for reference
       rmsMeasCm:0d,$ ; over the area [cm]
        x0F: 0D,$   ; x center ft
        y0F: 0D,$   ; y center ft
		hphM: lonarr(63),$; histogram of phase
		hlocM: dblarr(63),$; x of histogram of phase
    gainLoss:0d,$  ; linear
      npnts :0l}   ; included in area
;
	gIar=replicate(a,nlambda)
	giAr.lambdaCm=lambdaArcm
;   compute the radius about the center
    r= sqrt((xxF  - x0F)^2  + (yyF - y0F)^2)
;   select the points to use
	ii=where((imgCm ne noDataVal) and (r lt radF),npnts)
	if npnts eq 0 then begin
		print,"no points in img centered at x0,y0"
		return,-1
	endif
;    use this for the reference integral
	errCmRef=randomN(0.,gridLen,gridLen,/double)*rmsRefcm 
	a=rms(imgCm[ii],/quiet)
	rmsMeas=a[1]
	for  i=0,nlambda-1 do begin
		l=double(lambdaArcm[i])
;		generate reference phase using above noise
		phR=2*(errCmRef[ii]/l)*!dpi*2d
		Eref=dcomplex(cos(phR),sin(phR))
	    ErefSum=total(Eref)
;       compute E with errors
		phM=2d*(imgcm[ii]/l)*2d*!dpi
		EMeas=dcomplex(cos(phM),sin(phM))
		EMeasSum=total(Emeas)
		intenRef=(abs(ErefSum))^2
		intenMeas=(abs(EMeasSum))^2
;		print,"-->",ErefSum,EmeasSum
		y=phM mod (2d*!dpi)
		jj=where(y lt 0,cnt)
	    if cnt gt 0 then y[jj]+= 2d*!dpi
		giar[i].hphM=histogram(y,bin=.1d,min=0d,max=2d*!dpi,loc=hloc)
		giar[i].hlocM=hloc
		gIAr[i].rmsRefCm=rmsRefCm
		gIAr[i].rmsMeasCm=rmsMeas
		gIAr[i].x0F=x0F
		gIAr[i].y0F=y0F
		gIAr[i].gainLoss=(intenMeas/intenRef)
		gIar[i].npnts=npnts
	endfor 
	return,0
end