clear all; close all;


%=================== File list to read in ===================

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];

% 300 K and 350 K measurement at I0/50
filelist = [11939, 11940, 11941, 11943, 11944, 11945, 11946, ...
            11953, 11954, 11955, 11999, 12000, 12001, 12002, 12003, 12004];

x_vals   = [-2.71, -2.72, -2.71, -2.76, -2.71, -2.72, -2.78, ...
            -2.73, -2.76, -2.77, -2.77, -2.71, -2.71, -2.76, -2.72, -2.77];

dksc     = [-0.2871, -0.3888, -0.4886, -0.6934, -0.7897, -0.8903, -0.9938, ...
            -0.1882, -0.0904,  NaN, -0.3170, -0.4131, -0.5137, -0.6183, -0.7153, -0.8190];

 
% Initialize result array
Gamma = zeros(size(filelist));
dGamma = zeros(size(filelist));


for kk=1:length(filelist);
    
    %=================== Load .mat file and extract information ===================-0.2240   -0.3016   -0.7884
    filename=char(['dy0',num2str(filelist(kk),'%05d')]);
    disp(['Loading #' num2str(filelist(kk))]);
    load(['data\',filename,'.mat']);
    
    dK = meas.dK;
    numloops = meas.numloops;
    offspec = meas.endStatus.gamma-meas.endStatus.gammaSpecular;
    base_current = meas.ibase;
    Temp=str2num(meas.endStatus.tSample);
    E0=meas.beam.E0;
    alpha1=meas.tilt;
    
    polarisation=meas.mean.Preal+1i*meas.mean.Pimag;
    real_sig=meas.mean.Preal;
    imag_sig=meas.mean.Pimag;

    
    opt.forcesymm=0;
    opt.binning=1;
    [lambda_pos,energy,dE,spectrum,corrected_spectrum]=reconstruct_spectra_at677(base_current,polarisation,E0,alpha1,opt);
    
    data.current=base_current;
    data.polarisation=polarisation;
    data.ftpolarisation=spectrum;
    data.xaxisftspace=lambda_pos;
    
    xdat=dE;
    ydat=abs(corrected_spectrum);
    ydat=ydat./max(ydat);
    
    % Filter the region of interest
    mask = (xdat >= -3.5) & (xdat <= -2.2);
    x_fit = xdat(mask); y_fit = ydat(mask);
    
    [xData, yData] = prepareCurveData( x_fit, y_fit );
    % Set up fittype and options.
    ft = fittype( 'a*exp(-((x-b)^2)/(2*c^2))', 'independent', 'x', 'dependent', 'y' );
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.StartPoint = [0.1 -2.75 0.27];
    % Fit model to data.
    [fitresult, gof] = fit( xData, yData, ft, opts );
    c = fitresult.c;
    ci = confint(fitresult,0.9);   
    c_val = fitresult.c;       % Best-fit value of c
    c_low = ci(1,3);           % Lower bound of c
    c_high = ci(2,3);          % Upper bound of c
    
    Gam = 2.355 * abs(c);
    Gamma(kk)=sqrt(Gam.^2-0.06^2);
    FWHM_low = 2.355 * c_low;
    FWHM_high = 2.355 * c_high;
    dGamma(kk) = (FWHM_high - FWHM_low) / 2;

end

figure1=figure('units','centimeters','position',[5,4,16,12],'color','white','DefaultTextInterpreter','LaTex','DefaultLineLineWidth',2);
set(gca,'LineWidth',1.5,'FontSize',14,'Layer','top','TickLabelInterpreter','LaTex'); box(gca,'on'); hold(gca,'all');

%along GX plot
errbar(abs(dksc),Gamma,dGamma,'Color',cvec(1,:),'LineWidth',1.5);
p(2)= plot(abs(dksc),Gamma,'Marker','d','MarkerSize',6,'LineStyle','None','Color',cvec(1,:),'MarkerFaceColor',cvec2(1,:));
xlim([0,1.1]); ylim([0,2]);
xlabel('Wavecector (\AA$^{-1}$)'); ylabel('Linewidth $\Gamma$ of the mol. vibration (meV)');


%350 K measurement at I0/3
filelist = [011876,011875,11788, 11789, 11790, 11791, 11792, 11793, 11798, 11799, 11800, 11803, 11804];
x_vals   = [-2.72, -2.71 ,-2.75, -2.75, -2.74, -2.74, -2.741, -2.734, -2.67, -2.67, -2.77, -2.73, -2.85];
dksc     = [ -0.0878, -0.1371 ,-0.2715, -0.2915, -0.3157, -0.3408, -0.3659, -0.3904, -0.6356, -0.6856, -0.7433, -0.9392, -0.9967];


% Initialize result array
Gamma = zeros(size(filelist));
dGamma = zeros(size(filelist));


for kk=1:length(filelist);
    
    %=================== Load .mat file and extract information ===================
    filename=char(['dy0',num2str(filelist(kk),'%05d')]);
    disp(['Loading #' num2str(filelist(kk))]);
    load(['data\',filename,'.mat']);
    
    dK = meas.dK;
    numloops = meas.numloops;
    offspec = meas.endStatus.gamma-meas.endStatus.gammaSpecular;
    base_current = meas.ibase;
    Temp=str2num(meas.endStatus.tSample);
    E0=meas.beam.E0;
    alpha1=meas.tilt;
    
    polarisation=meas.mean.Preal+1i*meas.mean.Pimag;
    real_sig=meas.mean.Preal;
    imag_sig=meas.mean.Pimag;
    
    
    opt.forcesymm=0;
    opt.binning=1;
    [lambda_pos,energy,dE,spectrum,corrected_spectrum]=reconstruct_spectra_at677(base_current,polarisation,E0,alpha1,opt);
    
    data.current=base_current;
    data.polarisation=polarisation;
    data.ftpolarisation=spectrum;
    data.xaxisftspace=lambda_pos;
    
    xdat=dE;
    ydat=abs(corrected_spectrum);
    ydat=ydat./max(ydat);
    
    % Filter the region of interest
    mask = (xdat >= -3.5) & (xdat <= -2.2);
    x_fit = xdat(mask); y_fit = ydat(mask);
    
    [xData, yData] = prepareCurveData( x_fit, y_fit );
    % Set up fittype and options.
    ft = fittype( 'a*exp(-((x-b)^2)/(2*c^2))', 'independent', 'x', 'dependent', 'y' );
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.StartPoint = [0.1 -2.75 0.27];
    % Fit model to data.
    [fitresult, gof] = fit( xData, yData, ft, opts );
    c = fitresult.c;
    ci = confint(fitresult,0.9);   % 95% confidence interval
    c_val = fitresult.c;       % Best-fit value of c
    c_low = ci(1,3);           % Lower bound of c
    c_high = ci(2,3);          % Upper bound of c
    
    Gam = 2.355 * abs(c);
    Gamma(kk)=sqrt(Gam.^2-0.06^2);
    FWHM_low = 2.355 * c_low;
    FWHM_high = 2.355 * c_high;
    % Symmetric uncertainty (optional)
    dGamma(kk) = (FWHM_high - FWHM_low) / 2;

end

errbar(abs(dksc),Gamma,dGamma,'Color',cvec(4,:),'LineWidth',1.5);
p(1)=plot(abs(dksc),Gamma,'Marker','s','MarkerSize',6,'LineStyle','None','Color',cvec(4,:),'MarkerFaceColor',cvec2(4,:));
l1=legend(p,'0.2 ML CoPc','0.5 ML CoPc');
set(l1,'Interpreter','LaTeX');
xtickformat('%.1f'); ytickformat('%.1f');