function bs_alb = modis_bs_alb_mcd43b1(fiso,fvol,fgeo,sza)

% fiso,fvol,fgeo,sza need to have the same dimension

sza2 = sza.^2;
sza3 = sza2.*sza;
bs_alb = fiso+fvol.*(-0.007574-0.070987*sza2+0.307588*sza3)+...
    +fgeo.*(-1.284909-0.166314*sza2+0.041840*sza3);