function [y, Alpha, F] = late_reverb(To, L, f, alpha, Fs, N)
   // Late reverberation impulse response
   //
   // Usage:
   //        y = late_reverb(To, L, f, alpha, Fs, N)
   //
   // where 
   //        To:    Late reverberation duration (s)
   //        L:     Three component vector with dimensions
   //               of rectangular room (m)
   //        f:     Vector of frequencies where alfa is 
   //               specified (Hz)
   //        alpha: Vector of sound absorption coefficients
   //               at frequencies f
   //        Fs:    Sampling rate (Hz)
   //        N:     Size of filtering window
   //        y:     Impulse response  
   //
   // Fs is optional; its default value is 44100 Hz. N is
   // also optional; its default value is 2048. Typically
   // To is selected to be 1.5 or more times the largest
   // reverberation time so that the complete reverberation
   // tail is computed.
   //
   // Technical description
   // A filter with time-variable frequency response is
   // applied to a white noise sample of duration To. For
   // any given instant t the frequency response of this
   // filter is
   //
   //             -3*t/T60(f)
   // H(f, t) = 10
   //
   // where T60(f) is the reverberation time given by 
   // Sabine's formula
   //
   //                    V
   // T60(f) = 0.161 ----------
   //                S*alpha(f) 
   //
   // Since the filter is implemented with an FFT filter
   // technique, the response is updated only once every
   // N/2 samples. If not enough values of alpha are
   // given the rest are obtained by linear interpolation
   // and constant extrapolation. 
   // Often this function is used in combination with an
   // early response generator to simulate the impulse 
   // response of a natural reverberation.
   //
   // NOTE: timevarfilter.sci is required in Scilabs's 
   //       working directory 
   //
   // -----------------------
   // Author: Federico Miyara
   // Date:   2010-05-25
   //         2022-08-11
   //         2023-11-05

   // Default window size if omitted
   if nargin < 5
      N = 2048;
   end

   // Default sampling rate if omitted
   if nargin < 4
      Fs = 44100;
   end

   // Volume and area of room from its dimensions
   V = L(1)*L(2)*L(3);
   S = 2*(L(1)*L(2) + L(1)*L(3) + L(2)*L(3));

   // Complete alpha specification if too short for
   // subsequent interpolation (particularly, interp1 
   // does not support extrapolation, so it is necessary
   // to provide alpha for the minimum and maximum 
   // frequencies
   if f(1) > Fs/N
      f = [Fs/N, f];
      alpha = [alpha(1), alpha];
   end
   if f($) < Fs/2
      f = [f, Fs/2];
      alpha = [alpha, alpha($)];
   end

   // Linear interpolation of alpha extended to the N/2 
   // non-aliased frequencies representable with an FFT
   // of size N. Interpolation is performed with respect
   // to the logarithm of frequency
   F = [1:N/2]*Fs/N;
   Alpha = interp1(log(f), alpha, log([1:N/2]*Fs/N));

   // Sabine reverberation time
   // If the absorption coefficients were those according
   // to the physically meaningful definition, Eyring's 
   // equation should be used. However, most of the 
   // coefficients specified in the technical literature
   // have been measured using Sabine's approximation
   T60 = 0.161 * V/S./Alpha;

   // Vector of discrete times at which the frequency
   // response will be specified
   T = [1:N/2:Fs*To];

   // Frequency response. Column k contains the N/2 values
   // of the frequency response at time k*N/2/Fs. This
   // frequency response is kept constant up to time
   // ((k+1)*N/2 - 1)/Fs (i.e., for exactly N/2 samples)
   W = 10.^(-3 * kron(T/Fs, 1./T60'));

   // White noise generation
   x = rand(1, Fs*To, "normal");

   // Reverberation filtering
   exec morphfilter.sci;
   y = morphfilter(x, W, T, N);

endfunction
