// 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 volume et accord optimum

// Données d'entrées ----------------------------------------------------------

// Paramètres du haut-parleur monté dans l'enceinte
Re = 5 ;   // Résistance de la bobine mobile
Fso = 31;   // Frequence de résonnance (hz)
Qeo = 0.31 ;  // Facteur de qualité électrique
Qmo = 7.2 ;    // Facteur de qualité mécanique
Vas = 297;   // Volume equivalent à la raideur de la suspension (l)

// Paramètres de pertes de l'enceinte
Ql = 10;   // Facteur de perte par fuite
Qa = 100;   // Facteur de perte par absorption dans l'enceinte
Qp = 100;   // Facteur de perte par frottement dans l'évent

// Paramètres de l'évent (hz)
Fbmin = 15 ; 
Fbmax = 35 ;
Fbstep = 1 ;

// Paramètres de volume de l'enceinte (litres)
Vbmin = 100 ;
Vbmax = 1500 ;
Vbstep = 10 ;

// Depassement autorisé fonction de transfert
dbMax = 0 ;

// Paramétres des tracés 
Fmin = 10 ;    // Frequence minimum (hz)
Fmax = 200 ;   // Fréquence maximum (hz)
Nbp = 1000 ;   // Nombre de points

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

// Liste des frequences dont on veut la réponse
frq = linspace(Fmin,Fmax,Nbp)

// Variable de Laplace
s = poly(0,'s');

// Paramètres de l'enceinte
Qto = Qeo*Qmo/(Qeo+Qmo);

// Frequence de coupure min
f3min = 1000;

// Boucle sur le volume
for Vb = Vbmin:Vbstep:Vbmax
    // Boucle sur la frequence d'accord
    for Fb=Fbmin:Fbstep:Fbmax

      // Calcul parametre enceinte
      h = Fb/Fso;
      a = Vas/Vb ;

      // Calcul de la réponse en fréquence (amplitude db et phase deg)
      G = GBR(Fso,Qto,Ql,Qa,Qp,h,a);
      Gr = repfreq(G,frq);
      [Gd,Gp] = dbphi(Gr); 

      // Calcul frequence de coupure
      f = find(Gd<-3);
      i = f(length(f));
      f3 = frq(i);
      f3 = F3(Gd,f3);

      // Amplitude max
      [m,k] = max(Gd);
      
      // Test su meilleure réponse (f3 inférieure et pas de depassement de 3dB)
      if( f3<f3min & m<dbMax ) then
          f3min = f3 ;
          Vbopt = Vb;
          Fbopt = Fb;    
          dbopt = m ;     
      end 
      
      // Display
      disp("f3 "+string(f3)+" dbMax "+string(m)+" dB à "+string(frq(k))+" hz")
   end
end 

// Optimum
disp("f3 min "+string(f3min)+" hz");
disp("Vb optimal "+string(Vbopt)+" l");
disp("Fb optimal "+string(Fbopt)+" hz");
disp("db max "+string(dbopt)+" dB");
h = Fbopt/Fso;
a = Vas/Vbopt ;

// Tracés
scf(0);
plot2d(frq,Gd,style=11,logflag="ln",rect=[Fmin,-20,Fmax,10])
hg=gca();
hg.grid=[1,1]
hg.log_flags="lnn"; // x:log,y=log,z=lin
hg.sub_ticks=[0,0];
// X ticks
ticks = hg.x_ticks
ticks.labels = ["10";"20";"40";"60";"80";"100";"150";"200"];
ticks.locations = [10;20;40;60;80;100;150;200];
hg.x_ticks = ticks;











