// 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".

// Simulation réponse mode dans une pièce rectangulaire avec amortissement
// des modes

// Données d'entrées ----------------------------------------------------------
clear nx
clear ny
clear nz
clear tm;
clear resp;
clear dB;

// Color de la courbe
// 2 bleu 3 vert 4 bleu clair 5 rouge 6 mauve 7 jaune
clr = 2;

// Dimensions de la salle
//Dimension salle
lz = 6.5;
ly = lz*1.45;
lx = lz*2.10;

// Vitesse du son
C = 344;

// densité atmosphère
rho = 1.204;

// Pression de référence pour Lp
pref = 2e-5;

// Positions des sources
// x longueur
// y largeur
// z hauteur
xs = [0];
ys = [0];
zs = [0];

// Retard des sources (m)
tr = [0];

// Gain des sources
gs = [1];

// Position microphone
xm = lx ;
ym = ly ;
zm = lz ;

// Absorption des paroies
ax1 = 0.55 ;
ax2 = 0.55 ;

ay1 = 0.55 ;
ay2 = 0.55 ;

az1 = 0.55 ;
az2 = 0.55 ;

// Niveau ref (dB))
nref = 75-12 ;

// Echelle dB
dbMin = 35 ;
dbMax = 115 ;

dbMin = 35 ;
dbMax = 145 ;

// =========================
fmax = 100 // Frequence max de calcul des modes
nullMod = 0 ; // 1 : mode frequence nulle pris en compte
absorbtion = 0 ; // 0 : pas d'abosorbtion des paroies

// Functions
// Calcul fonction de green
// k nombre d'onde
// xm,ym,zm : position microphone
function green=greenfunction(k,xm,ym,zm)
    global fn ;
    global tm ;
    green = 0;

    // Boucle sur les modes
    for i=1:size(nx,"r"); 
      // Nombre d'onde des modes
      kx = nx(i)*%pi/lx;
      ky = ny(i)*%pi/ly;
      kz = nz(i)*%pi/lz;
      km = sqrt(kx^2+ky^2+kz^2);
      
      //disp("freq "+string(C * k / (2*%pi)))
      //disp("k "+string(k))
      
      // Calcul fonction psi au point microphone 
      psia = cos(kx*xm)*cos(ky*ym)*cos(kz*zm);
      // Boucle sur les source
      for j=1:size(xs,"r");  
        
        // Calcul facteur delay multiplié par gain
        Bn = gs(j)*exp(-%i*k*tr(j));
        
        // Calcul fonction psi au point source 
        psis = cos(kx*xs(j))*cos(ky*ys(j))*cos(kz*zs(j));
        
        // Calcul facteur normalisation
        if nx(i) == 0 then
            upx = 1;
        else
            upx = 2;
        end
        if ny(i) == 0 then
            upy = 1;
        else
            upy = 2;
        end
        if nz(i) == 0 then
            upz = 1;
        else
            upz = 2;
        end
        An = upx*upy*upz;
        
        //disp("mode "+string(i))
        //disp("psis "+string(psis))
        //disp("psia "+string(psia)) 
        //disp("km "+string(km))
        // disp("green "+string(-An*psis*psia/(k^2-km^2)/V))
        
        // Calcul extinction mode
        A = nx(i)/lx^2*log((1-ax1)*(1-ax2))+ny(i)/ly^2*log((1-ay1)*(1-ay2))+nz(i)/lz^2*log((1-az1)*(1-az2));
        fn(i) = C/2*sqrt((nx(i)/lx)^2+(ny(i)/ly)^2+(nz(i)/lz)^2);
        tm(i) = -55.3*fn(i)/C^2/A ;
                       
        // Somme 
        if absorbtion == 0 | nx(i)+ny(i)+nz(i) == 0  then
          green = green+Bn*An*psis*psia/(k^2-km^2);
        else
          green = green+Bn*An*psis*psia/(k^2-km^2-%i*km/(tm(i)*C/13.81));
        end
      end 
    end
    green = -green/V;
    //disp("green "+string(green))
    //disp("")
endfunction

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

// Volume de la salle
V = lx*ly*lz;

// Surface parois
S = 2*(lx*ly+lx*lz+ly*lz);

clear nx;
clear ny;
clear nz;
clear fn;
global fn ;
clear tm;
global tm ;

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

// 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( nullMod == 1 | ( nullMod <> 1 & i+j+k>0)) then
               // Calcul mode
               fmode = C/2*sqrt((i/lx)^2+(j/ly)^2+(k/lz)^2);
               
               // Test di mode dans les bornes
               if(fmode<=fmax) then
                 l=l+1;
                 nx(l) = i ;
                 ny(l) = j ;
                 nz(l) = k;
               end
             end
         end
     end
end

// Impression modes
disp("Modes")
for i=1:size(nx,"r")      
    kx = nx(i)*%pi/lx;
    ky = ny(i)*%pi/ly;
    kz = nz(i)*%pi/lz;
    km = sqrt(kx^2+ky^2+kz^2);
    fm = C * km / (2*%pi);
    disp("Mode ("+string(nx(i))+" "+string(ny(i))+" "+string(nz(i))+") "+string(fm)+" hz");
end
disp("Nombre de mode "+string(size(nx,"r")))

// Boucle sur les frequences
fq = linspace(5, 200, 2000);
for i = 1:size(fq,"c");
    // Nombre d'onde
    k = 2*%pi*fq(i)/C;   
    
    // Calcul green fonction
    pres(i)=abs(greenfunction(k,xm,ym,zm));
end

// Distance source microphone
rsm = sqrt((xs(1)-xm)^2+(ys(1)-ym)^2+(zs(1)-zm)^2) ;

// Débit source
// Niveau nref à la distance rsm
a = pref*4*%pi*rsm/rho*10^(nref/20) ;

// Pression reverbérée
pres = rho*a*pres ;

// Reponse reverbéré en dB
dB = 20*log10(pres/pref);

dBT = 10*log10(10^(nref/10)+10^(dB/10));


//plot2d(fq,dB,logflag="ln",style=clr,rect=[20,-50,80,70])
plot2d(fq,dBT,logflag="ln",style=clr,rect=[20,dbMin,fmax,dbMax])
//plot2d(fq,dB2,logflag="ln",style=clr,rect=[20,-20,100,40])
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";"40";"60";"80";"100";"200";"400";"800"];
ticks.locations = [20;30;40;60;80;100;200;400;800];
h.x_ticks = ticks;
title("Room h=6.5");
xlabel("Hz") ;
ylabel("dB") ;








