% Chirp Echo FT EPR-detected NMR
%
% 11/2017 Nino Wili: nino.wili@phys.chem.ethz.ch
% EPR @ ETH Zurich

% Comparison of Diagonal Slice of EPR-correlated chirp echo EPR-detected
% NMR with conventional EDNMR

clear, close all;


%% Diagonal Slice throuhg EPR correlated spectrum

%load data with function tailored to our AWG spectrometer
options.plot = 0; % 0 for do not plot by uwb_eval
options.evlen = 2000; % length of evaluation window in points
data=uwb_eval('20171122_0903_chirp_echo_fstep_EDNMR',options);

%extract needed data from struct
t=data.t_ax; %time axis
sig=data.dta_avg; %averaged echo transients

%No apodiazation, which would be unfair in a comparison with the
%conventional spectrum
window=chebwin(numel(t))*0+1;
sig_apo=sig.*window;

%Find echo maximum and extract phase for the reference spectrum 
[~,poi]=max(abs(sig_apo(:,1)));
phase0=angle(sig_apo(poi,1)); 
%loop over all spectra in data
for ii=1:size(sig,2)
   y=sig_apo(:,ii);%loop copy for convenience
   zero_filler=round(size(y)/2); %zero filler at start and end of echo
   y=[zeros(zero_filler); y; zeros(zero_filler)]; 
   y=circshift(y,-(poi+zero_filler)); %circular shift (=1st order phase correction)
   y=y.*exp(-1i*phase0); %oth order phase correction (usually irrelevant)
   %calculate spectra by FT
   spc(:,ii)=abs(fftshift(fft(y))); 
end
spc=spc/max(spc(:));


%calculate reference spectrum as the maximum of the spectra
%alternatively, one might take the median or actually record a reference
%spectrum w/o any HTA pulse
ref=max(spc')';

%calculate ednmr spectrum and frequency axis
center_frqs=data.dta_x{1}-data.det_frq;

for ii=1:size(spc,2)
    %take the difference spectrum
    diff_spc(:,ii)=(real(ref)-real(spc(:,ii)));
end

% projection
proj=diff_spc(round(size(spc,1))/2+1,:);
proj=proj-min(proj);
proj=proj/max(proj);

%% Conventional EDNMR spectrum
clear options, clear data

%load data with function tailored to our AWG spectrometer
options.plot = 0; % 0 for do not plot by uwb_eval
options.evlen = 2000; % length of evaluation window in points
data=uwb_eval('20171122_1015_chirp_echo_fstep_EDNMR',options);

%extract needed data from struct
f_HTA=data.dta_x{1};  %frequency axis of HTA pulse
spc=real(data.dta_ev); %averaged and integrated echo transients

%calculate ednmr spectrum and frequency axis
f_ednmr=f_HTA-data.exp.events{end}.det_frq;
ednmr=1-(spc/max(spc)); 
ednmr=ednmr-min(ednmr);
ednmr=ednmr/max(ednmr);


%% Plotting
h=figure(1);
clf
hold on
plot(center_frqs*1e3,proj,'k')
plot(center_frqs*1e3,ednmr+1,'k')
axis([-350 350 -0.04 2.04])


h=niwi_niceplot(h,[15 8]);
set(findall(h,'-property','LineWidth'),'LineWidth',1); % change linewidth
box off

%annotation
text(-330,0.3,'Diagonal slice')
text(-330,1.3,'Conventional')

xlabel('$( \nu_{\\text{HTA}}-\\nu  )$ / MHz')
ylabel('EDNMR Intensity / a.u.')
set(gca,'Ytick',[])
set(findall(h,'-property','Fontsize'),'Fontsize',12); % change fontsize

%Only needed for integration with Latex
% niwi_Plot2LaTex(h,'SI_MnCl2_comparison');