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;
 
figure1=figure('units','centimeters','position',[5,4,18,16],'color','white','DefaultTextInterpreter','LaTex','DefaultLineLineWidth',1.5);
splitIndex = ceil(length(Temp)/2);  % Index to split the plots



% ----------- LEFT PLOT (First half of temperatures) -----------
subplot(1,2,1);
hold on; box on;
xlim([-2.39,-1.47]); ylim([0,10]);
set(gca,'LineWidth',1.5,'FontSize',14,'TickLabelInterpreter','LaTeX');
grid on; set(gca,'GridLineStyle','--','GridAlpha',0.4); 
xlabel('Energy transfer (meV)');
ylabel('Normalised intensity (a.u.)');
xtickformat('%.1f');
set(gca,'ytick',[]);

for kk = 1:splitIndex
   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); 
     
    if kk==1
       yshift=0.2;
    else
       yshift=(kk-1)*0.73+0.2;
    end
    plot(-xfit(xfit<2300)/1e3,yfit(xfit<2300)+yshift,'-','Color',cvec2(1,:));
    plot(-xdata(xdata<2300)/1e3,ydata(xdata<2300)+yshift,'Marker','o','MarkerSize',5,'LineStyle','None','Color',cvec(1,:),'MarkerFaceColor',cvec2(1,:));
    if kk==length(Temp);
        text(-2.37,yshift-0.0,['T = ',num2str(Temp(kk),3),' K'],'BackgroundColor',[1 1 1],'VerticalAlignment','bottom','HorizontalAlignment','left','FontSize',14); 
    else
        text(-2.37,yshift-0.08,[num2str(Temp(kk),3),' K'],'BackgroundColor',[1 1 1],'VerticalAlignment','bottom','HorizontalAlignment','left','FontSize',14); 
    end
            
    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



% ----------- RIGHT PLOT (Second half of temperatures) -----------
subplot(1,2,2);
hold on; box on;
xlim([-2.39,-1.47]); ylim([0,10]);
set(gca,'LineWidth',1.5,'FontSize',14,'TickLabelInterpreter','LaTeX');
grid on; set(gca,'GridLineStyle','--','GridAlpha',0.4); 
xlabel('Energy transfer (meV)');
ylabel('Normalised intensity (a.u.)');
xtickformat('%.1f');
set(gca,'ytick',[]);

for kk = splitIndex+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); 
     
    if kk-splitIndex==1
       yshift=0.2;
    else
       yshift=(kk-1-splitIndex)*0.73+0.2;
    end
    plot(-xfit(xfit<2300)/1e3,yfit(xfit<2300)+yshift,'-','Color',cvec2(1,:));
    plot(-xdata(xdata<2300)/1e3,ydata(xdata<2300)+yshift,'Marker','o','MarkerSize',5,'LineStyle','None','Color',cvec(1,:),'MarkerFaceColor',cvec2(1,:));
    if kk==length(Temp);
        text(-2.37,yshift-0.0,['T = ',num2str(Temp(kk),3),' K'],'BackgroundColor',[1 1 1],'VerticalAlignment','bottom','HorizontalAlignment','left','FontSize',14); 
    else
        text(-2.37,yshift-0.08,[num2str(Temp(kk),3),' K'],'BackgroundColor',[1 1 1],'VerticalAlignment','bottom','HorizontalAlignment','left','FontSize',14); 
    end
            
    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