// Copyright J. Fourcade (volucres.fr)
// Ce logiciel est régi par la licence CeCILL soumise au droit français et
// respectant les principes de diffusion des logiciels libres. Vous pouvez
// utiliser, modifier et/ou redistribuer ce programme sous les conditions
// de la licence CeCILL| telle que diffusée par le CEA, le CNRS et l'INRIA 
// sur le site "http://www.cecill.info".

// Calcul de l'amortissement des modes
// Données d'entrées ----------------------------------------------------------
// Vitesse du son
c = 344;

//Dimension salle
lz = 4;
ly = lz*1.45;
lx = lz*2.10;

lz=9.42;
ly=lz;
lx=lz;

lz=5.5;
ly=11;
lx=15;

// Coef absorption general
alpha = 0.55 ;

// Coef absorbtion paroies
ax1 = 0.55 ;
ax2 = 0.55 ;

ay1 = 0.55 ;
ay2 = 0.55 ;

az1 = 0.55 ;
az2 = 0.55 ;

// Limitation des modes
fnmin = 0;
fnmax = 400;

// Affichage mode
nmodemax = 10 ;

// Calculs --------------------------------------------------------------------

// Frequence charnière 1/3 octave
fch = [13.9 ; 17.5 ; 22.1 ; 27.8 ; 35.1 ; 44.2 ; 55.7 ; 70.2 ; 88.4 ; 114.4 ; 140.3 ; 176.8 ; 222.7 ; 280.6 ; 353.6]
fcentrale = [16 ; 20 ; 25 ; 31.5 ; 40 ; 50 ; 63 ; 80 ; 100 ; 125 ; 160 ; 200 ; 250 ; 315]

// Coefficients modes et parcours
clear Ci;
clear Cj;
clear Ck;
clear res;
clear dns;
clear lmoy;
clear trmoy;

// Calcul valeurs max des coefficients i, j, k
imax = floor(fnmax*2*lx/c)+1 ;
jmax = floor(fnmax*2*ly/c)+1 ;
kmax = floor(fnmax*2*lz/c)+1 ;

disp("Indice max nx "+string(imax)) ;
disp("Indice max ny "+string(jmax)) ;
disp("Indice max nz "+string(kmax)) ;

// Parcours moyenné
lm2 = 0;

// Temps de réverberation moyenné
tr2 = 0 ;

// Boucle sur les coefficiens i,j,k
l=0
for i=0:imax
    for j=0:jmax
        for k=0:kmax
            
            // On supprime le mode nul
            if(i+j+k>0) then
               // Calcul mode
               fn = c/2*sqrt((i/lx)^2+(j/ly)^2+(k/lz)^2);
               
               // Test di mode dans les bornes
               if(fn>=fnmin & fn<=fnmax) then
                 l=l+1;
            
                 // Coefficients 
                 Ci(l)=i;
                 Cj(l)=j;
                 Ck(l)=k;
            
                 // Calcul mode
                 res(l,1) = fn;
               
                 // Calcul parcours moyen
                 res(l,2) = 2*fn/c/(i/lx^2+j/ly^2+k/lz^2) ;
                 lm2 = lm2 + res(l,2);
               
                 // Calcul temps de reverberation
                 A = i/lx^2*log((1-ax1)*(1-ax2))+j/ly^2*log((1-ay1)*(1-ay2))+k/lz^2*log((1-az1)*(1-az2));
                 res(l,3) = -55.3*fn/c^2/A ;
                 tr2 = tr2 + res(l,3) ;
               end
            end
        end
    end
end
        
// Tri des modes
[res2,index] = gsort(res(:,1),'g','i');
res3 = res(index,:);

// Calcul parcours moyen statistique
S = 2*(lx+ly)*lz + 2*lx*ly;
V = lx*ly*lz ;
disp("Volume : "+string(V)+" m3") ;
disp("Surface sol : "+string(lx*ly)+" m2");
disp("Surface parois : "+string(S)+" m2");

// Nombre de mode
disp("Nombre de mode: "+string(l)) ;

// Parcours moyen
lm2 = lm2/l;
disp("Parcours moyen sur la bande de fréquence : "+string(lm2)+" m");
lms = 4*V/S;
disp("Parcours moyen statistique : "+string(lms)+" m");
disp("Ratio parcours : "+string(lm2/lms));

// Tr Eyring
trs = -0.16*V/S/log(1-alpha) ;
disp("Tr Eyring  : "+string(trs)+" s");

// Tr moyen
tr2 = tr2/l;
disp("Tr moyen des modes : "+string(tr2)+" s");
tr2norm = tr2 /trs ;
disp("Ratio Tr: "+string(tr2norm));

scf(0);
plot2d(res3(:,1),res3(:,2),-3);
xtitle("Parcours (h=2.8)")  
xlabel("hz")
ylabel("mètre")
scf(1);
plot2d(res3(:,1),res3(:,3),-3);
xtitle("Temps de décroissance des modes (h=4)")  
xlabel("hz")
ylabel("sec")

// Calcul ecarts des modes
// Bande des modes
Bm = 2.2/trs/3 ;
disp("Bande mi puissance mode  "+string(Bm)+"hz") ;
fSch = 0 ;
for i=1:l-1
    df = res3(i+1,1)-res3(i,1) ; 
    //disp("fn "+string(res3(i,1))+" hz df "+string(df)+" hz")
    if df <= Bm 
        // Cas critère respecté
        if fSch == 0
          fSch = res3(i,1) ;
        end  
    else
        fSch = 0;
    end 
end
fSchT = 2000*sqrt(trs/V) ;
disp("Frequence de Schroeder théorique "+string(fSchT)+" hz")
disp("Frequence de Schroeder calculer sur les modes "+string(fSch)+" hz")

// Calcul courbe de densite
for i=1:size(fch,1)-1
    nbmod = 0 ;
    lmoy(i) = 0;
    trmoy(i) = 0;
    for j=1:l
       // On compte les modes dans l'intervale
       if res3(j,1)>=fch(i) & res3(j,1)<=fch(i+1)
          nbmod = nbmod + 1;
          
          // Calcul libre parcours moyen
          lmoy(i) = lmoy(i)+res3(j,2);
          
          // Calcul libre parcours moyen
          trmoy(i) = trmoy(i)+res3(j,3);
       end
    end
    dns(i) = nbmod;
    if nbmod <> 0 
        lmoy(i) = lmoy(i)/nbmod ;
        trmoy(i)= trmoy(i)/nbmod;
    end
end


scf(2);
plot2d(fcentrale,dns,-3,logflag="ln",rect=[10,0,50,12]);
xtitle("Densité de modes / tiers octave (h=6.5)")  
h=gca();
h.grid=[1,1]
h.log_flags="lnn"; // x:log,y=log,z=lin
h.sub_ticks=[0,0];
// X ticks
ticks = h.x_ticks
ticks.labels = ["16";"20";"30";"50";"80";"125";"200";"300";"400";"800"];
ticks.locations = [16;20;30;50;80;125;200;300;400;800];
h.x_ticks = ticks;
xlabel("hz")

scf(3);
plot2d(fcentrale,trmoy,-3,logflag="ln",rect=[20,0,400,0.6]);
xtitle("Temps de décroissance des modes / tiers octave (h=6.5)")  
h=gca();
h.grid=[1,1]
h.log_flags="lnn"; // x:log,y=log,z=lin
h.sub_ticks=[0,0];
// X ticks
ticks = h.x_ticks
ticks.labels = ["20";"30";"50";"80";"125";"200";"300";"400";"800"];
ticks.locations = [20;30;50;80;125;200;300;400;800];
h.x_ticks = ticks;
xlabel("hz")
xlabel("s")



