; =======================================================================================
  pro calc_xy,sva,ssa,sza,alb_gary,x,y,sigma,cden_o3,flag,  $
              slope=slope,intercept=intercept,fix_sigma=fix_sigma
; =======================================================================================
csec_ray  = 9.708e-26  ; cm2
csec_o3   = 9.261e-3   ; cm^2/peta 

alb=alb_gary*1e-6
npts=n_elements(sza)
alt=sza_alt(sza)

oomu2=chap(alt,sza,5) 

alt0=[40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0, 80.0, 85.0, 90.0]
cdtot=[6.65061e+22,3.53402e+22,1.90120e+22,1.01581e+22,5.24774e+21,  $
       2.53401e+21,1.09209e+21,3.96918e+20,1.21292e+20,4.13237e+19,1.72132e+19]
cdtot=alog(cdtot)

z0=alt < 90
cden_ref=exp(interpol(cdtot,alt0,z0))

rp=reform(rayleigh_phase(ssa))
mu1 = cos(sva / !radeg)
mucoeff=(1./mu1 + oomu2)

; Initialize output

flag=0
cden_o3=0.
sigma=0.
x=fltarr(npts)
y=fltarr(npts)

if keyword_set(fix_sigma) then begin  ; no need to do X/Y fit

  sigma=fix_sigma 
  cden=(alb-alb)/0.

  ok=where(finite(alb) and finite(mu1),n_ok)

  if n_ok gt 1 then begin
    cden(ok)=((rp(ok)*gamma(sigma+1)*csec_ray*cden_ref(ok))/(alb(ok)*(mu1(ok)*mucoeff(ok)^sigma)*(csec_o3^sigma)))^(1./sigma)
    crit=where(finite(cden) eq 1,nn)
    cden_o3=total(cden[crit])/float(nn)
  endif

endif else begin

  x = alog(mucoeff)
  y = alog(mu1*alb/rp)

  ok=where((finite(x)) and (finite(y)),n_ok)
  if n_ok gt 1 then begin
    lf,x(ok),y(ok),slope,intercept,r_corr
    sigma=-1.*slope
    if(sigma lt 0.2) then flag=1

    avg_cd=total(cden_ref[ok])/float(n_ok)
    cden_o3=((exp(-intercept)*gamma(sigma+1.)*csec_ray*avg_cd)^(1./sigma))/csec_o3
  endif

endelse

if(finite(cden_o3) ne 1 or cden_o3 gt 50) then flag=1

return
end