function out = calc_spec_cep(Espec, cep_list, traj_mat, Theta_Fit_list_dir, Theta_Fit_list_rec, input, options)
% calculate eSMM spectra with broadening and phases
% Input: energy axis (Espec [eV]), CEP axis (cep_list [pi]), trajectory simulation results (traj_mat), list of
% effective phases for all CEPs of direct and recoll electrons (Theta_Fit_list_dir, Theta_Fit_list_rec), relevant variables (input)
% Output: spectrum of all, direct and recollision electrons (complex amplitudes)

arguments
    Espec
    cep_list
    traj_mat
    Theta_Fit_list_dir
    Theta_Fit_list_rec
    input
    options.scale (1,1) double {mustBeNumeric} = 0.5
end
load_const;

% calculate energy widths
k = sqrt(input.W/(2*input.Up));
c = asinh(k) - k/sqrt(1+k^2);

sigma_d = options.scale * 1 * sqrt(input.Ephot / c);
sigma_r = options.scale * 0.425 * sqrt(input.Ephot / c); %0.425

fprintf("sigma_d^2 = %.2f eV, sigma_b^2 = %.2f eV \n", sigma_d^2, sigma_r^2);

N_cycle = numel(traj_mat(1).cycle.t_min);

% set matrices
out.all = nan(numel(Espec), numel(cep_list));
out.direct = nan(numel(Espec), numel(cep_list));
out.recoll = nan(numel(Espec), numel(cep_list));

% CEP loop 
for cep_i = 1:numel(cep_list)

    traj = traj_mat(cep_i);
    roi_x = ~isnan(traj.x) & traj.x >= 0;

    dir = zeros(size(Espec));
    recoll = zeros(size(Espec));

    % direct electrons
    if traj.use_dir
        for ci= 1:N_cycle
            roi_t = traj.t_birth >= traj.cycle.t_min(ci) & traj.t_birth <= traj.cycle.t_max(ci);
            roi_d = traj.n_recoll == 0 & roi_x & roi_t;
            rate_curr = traj.rate(roi_d);

            if ~isempty(rate_curr) & max(rate_curr) > 1e-4*max(traj.rate)
                Theta_fitted = Theta_Fit_list_dir(:, cep_i, ci).';

                for i = find(roi_d == 1)

                    E_curr = traj.Ekin(i);
                    rate_curr = traj.rate(i);

                    % gaussian in momentum space
                    psi_d = 1/sqrt(sqrt(pi)*sigma_d)*exp(-(sqrt(Espec) - sqrt(E_curr)).^2 / sigma_d^2);

                    % coherent summation including phase term 
                    dir = dir + sqrt(rate_curr) * psi_d .* exp(1i*Theta_fitted);

                end
            end
        end
    end

    % recollision electrons
    if traj.use_recoll
        for ci = 1:N_cycle
            roi_t = traj.t_birth >= traj.cycle.t_min(ci) & traj.t_birth <= traj.cycle.t_max(ci);
            roi_r = traj.n_recoll > 0 & roi_x & roi_t;
            rate_curr = traj.rate(roi_r);

            if ~isempty(rate_curr) & max(rate_curr) > 1e-4*max(traj.rate)

                Theta_fitted = Theta_Fit_list_rec(:, cep_i, ci).';

                for i = find(roi_r == 1)
                    E_curr = traj.Ekin(i);
                    rate_curr = traj.rate(i);

                    % gaussian in momentum space
                    psi_r = 1/sqrt(sqrt(pi)*sigma_r)*exp(-(sqrt(Espec) - sqrt(E_curr)).^2 / sigma_r^2);

                    % coherent summation including phase term
                    recoll = recoll + sqrt(rate_curr) * psi_r .* exp(1i*Theta_fitted);
                end

            end
        end
    end

    out.all(:, cep_i) = dir + recoll;
    out.direct(:, cep_i) = dir;
    out.recoll(:, cep_i) = recoll;
    out.sigma_d = sigma_d;
    out.sigma_r = sigma_r;
    out.psi = @(E, E_curr, sigma) 1/sqrt(sqrt(pi)*sigma)*exp(-(sqrt(E) - sqrt(E_curr)).^2 / sigma^2);

end
end
