function  gamma_s  = scat( model, g )
%Function computing the scattering angle

q=rand;

if model == 1 %isotropic
    gamma_s=1-2*q;
end
if model == 2 %H-G
    if g==0;
        gamma_s=1;
    else
        gamma_s=(2*g+1-2*q*g-2*q+2*q^2*g+g^2-g^3*q-2*q*g^2+2*g^3*q^2)/(1+g-2*q*g)^2;
    end
end


if model == 3 %Petzold
    %Petzold commulative distribution function, angle, anc cosine of the anglem:
    P=[0	0.517	0.666	0.751	0.811	0.851	0.879	0.901	0.918	0.931	0.942	0.95	0.9832	0.9853	0.987	0.9888	0.9902	0.9917	0.9928	0.9941	0.995	0.996	0.9968	0.957	0.963	0.968	0.972	0.975	0.978	0.981	0.9976	0.9982	0.9988	0.9992	0.9997	0.9998	1];
    theta_s=[0	2.5	7.5	12.5	17.5	22.5	27.5	32.5	37.5	42.5	47.5	52.5	92.5	97.5	102.5	107.5	112.5	117.5	122.5	127.5	132.5	137.5	142.5	57.5	62.5	67.5	72.5	77.5	82.5	87.5	147.5	152.5	157.5	162.5	167.5	172.5	180];
    mu_s=[1	0.999048	0.991445	0.976296	0.953717	0.92388	0.887011	0.843391	0.793353	0.737277	0.67559	0.608761	-0.04362	-0.13053	-0.21644	-0.30071	-0.38268	-0.46175	-0.5373	-0.60876	-0.67559	-0.73728	-0.79335	0.5373	0.461749	0.382683	0.300706	0.21644	0.130526	0.043619	-0.84339	-0.88701	-0.92388	-0.95372	-0.9763	-0.99144	-1];
    gamma_s=interp1(P,mu_s,q,'cubic');
end


end

