clear all; close all;


%=================== File list to read in ===================
addpath /home/shared/SpinEcho/Spinecho2Icon/data/Ag100/

filelist=[012106:012113,012093:012105,012116:012117]; 

lambda=zeros(2050,length(filelist))*NaN;
delE=zeros(2050,length(filelist))*NaN;
speclambda=zeros(2050,length(filelist))*NaN;
specenergy=zeros(2050,length(filelist))*NaN;
lambdatilted=zeros(2050,length(filelist))*NaN;
spectilted=zeros(2050,length(filelist))*NaN;

figure1=figure('units','centimeters','position',[5,4,16,22],'color','white','DefaultTextInterpreter','LaTex','DefaultLineLineWidth',2);
axes1 = axes('Parent',figure1,'LineWidth',1.5,'FontSize',16,'Layer','top'); box(axes1,'on'); hold(axes1,'all');

%=================== Load .mat file and extract information ===================
for kk=1:length(filelist);
    
    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(kk)=str2num(meas.endStatus.tSample);
    E0=meas.beam.E0;
    alpha1=meas.tilt;
    alpha(kk)=alpha1;

    base_current=[-1*fliplr(base_current(2:end)),base_current];
    real_sig=[fliplr(meas.mean.Preal(2:end)),meas.mean.Preal];
    imag_sig=[-1*fliplr(meas.mean.Pimag(2:end)),meas.mean.Pimag];
    polarisation=real_sig(:)'+1i*imag_sig(:)';

    opt.forcesymm=0;
    opt.binning=0;

    [lambda_pos,tiltlambda,energy,dE,bsf,spectrum,corrected_spectrum]=reconstruct_spectra_tilted(base_current,polarisation,E0,alpha1,opt);
    

    lambda(1:length(lambda_pos),kk)=lambda_pos;
    delE(1:length(dE),kk)=dE;
    speclambda(1:length(spectrum),kk)=abs(spectrum);
    specenergy(1:length(corrected_spectrum),kk)=abs(corrected_spectrum);
    lambdatilted(1:length(tiltlambda),kk)=tiltlambda;
    spectilted(1:length(bsf),kk)=abs(bsf);
    
    xdata=abs(dE(-2.4<dE & dE<-1.4)*1000);
    ydata=corrected_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); 

    plot(xdata,ydata+(kk-1)*0.2,'o');
    plot(xfit,yfit+(kk-1)*0.2,'-');
        
    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))/2;
    
        
end

data.xdata=delE;
data.ydata=specenergy;
data.Temp=Temp;
%save('LWdata.mat','data');

hold off; grid on;