function [s_out, Pulse] = calculate_smm_propagation(s_in, pulse_param, t_list, t_birth_param, options)
% simulation of electron propagation in electric field (SMM)
% Input: relevant parameters (s_in), relevant pulse parameters
% (pulse_param), time list (t_list [fs]), birth time list (t_birth_param
% [fs])
arguments
    s_in
    pulse_param
    t_list
    t_birth_param
    options.use_input_field (1,1) {mustBeNumericOrLogical} = 0              % use a defined field E(t) instead of pulse parameters
    options.Efield (1,:) double {mustBeNumeric} = zeros(size(t_list))       % defined field E(t), used when use_input_field == 1 
    options.nonlinearity (1,1) double {mustBeNumeric} = 0                   % for ionization rate: p = E^nonlinearity. If 0: WKB ionization rate
    options.near_field_power_law (1,1) double {mustBeNumericOrLogical} = 0  % use power law or exponential decay for near-field decay
end

load_const;
load_cm;
load_color;

% generate electric field from laser parameters
Pulse = get_laserfield("intensity",pulse_param.intensity, "wavelength", pulse_param.wavelength, "tlist", t_list, "t0", pulse_param.t0, "phase", pulse_param.phi, "width", pulse_param.width);

% near-field enhancement and static DC field
if options.near_field_power_law
    R = s_in.lambda_nf / 0.4;
    w_foc = 8000;
    a = 2*log(2)/(2*w_foc)^2;

    gamma = @(x) (x>=0).*((s_in.gamma_0-1)*(R./(x+R)).^3 + exp(-2*log(2)*x.^2 /(2*w_foc)^2));
    gamma_int = @(x) (x>=0).*((s_in.gamma_0-1)*(-0.5)*R^3./(x+R).^2 + sqrt(pi)*erf(sqrt(a)*x)/(2*sqrt(a)) - (s_in.gamma_0-1)*(-0.5)*R^3./(R).^2);

else
    gamma = @(x) (x>=0).*(1 + (s_in.gamma_0-1)*exp(-x/s_in.lambda_nf));
    gamma_int = @(x) (x>=0).*(x-s_in.lambda_nf * (s_in.gamma_0 - 1) * (exp(-x/s_in.lambda_nf) -1));
end

Estat = @(x) s_in.U_0 ./ (s_in.r_tip * (1 + x/s_in.r_tip).^2);
Ustat = @(x) s_in.U_0 ./ (1 + x/s_in.r_tip);

if options.use_input_field
    if ~any(options.Efield)
        warning('Input field is selected, but no input field is transferred.')
    end
else
    if any(options.Efield)
        warning('Input field is not selected, but input field is transferred.')
    end
    options.Efield = Pulse.E;
end

% create birth time list
dt = t_list(2) - t_list(1);
s_in.roi_tbirth = false(size(t_list));
[~, closest_indices] = arrayfun(@(t) min(abs(t_list - t)), t_birth_param);
s_in.roi_tbirth(closest_indices) = true;
s_in.roi_tbirth = s_in.roi_tbirth & options.Efield < -0.01*max(abs(options.Efield));
s_in.t_birth = t_list(s_in.roi_tbirth);
s_in.Ebirth = options.Efield(s_in.roi_tbirth);

if isempty(s_in.Ebirth)
    disp('no trajectories');
    s_in.rate = []; s_in.x = []; s_in.v = []; s_in.Ekin = []; s_in.n_recoll = []; s_in.t_recoll = []; s_in.S.W = []; s_in.S.T = []; s_in.S.V = []; s_in.S.tot = []; s_in.S_rec = []; 
    s_in.cycle.tmin = []; s_in.cycle.tmax = []; s_in.S_korr = []; s_in.roi_recoll = []; s_in.roi_direct = [];
    s_in.Up = Pulse.Up * s_in.gamma_0^2;
    s_in.w = Pulse.w;
    s_in.Ephot = Pulse.Ephot;
    s_out = s_in; 
    return
end

% ionization rate (weight of trajectories)
if options.nonlinearity == 0
    s_in.rate = p_WKB(abs((s_in.Ebirth)*s_in.gamma_0), s_in.W) * dt;
else
    s_in.rate = (abs(s_in.Ebirth)*s_in.gamma_0).^options.nonlinearity*dt;
end

% phase contribution bound state (befor birth)
s_in.S.W = -s_in.W*(s_in.t_birth-t_list(1))*const.e*1e-15;                           % classical action [Js] via workfunction (W)

% separate birth time cycles/windows
t_min = find(abs(diff([-1000, s_in.t_birth])) > 0.15*pi/Pulse.w);
t_max = find(abs(diff([s_in.t_birth, 1000])) > 0.15*pi/Pulse.w);
s_in.cycle.t_min = s_in.t_birth(t_min);
if s_in.t_birth(1) == 0 && isempty(s_in.cycle.t_min)
    s_in.cycle.t_min = 0;
end
s_in.cycle.t_max = s_in.t_birth(t_max);

% set matrices
s_in.x = zeros(size(s_in.t_birth));                           % x-values for each electron [Å]
s_in.v = zeros(size(s_in.t_birth));                           % velocity (in x-dir) for each e [Å/fs]
s_in.Ekin = zeros(size(s_in.t_birth));                        % kinetic energy of each electron [eV]
s_in.n_recoll = zeros(size(s_in.t_birth));                    % number of recollisions for each electron (at x = 0)
s_in.t_recoll = nan(size(s_in.t_birth));                      % time of last recollision
s_in.S.T = zeros(size(s_in.t_birth));                           % classical action [Js] via integration of kinetic energy (int Ekin dt)
s_in.S.V = zeros(size(s_in.t_birth));                           % classical action [Js] via integration of potential energy (int V dt)
s_in.S_rec = zeros(size(s_in.t_birth));
SW_test = zeros(size(s_in.t_birth));

% set parameters
s_in.Up = Pulse.Up * s_in.gamma_0^2;
s_in.w = Pulse.w;
s_in.Ephot = Pulse.Ephot;
Epoint = options.Efield(1)*gamma(s_in.x) + Estat(s_in.x);                       % calc field for given time (incl. nearfield)

% time loop
for ti = 1:numel(t_list)
    t = t_list(ti);

    if mod(ti,5000) == 0
        fprintf('t = %0.2f fs at phi = %0.2fpi\n', t, Pulse.phase/pi);
    end

    % set all electrons not born yet to zero (or tunneling exit if field is negative)
    roi = s_in.t_birth > t;

    if options.Efield(ti) < 0
        s_in.x(roi) = -s_in.tunnel_exit * s_in.W/(options.Efield(ti)*s_in.gamma_0);
    else
        s_in.x(roi) = 0;
    end
    s_in.v(roi) = 0;
    s_in.S.T(roi) = 0;
    s_in.S.V(roi) = 0;
    s_in.S_rec(s_in.n_recoll == 0) = 0;
    SW_test(roi) = SW_test(roi) - s_in.W*dt*const.e*1e-15;

    % Velocity Verlet Algorithm
    xold = s_in.x;
    s_in.x = s_in.x + s_in.v*dt - 0.5*const.e/const.m_e * (Epoint*1e-10) * dt.^2;

    Eold = Epoint;
    Epoint = options.Efield(ti)*gamma(s_in.x) + Estat(s_in.x);

    s_in.v = s_in.v - 0.5*const.e/const.m_e * (Epoint + Eold)*1e-10 * dt;
    s_in.Ekin = 0.5*const.m_e*(s_in.v*1e5).^2 / const.e;  % [eV]

    % for recollisions:
    if(s_in.recoll_on)
        roi = xold .* s_in.x < 0;
        s_in.rate(roi) =  s_in.rate(roi).*s_in.recoll_prob(s_in.Ekin(roi));
        s_in.x(roi) = -s_in.x(roi);
        s_in.v(roi) = -s_in.v(roi);
        s_in.n_recoll(roi) = s_in.n_recoll(roi) + 1;
        s_in.t_recoll(roi) = t;

        % Reflection - phase shift of pi
        s_in.S.T(roi) = s_in.S.T(roi) + pi * const.hbar;
    end
    
    % phase contribution kinetic part
    s_in.S.T = s_in.S.T + s_in.Ekin*const.e * dt * 1e-15;  % [Js]
    % phase contribution near-field potential part
    V_nf = gamma_int(s_in.x) * options.Efield(ti) * const.e;
    s_in.S.V = s_in.S.V - V_nf * dt * 1e-15;

    if s_in.plot_trajectory
        roi_save = s_in.t_birth < t;
        s_in.x_save(roi_save,ti) = s_in.x(roi_save);
        s_in.x_save(~roi_save,ti) = nan;
        s_in.v_save(roi_save,ti) = s_in.v(roi_save);
        s_in.v_save(~roi_save,ti) = nan;
        s_in.S_save(roi_save,ti) = -s_in.S.W(roi_save)+s_in.S.T(roi_save)+s_in.S.V(roi_save);
        s_in.S_save(~roi_save,ti) = -SW_test(~roi_save);
        s_in.SW_save(roi_save,ti) = -s_in.S.W(roi_save);
        s_in.SW_save(~roi_save,ti) = -SW_test(~roi_save);
        s_in.ST_save(roi_save,ti) = s_in.S.T(roi_save);
        s_in.ST_save(~roi_save,ti) = nan;
        s_in.SV_save(roi_save,ti) = s_in.S.V(roi_save);
        s_in.SV_save(~roi_save,ti) = nan;
    end

end

s_in.S.Wtest = SW_test;

% fix recollision of slow electrons manually
roi_coll = s_in.v < 0;
dt_coll = -(s_in.x(roi_coll))./s_in.v(roi_coll); %fs
t_coll = t + dt_coll;
dx = s_in.v(roi_coll).*dt_coll;
s_in.x(roi_coll) = s_in.x(roi_coll) + dx;
s_in.v(roi_coll) = -s_in.v(roi_coll);
s_in.t_recoll(roi_coll) = t_coll;
s_in.n_recoll(roi_coll) = s_in.n_recoll(roi_coll) + 1;
s_in.S.T(roi_coll) = s_in.S.T(roi_coll) + pi * const.hbar;
s_in.S.T(roi_coll) = s_in.S.T(roi_coll) + s_in.Ekin(roi_coll)*const.e.*dt_coll*1e-15;

% unitize phases (not the same position at the end and propagation time)
s_in.Ekin = s_in.Ekin - Ustat(s_in.x);
omega_fin = s_in.Ekin*const.e/const.hbar;
k_fin = sqrt(2*const.m_e*s_in.Ekin*const.e)/const.hbar;
s_in.S_korr = (k_fin.*s_in.x*1e-10 - omega_fin*t*1e-15)*const.hbar;
s_in.S.T = s_in.S.T - s_in.S_korr;
s_in.S.tot = -s_in.S.W + s_in.S.T + s_in.S.V;

s_in.roi_direct = s_in.n_recoll == 0;
s_in.roi_recoll = s_in.n_recoll == 1;

% plots
if s_in.plot_trajectory
    Ekin = 0.5*const.m_e*(s_in.v_save*1e5).^2 / (const.e * s_in.Up);
    Ekin(isnan(Ekin)) = -s_in.W/s_in.Up;

    x_plot = -50:0.1:100;
    roi_tplot = t_list > -50 & t_list < 50 & mod(1:numel(t_list), 10) == 0;
    for i = 1:numel(x_plot)
        xi = x_plot(i);
        Eplot(:,i) = gamma(xi)*options.Efield(roi_tplot);
    end

    figure
    imagesc(t_list(roi_tplot), x_plot, Eplot.')
    colormap(cm.asym_bright)
    clim([-2 2])
    hold on
    if sum(s_in.roi_direct) > 0
        plot(t_list, s_in.x_save(s_in.roi_direct,:), 'LineWidth', 1.5, 'Color', my_B)
    end
    if sum(s_in.roi_recoll) > 0
        plot(t_list, s_in.x_save(s_in.roi_recoll,:), 'LineWidth', 1.5, 'Color', my_R)
    end
    if sum(~s_in.roi_recoll&~s_in.roi_direct) > 0
        plot(t_list, s_in.x_save(~s_in.roi_recoll&~s_in.roi_direct,:), 'LineWidth', 0.5, 'Color', my_G)
    end
    ylim([0 100])
    xlim([-2 10])
    ylabel 'x [Å]'
    xlabel 't [fs]'
    set(gca, 'FontSize', 16)

end

if s_in.make_plots
    figure
    subplot(3,1,1)
    plot(t_list, options.Efield, 'Color', my_B)
    hold on
    plot(s_in.t_birth, s_in.Ebirth, 'Color', my_R)
    plot(s_in.cycle.t_min, 0, 'o', s_in.cycle.t_max, 0, 'o')
    xlim([s_in.t_birth(1) s_in.t_birth(end)])

    subplot(3,1,2)
    plot(s_in.t_birth(s_in.roi_direct), s_in.Ekin(s_in.roi_direct)/s_in.Up, '.', 'Color', my_B)
    hold on
    plot(s_in.t_birth(s_in.roi_recoll), s_in.Ekin(s_in.roi_recoll)/s_in.Up, '.', 'Color', my_R)
    plot(s_in.t_birth(~s_in.roi_recoll&~s_in.roi_direct), s_in.Ekin(~s_in.roi_recoll&~s_in.roi_direct)/s_in.Up, '.', 'Color', my_G)
    xlim([s_in.t_birth(1) s_in.t_birth(end)])

    subplot(3,1,3)
    plot(s_in.t_birth(s_in.roi_direct), s_in.rate(s_in.roi_direct), '.', 'Color', my_B)
    hold on
    plot(s_in.t_birth(s_in.roi_recoll), s_in.rate(s_in.roi_recoll), '.', 'Color', my_R)
    plot(s_in.t_birth(~s_in.roi_recoll&~s_in.roi_direct), s_in.rate(~s_in.roi_recoll&~s_in.roi_direct), '.', 'Color', my_G)
    xlim([s_in.t_birth(1) s_in.t_birth(end)])
end


s_out = s_in;
end