function out = calc_spec_with_eff_phase(Espec, loop_list, traj_mat, options)

arguments
    Espec
    loop_list
    traj_mat
    options.weight_scale (1,1) double {mustBeNumeric} = 1
    options.scale (1,1) double {mustBeNumeric} = 1
end

load_const;
load_color;

N_spec = 200;
W_scale = 1; %0.75;
V_scale = 1;
T_scale = 1;

N_cycle = numel(traj_mat(1).cycle.t_min);
out.theta_eff_dir = nan(N_spec, numel(loop_list), 4*N_cycle);
out.Theta_Fit_list_dir = nan(numel(Espec), numel(loop_list), 4*N_cycle);
out.offset_dir = nan(numel(loop_list), 4*N_cycle);
out.Edir = nan(numel(loop_list), 10000, 4*N_cycle);
out.timedir = nan(numel(loop_list), 10000, 4*N_cycle);
out.ratedir = nan(numel(loop_list), 10000, 4*N_cycle);
out.Thetadir = nan(numel(loop_list), 10000, 4*N_cycle);
out.spec_dir_parts = nan(numel(Espec), numel(loop_list), 4*N_cycle);

out.theta_eff_rec = nan(N_spec, numel(loop_list), 4*N_cycle);
out.Theta_Fit_list_rec = nan(numel(Espec), numel(loop_list), 4*N_cycle);
out.offset_rec = nan(numel(loop_list), 4*N_cycle);
out.Erec = nan(numel(loop_list), 10000, 4*N_cycle);
out.timerec = nan(numel(loop_list), 10000, 4*N_cycle);
out.raterec = nan(numel(loop_list), 10000, 4*N_cycle);
out.Thetarec = nan(numel(loop_list), 10000, 4*N_cycle);
out.spec_rec_parts = nan(numel(Espec), numel(loop_list), 4*N_cycle);

weight_co_scale = 25e-4;

% for extracting local max and min
min_sep_rec = 150; % 550; %
min_sep_dir = 150;
min_prom = 1e-4; % 1e-2;%
min_t_abs = 0.1; % 0.1

% energy broadening
Up = traj_mat(1).Up; W = traj_mat(1).W; Ephot = traj_mat(1).Ephot;
k = sqrt(W/(2*Up));
c = asinh(k) - k/sqrt(1+k^2);

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

out.all = nan(numel(Espec), numel(loop_list));
out.direct = nan(numel(Espec), numel(loop_list));
out.recoll = nan(numel(Espec), numel(loop_list));

% traj_mat(1).make_plots = 1;
if traj_mat(1).make_plots
    figure(501); subplot(2,3,1); hold off; subplot(2,3,2); hold off; subplot(2,3,3); hold off; subplot(2,3,4); hold off; subplot(2,3,4); hold off; subplot(2,3,5); hold off;
end

for loop_i = 1:numel(loop_list) %1:numel(cep_list)
    disp(loop_list(loop_i))

    input = traj_mat(loop_i);
    % input.make_plots = 1;
    N_cycle = numel(input.cycle.t_min);

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

    if isempty(input.Ebirth)
        disp('nothing to see here')
        continue
    end

    roi_x = ~isnan(input.x(1,:)) & input.x(1,:) >= 0;

    max_rate_tot = max(input.rate);

    idx_dir = 1;
    if input.use_dir
        min_sep = min_sep_dir;
        for ci = 1:N_cycle
            roi_t = input.t_birth >= input.cycle.t_min(ci) & input.t_birth <= input.cycle.t_max(ci);
            roi_d = input.n_recoll == 0 & roi_x & roi_t;

            E_curr_cycle = input.Ekin(roi_d);
            t_curr_cycle = input.t_birth(roi_d);

            ind_max = islocalmax(E_curr_cycle, 'MinSeparation', min_sep, 'MinProminence', min_prom);
            ind_min = islocalmin(E_curr_cycle, 'MinSeparation', min_sep, 'MinProminence', min_prom);
            ind_tmin = find(abs(diff([-1000, t_curr_cycle])) > 0.15*pi/input.w);
            ind_tmax = find(abs(diff([t_curr_cycle, 1000])) > 0.15*pi/input.w);

            ind_start = 1;

            list = sort([find(ind_max|ind_min == 1), numel(t_curr_cycle), ind_tmin, ind_tmax]);
            

            for idx_list = 1:numel(list)
                ind = list(idx_list);
                if (idx_list < numel(list) && abs(t_curr_cycle(ind) - t_curr_cycle(list(idx_list+1))) < min_t_abs) || numel(t_curr_cycle) < 1
                    continue;
                end
                
                roi_dc = roi_d & input.t_birth > t_curr_cycle(ind_start) & input.t_birth < t_curr_cycle(ind);
                ind_start = ind;


                E_curr = input.Ekin(roi_dc);
                t_curr = input.t_birth(roi_dc);
                rate_curr = input.rate(roi_dc);
                S_curr = T_scale*input.S.T(roi_dc) + V_scale*input.S.V(roi_dc) - W_scale*input.S.W(roi_dc);

                if numel(t_curr)>2 & max(rate_curr) > 1e-4*max_rate_tot
                    weight = rate_curr / max(rate_curr);
                    weight = weight.^options.weight_scale;

                    cut_off_dir = max(E_curr);

                    Espec_hist = linspace(0, 1.1*cut_off_dir, N_spec);

                    A_dir = histw_phase(E_curr, rate_curr, S_curr, Espec_hist);
                    theta_dir = unwrap(angle(A_dir));

                    % figure
                    % plot(Espec_hist/input.Up, log10(abs(A_dir).^2))
                    % 
                    % pause()

                    roi_weight = weight > weight_co_scale;
                    Etestmax = max(E_curr(roi_weight));
                    Etestmin = min(E_curr(roi_weight));

                    roi_smm_dir = Espec_hist < Etestmax & Espec_hist > Etestmin;

                    if numel(Espec_hist(roi_smm_dir)) < 2 | input.use_phases == 0
                        Theta_fitted = zeros(size(Espec));
                    else
                        Theta_fitted = interp1(Espec_hist(roi_smm_dir), theta_dir(roi_smm_dir), Espec, "linear", "extrap");
                    end

                    [~, Eminidx_curr_all] = findpeaks(-abs(E_curr-Etestmax));
                    weight_help = zeros(size(weight)); weight_help(Eminidx_curr_all) = weight(Eminidx_curr_all);
                    [~, Eminidx_curr] = max(weight_help);
                    [~, Eminidx] = min(abs(Espec-Etestmax));

                    offset = round(S_curr(Eminidx_curr)/(const.hbar*pi) - Theta_fitted(Eminidx)/pi, 0);

                    dir_part = zeros(size(dir));
                    for i = find(roi_dc == 1)

                        E_i = input.Ekin(i);
                        rate_i = input.rate(i);

                        psi_d = 1/sqrt(sqrt(pi)*sigma_d)*exp(-(sqrt(Espec) - sqrt(E_i)).^2 / sigma_d^2);
                        dir_part = dir_part + sqrt(rate_i) * psi_d .* exp(1i*Theta_fitted);

                    end
                    dir = dir + dir_part;

                    if input.make_plots
                        subplot(2,3,1)
                        scatter(t_curr, E_curr/input.Up, 0.1+20*weight, my_B, 'filled')
                        hold on
                        ylim([0 2*cut_off_dir/input.Up])

                        subplot(2,3,2)
                        scatter(E_curr/input.Up, S_curr/const.hbar/pi, 0.1+20*weight, my_B, 'filled')
                        hold on
                        plot(Espec/input.Up, Theta_fitted/pi + offset, '--', 'LineWidth', 1.5, 'Color', my_G)
                        plot(Espec_hist/input.Up, theta_dir/pi + offset, 'LineWidth', 1, 'Color', 'c')
                        xlim([0 2*cut_off_dir/input.Up])
                        ylim([0 100])
                        hold off

                        subplot(2,3,3)
                        plot(Espec/Up, log10(abs(dir).^2), 'Color', my_B, 'LineWidth', 1.5)
                        xlim([0 2*cut_off_dir/input.Up])
                        % hold off
                        drawnow;
                        % pause()
                    end

                    out.theta_eff_dir(:, loop_i, idx_dir) = theta_dir; 
                    out.Theta_Fit_list_dir(:, loop_i, idx_dir) = Theta_fitted;
                    out.offset_dir(loop_i, idx_dir) = offset;
                    out.Edir(loop_i, 1:numel(E_curr), idx_dir) = E_curr;
                    out.timedir(loop_i, 1:numel(t_curr), idx_dir) = t_curr;
                    out.ratedir(loop_i, 1:numel(rate_curr), idx_dir) = rate_curr;
                    out.Thetadir(loop_i, 1:numel(S_curr), idx_dir) = S_curr;
                    out.Espechist(loop_i,:) = Espec_hist;
                    out.spec_dir_parts(:,loop_i,idx_dir) = dir_part;

                    idx_dir = idx_dir + 1;


                else
                    disp('nothing to fit here')
                end
            end
        end
    end


    idx_rec = 1;
    if input.use_recoll
        min_sep = min_sep_rec;
        for ci = 1:N_cycle
            roi_t = input.t_birth >= input.cycle.t_min(ci) & input.t_birth <= input.cycle.t_max(ci);
            roi_r = input.n_recoll == 1 & roi_x & roi_t;

            E_curr_cycle = input.Ekin(roi_r);
            t_curr_cycle = input.t_birth(roi_r);

            ind_max = islocalmax(E_curr_cycle, 'MinSeparation', min_sep, 'MinProminence', min_prom);
            ind_min = islocalmin(E_curr_cycle, 'MinSeparation', min_sep, 'MinProminence', min_prom);
            ind_tmin = find(abs(diff([-1000, t_curr_cycle])) > 0.15*pi/input.w);
            ind_tmax = find(abs(diff([t_curr_cycle, 1000])) > 0.15*pi/input.w);

            ind_start = 1;

            list = sort([find(ind_max|ind_min == 1), numel(t_curr_cycle), ind_tmin, ind_tmax]);
            

            for idx_list = 1:numel(list)
                ind = list(idx_list);
                if (idx_list < numel(list) && abs(t_curr_cycle(ind) - t_curr_cycle(list(idx_list+1))) < min_t_abs) || numel(t_curr_cycle) < 1
                    continue;
                end
                
                roi_rc = roi_r & input.t_birth > t_curr_cycle(ind_start) & input.t_birth < t_curr_cycle(ind);
                ind_start = ind;


                E_curr = input.Ekin(roi_rc);
                t_curr = input.t_birth(roi_rc);
                rate_curr = input.rate(roi_rc);
                S_curr = T_scale*input.S.T(roi_rc) + V_scale*input.S.V(roi_rc) - W_scale*input.S.W(roi_rc);

                if numel(t_curr)>2 & max(rate_curr) > 1e-4*max_rate_tot
                    weight = rate_curr / max(rate_curr);
                    weight = weight.^options.weight_scale;

                    cut_off_rec = max(E_curr);

                    Espec_hist = linspace(0, 1.1*cut_off_rec, N_spec);

                    A_rec = histw_phase(E_curr, rate_curr, S_curr, Espec_hist);
                    theta_rec = unwrap(angle(A_rec));

                    % figure
                    % plot(Espec_hist/input.Up, log10(abs(A_rec).^2))
                    % 
                    % pause()

                    roi_weight = weight > weight_co_scale;
                    Etestmax = max(E_curr(roi_weight));
                    Etestmin = min(E_curr(roi_weight));

                    roi_smm_rec = Espec_hist < Etestmax & Espec_hist > Etestmin;

                    if numel(Espec_hist(roi_smm_rec)) < 2 | input.use_phases == 0
                        Theta_fitted = zeros(size(Espec));
                    else
                        Theta_fitted = interp1(Espec_hist(roi_smm_rec), theta_rec(roi_smm_rec), Espec, "linear", "extrap");
                    end

                    [~, Eminidx_curr_all] = findpeaks(-abs(E_curr-Etestmax));
                    weight_help = zeros(size(weight)); weight_help(Eminidx_curr_all) = weight(Eminidx_curr_all);
                    [~, Eminidx_curr] = max(weight_help);
                    [~, Eminidx] = min(abs(Espec-Etestmax));

                    offset = round(S_curr(Eminidx_curr)/(const.hbar*pi) - Theta_fitted(Eminidx)/pi, 0);

                    rec_part = zeros(size(rec));
                    for i = find(roi_rc == 1)

                        E_i = input.Ekin(i);
                        rate_i = input.rate(i);

                        psi_r = 1/sqrt(sqrt(pi)*sigma_r)*exp(-(sqrt(Espec) - sqrt(E_i)).^2 / sigma_r^2);

                        rec_part = rec_part + sqrt(rate_i) * psi_r .* exp(1i*Theta_fitted);
                        
                    end
                    rec = rec + rec_part;

                    if input.make_plots
                        subplot(2,3,4)
                        scatter(t_curr, E_curr/input.Up, 0.1+20*weight, my_R, 'filled')
                        hold on
                        ylim([0 2*cut_off_rec/input.Up])

                        subplot(2,3,5)
                        scatter(E_curr/input.Up, S_curr/const.hbar/pi, 0.1+20*weight, my_R, 'filled')
                        % scatter(E_curr/input.Up, S_curr/const.hbar/pi, 10, my_R, 'filled')
                        hold on
                        plot(Espec/input.Up, Theta_fitted/pi + offset, '--', 'LineWidth', 1.5, 'Color', my_G)
                        plot(Espec_hist/input.Up, theta_rec/pi + offset, 'LineWidth', 1, 'Color', 'm')
                        xlim([0 2*cut_off_rec/input.Up])
                        ylim([0 200])
                        hold off

                        subplot(2,3,6)
                        plot(Espec/Up, log10(abs(rec).^2), 'Color', my_R, 'LineWidth', 1.5)
                        xlim([0 2*cut_off_rec/input.Up])
                        drawnow;
                        % pause()
                    end

                    out.theta_eff_rec(:, loop_i, idx_rec) = theta_rec; 
                    out.Theta_Fit_list_rec(:, loop_i, idx_rec) = Theta_fitted;
                    out.offset_rec(loop_i, idx_rec) = offset;
                    out.Erec(loop_i, 1:numel(E_curr), idx_rec) = E_curr;
                    out.timerec(loop_i, 1:numel(t_curr), idx_rec) = t_curr;
                    out.raterec(loop_i, 1:numel(rate_curr), idx_rec) = rate_curr;
                    out.Thetarec(loop_i, 1:numel(S_curr), idx_rec) = S_curr;
                    out.spec_rec_parts(:, loop_i, idx_rec) = rec_part;

                    idx_rec = idx_rec + 1;

                else
                    disp('nothing to fit here')
                end
            end
        end
    end

    out.all(:, loop_i) = dir + rec;
    out.direct(:, loop_i) = dir;
    out.recoll(:, loop_i) = rec;
    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
