;
; compute pure seawater scattering coefficient
; 
; based on:
; Buiteveld et al. 1994, SPIE Ocean Optics XII 2258, 174-183
; Twardowski et al. 2007, Biogeosciences 4, 1041-1058
; references listed below expanded upon in above documents
;
; bsw = b_sw(wl, S, T, betaw=betaw, bw=bw, betasw=betasw)
; where wl is wavelength (nm), S is salinity (PSU), and T is temperature (C)
; and keywords betaw, bw, betasw provide optional output products
;
; jeremy.werdell@nasa.gov, 7 Aug 2008
;

function b_sw, wl, S, T, betaw=betaw, bw=bw, betasw=betasw

rad = !dtor * findgen(18000) * 0.01

; Boltzmann constant
;
k = 1.38054e-23

; depolarization ratio from Farinato and Roswell 1975
;
dr = 0.051  

; refractive index from Mcneil 1977
; salinity dependent term removed (Morel 1974 enhancement for salinity applied below)
;
nw = 1.3247 + 3.3e3*wl^(-2) - 3.2e7*wl^(-4) - 2.5e-6*T^2

; isothermal compression from Lepple and Millero 1971
;
ic = (5.062271 - 0.03179*T + 0.000407*T^2) * 1e-10

; pressure derivative of refractive index
; from O'Conner and Schlupf 1967 and Evtyushenkov and Kiyachenko 1982
;
dndp_633_20 = (-0.000156 * 633. + 1.5989) * 1e-10
dndp1 = (-0.000156 * wl + 1.5989) * 1e-10
dndp2 = (1.61857 - 0.005785 * T) * 1e-10
npd = (dndp1 * dndp2) / dndp_633_20

; Einstein-Smoluchowski theory
;
betaw90 = 2 * !dpi^2 * k * (273+T) * nw^2 * npd^2 / (wl*1e-9)^4 / ic * (6+6*dr) / (6-7*dr)

betaw = betaw90 * (1 + cos(rad)^2 * (1-dr) / (1+dr))
bw = 8 * !dpi / 3 * betaw90 * (2+dr) / (1+dr)

; Morel 1974 enhancement for 37 ppt seawater
;
betasw = betaw * (1 + 0.3 *S / 37.)
bsw = bw * (1 + 0.3 * S / 37.) 

return, bsw
end