; =========================================================================================
  pro retcsigma,iter,n_trials,doy,orbit,high_sza,high_sza_sigma,sza_l1b,sva_l1b,ssa_l1b,    $
                alb_l1b,bin_sza,n_sza,dsza,cden_clim,hratio_clim,scale_cden,       $
                scale_hratio,cden_fg,hratio_fg,nadir_fg,sva60_fg,                  $
                n_all,cden_all,hratio_all,n_back,cden_back,hratio_back,            $
                min_back=min_back,verbose=verbose,nofit=nofit,backonly=backonly
; =========================================================================================
; This procedure obtains the first guess at C and sigma from an orbits worth 
; of observations. This is done by looking at the backward scattered
; observations and doing an XY analysis on them.
;
; SMB 8/2009
; SMB 9/10/2009 - converted to use backward and forward observations
; if you get a different answer there, you know there are clouds and
; should get C & sigma through interpolation or other means
; Kim Nielson 11/2009 brought in climatology for firstguess
; SMB 12/21/2009 removed fit based on keyword parameter
; ------------------------------------------------------------------------------------------

ind_nohi=where(bin_sza lt high_sza)
dum=max(ind_nohi,index_maxnohi)

if not keyword_set(min_back) then min_back=110.

cden_fg=fltarr(n_sza)   ; first guess at C
hratio_fg=fltarr(n_sza) ; first guess at sigma
cden_back=fltarr(n_sza)   ; first guess at C, back scatter only
hratio_back=fltarr(n_sza) ; first guess at sigma, back scatter only
cden_all=fltarr(n_sza)   ; first guess at C, all obs.
hratio_all=fltarr(n_sza) ; first guess at sigma, all obs.
n_back=fltarr(n_sza) ; number of back scattered observations
n_all=fltarr(n_sza) ; number of observations, all points

del_sza_norm=dsza*1.0
del_sza_highsza=dsza*1.0

for i=0,n_sza-1 do begin

; get all the backscattered and forward scattered observations 

  del_sza=del_sza_norm
  if bin_sza(i) gt high_sza then begin
     del_sza=del_sza_highsza
  endif

  crit_all = abs(sza_l1b - bin_sza(i)) le (del_sza/2.0) and $
             finite(sva_l1b) and $
             finite(ssa_l1b) and $
             finite(sza_l1b) and $
             finite(alb_l1b) and $
             alb_l1b gt 0.

  crit_backward = ssa_l1b ge min_back
  crit_back = crit_all and crit_backward

  ind_back=where(crit_back,n_ind_back)
  ind_all =where(crit_all, n_ind_all)

  n_back(i)=n_ind_back
  n_all(i)=n_ind_all

  if n_ind_back gt 1 then begin

     t_sva=sva_l1b(ind_back)
     t_ssa=ssa_l1b(ind_back)
     t_sza=sza_l1b(ind_back)
     t_alb=alb_l1b(ind_back)

     if bin_sza(i) le high_sza then begin
        calc_xy,t_sva,t_ssa,t_sza,t_alb,x,y,hratio_l1b,cden_l1b,flag, $
                    slope=xym,intercept=xyb 
     endif else begin
        arr_hisza_hratio=hratio_back(index_maxnohi-20:index_maxnohi)
        ind_arr_hisza_hratio=where(arr_hisza_hratio gt 0.,n_ind_arr_hisza_hratio)
        if n_ind_arr_hisza_hratio gt 2 then mn_hisza_hratio=median(arr_hisza_hratio[ind_arr_hisza_hratio]) $
           else mn_hisza_hratio=0.7
        calc_xy,t_sva,t_ssa,t_sza,t_alb,x,y,hratio_l1b,cden_l1b,flag, $
                    slope=xym,intercept=xy,fix_sigma=mn_hisza_hratio
     endelse
  endif else begin ; if i=0 (1st point) use climatology, otherwise make continuous
     if (i eq 0) then cden_l1b=cden_clim[i] else cden_l1b=cden_back[i-1]
     if (i eq 0) then hratio_l1b=hratio_clim[i] else hratio_l1b=hratio_back[i-1]
  endelse
  cden_back(i)=cden_l1b
  hratio_back(i)=hratio_l1b

  if n_ind_all gt 1 then begin

     t_sva=sva_l1b(ind_all)
     t_ssa=ssa_l1b(ind_all)
     t_sza=sza_l1b(ind_all)
     t_alb=alb_l1b(ind_all)

     if bin_sza(i) le high_sza then begin
        calc_xy,t_sva,t_ssa,t_sza,t_alb,x,y,hratio_l1b,cden_l1b,flag, $
                slope=xym,intercept=xyb 
     endif else begin
        if total(size(mn_hisza_hratio)) eq 0 then mn_hisza_hratio=0.7
        calc_xy,t_sva,t_ssa,t_sza,t_alb,x,y,hratio_l1b,cden_l1b,flag, $
                slope=xym,intercept=xy,fix_sigma=mn_hisza_hratio  ; use same values for all and back points
     endelse
  endif else begin ; if i=0 (1st point) use climatology, otherwise make continuous
     if (i eq 0) then cden_l1b=cden_clim[i] else cden_l1b=cden_all[i-1]
     if (i eq 0) then hratio_l1b=hratio_clim[i] else hratio_l1b=hratio_all[i-1]
  endelse
  cden_all(i)=cden_l1b
  hratio_all(i)=hratio_l1b
 
endfor   ; Loop over SZA

; Use back/all points to screen for clouds. 
; If CDEN from two methods agrees to within 10% - take result from BACK points fit.
; If there are no backscattered observations - take result from ALL points fit.
; if CDEN are different by more than 10% - use climatology
; if SZA greater than HIGHSZA, fix sigma/hratio to a specific value
; after this is all puttogether, fit the various points to a 
; fourth order polynomial to make sure it's smooth
; fit is done to SZA=HIGHSZA, then the value of sigma past that
; is fixed to the fit value at SZA=HIGHSZA

diff=abs((cden_back-cden_all)/cden_back)

; if the BACK and ALL calculations agree, take the BACK values

; testing.... check to see how things change if using 10% for 1st iteration only (not first 2)
;IF iter LE 1 THEN ind=where(diff le 0.1 and n_back gt 0,n_ind) ELSE $
;	ind=where(diff le 0.2 and n_back gt 0,n_ind)
crit0 = diff le 0.1 and n_back gt 0 AND bin_sza LT 92.
crit1 = diff le 0.2 AND n_back GT 0 AND bin_sza GE 92.
crit2 = diff le 0.2 and n_back gt 0
IF iter LE 1 THEN ind=where(crit0 OR crit1,n_ind) ELSE ind=where(crit2,n_ind)

if n_ind gt 0 then begin
   cden_fg[ind]=cden_back[ind]
   cden_back_scale = cden_back[ind]   ;we use this to scale clim...
   ind_scale = ind                   ;keep track of which indices we are using to calculate scaling factor.
   hratio_fg[ind]=hratio_back[ind]
   hratio_back_scale = hratio_back[ind]
endif

if keyword_set(backonly) then begin
   cden_fg=cden_back
   hratio_fg=hratio_back
endif

ind=where(n_back eq 0,n_ind)
if n_ind gt 0 then begin
   cden_fg(ind)=cden_all(ind)
   hratio_fg(ind)=hratio_all(ind)
endif

;+++++++++++++++++++++++++++++++++++++
;+++++++++++++++++++++++++++++++++++++
;K. Nielsen 11 Oct, 2010
;using interpolation in the lastr iteration is troublesome for many orbits
;with high roll as the number of bins (and pixels) we apply the interpolation
;is large (~>50 as compared to ~5 in NH2007). Instead, we should look at the 
;number of bins we are applying interpolated values to. If this number exceed
;30 (based on looking at 10 orbits in NH2007 and NH2009) we should use a
;polynomial fit instead of an interpolation
;++++++++++++++++++++++++++++++++++++++
;PRINT, iter+1, n_trials
IF iter EQ n_trials-1 THEN BEGIN   ;if last iteration, then interpolate (if n_ind0 LT 30)
	ind0=where(cden_fg eq 0,n_ind0)
	;PRINT, n_ind0
  	if (n_ind0 gt 0 AND n_ind0 LT 30.) then begin
    	indnot0=where(cden_fg gt 0.,n_indnot0)
    	if n_indnot0 ge 2 then begin
      		cden_fg(ind0)=interpol(cden_fg(indnot0),bin_sza(indnot0),bin_sza(ind0))
      		hratio_fg(ind0)=interpol(hratio_fg(indnot0),bin_sza(indnot0),bin_sza(ind0))
    	endif
    	;we interpolated, now skip polynomial fit
    	;PRINT, 'Applied interpolation...'
    	GOTO, skip_polyfit
  	endif
ENDIF

; ----------------------------------------------------------------------
; Now insert climatology where there is missing information and do a
; final polynomial fit to entire orbit to smooth out.
; NOTE: Not implemented if either keyword "nofit" or "backonly" is set
; ----------------------------------------------------------------------

;if (not keyword_set(nofit)) and (not keyword_set(backonly)) then begin

; -------------------------------------------------------------------
; Put the climatology anywhere that we don't yet have a value.
; The input climatology is scaled to the current C/sigma values 
; in the center of the orbit before being incorporated.
; -------------------------------------------------------------------

scale_hratio=0.
scale_cden=0.
ind_eq0=where(hratio_fg eq 0.,n_ind_eq0)
if n_ind_eq0 gt 0 then begin

;     Determine portion of orbit to scale to
  tmp_bin = bin_sza[ind_scale]
  ind_good = WHERE(tmp_bin GT 40 AND tmp_bin LT 70.,n_ind_good)

;     Scale and apply climatology for C
  cden_clim_tmp = cden_clim[ind_scale]   ;match the elements in cden_back
  ratio_clim_cden = cden_clim_tmp[ind_good] / cden_back_scale[ind_good]
  scale_cden = total(ratio_clim_cden)/float(n_ind_good)
  cden_clim = cden_clim / scale_cden
       
;     Scale and apply climatology for Sigma
  hratio_clim_tmp = hratio_clim[ind_scale]
  ratio_clim_hratio = hratio_clim_tmp[ind_good] / hratio_back_scale[ind_good]
  scale_hratio = total(ratio_clim_hratio)/float(n_ind_good)
  hratio_clim = hratio_clim / scale_hratio
       
;     Insert into missing bins of orbit
  cden_fg(ind_eq0)=cden_clim(ind_eq0)
  hratio_fg(ind_eq0)=hratio_clim(ind_eq0)
endif else ind_eq0=indgen(3) ; dummy values for plotting

; ------------------------------------------------------------------------
; Now do a polynomial fit to the entire orbit to make smooth & continuous
; ------------------------------------------------------------------------

crit_change=bin_sza gt 30. and bin_sza le high_sza

crit_fit=crit_change       and $
       hratio_fg gt 0.   and $
       finite(hratio_fg) and $
       cden_fg gt 0.     and $
       finite(cden_fg)   

ind_fit=where(crit_fit,n_ind_fit)
ind_change=where(crit_change)
nfit=4  ; order of polynomial fit
if n_ind_fit ge nfit then begin
	hratio_coef=poly_fit(bin_sza(ind_fit),hratio_fg(ind_fit),nfit,yfit=yfit)
  	cden_coef=poly_fit(bin_sza(ind_fit),cden_fg(ind_fit),nfit,yfit=yfit)

	hratio_fit=hratio_coef(nfit)
	cden_fit=cden_coef(nfit)
	for mm=nfit-1,0,-1 do begin
   		hratio_fit=hratio_fit*bin_sza+hratio_coef(mm)
   		cden_fit=cden_fit*bin_sza+cden_coef(mm)
	endfor
	hratio_fg(ind_change)=hratio_fit(ind_change)
	cden_fg(ind_change)=cden_fit(ind_change)
	;PRINT, 'Performed Polynomial fit...'
endif

skip_polyfit:

;above high_sza we want hratio_fg to be set to hratio_back (=hratio_all)
;K.Nielsen 08/19/10
ind = WHERE(bin_sza GT high_sza)
hratio_fg[ind] = hratio_back[ind]

; ------------------------------------------------------------------
; Fix C/sigma at low sza (<25). (This is not currently implemented
; but may come into play later)
; ------------------------------------------------------------------

ind_lsza = WHERE(bin_sza LT 25.,count)
nn = MAX(ind_lsza)
IF count GT 1 THEN BEGIN
    cden_fg[ind_lsza] = cden_fg(nn+1)
    hratio_fg[ind_lsza] = hratio_fg(nn+1)
ENDIF    

; ------------------------------------------------------------------
; Finally, calculate the nadir albedo and 60 deg albedo to pass out
; ------------------------------------------------------------------

rp=rayleigh_phase(180.-bin_sza)
nadir_fg=calc_albedo(cden_fg,hratio_fg,rp,0.,bin_sza)
rp=rayleigh_phase(180.-bin_sza-60.)
sva60_fg = calc_albedo(cden_fg,hratio_fg,rp,60.,bin_sza)

return
end