
function [d, f, fm] = bonello(Lx,Ly,Lz)
   // Criterio de Bonello para salas rectangulares
   // Obtiene la densidad de modos normales en las bandas
   // de tercio de octava desde 20 Hz a 200 Hz
   //
   // [d, f, fm] = bonello(Lx, Ly, Lz)
   // 
   //     [d, f] = bonello(Lx, Ly, Lz)
   //
   //          d = bonello(Lx, Ly, Lz)
   //
   // Lx, Ly, Lz: dimensiones de la habitación
   // fm:         frecuencias de los modos normales
   // f:          frecuencias de las bandas de tercio de octava
   // d:          densidad de modos en modos por octava
   //
   // ---------------------------------------------------------
   // Federico Miyara
   // 2007-03-19
   // 2007-03-22
   // 2007-04-01
   // 2024-08-01

   // Velocidad del sonido
   c = 345;

   // Frecuencias normalizadas
   f = [20, 25, 31.5, 40, 50, 63, 80, 100, 125, 160, 200];

   // Frecuencias exactas
   // NOTA: Si se trabajara con las frecuencias normalizadas los 
   //       límites comunes a dos bandas contiguas calculados como
   //       tercios de octava exactos no coincidirían
   ff = 1000*10.^([-17:-7]/10);

   // Máxima frecuencia a considerar
   maxfrec = (20/Lx)^2 + (20/Ly)^2 + (20/Lz)^2;

   // Máximo orden de armónico correspondiente a ondas axiales
   // de modo de incluir hasta la banda de tercio de octava de 
   // 200 Hz  
   Nx = ceil(200*2^(1/6) * 2 * Lx / c);
   Ny = ceil(200*2^(1/6) * 2 * Ly / c);
   Nz = ceil(200*2^(1/6) * 2 * Lz / c);

   // Frecuencias de los modos normales 
   fm = [];
   for l = 0:Nx
      for m = 0:Ny
         for n = 0:Nz
            nextf =   sqrt((l/Lx)^2 + (m/Ly)^2 + (n/Lz)^2)*c/2
            if nextf <= 200*2^(1/6)
                fm = [fm, nextf];
            end
         end
      end
   end

   // Reordenamiento de menor a mayor
   fm = gsort(fm, "g", "i");

   // Clasificación en casilleros de un tercio de octava
   // NOTA: un tercio de octava es aproximadamente un décimo de
   //       década
   d = zeros(1,11);
   ko = 2;
   for i = 1:11
      for k = ko:length(fm)
         if fm(k) <= ff(i)*10^(1/20)
            d(i) = d(i) + 1;
         else
            ko = k;
            break
         end
      end
   end

endfunction


