function [Kd_layers, R0] = MC_3D_2024(a,b,bw,theta_s,N,gamma_cutoff)
% a = .1;
% b = 0.5;
% theta_s = 0; % sun position
% N = 100000; % photons
% gamma_cutoff = 10; % angle for transport approximation

w = b/(a+b);
n_refr = 1.34;
theta_0 = asind(sind(theta_s)/n_refr); % underwater incident angle. After refraction theta_0 = asind(sind(theta_s)/1.34); 60->~40

% perallocating 
r = zeros(3,N); % photons coords
r0 = [0 0 0]'; % initial position
dr0 = [sind(theta_0) 0 cosd(theta_0)]'; % initial direction
alive = zeros(1,N); 
theta = zeros(1,N); 

%% Fresnel reflection of the incident radiation 
if theta_0==0
    rf_sun = ((n_refr-1)/(n_refr+1))^2;
else
    if theta_s == 60
        rf_sun = 0.063226; % The formula below is correct, Mobley uses the same (surf_refl.f95). However, Hydrolight for some reason uses this value...
    else
        rf_sun = 0.5*((tand(theta_s-theta_0)/tand(theta_s+theta_0))^2 + (sind(theta_s-theta_0)/sind(theta_s+theta_0))^2);
    end
end
N = round(N*(1-rf_sun),0); % N inside

%% Phase Functions
P_water = @(mu) 0.06083*(1+0.92493*(mu).^2);
g = 0.95; % 0 - isotropic 
P = @(mu) (1-g.^2)./(1+g.^2-2*g*mu).^(3/2)/(4*pi);

dmu = -0.0001;
mu = 1:dmu:-1;
beta = (P_water(mu)*bw + P(mu)*(b-bw))./b; % VSF_integral = 1
% beta = P_water(mu); % water only 
% beta = P(mu); % HG only

mu_int = [1+dmu/2:dmu:-1-dmu/2];
beta_int = interp1(mu, beta, mu_int)';

D = -2*pi*dmu*cumsum(beta_int);
D = D/max(D);

% reinterpolation on the regular grid
dD = (D(end)-D(1))/(size(D,1)-1);
mu_int_reg = interp1(D, mu_int, [D(1):dD:D(end)]);
gamma_int_reg = acosd(mu_int_reg);

% plot(mu_int_reg, [D(1):dD:D(end)],'o')
% gamma_values = acosd(mu_int);
% mu_int_reg(round(1+rand*(size(mu_int_reg,2)-1)))

% tests -> try for the new VSF
% VSF_integral = -2*pi*dmu*sum(beta_int) % should be equal to 1. '-' from inversion of limits 
% g_test = -2*pi*dmu*sum(beta_int.*mu_int') % % should be equal to g
%
% ind = mu_int<0;
% bb_b_ratio = -2*pi*dmu*(sum(beta_int(ind)))'
% bb = b.*bb_b_ratio

%% main part
n_out = 0; % out of water 
interactions = zeros(1,N);
bottom = 100;
down = zeros(1,bottom);
% up = zeros(1,bottom);

for n=1:N % runs through the entire array several times while there are photons alive

% r_prev = [(1-2*rand) (1-2*rand) 0]'; % random input position
dr = dr0;
r_prev = r0;
l = -log(rand())/(a+b);

% moving down, count for Ed
ind_d = [ceil(r_prev(3)):floor(r_prev(3)+l*dr(3))]+1; % +1 to take into account zero level
down(ind_d) = down(ind_d) +1 ;

r_prev = r0 + l*dr;

while ~isnan(alive(n))
    if w<rand % absorption (L&W 6.1)
        theta(n) = acosd(dr(3)); % = gamma;  angle before abs-n: depends on last scattering direction
        alive(n) = NaN;
        
    else % scattering
        l = -log(rand())/(a+b); % random path length, Mobley, 2021, p 775
        gamma = gamma_int_reg(round(1+rand*(size(gamma_int_reg,2)-1))); % fast way to get angles
%         gamma = acosd(1-2*rand()); % isotropic VSF without water
        phi = 360*rand();
        
        % get path
        if gamma>gamma_cutoff % transport approximation (check its value for possible difference!)
            dr = new_direction_3D(dr,gamma,phi)'; %!!!
            % dr = [sind(gamma)*cosd(phi); sind(gamma)*sind(phi); cosd(gamma)]; % without rotation for isotropic

        %     % faster way        
        %     dr =dr/norm(dr);
        %     if dr(3)==1
        %         dr = [sqrt(1-mu^2)*cos(phi); sqrt(1-mu^2)*sin(phi); mu];
        %     else
        %         dr(1) = dr(1)*mu-(dr(2)*sin(phi)+dr(1)*dr(3)*cos(phi))*sqrt((1-mu^2)/(1-dr(3)^2));
        %         dr(2) = dr(2)*mu+(dr(1)*sin(phi)-dr(2)*dr(3)*cos(phi))*sqrt((1-mu^2)/(1-dr(3)^2));
        %         dr(3) = dr(3)*mu+(1-dr(3)*dr(3))*cos(phi)*sqrt((1-mu^2)/(1-dr(3)^2));
        %     end
          
        end
        
        if dr(3) > 0 % moving down, count for Ed
            ind_d = [ceil(r_prev(3)):floor(r_prev(3)+l*dr(3))]+1; % +1 to take into account zero level 
            ind_d(ind_d>bottom) = bottom; % everything that flies below the bottom stays at the bottom
            down(ind_d) = down(ind_d) +1 ;
        end
        
        r(:,n) = r_prev +  l*dr; % new position
        if r(3)>bottom % absorption at the bottom
            theta(n) = acosd(dr(3)); % = gamma;  angle before abs-n: depends on last scattering direction
            alive(n) = NaN;
        end
        
        r_prev = r(:,n);

        if r(3, n) < 0 % across the boundary
        theta(n) = acosd(dr(3));
        theta_suppl = 180-theta(n); % supplementary angle
        
        if theta(n)<(180-asind(1/n_refr)) %  total internal reflection, critical value for n=1.34 is 131.73°
            %             theta(n) = theta_suppl;
            
            l = -log(rand())/(a+b); % new jump
            r(3,n) = 0;
            r_prev = r(:,n);
            dr(3) = -dr(3);
            
%             if dr(3) > 0 % moving down, count for Ed <- always true here
                ind_d = [ceil(r_prev(3)):floor(r_prev(3)+l*dr(3))]+1; % +1 to take into account zero level
%                 ind_d(ind_d>bottom) = bottom; % everything that flies below the bottom stays at the bottom <- always false
                down(ind_d) = down(ind_d) +1 ;
%             end
            
            r(:,n) = r_prev +  l*dr;
            r_prev = r(:,n);

        else
            theta_refr = asind(sind(theta_suppl)*n_refr);
            
            if theta_suppl==0
                rf = ((n_refr-1)/(n_refr+1))^2;
            else
                rf = 0.5*((tand(theta_suppl-theta_refr)/tand(theta_suppl+theta_refr))^2 + (sind(theta_suppl-theta_refr)/sind(theta_suppl+theta_refr))^2);
            end

            n_out = n_out+1*(1-rf); % alternative way for R / works good 
            
            if rand>rf % Refraction to the air
                r(3,n) = r(3,n)/n_refr; % geometrical simplification
                theta(n) = 180 - theta_refr;
                alive(n) = NaN;
                
            else % Reflection inside
            l = -log(rand())/(a+b); % new jump
            r(3,n) = 0;
            r_prev = r(:,n);
            dr(3) = -dr(3);
                        
            % if dr(3) > 0 % moving down, count for Ed <- always true
            ind_d = [ceil(r_prev(3)):floor(r_prev(3)+l*dr(3))]+1; % +1 to take into account zero level
            % ind_d(ind_d>bottom) = bottom; % everything that flies below the bottom stays at the bottom <- always false
            down(ind_d) = down(ind_d) +1 ;
            % end
            
            r(:,n) = r_prev +  l*dr;
            r_prev = r(:,n);
            end
        end
    end
    end
    interactions(n) = interactions(n) + 1; % number of scaterrings for every photon
end
end
% clearvars theta_refr theta_suppl rf

if any(imag(theta) ~= 0)
    'theta is complex, find a reason above'
    return
end
% mean_interactions = mean(interactions) % ~(a+b)/a

%% R
ind = r(3,:)<0;
Eu0 = sum(ind) + N/(1-rf_sun)*rf_sun; % emerged and reflected parts 
% Eu0 = n_out + N/(1-rf_sun)*rf_sun; % emerged and reflected parts // alternative
R0 = Eu0/(N/(1-rf_sun));

%% rigorous Kd calculation method taking into account passage through layers
down2 = down(down>down(1)/10); % in z_0.1 layer 
% cftool([0:2:length(down2)-1], log(down2(1:2:end)))
Ed_fit = robustfit([0:length(down2)-1], log(down2)); % bisquare
Kd_layers = abs(Ed_fit(2));

end
