function phi_0=fluence(mua,mus,g,N)
%***************************************
r1 = 0.01;                                 %start position
r2 = 3;                                    %end position
dr = 0.01;                                 %step size
n_end = round((r2-r1)/dr);                 %number of detectors
%***************************************
l=0:N;
mua(l+1) = mua + (1-g.^l).*mus./(1-g);     %absorption coefficients of order l
%***************************************
k=1:N;
A=diag(k,1)+diag(k,-1);                   
B=diag(1./sqrt((2*l+1).*mua(l+1)),0);
[V,D] = eig(B*A*B);
%***************************************
s = 3;                                     %dimension
for m=0:n_end
    r(m+1) =r1 + m*dr;
    phi_0=0;
    for k=1:N+1
        if (D(k,k) > 0)
             %three-dimensional fluence
             phi_0 = phi_0 + 1./(2*pi*mua(1)*r).*exp(-r./D(k,k))./D(k,k).^2.*V(1,k).^2;
             
             %s-dimensional fluence
             %phi_0 = phi_0 + V(1,k).^2./(mua(1)*pi*D(k,k).^2)./(2*pi*D(k,k).*r).^(s/2-1).*besselk(s/2-1,r/D(k,k));
        end
    end
end
semilogy(r,phi_0,'black');