

function [beta, bb_tilde] = ff(mu,n);
%function returns the Fournier Forand phase function every degree. 
%n is the real part of the index of refraction, mu is the power-law differential number distribution slope. 


nu=(3-mu)/2;
phi=[0.01 1:1:179 179.99]*pi/180;
mm=length(phi);
u=2*sin(phi/2);
delta=u.^2/(3*(n-1)^2);
delta_pi=(2*sin(pi/2))^2/(3*(n-1)^2);

beta=1./(4*pi*((1-delta).^2).*(delta.^nu)).*(nu*(1-delta)-(1-delta.^nu)+(4./(u.^2)).*(delta.*(1-delta.^nu)-nu*(1-delta)));
beta=beta+beta(mm)*(3*cos(phi).^2-1)/4;


%B=1-((1-delta_pi^(nu+1))-0.5*(1-delta_pi^nu))/((1-delta_pi)*delta_pi^nu);

sin_factor=sin(phi);

dtheta=phi(3)-phi(2);
b_back=2*pi*sum(beta(mm/2+1:mm).*sin_factor(mm/2+1:mm))*dtheta;
b=2*pi*dtheta*(sum(beta(1:mm).*sin_factor(1:mm))-0.5*(beta(1)*sin_factor(1)+beta(181)*sin_factor(181)));
bb_tilde=b_back./b;

