function slbaverage
% slbaverage permet de calculer des valeurs moyennes et incertitudes (à 68%) à partir de plusieurs
% simulations. Slbaverage se base sur les noms de composantes et sur la nature de la 1ere simul selectionnee (frozen, slow motion, couplage hyperfins...).
evalin('base','clear'); %clear workspace de base
clc

[filename,path2simul]=uigetfile('*.slb','Select Spectra to Average','MultiSelect','on');

%coefficient de student pour 1 niveau de confiance d'1 ecart type (68%)
nbmesure=[2 3 4 5 6 7 8 9 10 20 40];
student=[1.84 1.32 1.20 1.14 1.11 1.09 1.08 1.07 1.06 1.03 1.01]; %d'après cours ENS

Syssimul={};
Exp=[];
%datasimul

if ~isreal(filename) %0 si cancel ou 1! simul selectionnee
    numslb=numel(filename); %nb de simul à moyenner
    for k=1:numslb
        disp(filename{k})
        load (fullfile(path2simul,filename{k}),'-mat');
        data2av{k}=datasimul;
        clear datasimul %on garde dernier Exp et Syssimul en mémoire (on suppose que meme Exp pour toutes les simuls, meme Nucs pour toutes les simuls)
    end
    numcmp=size(data2av{1},2); %nb de composante(s)
    
    simultype=data2av{1}(1).type; %Slow Motion ou Frozen Solution
    switch simultype
        case 'Frozen Solution'
            easyspinfunc=@pepper;
        case 'Slow Motion'
            easyspinfunc=@chili;
    end

    cmpname={data2av{1}(:).name}; %cellarray contenant le noms des composantes de la 1ere simul
    %%% preparation des donnees; si besoin remise dans l'ordre des composantes d'après leur nom (/ref)%%%
    for k=2:numslb
        datatemp=data2av{1,k};
        for ii=1:numcmp
            index=contains({datatemp(:).name},cmpname{ii});
            data2av{1,k}(ii)=datatemp(index);
        end
    end
    
    %%% preparation des donnees: weight de la 1ere composante=1 %%%
    for k=1:numslb
        datatemp=data2av{k};
%         weighttot=sum([datatemp.weightval]);
        weightnorm=data2av{k}(1).weightval;
        for ii=1:numcmp
            data2av{k}(ii).weightval=datatemp(ii).weightval/weightnorm;
        end
    end
    
    %%% ref=1ere simul %%%
    dataref=data2av{1};
    %datafieldname=fieldnames(dataref); %noms des champs

    
    %%% moyennage de chaque champ  de chaque composante %%%
    dataout=dataref;
    dataoutincert=dataref;
    for k=1:numcmp
        
        %weight gaussian lorentzian
        tempw=zeros(1,numslb);
        tempgaus=zeros(1,numslb);
        templor=zeros(1,numslb);
        
        for ii=1:numslb
            tempw(ii)=data2av{ii}(k).weightval;
            tempgaus(ii)=data2av{ii}(k).gaussianval;
            templor(ii)=data2av{ii}(k).lorentzianval;
        end
        %tempw
        [dataout(k).weightval,dataoutincert(k).weightval]=mean_err(tempw);
                %toto=dataoutincert(k).weightval
                %tata=dataoutincert(k).weightval
        Syssimul{k}.weight=dataout(k).weightval;
        [dataout(k).gaussianval,dataoutincert(k).gaussianval]=mean_err(tempgaus);
        [dataout(k).lorentzianval,dataoutincert(k).lorentzianval]=mean_err(templor);
        Syssimul{k}.lw=[dataout(k).gaussianval dataout(k).lorentzianval];
        
        %tcorr
        for jj=1:3
            tempt=zeros(1,numslb);
            tempg=zeros(1,numslb);
            for ii=1:numslb
                tempt(ii)=data2av{ii}(k).tcorrval(jj);
                tempg(ii)=data2av{ii}(k).gval(jj);
            end
            [dataout(k).tcorrval(jj),dataoutincert(k).tcorrval(jj)]=mean_err(tempt);
            [dataout(k).gval(jj),dataoutincert(k).gval(jj)]=mean_err(tempg);
        end
        Syssimul{k}.g=dataout(k).gval;
        Syssimul{k}.logtcorr=log10(dataout(k).tcorrval);
        
        %A
        for jj=1:2
            for kk=1:3
                tempA=zeros(1,numslb);
                for ii=1:numslb
                    tempA(ii)=data2av{ii}(k).Aval(jj,kk);
                end
                [dataout(k).Aval(jj,kk),dataoutincert(k).Aval(jj,kk)]=mean_err(tempA);
            end
        end
        Ainit=dataout(k).Aval;
        Ainit(isnan(Ainit))=[];
        Syssimul{k}.A=mt2mhz(Ainit); %en g=2
        clear Ainit
        
        %specsimul
        dataout(k).specsimul=easyspinfunc(Syssimul{k},Exp);
        dataoutincert(k).specsimul=zeros(1,1024);
        dataoutincert(k).specsimul(end)=1; %pour affichage dans simlabel
        
    end
    
    %sauvegarde
    [averagefile,path2averagefile]=uiputfile('*.slb','slb Average result',fullfile(path2simul,'*.slb'));
    if averagefile~=0
     datasimul=dataout;
     save(fullfile(path2averagefile,averagefile),'Syssimul','datasimul','Exp');
     clear datasimul
     datasimul=dataoutincert;
     save(fullfile(path2averagefile,[averagefile(1:end-4),'_incert.slb']),'Syssimul','datasimul','Exp');
     disp(['averaged and saved in ',averagefile,' and ',averagefile(1:end-4),'_incert.slb']);
    else
        disp('No result was saved.. Try again...');
    end
    
else
    disp('Try again...');
end


function [moy,incert]=mean_err(u)%,method)
% mean_err(u) retourne la valeur moyenne (moy) et l'incertitude (+-)
% (incert=t*std(u)/sqrt(nb de mesure) où t(coefficient de Student)) 
% d'un ensemble de plus de 2 mesures stockées dans le vecteur u.
% Affichage des données sous forme d'histogramme.
% Option: Si method=1 : recherche de mesures aberrantes
% les mesures ayant une occurence prévue < 0.5 sont éliminées. 
% (occurence prévue = nb de mesure * Proba(loi gaussienne normalisée))





nbelu=numel(u);

if numel(u)==2
    moy=mean(u);
    incert=mean([max(u)-moy,moy-min(u)]);

else
    
    

    if nbelu>40
        t=1;
    else
        t=interp1(nbmesure,student,nbelu);
    end

    moy=mean(u);
    devstd=std(u);
    %log10(devstd)


%     %%%%%%%%%% Recherche de mesure aberrantes si method=1? %%%%%%%%%%
%     if exist('method','var') && method==1
% 
%        %definition loi normale de probabilité (gaussienne)
%     demiplage=30*devstd;
%     pas=devstd/10;
%     x=-1*demiplage:pas:demiplage;
%     probanorm=gauss_distribution(x,0,devstd);
%        %recherche du nb alpha / nb occurence d'écart à la moyenne(=(1-P([x-alpha*stdev,x+alpha*stdev])*nbelu)<0.5
%     alpha=1;
%     incralpha=0.05;
%     while(nbelu*(1-integrale_a_b(x,probanorm,-alpha*devstd,alpha*devstd))>0.5)
%     alpha=alpha+incralpha;
%     end
%     alpha=alpha-incralpha;
%         %recherche des mesures/ écart à la moyenne>=90%*alpha*devstd
%     testaberr=abs(u-moy)./devstd/alpha; %test
%     isaberr= find(testaberr>0.9); %  indices des mesures aberrantes
% 
%     % supression de(s) mesure(s) aberrante(s) s'il en existe
%     if not(isempty(isaberr))
%         disp('Les moyenne et incertitude sont calculées en éliminant:');
%        for i=1:numel(isaberr)
%         disp([inputname(1),'(',num2str(isaberr(i)),')=',num2str(u(isaberr(i)))]);
%         plot(u(isaberr(i)),0.5,'X','MarkerSize',20,'LineWidth',4,'Color','g')
%        end
%        disp(' ')
%     u(isaberr)=[];
%     moy=mean(u);
%     devstd=std(u);
%     nbelu=numel(u);
%     end
% %%%%%%%%%%% fin de Recherche de mesure aberrantes %%%%%%%%%%
    
incert=t*devstd/sqrt(nbelu);

% if moy>1000
% text2disp=[num2str(moy,'%e'),177,num2str(incert,'%e')];
% else
% text2disp=[num2str(moy,'%f'),177,num2str(incert,'%f')];
% end
% disp([inputname(1),' = ',text2disp]);
% end

end
end



end