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