%function res=BRDF(teta_v,teta_s,az)

%   Hapke BRDF reflectivities
%   teta_v - view zenith angle from surface to sensor
%   teta_s - is yhe source zenith angle at the surface
%   az - is the view to source relative azimuth angle from the surface

% parameters:
% omega - single scatter albedo
% g - asymetry parameter
% h - 
% Sh - 
omega=0.99;
g=0.6;
h=0.0;
Sh=0.995;

teta_v=0;
teta_s=0*pi/180;
az=0;

cos_r=cos(teta_s).*cos(teta_v)+sin(teta_s).*sin(teta_v).*cos(az);

P_HG=(1-g^2)/(1+g^2+2*g*cos_r).^1.5;    % Henyey-Greenstain function

B=(1-g)/(1+g^2)*(1+sqrt((1+cos_r)*1-cos_r)/h);

x=cos(teta_v);
Hv=(1+2*x)/(1+2*x*sqrt(1-omega));

x=cos(teta_s);
Hs=(1+2*x)/(1+2*x*sqrt(1-omega));

res=0.25*omega/(cos(teta_v)+cos(teta_s))*((1+Sh/(omega*B))*P_HG+Hv*Hs-1);

