clear all; close all;

%% parameters to be set
r_p_to_lambda = 5;
L_p_to_r_p = 10;
zeta_0_factor = 1;

%% fixed parameters and constants
lambda = 1e-9;
r_blunt = 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;
t = tau_l;

%% analytical solution
r = linspace(0,r_p,1000);

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

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

% i_z(r,z,t)
i_z = rho_v .* w;

%% export the analytical solution
r_analyt = r/r_p;
iz_analyt = i_z;

analyt = [r_analyt.' iz_analyt.'];

writematrix(analyt, "analyt_r_over_t.txt", 'Delimiter','tab');