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

% HYSCORE Type Correlations in EDNMR on Hexaquamanganese

clear, close all;

%% The 1D Spectrum as reference

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

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

%apodization
window=chebwin(numel(t),200);
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 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 ednmr 
ednmr=spc(:,1)-spc(:,2);
ednmr=ednmr/max(ednmr);

% plot 1d as reference on top
h=figure(1);
clf
hold on
subplot(4,1,[1 2])
hold on
plot(freq*1e3,ednmr,'k')
set(gca,'Xtick',[],'Ytick',[])
axis([180*[-1 1] -0.04 1.54])



%% The Hyscore type spectra Number 1)

%load data with function tailored to our AWG spectrometer
data=uwb_eval('20171023_1620_chirp_echo_inv_n_HTA_rescued',options);

%extract needed data from struct
sig=data.dta_avg; %averaged echo transients

%apodization
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 hys_ednmr of the first trace
hys_ednmr1=spc(:,1)-spc(:,2);


%% The Hyscore type spectra Number 2)

%load data with function tailored to our AWG spectrometer
data=uwb_eval('20171023_1624_chirp_echo_inv_n_HTA',options);

%extract needed data from struct
sig=data.dta_avg; %averaged echo transients

%apodization
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 hys_ednmr 
hys_ednmr2=spc(:,1)-spc(:,2);


%% Plotting of the HYS-type spectra 1) and 2)
xlimtop=xlim;

subplot(4,1,[ 3 4])
hold on

plot(freq*1e3,hys_ednmr1,'k') 
plot(freq*1e3,hys_ednmr2+0.1,'k') 

ylim1=[-0.1 0.2];
axis([xlimtop -0.1 0.2])
set(gca,'YTick',[])
xlabel('$( \nu - \\nu_{\\text{HTA}} )$ / MHz')


h=niwi_niceplot(h,[15 12]*1.5);

%Annotations
%assignment
peaklist=[ -157.2 -113.2  114.3 156.7];

subplot(4,1,[1 2])

%asterisks
starsx=[-167.2 -104.5 106.5 165.7];
starsy=[0.25 0.28 0.25  0.2];
plot(starsx,starsy,'*k')
%protons
protx=[-52 52];
proty=[0.65 0.65];
text(protx,proty,'\tiny{H}','horizontalalignment','center')
%manganese
Mnx=[-135 136];
Mny=[0.45 0.45];
text(Mnx,Mny,'\tiny{Mn}','horizontalalignment','center')

subplot(4,1,[1 2])
plot([peaklist(1)*[1 1] peaklist(3)*[1 1]],[0.6 1.1 1.1 0.5],'k')
text(peaklist(1),1.2,'\tiny{1)}','horizontalalignment','center')
plot([peaklist(2)*[1 1] peaklist(4)*[1 1]],[0.9 1.2 1.2 0.8]+0.1,'k')
text(peaklist(2),1.4,'\tiny{2)}','horizontalalignment','center')

subplot(4,1,[3 4])
text(-137,0.03,'\tiny{1)}','horizontalalignment','center')
text(-137,0.13,'\tiny{2)}','horizontalalignment','center')

%only needed for latex integration
% niwi_Plot2LaTex(h,'hyscore_type_Mn');
