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

% EPR-correlated chirp echo FT EPR spectrum of protonated nitroxide

clear, close all;

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

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

%apodization
window=chebwin(numel(t));
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)=real(fftshift(fft(y))); 
end
spc=spc/max(spc(:));


%calculate frequency axis
dt=t(2)-t(1); %time step
nyq=1/(2*dt); %nyquist frequency from time step
freq=linspace(-nyq,nyq,size(spc,1));

%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 and shift such that the HTA corresponds to 0 frequency
%of the hyperfine spectrum
center_frqs=data.dta_x{1}-data.det_frq;
df=freq(2)-freq(1);
f_range=100e-3;  %range that is cut out from the data
Nf=ceil(f_range/df);
f_nmr=(-Nf:Nf)*df;

for ii=1:size(spc,2)
    %take the difference spectrum
    diff_spc(:,ii)=(real(ref)-real(spc(:,ii)));
    
    %cut out range symmetrically around the HTA frequency
    [~,poi]=min(abs(freq-center_frqs(ii)));
    ednmr(:,ii) = diff_spc(poi-Nf:poi+Nf,ii);
    ednmr(:,ii)=ednmr(:,ii)-min(ednmr(:,ii));
end



%Plot
h=figure(14);
clf
colormap(jet)

%gaussian "blur" filter, just makes the contour plot look nicer
filt_width=5;
x_filt=-filt_width:filt_width;
y_filt=-filt_width:filt_width;
[x_filt,y_filt]=meshgrid(x_filt,y_filt);
H=normpdf(x_filt,0,1).*normpdf(y_filt,0,1);
H=H/sum(H(:));
ednmr_filtered=filter2(H,ednmr);

%set axis of 2D plot
LO=data.exp.LO;
nu_obs=data.exp.events{end}.det_frq;
x_2=center_frqs+data.det_frq+LO;
y_2=linspace(-f_range,f_range,size(ednmr,1))*1e3;
[X_2,Y_2]=meshgrid(x_2,y_2);
ax_x=[34.52 34.83];

%plot reference spectrum on top
subplot(4,1,1)
plot(freq+LO+nu_obs,ref*0.95,'k','linewidth',2)
axis([ax_x ylim-[0.04 0]])
set(gca,'YTick',[],'XTick',[])

%contour plot (normalized to max of 2d spectrum)
ednmr_filtered=ednmr_filtered/max(ednmr_filtered(:));
levels=linspace(0.01,0.06,20);
subplot(4,1,[2 3 4])
contour(X_2,Y_2,(ednmr_filtered),'levellist',(levels))
axis([ax_x ylim])


%set axis labels and make look nice
xlabel('$\nu_{\\text{HTA}} $ / GHZ')
ylabel('$( \nu - \\nu_{\\text{HTA}}  )$ / MHz')
h=niwi_niceplot(h,[15 12]*1.5);


%annotations
hold on
p13C=[34.65 16];
plot([34.675 p13C(1)],[13 p13C(2)],'k'); plot([34.682 p13C(1)],[-11 p13C(2)],'k');
text(p13C(1)-0.005,p13C(2),'$^{13}$C','Horizontalalignment','right')

text(34.76,52,'$^{1}$H','Horizontalalignment','left')
text(34.605,-62,'$^{1}$H','Horizontalalignment','right')

pcombi1=[34.77 85];
plot([34.72 pcombi1(1)],[73 pcombi1(2)],'k'); 
text(pcombi1(1)-0.005,pcombi1(2),'$^{14}$N+$^{1}$H','Horizontalalignment','left','verticalalignment','bottom')

pcombi2=[34.73 -80];
plot([34.75 pcombi2(1)],[-73 pcombi2(2)],'k'); 
text(pcombi2(1)-0.005,pcombi2(2),'$^{14}$N+$^{1}$H','Horizontalalignment','right','verticalalignment','top')

plot(34.58*[1 1 ],95/2*[ -1 1],'k'); 
plot([34.58 34.59],95/2*[ 1 1],'k');
plot([34.58 34.59],-95/2*[ 1 1],'k')
text(34.58,25,'$A_{zz}(^{14}$N) ','Horizontalalignment','right')

plot(34.82*[1 1 ],11.5/2*[ -1 1],'k'); 
plot([34.82 34.81],11.5/2*[ 1 1],'k');
plot([34.82 34.81],-11.5/2*[ 1 1],'k')
text(34.81,20,'$A_{xx,yy}$','Horizontalalignment','center')

%only needed for integration with Latex
h.Renderer='painters';
% niwi_Plot2LaTex(h,'2d_protonated_nitroxide');
