Расчет энергетического спектра электронов

%% Очистка памяти и экрана
clear; clc;
%% Начальные условия
N = 10;                     % кол-во частиц
Theta0 = 0;                 % начальный полярный угол падения в градусах
l_tot = 1;                      % средняя длина свободного пробега
V = 0.9;                    % однократное альбедо
tmax = 100;                 % максимальный пробег
q = 0.9;                    % параметр ХГ
mu_HG = @(x,q) (1+q^2-(1-q^2)/(1+q-2*q*x)^2)/(2*q);% функция розыгрыша косинуса угла рассеяния по ХГ 
alf = 0;
J = 1;
mu_Li = @(delta,alf) ((1+alf)*J^(1+alf))/delta^(2+alf); %инд Лиликвиста
mu_obrLi= @(delta) alf^(1/(-delta - 2)); %обратная функция 
Emin = 10; %Минимальная 
%% Инициализация
u = 0;          % пробег частицы
r = zeros(3,N);         % начальные координаты частиц
Theta0(Theta0 == 0)=10^-5;
d0=[sind(Theta0) 0 cosd(Theta0)];
            % начальные вектор направления
%% Расчет траектории
for i = 1:N
    l = -l_tot*log(rand(1,1));
    tt=l; dd=d0; rr=l*dd;
    iFlag=1;
    while iFlag
        alfa=rand(1,4);
        if alfa(1) < V
    mu = mu_HG(alfa(2),q);
                phi=2*pi*alfa(3);
                cosP = cos(phi); sinP = sin(phi); sinS = sqrt((1-mu^2)/(1-dd(3)^2)); 
                u = -dd(1)*dd(3)*cosP+dd(2)*sinP;  v = -dd(2)*dd(3)*cosP+dd(1)*sinP; 
                w = (1-dd(3)^2)*cosP;  dd = dd* mu + [u v w]*sinS;
                l = -l_tot*log(alfa(4));
                rr = rr + l*dd;  tt=tt+l;
                iFlag = (tt<tmax) & ((rr(3)>0));
        else
                iFlag = 0;
        end
    end
end