; --------------------------------------------------------------------------
; Function calc_albedo
; --------------------------------------------------------------------------
; PURPOSE:
;    This function calculates atmospheric albedo for 265 nm given ozone column density
;    at 50 km, and the ratio of the ozone and air scale heights
;
; INPUT PARAMETERS:
;   ozone_cden - Ozone column density above peak alt, Peta/cm^2 (1Peta=1e15 molecules)
;   hratio     - Ratio of ozone number density scale height to neutral atmosphere scale height
;   rp         - Rayleigh phase
;   cva        - view angle at rayleigh peak, degrees
;   sza        - Solar zenith angle, degrees
;
; RETURN VALUE (if applicable):
;   Albedo in Garys (1G=albedo of 1e-6)
;
; REFERENCES:
;   CIPS Algorithm Theoretical Basis Document
;
; TODO:
;   comment
;   add error handling
;   under what circumstances will this fail?  how will that affect calling procedure?
;
function calc_albedo,ozone_cden,hratio,rp,cva,sza

csec_ray  = 9.708e-26  ; Rayleigh scattering cross section, cm^2/molecule
csec_o3   = 9.261e-3   ; Ozone cross section, cm^2/peta

alt=sza_alt(sza)
npts=n_elements(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))

;This depends solely on C and sigma, and the above constants, not on geometry
ratio = csec_ray * cden_ref / ((csec_o3 * ozone_cden) ^ hratio)

mu1 = cos(cva*!const.dtor)
mucoeff=(1./mu1+oomu2)

;Set up albedo to be an array of the correct size
albedo=mucoeff*0.

;Handle bad results for SZA>90 (really, mucoeff<0)
whereok=where(mucoeff gt 0, okcount,complement=wherenotok,ncomplement=notokcount)
if n_elements(ratio) gt 1 then begin
  w_ratio=whereok
end else begin
  w_ratio=0
end
if n_elements(ozone_cden) gt 1 then begin
  w_o3cden=whereok
end else begin
  w_o3cden=0
end
if n_elements(hratio) gt 1 then begin
  w_hratio=whereok
end else begin
  w_hratio=0
end

if okcount gt 0 then begin
  albedo[whereok] = rp[whereok]*gamma(hratio[w_hratio]+1.0)*ratio[w_ratio]/ $
                  (mu1[whereok]*(mucoeff[whereok]^hratio[w_hratio]))
end
if notokcount gt 0 then begin
  albedo[wherenotok]=0
end

return, albedo*1e6
end