clear
clc

filename1 = load('filename.mat');

for atom_index = 1:length(filename1)
    clearvars -except atom_index filename1
    load(filename1{atom_index});

    transitionRate = transitionRate*1e2;  % fs^-1 -> 1e2 fs^-2
    finalTime = 1e3;
    steplength = 1e-3;
    opts.MaxStep = steplength;
    time = (0:steplength:finalTime).';
    
    tr1 = bsxfun(@plus,transitionRate(:,1:end-1),-transitionRate(:,end));
    tr2 = transitionRate(:,end);
    hfunc = @(x,y)(tr1*y(1:end-1,:)+tr2);
    
    corrs = cell(length(dipoleMoment),1);
    for kk = 1:length(dipoleMoment)
        corr_t0 = zeros(length(eigenEnergy),1);
        corr_t0(kk) = 1;
        corrs{kk} = ode45(hfunc,[0,finalTime],corr_t0,opts);
    end
    
    %%
    corr_dipole = 0;
    for kk = 1:length(dipoleMoment)
        corr_rho = spline(corrs{kk}.x,corrs{kk}.y,time);
        corr_dipole = corr_dipole +dipoleMoment(kk)*sum(dipoleMoment.*corr_rho,1)*corrs{1}.y(kk,end);
    end
    corr_dipole = corr_dipole.';
    
    %%
    freq = (0:1/time(end):1/steplength).';
    freq = freq*1e13;
    
    sp_0 = trapz(time,corr_dipole-corr_dipole(end));
    sp = ifft(corr_dipole-corr_dipole(end));
    sp_after = sp/real(sp(1))*sp_0;
    
    %%
    coverage = 1e18;  % m^-2
    epsilon_0 = 8.854187817e-12;  % C/V/m
    distance = 40e-6;  % m
    sp_after = sp_after*1e-13*(3.33564e-30)^2;  % Debye^2*fs > (C*m)^2*s
    sp_after = sp_after*3*pi*coverage/(4*pi*epsilon_0)^2/distance^4;  % V^2*m^-2*s
    
    freq_tm = [log10(freq(2)),log10(freq(end))];
    freq_new = power(10,linspace(freq_tm(1),freq_tm(2),1e4)).';
    sp = spline(freq,sp_after,freq_new);
    sp = [freq_new, sp];
    save(['sp_',filename1{atom_index}],'sp');
    
    time_new = linspace(time(1),time(end),1e4).';
    corr = spline(time,corr_dipole,time_new);
    corr = [time_new, corr];
    save(['corr_',filename1{atom_index}],'corr');

end
clear