ang_dist [Biomedical Optics]

clear all; close all; clc;
%%
lambda = 400*1e-9;
rayleghCondition = 0.5*lambda;
n_s = 1.57;
n_b = 1.33;
w_s = 1.05*1e3;
w_b = 1.0*1e3;
concentration = 0.002;

dSphere = 50:100:6500;
diameter = zeros(1,length(dSphere));
itera = 0;

for iter = 1:length(dSphere)
    diameter(iter) = dSphere(iter)*1e-9;
    radius(iter) = diameter(iter)/2;
    %%
    k = 2*pi*n_b/lambda;
    x(iter) = k*radius(iter);
    n_rel = n_s/n_b;
    y = n_rel*x(iter);
    % Calculate the summations
    err = 1e-8;
    Qs = 0;
    gQs = 0;
    %%

    for n = 1:100000
        Snx = sqrt(pi*x(iter)/2)*besselj(n+0.5,x(iter));
        Sny = sqrt(pi*y/2)*besselj(n+0.5,y);
        Cnx = -sqrt(pi*x(iter)/2)*bessely(n+0.5,x(iter));
        Zetax = Snx+i*Cnx;
        % Calculate the first-order derivatives
        Snx_prime = - (n/x(iter))*Snx+sqrt(pi*x(iter)/2)*besselj(n-0.5,x(iter));
        Sny_prime = - (n/y)*Sny+sqrt(pi*y/2)*besselj(n-0.5,y);
        Cnx_prime = - (n/x(iter))*Cnx-sqrt(pi*x(iter)/2)*bessely(n-0.5,x(iter));
        Zetax_prime = Snx_prime + i*Cnx_prime;
        an_num = Sny_prime*Snx-n_rel*Sny*Snx_prime;
        an_den = Sny_prime*Zetax-n_rel*Sny*Zetax_prime;
        an = an_num/an_den;
        bn_num = n_rel*Sny_prime*Snx-Sny*Snx_prime;
        bn_den = n_rel*Sny_prime*Zetax-Sny*Zetax_prime;
        bn = bn_num/bn_den;
        Qs1 = (2*n+1)*(abs(an)^2+abs(bn)^2);
        Qs = Qs+Qs1;
        if n > 1
            gQs1 = (n-1)*(n+1)/n*real(an_1*conj(an)+bn_1*conj(bn))+(2*n-1)/((n-1)*n)*real(an_1*conj(bn_1));
            gQs = gQs+gQs1;
        end
        an_1 = an;
        bn_1 = bn;
        if abs(Qs1)<(err*Qs) & abs(gQs1)<(err*gQs)
            break;
        end
    end
    %%
    QsArr(iter) = (2/x(iter)^2)*Qs;

    if diameter(iter) < rayleghCondition
        itera = itera + 1;
        QsRay(itera) = 8*x(itera)^4/3*abs((n_rel^2 - 1)/(n_rel^2 + 2))^2;
    end

    gQsArr(iter) = (4/x(iter)^2)*gQs;
    g(iter) = gQsArr(iter)/QsArr(iter);
    alpha(iter) = acosd(g(iter));
    vol_s(iter) = 4*pi/3*radius(iter)^3;
    N_s = concentration*w_b/(vol_s(iter)*w_s);
    sigma_s(iter) = QsArr(iter)*pi*radius(iter)^2;
    mu_s(iter) = N_s*sigma_s(iter);
    mu_s_prime(iter) = mu_s(iter)*(1-g(iter));
    fprintf('%d\n',iter);

     dndomega1(iter) = 0.06225*(1+0.835*g(iter)*g(iter));
     dndomega2(iter) = (1-g(iter)*g(iter))/(4*pi*(1+g(iter)*g(iter)-2*g(iter)*g(iter))^1.5);
end

figure()
loglog(x,QsArr);
hold on
loglog(x(1:itera),QsRay);
title('Scattering efficiency versus x=ka')
xlabel('x = ka');
ylabel('Scattering efficiency');
grid on
hold off

figure()
semilogx(x,g);
title('Anisotropy g versus x=ka')
xlabel('x = ka');
ylabel('Anisotropy g');
grid on

figure()
semilogx(x,sigma_s);
title('Scattering cross section versus x=ka')
xlabel('x = ka');
ylabel('Scattering cross section');
grid on


figure()
semilogy(alpha,dndomega1);
title('Scattering in pure sea water')
xlabel('Angle alpha');
ylabel('dN/dOmega');
grid on

figure()
semilogy(alpha,dndomega2);
title('Henyey-Greenstein function :')
xlabel('Angle alpha');
ylabel('dN/dOmega');
grid on


% hold on
% loglog(x(1:itera),QsRay);




% Output results
% {'wavelength(nm)','Qs ( - )','g (-)','mus (/cm)','mus_prime(/cm)';lambda*1e9,Qs,g,mu_s*1e-2,mu_s_prime*1e-2}