clear;
clc;
load_au;
load_cm;

list = 0.0:0.01:2;
list = round(list,2,TieBreaker="minusinf");

% Read TDSE results calculated with C++
for j = 1:numel(list)
    l = list(j);
    disp(l);
    path = append('800nm_1e12_SC_cep\phase_', num2str(l, '%.2f'), '\');
    out = read_tdse_run(path);

    spec(j,:) = out.spec;
    
end
specr = real(spec);
roi = specr < 0;
specr(roi) = 0;

map = cm.vlasov;
Up = out.Up;

%% Plot spectra 

figure
imagesc(out.E_list/Up, list, log10(specr))
colormap(map)
clim([-11 -1])
xlim([0 15])
xlabel('E [U_p]')
ylabel('phase [pi]')
set(gca, 'FontSize', 14)
