% 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 deuterated 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 = 6000; % length of evaluation window in points
data=uwb_eval('20170828_1711_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),500);
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.346 34.654];

%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.007,0.045,40);

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

p2H=[34.37 0];
plot([34.375 p2H(1)],[15.5 p2H(2)],'k'); 
plot([34.38 p2H(1)],[7 p2H(2)],'k');
plot([34.395 p2H(1)],[-7.333 p2H(2)],'k');
plot([34.405 p2H(1)],[-16.7 p2H(2)],'k');

text(p2H(1)-0.005,p2H(2),'$^{2}$H','Horizontalalignment','right')

text(34.59,54,'$^{1}$H','Horizontalalignment','left')
text(34.45,-62,'$^{1}$H','Horizontalalignment','right')

text(34.435,-40,'$^{14}$N$_{\text{SQ}}$','Horizontalalignment','right')
text(34.435,+40,'$^{14}$N$_{\text{SQ}}$','Horizontalalignment','right')

text(34.54,+86,'$^{14}$N$_{\text{DQ}}$','Horizontalalignment','left')

starsx=[34.6];
starsy=[-75];
plot(starsx,starsy,'*k','Markersize',15)

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

