;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;+
; NAME:
;   LAMBDA_RICHNESS
; PURPOSE:
;   Optimized red-sequence redshift estimator for galaxy clusters
; EXPLANATION:
;   Calculate lambda richness from Rykoff et al. (2011).  Lambda is an
;   optimized richness estimator, designed to minimize the scatter in
;   LX at fixed richness with extensive tests on the maxBCG cluster
;   catalog (Koester et al. 2007).  The probabilistic methodology is
;   very robust to various perturbations, including changes in the
;   color used to calculate the red sequence, background shifts, and
;   zero-point uncertainties.
;   
;   The specific implementation here is appropriate for SDSS data and
;   galaxy clusters at 0.1<z<0.35.  If modifying, please be sure to
;   update the red sequence model, estimation of m*, and background structure.
;
;   The aperture used to measure each cluster is optimized based on the
;   richness of the cluster, such that the cutoff radius is given by:
;     rc = r0*(lambda/100)^beta
;
; CALLING SEQUENCE:
;   lambda = lambda_richness(z,gmr,gmr_err,imag,r,r0=,beta=,lstarcut=,h0=,bkg=,
;                            rlambda=,lambda_err=,wtvals=,pvals=)
;
; INPUTS:
;   z - the redshift of the cluster
;   gmr - the g-r color (dereddened MODEL_MAG) for each galaxy within the
;          region of interest (at least 1.5 h^-1 Mpc from the center)
;   gmr_err - the g-r color error for each galaxy
;   imag - the i-band magnitude (dereddended CMODEL_MAG) for each galaxy
;   r - the radius from the center (h^-1 Mpc) for each galaxy
;
; OPTIONAL INPUTS:
;   r0 - the constant in the radial scaling (h^-1 Mpc) [default 1.0]
;   beta - the exponent in the radial scaling [default 0.2]
;   lstarcut - the fraction of L* to use as a cut [default 0.2]
;   h0 - value for H0 used to define R0, r [default 100.0]
;   bkg - the background structure, consisting of the following fields:
;         .gmr - the gmr color bins [findgen(40)*gmrbinsize-1.0]
;         .imag - the imag magnitude bins [findgen(100)*imagbinsize+12.0]
;         .gmrbinsize - the size of the gmr bins [0.1 mag]
;         .imagbinsize - the size of the imag bins [0.1 mag]
;         .sigma_g - background density in N/deg.sq/mag/mag
;
; RETURN VALUE:
;   lambda - the lambda richness for the cluster
;
; OPTIONAL OUTPUTS:
;   rlambda - the cutoff radius used: rlambda = r0*(lambda/100)^beta
;   lambda_err - error in lambda estimated from Eqn. 3 in Rykoff et al. 2011
;   wtvals - weights for individual galaxies brighter than the luminosity
;             cutoff and within rlambda.  lambda = total(wtvals)
;             identically 0 for all galaxies outside the luminosity/radial
;             range
;   pvals - probabilities for all galaxies input.  For galaxies brighter than
;             the luminosity cutoff and within rlambda wtvals=pvals.
;             Outside the range, we can still calculate a probability for the
;             galaxy, with the caveat that the sum of probabilities will not
;             be properly normalized.
;
; CALLED ROUTINES:
;   LAMBDA_NFW_WEIGHT, LAMBDA_LUM_FUNC, LAMBDA_ANGDIST, LAMBDA_BCOUNTS,
;    LAMBDA_WEIGHTS, LAMBDA_BISECTION_ZERO, LAMBDA_SDSS_BACKGROUND
;
; REFERENCES:
;   Rykoff, E. S. et al 2011 (submitted)
; 
; REVISION HISTORY
;   Author: Eli S. Rykoff, LBNL, March 2011
;-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;-
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; routine to compute the NFW filter weight
;   (see Eqns 6-8 in Rykoff et al. 2011)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


function lambda_nfw_weight,r,rscale=rscale

if n_params() eq 0 then begin
    print,'syntax - weights = lambda_nfw_weight(r,rscale=rscale)'
    return,-1
endif

if n_elements(rscale) eq 0 then rscale=0.15

x=r/rscale

nx=n_elements(x)
sigx=fltarr(nx)

;; core is 100 kpc
corer=0.1
corex=corer/rscale

low=where(x lt corex,nlow)
mid=where(x ge corex and x lt 1.0,nmid)
high=where(x gt 1.0 and x lt 10.0/rscale,nhigh)
other=where(x gt 0.999 and x lt 1.001,nother)

if (nlow gt 0) then begin
    arg=sqrt((1.-corex)/(1.+corex))
    pre=2./(sqrt(1.-corex^2.))
    front=1./(corex^2.-1)
    sigx[low]=front*(1.-pre*0.5*alog((1.+arg)/(1.-arg)))
endif

if (nmid gt 0) then begin
    arg=sqrt((1.-x[mid])/(1.+x[mid]))
    pre=2./(sqrt(1.-x[mid]^2.))
    front=1./(x[mid]^2.-1)
    sigx[mid]=front*(1.-pre*0.5*alog((1+arg)/(1.-arg)))
endif

if (nhigh gt 0) then begin
    arg=sqrt((x[high]-1.)/(x[high]+1.))
    pre=2./(sqrt(x[high]^2.-1.))
    front=1./(x[high]^2.-1.)
    sigx[high]=front*(1.-pre*atan(arg))
endif

if (nother gt 0) then begin
    xlo=0.999
    xhi=1.001
    arglo=sqrt((1.-xlo)/(1.+xlo))
    prelo=2./(sqrt(1.-xlo^2.))
    frontlo=1./(xlo^2.-1)
    testlo=frontlo*(1.-prelo*0.5*alog((1+arglo)/(1.-arglo)))
    arghi=sqrt((xhi-1.)/(xhi+1.))
    prehi=2./(sqrt(xhi^2.-1.))
    fronthi=1./(xhi^2.-1.)
    testhi=fronthi*(1.-prehi*atan(arghi))
    sigx[other]=(testlo+testhi)/2.0

endif

wt=2.*!pi*r*sigx

;; we need to make sure that things aren't blowing up at 0:
smallr = where(r lt 1e-6,nsmall)
if (nsmall gt 0) then wt[smallr] = 2.*!pi*(1e-6)*sigx[smallr]

return,wt
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; routine to compute luminosity filter weight
;   (see Eqn 10 in Rykoff et al. 2011)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


function lambda_lum_func,imag,mstar,limmag,alpha=alpha

if n_params() eq 0 then begin
    print,'syntax - weights = lambda_lum_func(imag,mstar,limmag,alpha=alpha)'
    return,-1
endif

if n_elements(alpha) eq 0 then alpha=-0.8

;; do the numerical integration
stepsize=0.01
range=[10,limmag]
nsteps=(range[1]-range[0])/stepsize + 1

;; calculate the binned values for the normalization, as well as the values for
;; the specific magnitudes
steps=[findgen(nsteps)*stepsize+range[0],imag]

;; Schechter function
f=10.^(0.4*(alpha+1.0)*(mstar-steps))*exp(-10.^(0.4*(mstar-steps)))

;; integrate
n=int_tabulated(steps[0:nsteps-1],f[0:nsteps-1])

;; and the end of the array are the values we want, after normalizing
wts=f[nsteps:n_elements(f)-1]/n

return,wts
end


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; routine to compute 1/E(z) for angular diameter distance calculation
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function lambda_oneovere,z,omega_l=omega_l

if n_elements(omega_k) eq 0 then omega_l = 0.7
omega_m = 1.0 - omega_l

return,1./sqrt(omega_m*(1+z)^3.+omega_l)
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; routine to compute angular diameter distance
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function lambda_angdist,z,h0=h0,omega_l=omega_l

if n_elements(h0) eq 0 then h0=100.0
if n_elements(omega_l) eq 0 then omega_l=0.7

c=2.99792458d5

dh=c/h0

dm=fltarr(n_elements(z))
for i=0l,n_elements(z)-1 do begin
    qsimp,'lambda_oneovere',0,z[i],dc,omega_l=omega_l,/double
    dm[i] = dc
endfor

if n_elements(dm) eq 1 then dm=dm[0]


return,dh*dm/(1+z)
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; routine to compute b(x) = 2*pi*r*Sigma_g(m,c)
;   (see S3.4 in Rykoff et al. 2011)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


function lambda_bcounts,bkg,z,r,gmr,imag,h0=h0

if n_params() eq 0 then begin
    print,'syntax- bcounts = lambda_bcounts(bkg,z,r,gmr,imag,h0=h0)'
    return,-1
endif

if n_elements(h0) eq 0 then h0=100.0

ngmr = n_elements(bkg.gmr)
gmrindex = fix((gmr-bkg.gmr[0])*ngmr/((bkg.gmr[ngmr-1]+bkg.gmrbinsize)-bkg.gmr[0]))
nimag = n_elements(bkg.imag)
imagindex = fix((imag-bkg.imag[0])*nimag/((bkg.imag[nimag-1]+bkg.imagbinsize)-bkg.imag[0]))

;; check for overruns
badgmr=where(gmrindex lt 0 or gmrindex ge ngmr,nbadgmr)
if (nbadgmr gt 0) then $
  gmrindex[badgmr] = 0
badimag=where(imagindex lt 0 or imagindex ge nimag,nbadimag)
if (nbadimag gt 0) then $
  imagindex[badimag] = 0

sigma_g=bkg.sigma_g[gmrindex,imagindex]

;; and set all bad elememts to infinity: these galaxies should not be given any
;; weight
if (nbadgmr gt 0) then sigma_g[badgmr]=1./0.0
if (nbadimag gt 0) then sigma_g[badimag]=1./0.0
badcombination = where(sigma_g eq 0.0,nbadcombination)
if (nbadcombination gt 0) then sigma_g[badcombination] = 1./0.0

;; get the conversion factor to go from deg. to Mpc
angdist = lambda_angdist(z,h0=h0)

c=(angdist*1.0)*(!dpi/180.)


;; calculate b_i
bcounts=2.*!pi*r*(sigma_g/c^2.)

return,bcounts
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; routine to calculate p_i = lambda*u/(lambda*u+b)
;   (see Eqn. 1 in Rykoff et al. 2011)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


function lambda_weights,inlambda,r0,beta,ucounts,bcounts,r,theta_i

if n_params() eq 0 then begin
    print,'syntax- pwt=lambda_weights(lambda,r0,beta,ucounts,bcounts,r,theta_i)'
    return,-1
endif

;; figure out the radial cut
rc=r0*(inlambda/100.)^beta

;; calculate the normalization for the NFW filter for this radial cut
;;  (see Eqns. 8-9 in Rykoff et al. 2011)
nfwnorm = exp(1.65169 -0.547850*alog(rc)+0.138202*alog(rc)^2.-0.0719021*alog(rc)^3.-0.0158241*alog(rc)^4-0.000854985*alog(rc)^5.)

;; we only want to weight those that are inside...
;; our weight function is theta
theta_r = fltarr(n_elements(r))
in=where(r le rc,nin)
if (nin gt 0) then theta_r[in] = 1.0 ;; 0.0 otherwise

;; store both "p" the probabilities for all galaxies in the region, and "wt"
;; the weights used to calculate lambda
pwt=fltarr(2,n_elements(ucounts))

pwt[0,*] = (inlambda*ucounts*nfwnorm)/(inlambda*ucounts*nfwnorm + bcounts)

;; fix any NaNs (if both ucounts and bcounts are 0.0)
bad=where(finite(pwt[0,*]) eq 0,nbad)
if (nbad gt 0) then pwt[0,bad]=0.0

pwt[1,*] = pwt[0,*] * theta_r * theta_i

return,pwt
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; routine to calculate the zero in the function lambda - Sum(p)=0
;   (see Eqn. 2 in Rykoff et al. 2011)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


function lambda_bisection_zero,r0,beta,ucounts,bcounts,r,theta_i,wtvals=wtvals,pvals=pvals,tol=tol

if n_params() eq 0 then begin
    print,'syntax- lambda_bisection_zero,r0,beta,ucounts,bcounts,r,theta_i,wtvals=wtvals,pvals=pvals,tol=tol'
    return,-1
endif

if n_elements(tol) eq 0 then tol=1e-3

;; simple bisector to find the zero of the function

lamhi=800.0
lamlo=0.5
outlo=-1
while (abs(lamhi-lamlo) gt 2*tol) do begin
    mid=(lamhi+lamlo)/2.
    if (outlo lt 0) then begin
        pwt=lambda_weights(lamlo,r0,beta,ucounts,bcounts,r,theta_i)
        outlo=total(pwt[1,*])
    endif
    pwt=lambda_weights(mid,r0,beta,ucounts,bcounts,r,theta_i)
    outmid = total(pwt[1,*])

    if (outlo lt 1.0) then outlo=0.9 ;; provides stability at low end
    if ((outlo-lamlo)*(outmid-mid) gt 0) then begin
        ;; we need the upper bracket
        lamlo=mid
        outlo=-1
    endif else lamhi=mid ;; we need the lower bracket
endwhile

lambda=(lamlo+lamhi)/2.0
pwt=lambda_weights(lambda,r0,beta,ucounts,bcounts,r,theta_i)
pvals=pwt[0,*]
wtvals=pwt[1,*]

if (lambda lt 1.0) then begin
    ;; failure to find any galaxies
    lambda = -1
endif

return,lambda
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; routine to generate the background structure, obtained from 
;  a counts-in-cells analysis of all SDSS DR7 galaxies
;  (see S3.4 in Rykoff et al. 2011)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


function lambda_sdss_background

;; load in the mean background structure

;; This has been generated from the full SDSS catalog, via a Counts-In-Cells
;; analysis.  The units of sigma_g are n/sq.deg./mag/mag

ngmrbins=40
gmrbinsize=0.1
gmrzero=-1.0
nimagbins=100
imagbinsize=0.1
imagzero=12.0

bkg=create_struct('gmr',findgen(ngmrbins)*gmrbinsize+gmrzero,$
                  'imag',findgen(nimagbins)*imagbinsize+imagzero,$
                  'gmrbinsize',gmrbinsize,$
                  'imagbinsize',imagbinsize,$
                  'sigma_g',fltarr(ngmrbins,nimagbins))

bkg.sigma_g[*,0]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.12,0.25,0.26,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,1]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.06,0.11,0.38,0.31,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,2]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.07,0.19,0.39,0.34,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,3]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.08,0.19,0.42,0.43,0.06,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,4]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.09,0.27,0.56,0.56,0.06,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,5]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.06,0.12,0.28,0.68,0.76,0.08,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,6]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.10,0.15,0.27,0.77,0.91,0.11,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,7]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.07,0.12,0.15,0.32,0.84,0.97,0.12,0.00,0.00,0.00,0.05,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,8]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.10,0.15,0.42,1.02,1.18,0.20,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,9]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.09,0.20,0.32,0.55,1.22,1.19,0.16,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,10]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.15,0.24,0.25,0.50,1.34,1.56,0.28,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,11]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.05,0.09,0.19,0.33,0.62,1.56,1.61,0.27,0.05,0.00,0.00,0.06,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,12]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.11,0.24,0.60,0.84,1.89,2.07,0.29,0.10,0.00,0.00,0.07,0.06,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,13]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.06,0.15,0.32,0.63,1.04,2.23,2.39,0.38,0.07,0.00,0.00,0.10,0.09,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,14]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.05,0.25,0.48,0.66,1.07,2.39,2.70,0.51,0.08,0.00,0.00,0.13,0.11,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,15]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.08,0.30,0.60,0.81,1.25,2.79,3.14,0.57,0.10,0.00,0.07,0.18,0.15,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,16]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.11,0.32,0.61,0.88,1.41,3.08,3.86,0.84,0.15,0.06,0.07,0.23,0.20,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,17]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.10,0.26,0.62,1.05,1.51,3.43,4.32,0.96,0.11,0.00,0.06,0.24,0.17,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,18]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.11,0.40,0.96,1.30,1.85,3.81,4.66,1.19,0.16,0.07,0.10,0.28,0.24,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,19]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.18,0.61,1.06,1.46,2.14,4.28,5.38,1.39,0.25,0.11,0.17,0.47,0.34,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,20]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.07,0.17,0.59,1.26,1.84,2.60,5.09,6.22,1.67,0.27,0.16,0.17,0.44,0.27,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,21]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.06,0.15,0.59,1.44,2.20,2.86,5.77,7.07,2.05,0.28,0.12,0.24,0.50,0.39,0.07,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,22]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.05,0.25,0.81,1.60,2.43,3.31,6.14,8.00,2.50,0.37,0.20,0.23,0.54,0.45,0.07,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,23]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.07,0.29,0.83,1.74,3.05,4.15,6.99,8.50,3.08,0.49,0.20,0.22,0.59,0.52,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,24]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.08,0.26,0.95,2.15,3.34,4.71,7.80,10.10,3.92,0.58,0.21,0.24,0.69,0.56,0.05,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,25]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.09,0.21,1.07,2.50,3.42,4.81,9.03,11.97,4.71,0.66,0.27,0.31,0.74,0.62,0.10,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,26]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.09,0.36,1.37,2.65,3.68,5.16,9.43,13.16,5.60,0.91,0.27,0.34,0.88,0.68,0.10,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,27]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.07,0.11,0.45,1.59,2.99,4.54,6.18,10.50,15.04,7.06,1.27,0.35,0.38,0.99,0.81,0.12,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,28]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.08,0.13,0.54,1.86,3.71,5.35,6.99,11.73,17.31,8.57,1.71,0.51,0.49,1.08,0.84,0.11,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,29]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.05,0.14,0.62,1.91,3.97,6.12,7.55,13.07,18.98,10.04,2.26,0.56,0.43,1.18,0.87,0.10,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,30]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.07,0.23,0.68,2.00,4.38,6.84,8.96,14.37,21.03,11.70,2.83,0.70,0.52,1.25,0.86,0.13,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,31]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.10,0.22,0.88,2.55,5.16,7.23,9.74,15.78,23.21,13.80,3.65,1.07,0.58,1.28,0.95,0.12,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,32]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.06,0.10,0.23,0.97,3.09,5.67,8.28,10.77,16.74,25.98,16.35,4.44,1.29,0.74,1.44,1.13,0.15,0.00,0.06,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,33]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.05,0.10,0.27,0.98,3.30,6.40,9.20,12.31,18.65,28.51,18.83,5.50,1.61,0.76,1.52,1.19,0.17,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,34]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.05,0.14,0.37,1.18,4.02,7.72,10.53,13.27,20.20,31.14,22.58,7.53,2.11,0.93,1.72,1.35,0.23,0.00,0.00,0.06,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,35]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.05,0.05,0.07,0.15,0.39,1.37,4.55,9.05,11.89,15.27,21.61,33.52,25.38,9.89,2.92,1.23,2.04,1.58,0.24,0.05,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,36]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.07,0.07,0.24,0.47,1.79,5.43,10.13,13.09,17.19,24.84,37.45,30.14,12.38,3.89,1.87,2.46,1.76,0.23,0.05,0.00,0.06,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,37]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.08,0.08,0.25,0.63,2.15,5.93,10.72,14.74,18.73,27.51,41.14,34.62,15.41,5.19,2.55,2.75,1.89,0.31,0.13,0.08,0.05,0.06,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,38]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.05,0.13,0.24,0.69,2.45,6.74,12.34,16.63,21.08,29.48,44.56,38.37,19.32,6.86,3.34,3.12,2.21,0.38,0.11,0.10,0.06,0.00,0.06,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,39]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.11,0.26,0.80,2.88,8.20,13.89,18.29,22.74,32.57,47.97,43.05,22.45,9.05,4.10,3.93,2.73,0.43,0.09,0.09,0.08,0.06,0.06,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,40]=[0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.08,0.16,0.38,0.90,3.30,9.17,15.44,20.94,25.98,35.43,53.25,49.14,26.47,11.47,5.22,4.45,2.96,0.48,0.16,0.09,0.07,0.09,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,41]=[0.00,0.00,0.00,0.00,0.00,0.00,0.07,0.10,0.10,0.26,0.42,1.00,3.80,10.47,17.72,22.96,28.75,39.74,57.60,54.30,31.68,14.67,6.61,5.47,3.49,0.76,0.31,0.15,0.08,0.00,0.00,0.05,0.00,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,42]=[0.00,0.00,0.05,0.00,0.00,0.00,0.09,0.13,0.16,0.24,0.42,1.27,5.05,12.47,19.62,25.34,32.04,42.71,62.12,60.84,38.84,18.89,9.15,6.63,4.30,1.00,0.35,0.20,0.10,0.09,0.07,0.06,0.05,0.00,0.00,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,43]=[0.00,0.00,0.06,0.00,0.05,0.07,0.09,0.09,0.13,0.30,0.49,1.43,5.47,14.55,22.71,29.04,35.80,46.72,66.51,68.28,45.07,23.11,11.73,8.53,5.19,1.21,0.32,0.24,0.12,0.11,0.06,0.00,0.06,0.06,0.06,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,44]=[0.00,0.06,0.06,0.05,0.00,0.08,0.13,0.16,0.21,0.36,0.59,1.58,5.67,15.92,24.99,31.79,39.47,51.45,72.90,76.09,53.88,28.94,15.00,10.42,6.24,1.58,0.51,0.29,0.17,0.13,0.06,0.00,0.06,0.11,0.07,0.05,0.05,0.00,0.00,0.00]
bkg.sigma_g[*,45]=[0.00,0.00,0.06,0.06,0.07,0.10,0.12,0.19,0.24,0.33,0.69,2.08,7.01,19.04,29.06,35.65,42.69,55.66,78.42,84.28,62.88,35.06,19.55,13.62,7.91,2.12,0.68,0.31,0.19,0.17,0.08,0.08,0.09,0.06,0.05,0.07,0.00,0.00,0.00,0.07]
bkg.sigma_g[*,46]=[0.05,0.00,0.00,0.06,0.12,0.14,0.14,0.18,0.25,0.38,0.80,2.41,8.40,21.04,31.98,39.78,48.31,61.67,84.58,92.22,71.86,42.35,24.57,16.89,10.36,2.95,1.00,0.44,0.27,0.17,0.12,0.10,0.07,0.08,0.08,0.00,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,47]=[0.05,0.12,0.12,0.11,0.10,0.12,0.16,0.19,0.27,0.48,1.02,2.59,9.48,22.93,35.29,44.10,54.32,67.55,92.68,101.16,81.60,50.80,31.31,21.72,13.28,4.25,1.31,0.67,0.34,0.20,0.14,0.14,0.09,0.12,0.09,0.06,0.00,0.00,0.00,0.00]
bkg.sigma_g[*,48]=[0.07,0.11,0.11,0.10,0.11,0.14,0.17,0.26,0.33,0.64,1.28,3.09,11.00,26.32,39.92,50.11,59.52,73.16,98.04,110.37,94.91,61.61,38.43,27.50,16.45,5.66,1.85,0.85,0.39,0.31,0.26,0.20,0.12,0.13,0.08,0.11,0.08,0.00,0.00,0.00]
bkg.sigma_g[*,49]=[0.06,0.10,0.13,0.14,0.20,0.19,0.18,0.27,0.43,0.71,1.24,3.71,12.79,30.16,44.85,55.71,64.82,79.88,107.12,120.43,105.51,72.68,47.52,34.24,20.74,7.55,2.76,1.48,0.60,0.33,0.29,0.22,0.14,0.16,0.13,0.08,0.11,0.07,0.11,0.06]
bkg.sigma_g[*,50]=[0.06,0.12,0.14,0.20,0.26,0.15,0.19,0.27,0.39,0.83,1.63,4.31,14.04,33.49,49.95,60.61,71.98,89.02,116.29,129.70,116.51,83.62,59.02,42.61,26.58,10.69,4.00,2.21,0.96,0.41,0.25,0.23,0.18,0.23,0.16,0.09,0.08,0.08,0.11,0.00]
bkg.sigma_g[*,51]=[0.08,0.14,0.16,0.24,0.22,0.22,0.20,0.30,0.53,0.99,1.76,4.72,16.03,37.31,54.99,66.45,79.30,95.91,123.35,142.27,128.83,96.92,69.84,52.13,33.98,13.70,5.79,3.27,1.51,0.52,0.35,0.30,0.18,0.17,0.13,0.08,0.09,0.10,0.08,0.08]
bkg.sigma_g[*,52]=[0.10,0.10,0.16,0.21,0.24,0.26,0.34,0.44,0.72,1.00,1.78,5.25,18.12,42.28,60.98,72.63,86.79,105.91,133.43,153.40,142.88,111.08,84.03,63.51,42.05,18.14,7.92,4.49,1.71,0.68,0.42,0.33,0.23,0.22,0.17,0.16,0.16,0.09,0.07,0.08]
bkg.sigma_g[*,53]=[0.15,0.15,0.15,0.22,0.31,0.35,0.49,0.57,0.77,1.15,2.21,6.33,20.33,47.35,67.92,81.30,98.10,117.58,144.12,164.81,159.37,128.69,98.52,78.00,52.12,24.45,11.28,6.48,2.18,0.61,0.49,0.36,0.26,0.26,0.17,0.20,0.20,0.10,0.12,0.09]
bkg.sigma_g[*,54]=[0.16,0.19,0.21,0.26,0.24,0.30,0.48,0.58,0.87,1.37,2.42,7.13,22.93,53.09,75.34,90.24,108.31,126.34,153.48,176.81,175.67,146.53,114.62,93.86,63.20,30.49,15.63,9.48,3.25,0.90,0.64,0.41,0.31,0.26,0.18,0.17,0.24,0.18,0.15,0.07]
bkg.sigma_g[*,55]=[0.18,0.23,0.28,0.28,0.27,0.37,0.50,0.66,0.84,1.38,2.88,7.98,25.88,58.77,84.62,99.80,117.62,139.61,166.38,189.73,189.47,163.36,134.87,111.84,77.56,39.19,21.13,13.31,4.77,1.18,0.73,0.41,0.37,0.31,0.25,0.20,0.21,0.21,0.18,0.13]
bkg.sigma_g[*,56]=[0.18,0.30,0.29,0.30,0.33,0.53,0.61,0.78,0.86,1.42,3.23,8.68,28.61,64.84,92.93,108.10,128.08,152.70,179.04,204.76,208.35,183.54,155.77,132.48,93.56,50.09,28.06,17.94,6.26,1.35,0.72,0.56,0.45,0.39,0.32,0.23,0.17,0.17,0.14,0.13]
bkg.sigma_g[*,57]=[0.20,0.25,0.32,0.41,0.41,0.50,0.68,0.87,1.12,1.66,3.68,10.15,32.42,72.41,103.51,121.30,143.12,166.50,197.10,219.73,225.32,203.86,177.40,153.92,111.41,61.55,37.03,25.06,8.24,1.91,0.87,0.67,0.53,0.49,0.46,0.29,0.25,0.19,0.16,0.17]
bkg.sigma_g[*,58]=[0.32,0.36,0.36,0.42,0.41,0.68,0.87,1.02,1.11,1.90,3.87,11.45,36.49,81.48,114.78,134.79,156.27,182.99,212.63,239.53,247.25,228.74,201.82,178.56,132.62,76.43,49.31,32.49,10.72,2.51,1.07,0.68,0.60,0.48,0.46,0.29,0.33,0.23,0.12,0.16]
bkg.sigma_g[*,59]=[0.29,0.36,0.45,0.46,0.55,0.76,0.80,0.95,1.45,2.26,4.35,12.84,41.06,91.57,128.63,149.44,171.91,200.24,230.47,256.78,264.86,248.44,225.70,203.20,154.96,94.07,64.94,43.97,15.08,3.42,1.31,0.86,0.69,0.46,0.42,0.38,0.47,0.30,0.20,0.20]
bkg.sigma_g[*,60]=[0.31,0.39,0.55,0.58,0.61,0.75,0.87,1.11,1.64,2.43,5.19,14.31,45.21,99.62,141.05,164.13,188.60,219.45,249.36,275.43,284.59,275.43,255.79,232.76,180.16,114.67,82.85,57.70,19.88,4.30,1.57,1.15,0.76,0.59,0.42,0.38,0.44,0.44,0.28,0.24]
bkg.sigma_g[*,61]=[0.28,0.48,0.55,0.54,0.77,1.01,1.15,1.30,1.86,2.84,5.68,15.75,49.36,110.44,157.61,182.24,209.32,239.14,268.59,296.57,306.04,301.85,286.80,265.43,207.22,138.09,106.79,74.01,25.60,6.04,2.14,1.32,0.84,0.72,0.59,0.51,0.45,0.42,0.37,0.32]
bkg.sigma_g[*,62]=[0.43,0.58,0.59,0.67,0.84,1.16,1.57,1.64,2.28,3.45,6.96,18.54,55.81,124.25,176.98,201.32,227.59,260.27,292.31,321.49,331.86,326.59,319.37,300.90,238.54,166.23,135.91,94.53,33.10,7.64,2.81,1.61,1.14,0.86,0.76,0.70,0.62,0.52,0.43,0.34]
bkg.sigma_g[*,63]=[0.52,0.73,0.88,0.88,1.02,1.46,1.60,2.01,2.58,4.22,7.94,20.68,61.98,136.03,195.98,223.68,251.80,285.89,317.28,345.13,359.10,357.31,354.36,336.67,274.71,199.86,166.56,116.43,42.23,9.99,3.51,1.81,1.61,1.22,0.95,0.72,0.59,0.52,0.43,0.45]
bkg.sigma_g[*,64]=[0.68,0.96,0.97,1.27,1.43,1.67,1.98,2.59,3.17,4.57,9.28,23.26,68.51,148.54,215.07,247.69,277.08,312.80,344.47,371.11,386.88,388.88,391.09,374.09,310.75,237.60,207.35,142.27,53.26,13.57,4.53,2.59,1.98,1.42,1.15,1.15,0.90,0.67,0.60,0.55]
bkg.sigma_g[*,65]=[0.71,0.94,1.05,1.42,1.67,1.93,2.49,3.01,3.71,5.49,10.57,26.70,75.97,162.65,237.90,274.26,307.24,341.12,373.22,406.09,420.21,425.97,427.69,412.69,351.81,279.93,250.14,172.91,66.74,17.37,6.25,3.32,2.44,1.90,1.50,1.34,1.11,1.05,0.84,0.67]
bkg.sigma_g[*,66]=[0.87,1.16,1.46,1.61,1.93,2.31,2.75,3.51,4.39,6.56,11.94,30.12,83.59,178.67,261.12,302.38,334.97,374.05,406.92,436.76,456.19,467.50,466.57,457.34,396.62,328.46,298.21,206.98,81.76,22.58,7.47,4.12,2.90,2.20,1.93,1.48,1.22,1.20,1.18,0.82]
bkg.sigma_g[*,67]=[1.28,1.50,1.69,1.94,2.45,2.88,3.33,3.96,5.27,7.58,13.13,33.40,91.62,195.77,287.91,335.64,371.41,409.09,444.13,475.50,495.38,505.76,512.84,506.50,444.79,383.54,350.47,242.60,99.17,29.19,9.98,5.25,3.58,2.78,2.36,1.88,1.35,1.32,1.39,1.14]
bkg.sigma_g[*,68]=[1.26,1.83,1.86,2.26,2.71,3.42,3.93,4.75,6.36,8.84,15.51,38.64,102.09,215.12,317.22,368.29,406.17,449.26,489.71,521.07,544.17,553.06,561.36,556.41,497.51,442.46,407.04,279.76,118.58,38.06,13.39,6.90,4.60,3.46,2.80,2.23,1.72,1.68,1.61,1.39]
bkg.sigma_g[*,69]=[1.35,2.10,2.12,2.50,3.01,3.70,4.36,5.40,6.96,9.99,18.17,41.56,108.63,232.91,345.59,404.92,447.91,489.48,534.34,571.20,592.47,604.67,613.67,606.10,553.70,505.21,468.89,319.86,141.28,47.51,17.10,8.62,5.47,3.96,3.11,2.59,2.05,1.95,1.81,1.41]
bkg.sigma_g[*,70]=[1.63,2.23,2.39,2.91,3.32,4.00,4.90,5.95,7.59,10.79,19.16,45.73,118.84,251.84,379.97,450.52,491.54,534.55,584.38,625.14,648.44,662.92,667.05,654.80,611.21,579.42,537.95,363.90,165.16,59.24,21.99,10.75,6.80,5.18,3.73,3.03,2.36,2.22,2.10,1.51]
bkg.sigma_g[*,71]=[1.95,2.48,2.73,3.41,3.68,4.65,5.37,6.64,8.53,12.19,21.72,49.49,128.20,271.43,412.40,490.16,538.75,588.23,638.13,685.00,712.82,725.94,725.29,710.67,675.07,657.21,604.42,412.58,191.13,72.85,28.57,13.53,8.55,5.63,4.41,3.70,2.82,2.50,2.32,1.80]
bkg.sigma_g[*,72]=[2.17,2.75,3.17,3.83,4.39,5.30,6.08,7.18,9.76,13.78,24.02,54.60,139.41,292.18,449.03,538.39,598.36,647.66,698.97,755.00,790.47,797.57,792.33,774.53,747.53,731.84,673.59,459.30,223.11,89.02,36.69,17.91,10.22,6.69,5.18,3.90,3.33,3.02,2.71,1.92]
bkg.sigma_g[*,73]=[2.40,3.03,3.30,4.07,4.86,5.94,6.76,8.19,10.53,15.19,26.44,60.81,152.46,317.65,484.88,589.36,656.77,711.67,772.61,838.01,875.30,876.44,862.30,841.72,818.60,815.02,745.05,509.06,255.29,106.92,46.32,23.46,13.36,8.77,6.09,4.97,4.16,3.42,2.69,2.23]
bkg.sigma_g[*,74]=[2.45,3.07,3.72,4.39,5.20,6.17,7.20,8.86,11.31,16.28,28.79,65.27,165.29,348.69,534.56,650.15,721.84,786.01,857.18,936.90,975.55,964.65,940.69,911.77,891.04,895.61,807.76,554.90,290.54,129.02,58.28,29.24,17.14,11.25,8.04,6.21,4.77,4.18,3.03,2.60]
bkg.sigma_g[*,75]=[2.89,3.69,4.17,4.89,6.03,6.74,7.65,9.59,11.92,18.37,31.69,72.96,180.12,374.11,581.80,715.13,796.19,872.11,958.91,1046.02,1080.25,1062.33,1022.63,982.96,968.02,978.61,874.28,604.57,328.20,152.99,71.83,38.53,22.55,14.50,9.83,7.90,5.97,4.90,4.15,3.02]
bkg.sigma_g[*,76]=[3.24,4.03,4.56,5.22,6.43,7.36,8.54,10.62,13.30,20.55,35.44,80.68,195.06,400.66,631.21,783.94,877.35,968.64,1069.85,1162.49,1200.09,1163.01,1108.95,1064.69,1058.74,1058.97,939.30,656.29,367.46,182.92,91.49,49.12,29.22,18.36,12.60,9.56,7.15,5.93,5.07,3.78]
bkg.sigma_g[*,77]=[3.54,4.17,5.10,6.01,7.18,8.43,10.12,12.31,15.41,22.24,39.78,88.58,212.53,430.74,679.98,855.78,970.68,1080.62,1203.45,1300.27,1325.12,1274.32,1192.93,1142.33,1142.00,1136.07,988.38,701.70,408.66,215.00,111.79,62.11,37.94,23.91,16.48,12.34,9.02,7.39,6.03,4.36]
bkg.sigma_g[*,78]=[3.62,4.68,5.82,6.56,7.97,9.53,11.58,14.48,17.54,25.37,46.04,99.71,234.18,466.38,739.14,940.22,1077.27,1210.13,1354.76,1457.41,1462.95,1388.58,1295.40,1231.91,1227.92,1205.89,1040.28,745.89,447.73,245.91,136.75,79.72,48.15,31.33,22.29,16.32,11.75,9.37,7.45,5.04]
bkg.sigma_g[*,79]=[4.33,5.44,6.65,7.48,9.20,11.08,13.30,15.81,20.18,29.84,51.98,109.84,254.85,506.17,802.18,1031.37,1200.75,1361.32,1521.91,1625.20,1614.81,1508.80,1394.61,1328.32,1311.91,1268.96,1087.67,780.30,487.39,281.73,163.66,98.52,62.04,40.52,28.98,21.24,15.79,11.58,9.43,6.96]
bkg.sigma_g[*,80]=[4.92,6.57,7.71,9.05,10.81,12.74,15.18,18.26,23.42,33.44,57.43,123.14,281.46,550.72,872.79,1132.19,1334.70,1530.81,1712.79,1805.09,1767.03,1644.73,1501.32,1422.69,1389.59,1322.02,1129.08,823.27,531.89,322.23,195.49,122.51,78.94,52.40,37.45,26.88,19.29,15.80,12.31,8.75]
bkg.sigma_g[*,81]=[5.18,7.45,8.74,10.60,12.39,14.59,17.63,21.56,27.28,38.03,64.66,137.94,310.35,596.45,943.39,1242.09,1488.79,1727.87,1920.56,2003.62,1934.05,1777.87,1618.45,1520.61,1463.39,1368.69,1164.05,859.62,574.72,365.62,229.86,148.22,98.36,67.03,47.64,35.10,25.50,19.92,15.98,11.44]
bkg.sigma_g[*,82]=[6.20,8.53,10.07,12.27,14.09,17.56,20.39,25.04,32.10,44.65,75.49,156.13,339.96,652.56,1027.06,1370.32,1664.74,1941.53,2147.80,2218.59,2114.99,1917.54,1731.57,1617.73,1537.78,1411.48,1186.70,892.87,612.17,405.70,266.05,175.98,121.50,84.90,60.81,44.41,33.82,25.99,21.11,15.07]
bkg.sigma_g[*,83]=[7.51,10.40,11.67,13.94,16.02,20.06,24.60,29.25,37.57,52.24,88.59,179.04,380.76,718.64,1122.57,1513.58,1862.63,2180.36,2399.92,2438.78,2293.70,2064.97,1852.73,1715.74,1602.82,1452.27,1213.27,927.57,656.09,451.81,305.05,207.05,143.15,103.35,74.60,56.21,42.79,33.50,26.72,19.39]
bkg.sigma_g[*,84]=[9.03,12.15,14.07,16.31,20.10,24.53,29.25,35.13,44.72,62.01,104.57,207.12,427.45,788.30,1236.00,1686.75,2094.14,2446.29,2676.25,2679.24,2477.68,2212.61,1975.97,1799.55,1656.47,1480.45,1238.47,956.02,696.77,495.69,346.97,245.46,172.12,124.35,92.00,70.08,53.82,41.87,33.79,24.24]
bkg.sigma_g[*,85]=[10.15,13.53,15.82,18.72,23.22,27.96,33.59,42.57,53.78,73.68,124.26,241.00,484.08,873.59,1368.38,1876.78,2346.26,2729.41,2941.00,2909.53,2667.74,2359.92,2089.67,1887.21,1710.13,1512.84,1261.65,986.50,739.21,544.21,393.67,281.67,205.11,153.59,115.05,87.08,67.08,52.39,42.72,30.73]
bkg.sigma_g[*,86]=[11.61,16.03,18.93,22.20,26.54,32.42,39.71,50.59,64.37,91.64,149.25,283.46,554.17,982.52,1528.69,2100.62,2621.45,3027.34,3223.67,3145.50,2856.17,2518.03,2211.52,1981.02,1768.19,1540.67,1281.56,1022.08,787.06,587.82,434.99,319.24,237.98,178.93,137.29,105.55,81.89,65.24,53.18,38.48]
bkg.sigma_g[*,87]=[13.31,17.87,21.45,25.98,31.28,38.59,48.76,60.12,80.38,113.27,184.19,336.66,644.24,1116.87,1707.19,2338.11,2906.89,3340.57,3506.59,3390.81,3043.21,2669.23,2334.81,2063.25,1817.27,1575.71,1313.72,1055.42,827.30,629.31,476.45,359.67,270.92,208.21,161.11,125.48,100.21,79.21,64.29,47.15]
bkg.sigma_g[*,88]=[15.61,21.41,24.61,30.12,37.36,46.47,58.57,73.85,100.28,143.60,230.16,408.14,747.57,1262.13,1915.96,2610.47,3224.59,3651.00,3790.97,3617.27,3236.48,2821.09,2442.95,2145.58,1871.51,1607.93,1345.29,1094.54,865.76,675.31,517.95,396.56,304.32,237.60,184.00,147.26,117.51,94.32,78.06,56.35]
bkg.sigma_g[*,89]=[18.47,24.51,30.00,36.41,45.53,56.14,72.08,94.23,127.67,185.28,290.19,497.59,876.41,1434.04,2143.65,2890.54,3541.32,3968.04,4064.04,3835.99,3411.87,2963.69,2557.50,2227.94,1923.39,1642.58,1382.41,1133.63,898.83,709.90,553.66,433.49,337.57,265.60,208.85,169.24,136.45,111.93,91.64,67.48]
bkg.sigma_g[*,90]=[21.90,30.01,35.63,44.08,55.90,68.96,91.50,120.70,166.03,238.98,370.00,610.92,1029.65,1639.47,2398.41,3191.76,3859.30,4264.00,4313.78,4048.59,3584.00,3092.42,2665.28,2297.82,1972.34,1681.76,1411.28,1161.87,936.89,745.26,586.60,465.99,367.47,294.35,234.08,190.80,154.92,128.10,104.75,77.20]
bkg.sigma_g[*,91]=[25.54,35.58,42.83,53.84,68.68,87.69,114.12,153.88,214.00,307.60,470.61,750.17,1216.39,1871.68,2679.03,3496.49,4164.16,4524.28,4537.28,4226.94,3732.54,3214.33,2758.20,2362.80,2025.85,1719.48,1434.82,1178.24,959.85,770.50,616.91,496.50,397.52,319.00,258.96,209.90,170.52,140.48,116.73,87.18]
bkg.sigma_g[*,92]=[31.34,43.38,51.91,64.82,84.38,111.68,147.76,197.60,274.69,394.73,588.63,905.28,1409.02,2116.15,2949.80,3777.01,4429.06,4764.81,4711.15,4363.62,3851.84,3314.01,2825.82,2413.38,2051.61,1739.27,1442.58,1188.51,967.83,787.93,634.89,516.06,415.52,336.31,274.80,225.28,186.36,154.45,127.86,92.25]
bkg.sigma_g[*,93]=[37.22,52.18,64.86,82.55,105.74,139.83,188.41,253.90,349.64,495.32,726.63,1077.46,1617.29,2351.76,3193.87,4014.80,4630.19,4917.21,4821.96,4450.01,3932.00,3378.03,2876.23,2450.37,2073.55,1743.45,1452.11,1193.35,973.55,795.84,645.57,523.90,426.11,347.73,285.46,233.06,192.94,160.49,132.82,97.26]
bkg.sigma_g[*,94]=[45.01,62.48,77.98,100.95,132.64,173.42,233.09,318.66,436.23,608.50,862.80,1254.25,1816.11,2548.05,3380.91,4168.62,4741.83,4971.85,4850.08,4456.58,3952.54,3392.50,2890.23,2444.17,2054.54,1723.33,1428.55,1172.42,960.74,782.58,641.11,523.15,428.59,350.56,287.06,238.76,194.80,161.40,131.60,96.91]
bkg.sigma_g[*,95]=[51.62,72.61,93.59,123.30,160.56,210.04,283.07,388.62,522.56,721.55,1001.49,1406.23,1967.97,2693.26,3487.72,4212.85,4726.54,4908.52,4759.31,4366.92,3860.01,3328.83,2835.97,2391.43,1997.73,1667.43,1376.68,1129.66,926.90,756.56,621.74,508.14,417.46,337.94,277.84,231.03,190.19,153.72,126.77,91.95]
bkg.sigma_g[*,96]=[59.67,85.50,107.97,141.00,185.97,249.00,337.39,452.27,605.01,816.65,1110.16,1509.24,2056.59,2740.57,3473.23,4122.07,4551.33,4690.05,4522.33,4144.13,3664.98,3174.43,2697.98,2268.90,1890.21,1566.49,1289.19,1055.42,869.50,708.65,576.27,471.49,385.43,314.32,258.58,211.86,174.57,142.03,116.55,82.54]
bkg.sigma_g[*,97]=[65.96,96.16,121.41,157.55,209.61,278.98,374.39,498.67,666.32,881.53,1170.16,1554.63,2051.88,2653.79,3299.98,3873.26,4218.46,4296.49,4127.85,3794.14,3372.53,2916.64,2472.27,2078.34,1719.77,1422.42,1168.44,951.16,775.78,638.02,517.82,421.07,341.71,279.72,230.43,185.06,152.05,124.43,99.00,71.00]
bkg.sigma_g[*,98]=[70.27,101.49,129.99,170.43,225.84,298.83,396.87,527.00,692.78,892.44,1154.01,1493.30,1931.92,2439.44,2974.30,3442.94,3709.16,3759.63,3618.24,3331.09,2962.16,2567.00,2174.11,1820.20,1501.91,1233.38,1005.36,817.98,659.99,539.63,438.04,356.58,286.72,231.73,189.60,154.12,123.89,100.35,79.30,57.55]
bkg.sigma_g[*,99]=[62.05,88.88,116.81,152.85,200.06,264.37,347.15,463.04,591.13,753.44,951.88,1196.64,1506.25,1868.77,2244.71,2564.94,2757.90,2787.08,2678.81,2475.79,2210.10,1920.76,1623.49,1353.41,1118.99,908.42,735.27,596.85,483.52,390.02,312.99,253.33,203.83,164.99,132.70,106.00,85.69,69.84,54.61,39.18]

return,bkg
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; MAIN ROUTINE
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


function lambda_richness,z,gmr,gmr_err,imag,r,r0=r0,beta=beta,lstarcut=lstarcut,wtvals=wtvals,pvals=pvals,bkg=bkg,h0=h0,rlambda=rlambda,lambda_err=lambda_err

if n_params() eq 0 then begin
    print,'syntax- lambda = lambda_richness(z,gmr,gmr_err,imag,r,r0=r0,beta=beta,lstarcut=lstarcut,wtvals=wtvals,pvals=pvals,bkg=bkg,h0=h0,rlambda=rlambda,lambda_err=lambda_err)'
    return,-1
endif

if n_elements(r0) eq 0 then r0=1.0
if n_elements(beta) eq 0 then beta=0.2
if n_elements(lstarcut) eq 0 then lstarcut=0.2
if n_elements(bkg) eq 0 then bkg=lambda_sdss_background()
if n_elements(h0) eq 0 then h0=100.0

;; z is redshift
;; gmr, gmr_err are MODEL colors and errors
;; imag is CMODEL magnitude (although 
;; r is radius from center in h^-1 Mpc

if (n_elements(z) ne 1) then begin
    print,'ERROR: can only specify one redshift z.'
    return,-1
endif

if ((n_elements(gmr) ne n_elements(gmr_err)) or $
    (n_elements(gmr) ne n_elements(imag)) or $
    (n_elements(gmr) ne n_elements(r))) then begin
    print,'ERROR: number of elements of gmr, gmr_err, imag, r must be identical.'
    return,-1
endif

;; First, calculate the color filter G(c)

;; define the red sequence tilt and intercept as a function of redshift
m17=-0.008969-0.0701*z
b17=0.58907+3.2982*z
sigma_int = 0.05

;; calculate G(c) filter
d=gmr - (m17*(imag-17.0)+b17)
d_err = sqrt(gmr_err^2. + sigma_int^2.)
d_wt = (1./(sqrt(2.*!pi)*d_err))*exp(-(d^2.)/(2.*d_err^2.))

;; Second, calculate the NFW filter

nfw_wt = lambda_nfw_weight(r)

;; Third, the luminosity function

;; we need mstar...
;; this is valid for 0.05<z<0.35
mstar=12.27+62.36*z-289.79*z^2.+729.69*z^3.-709.42*z^4

;; and the limiting mag
limmag = mstar-2.5*alog10(lstarcut)

lum_wt = lambda_lum_func(imag,mstar,limmag)

;; We can now calculate the filter u:
ucounts = d_wt*nfw_wt*lum_wt

;; make sure we don't have any infinities
bad=where(finite(ucounts) eq 0,nbad)
if (nbad gt 0) then ucounts[bad]=0.0

;; and now we need the background
bcounts = lambda_bcounts(bkg,z,r,gmr,imag)

;; finally, the weight function for the luminosity cutoff
;; 1.0 if it's above the cut, 0.0 below
theta_i = fltarr(n_elements(ucounts))
br=where(imag lt limmag,nbr)
if (nbr gt 0) then theta_i[br]=1.0

lambda = lambda_bisection_zero(r0,beta,ucounts,bcounts,r,theta_i,wtvals=wtvals,pvals=pvals)
rlambda = r0*(lambda/100.)^beta

;; and calculate the error on lambda

bar_p = total(wtvals^2.)/total(wtvals)
lambda_err = sqrt(lambda*(1.-bar_p))

return,lambda
end