function [T1,T2,ir12,ir11] = coupled_reverb(V,S,S12,alpha,T,Fs)
   // This function computes the reverberation times
   // and impulse responses of two acoustically coupled
   // rooms connected by an opening, at 500 Hz, 60 % RAH 
   // and 25 °C
   //
   // Usage:
   //       [T1,T2,ir12,ir11] = coupled_reverb(V,alpha,S,S12,T)
   // where
   //       T1, T2: Reverberation times in s
   //       ir12:   Impulse response envelope from source room
   //               power to receiver room energy density 
   //       ir11:   Impulse response envelope from source room 
   //               to source room energy density
   //       V:      Volume of rooms, [V1, V2], in m3
   //       S:      Inner area of rooms, [S1, S2], including
   //               common opening, in m2
   //       S12:    Opening area, in m2
   //       alpha:  Average absorption coefficients of rooms,
   //               [alpha1, alpha2]
   //       T:      Duration of impulse response, in s 
   //       Fs:     Sample rate, in Hz
   //
   // The formulas used are derived from an energy balance
   // where the sound power radiated on each room from the
   // other must equate the power that is absorbed by its 
   // surfaces. Only the reverberant field is considered.
   // If T is not provided it is taken as the largest of the 
   // reverberation times. If Fs is not provided its default
   // value is 44100. 
   // NOTE: Reverberation times must not be interpreted as 
   //       the reverberation times of each room but a joint
   //       dynamic property of the reverberation tail as
   //       a whole
   //
   // -----------------------------------------
   // Author: Federico Miyara
   // Date:   2016-05-22
   //         2022-04-19
   //         2022-07-19
   //         2022-08-05

   // Atmospheric absorption at f = 500 Hz, RAH = 60 %, 
   // T = 25 °C in Np/m
   gamma_m = 0.318/434;

   // Sound velocity in m/s 
   c = 345;
   
   // Average absorption coefficients
   alpha1_av = (alpha(1)*(S(1) - S12) + S12)/S(1);
   alpha2_av = (alpha(2)*(S(2) - S12) + S12)/S(2);
   
   // Parameters of the equation system 
   K1 = (gamma_m - S(1)/4/V(1) * log(1 - alpha1_av)) * c*V(1);
   K2 = (gamma_m - S(2)/4/V(2) * log(1 - alpha2_av)) * c*V(2);
   K12 = S12 / (alpha(1)*(S(1) - S12) + S12);
   K21 = S12 / (alpha(2)*(S(2) - S12) + S12); 

   aa = V(1)*V(2);
   bb = V(1)*K2 + V(2)*K1;
   cc = K1*K2*(1 - K12*K21);

   // System poles
   s1 = (-bb + sqrt(bb^2 - 4*aa*cc)) / (2*aa);
   s2 = (-bb - sqrt(bb^2 - 4*aa*cc)) / (2*aa);

   // Time constants
   tau1 = -1/s1;
   tau2 = -1/s2;

   // Reverberation times
   T1 = tau1 * log(10^6);
   T2 = tau2 * log(10^6);
   
   // Missing arguments
   if argn(2)<6
      // Default sample rate
      Fs = 44100;
   end
   if argn(2)<5
      // Default duration
      T = max(T1,T2)
   end
   
   // Time vector
   t = 0:1/Fs:T;
    
   // Impulse response envelope from acoustic power source
   // at room 1 to reverberant field at room 2 
   A = K12*K1/V(1)/V(2);
   A1 =  A / (s1 - s2);
   A2 = -A / (s1 - s2);
   ir12 = A1*exp(s1*t) + A2*exp(s2*t);
   
   // Impulse response envelope from acoustic power source
   // at room 1 to reverberant field at room 1
   B1 = (K2/V(1)/V(2) + s1/V(1)) / (s1 - s2);
   B2 = -(K2/V(1)/V(2) + s2/V(1)) / (s1 - s2);
   ir11 = B1*exp(s1*t) + B2*exp(s2*t); 

endfunction
