clear all; close all;

addpath 'F:\HeSurfer\home\at677\matlabwork'

%=================== File list to read in ===================
%addpath /home/shared/SpinEcho/Spinecho2Icon/data/Ag100/
%filelist=[012106:012113,012093:012105,012116:012117]; 

cvec =[ 0,0,0.85; 0,0.5,0; 0.9,0.5,0; 0.9,0,0; 0.2 0.2 0.2];
cvec2=[ 0,0.7,1;  0,0.8,0; 1,0.7,0;   1.0,0.4,0; 0.2 0.2 0.2];
 
load('LWdata.mat');
Temp=data.Temp;
 

%=================== Load .mat file and extract information ===================
for kk=1:length(Temp);
     
    dE=data.xdata(:,kk);
    spectrum=data.ydata(:,kk);
     
    xdata=abs(dE(-2.4<dE & dE<-1.4)*1000);
    ydata=spectrum(-2.4<dE & dE<-1.4);
     
    ydata=abs(ydata)./max(abs(ydata)); %normalised to 1
    [cmax,ind]=max(ydata); maxpos=xdata(ind);
     
    xmeas(kk,:)=xdata;
    ymeas(kk,:)=ydata;
         
     
    %% Fit: 'untitled fit 1'.
    [xData, yData] = prepareCurveData( xdata, ydata );
 
    ft = fittype( 'a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2)+d', 'independent', 'x', 'dependent', 'y' );
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.Lower = [0 0 0 0 0 0 -1];
    opts.StartPoint = [1 1 maxpos maxpos+200 50 100 0.2];
 
    % Fit model to data.
    [fitresult, gof] = fit( xData, yData, ft, opts );
     
    xfit=linspace(min(xdata),max(xdata),201);
    yfit = feval(fitresult,xfit); 
    
            
    Phen(kk) = fitresult.b1;
    Phsig(kk) = fitresult.c1;
    ci = confint(fitresult,0.66);
     
    dPhen(kk) = 3*abs(ci(1,3) - ci(2,3))/2;
    dPhsig(kk) =abs(ci(1,5)-ci(2,5))/1.5;
     
         
end

%plot of phonon energy vs. Temperature
figure1=figure('units','centimeters','position',[2,4,16,12],'Color','white','DefaultTextInterpreter','LaTex','DefaultLineLineWidth',1.5);
set(gca,'LineWidth',1.5,'FontSize',14,'Layer','top','TickLabelInterpreter','LaTex'); box(gca,'on'); hold(gca,'all');

plot(Temp,Phen*1e-3,'Marker','o','MarkerSize',6,'LineStyle','None','Color',cvec(1,:),'MarkerFaceColor',cvec2(1,:));
errbar(Temp,Phen*1e-3,dPhen*1e-3,'Color',cvec(1,:),'LineWidth',1.5);
ylim([1680,1860]*1e-3); xlim([0,800]);
xlabel('Surface temperature $T$ (K)'); ylabel('Phonon energy (meV)');
 
fo_ = fitoptions('method','NonlinearLeastSquares');
startingval=[2000,0.01];
set(fo_,'Startpoint',startingval);
ft_ = fittype('Ez*(1-xi*(2./(exp(Ez./(86.173324*x))-1 )+1 ))','dependent',{'y'},'independent',{'x'},'coefficients',{'Ez','xi'});
fitresult = fit(Temp',Phen',ft_,fo_);
 
E0 = fitresult.Ez; Xi = fitresult.xi;
ci = confint(fitresult,0.66);
dE0 = abs(ci(1,1) - ci(2,1))/2;
delXi=abs(ci(1,2)-ci(2,2))/2;
 
xfit=linspace(0,1.1*max(Temp),500);
yfit = feval(fitresult,xfit); 
plot(xfit,yfit*1e-3,'-.','Color',cvec(4,:));
ydiffEn=feval(fitresult,Temp)-Phen';
ytickformat('%.2f')
 

%plot of linewidth versus temperature
figure2=figure('units','centimeters','position',[18,4,16,12],'Color','white','DefaultTextInterpreter','LaTex','DefaultLineLineWidth',1.5);
set(gca,'LineWidth',1.5,'FontSize',14,'Layer','top','TickLabelInterpreter','LaTex'); box(gca,'on'); hold(gca,'all');
 
fwhmconv=2*sqrt(log(2));
Phsig=fwhmconv*Phsig;
Gamint=sqrt(Phsig.^2-60^2); %subtract angular resolution

for kk=1:length(Phsig); %solve scan curve equation
   
   syms x
   sol=vpasolve(Phsig(kk)==0.5346*x+sqrt( 0.2166*x^2 + 60^2 ),x);
   Gamnat(kk)=double(sol);

end
Gamma=Gamint./2;
plot(Temp,Gamma,'Marker','o','MarkerSize',6,'LineStyle','None','Color',cvec(1,:),'MarkerFaceColor',cvec2(1,:));
errbar(Temp,Gamma,dPhsig,'Color',cvec(1,:),'LineWidth',1.5);
 
ylim([40,100]./2);
xlabel('Surface temperature $T$ (K)'); ylabel('Linewidth $\Gamma_{RW}$ of the RW mode ($\mu$eV)');
 
fo_ = fitoptions('method','NonlinearLeastSquares');
startingval=[0.04,40];
set(fo_,'Startpoint',startingval);
ft_ = fittype('k*x+d','dependent',{'y'},'independent',{'x'},'coefficients',{'k','d'});
fitresult = fit(Temp',Gamma',ft_,fo_);
ydiffLw=feval(fitresult,Temp)-Gamma';
 
k = fitresult.k; d = fitresult.d;
ci = confint(fitresult,0.66);
dk = abs(ci(1,1) - ci(2,1))/2;
deld=abs(ci(1,2)-ci(2,2))/2;
 
xfit=linspace(0,1.03*max(Temp),100);
yfit = feval(fitresult,xfit); 
plot(xfit,yfit,'-.','Color',cvec(4,:));
xlim([0,800]); ylim([20,50]);

tau=6.58211928E-16/(d*1e-6)*1e12; %in ps,
dtau=(6.58211928E-16/((d-deld)*1e-6)*1e12)-tau; %in ps,

%plot of deviation from linear model
figure3=figure('units','centimeters','position',[18,4,16,12],'Color','white','DefaultTextInterpreter','LaTex','DefaultLineLineWidth',1.5);
set(gca,'LineWidth',1.5,'FontSize',14,'Layer','top','TickLabelInterpreter','LaTex'); box(gca,'on'); hold(gca,'all');

p(1)=plot(Temp,ydiffEn,'Marker','o','MarkerSize',6,'LineStyle','None','Color',cvec(1,:),'MarkerFaceColor',cvec2(1,:));
errbar(Temp,ydiffEn,dPhen,'Color',cvec(1,:),'LineWidth',1.5);

p(2)=plot(Temp,ydiffLw,'Marker','o','MarkerSize',6,'LineStyle','None','Color',cvec(2,:),'MarkerFaceColor',cvec2(2,:));
errbar(Temp,ydiffLw,dPhsig,'Color',cvec(2,:),'LineWidth',1.5);

ylim([-10,12]); 
line([0,800],[0,0],'Color',[0.4 0.4 0.4]);
l1=legend(p,'Phonon energy deviation','Line width deviation'); set(l1,'Interpreter','LaTeX');
xlabel('Surface temperature $T_S$ (K)'); ylabel({'Deviation from','linear behaviour ($\mu$eV)'});