clear all; close all;

%% parameters to be set
n = 100;
r_p_to_lambda = linspace(5, 20, n);
L_p_to_r_p = 5;
zeta_0_factor = 1;

%% fixed parameters and constants
lambda = 1e-9;
D = 1e-9;
mu = 1e-3;
T = 295;
k_B_const = physconst('Boltzmann');
e_const = 1.60217663e-19;
epsilon0_const = 8.854187e-12;
epsilon_r = 80;

C = 1/6;

%% derived parameters
r_p = r_p_to_lambda*lambda;
L_p = L_p_to_r_p*r_p;
Phi_th = k_B_const*T/e_const;
epsilon = epsilon0_const*epsilon_r;
zeta_0 = zeta_0_factor*Phi_th;
tau_s = lambda*r_p/D;
tau_l = 8/pi^2*L_p.^2*lambda./(r_p*D).*(1+pi^3/16*r_p./L_p);
tau_TLM = 0.5*L_p.^2*lambda./(r_p*D);

%% position to be evaluated
z = r_p;

%% find the max. velocity
for i=1:n
    r = linspace(0, r_p(i), 1000);
    t = tau_s(i);     % set t to some value >0 at the beginning, since the radial position of the maximum is independent of t

    % analytical solution
    Ez = -C*zeta_0/r_p(i);

    w_r = Ez*epsilon*zeta_0/mu.*(besseli(0,r/lambda)/besseli(0,r_p(i)/lambda)-1-4*(lambda/r_p(i)*besseli(1,r_p(i)/lambda)/besseli(0,r_p(i)/lambda)-0.5).*(1-(r/r_p(i)).^2));
    w = w_r * (1-erf(z/L_p*sqrt(tau_TLM(i)/t))) * (1-exp(-t/tau_s(i))) * exp(-t/tau_l(i));

    rho_v = -epsilon*zeta_0/lambda^2.*besseli(0,r/lambda)/besseli(0,r_p(i)/lambda) * (1-erf(z(i)/L_p(i)*sqrt(tau_TLM(i)/t))) * (1-exp(-t/tau_s(i)));

    iz = rho_v .* w;

    % radial position of the minimum
    [~, r_idx] = min(iz);
    r_min(i) = r(r_idx);

    % find the position in time of the minimum of iz(r=r_min, z=r_p, t)
    t = linspace(0, 2*tau_l(i), 1000);

    w_r = Ez*epsilon*zeta_0/mu*(besseli(0,r_min(i)/lambda)/besseli(0,r_p(i)/lambda)-1-4*(lambda/r_p(i)*besseli(1,r_p(i)/lambda)/besseli(0,r_p(i)/lambda)-0.5)*(1-(r_min(i)/r_p(i))^2));
    w = w_r .* (1-erf(z/L_p*sqrt(tau_TLM(i)./t))) .* (1-exp(-t/tau_s(i))) .* exp(-t/tau_l(i));

    rho_v = -epsilon*zeta_0/lambda^2*besseli(0,r_min(i)/lambda)/besseli(0,r_p(i)/lambda) .* (1-erf(z(i)/L_p(i).*sqrt(tau_TLM(i)./t))) .* (1-exp(-t/tau_s(i)));

    iz = rho_v .* w;

    [iz_min(i), t_idx] = min(iz);
    t_min(i) = t(t_idx);
end

%% export
result = [r_p_to_lambda.' iz_min.'];
plot(r_p_to_lambda, iz_min)
writematrix(result, "analyt_min(iz)_over_rp_to_lambda.txt", 'Delimiter','tab');