function [y, env] = cochlea_travelling_wave(x, t, f, alpha, e)
   // Travelling waves in the cochlea
   //
   // Computes the shapes of the cochlea at distances
   // x from the oval window, at instants t when  
   // excited by sine waves of frequencies f and 
   // phases alpha 
   //
   // Usage 
   //       y = cochlea_travelling_wave(x, t, f, alpha, e)
   // where 
   //       x:     Vector of positions along the 
   //              cochlea in m
   //       t:     Vector of time instants in s
   //       f:     Vector of frequencies in Hz
   //       alpha: Vector of phases, in radians. Must 
   //              have the same length as f
   //       e:     Resonance enhancement exponent  
   //       y:     Array of dimensions length(x), 
   //              length(t), length(f) containing 
   //              the displacement of the cochles 
   //              from the equilibrium at the 
   //              corresponding values of x and t
   //              for the corresponding f, in 
   //              arbitrary units.  
   //       env:   Array of dimensions length(f),
   //              length(x) that contains the 
   //              envelopes for each frequency 
   //
   // The algorithm uses a transmission line electric
   // circuit model given in the references. In 
   // addition, resonant envelopes may be enhanced
   // (i. ., improve selectivity) by applying an
   // exponent to the envelopes normalized by their
   // peak value, so that the peak ei kept but the
   // neigbouring values decrease much faster than 
   // their original slopes.
   //
   // References
   // Deutsch, Sid; Micheli-Tzanakou, Evangelina (1987)
   //    "Neuroelectric Systems". New york University
   //    Press. New York 
   // Deutsch, Sid; Deutsch, Alice (1993) "Understanding 
   //    the nervous system. An engineeering perspective".
   //    IEEE Press.New York
   // Miyara, Federico (2024) "Análisis de un modelo
   //    circuital de la cóclea" 
   //
   // See also: cochlea_model()
   //
   // ----------------------------------
   // Author: Federico Miyara
   // Date:   2024-10-94
   
   // Missing arguments
   rhs = argn(2);
   if rhs<2
      t = 0;
   end
   if rhs<3
      f = 1000;
   end
   if rhs<4
      alfa = zeros(f);
   end
   if rhs<5
      e = 1;
   end
  
   // testing values --DELETE--
   if 1==2
      // Distances to the oval window
      x = 0:0.00001:0.035;
      // Frequencies
      f = [400 1000 2500];
      // Phases
      alpha = [0 0 0];
      // Times 
      t = 0:0.00005:0.005;
      // Resonance enhancement exponent
      e = 1;
   end
   
   // Transfer function of the cochlea
   H = cochlea_model(x,f);
   
   // Normalized envelopes (respect to the 
   // maximum of each ebvelope)
   env = abs(H)./ (max(abs(H),"c") * ones(1,size(H,2)));
   
   // Transfer function phase in radians
   phi = phasemag(H, "c")*%pi/180;
   
   // Initialize y
   y = zeros(length(t), length(x), length(f));
   for i=1:length(f)
      // Phases due to the transfer function and
      // the position on the cóchlea, replicated
      // once per time instant
      phis = repmat(phi(i,:), length(t), 1);
      // Phases due to the oval window input sine
      // waves (includes initial phase and time 
      // dependence for current frequency)
      psi = 2*%pi*f(i)*t(:) + alpha(i);
      // Phases replicated once per position
      psis = repmat(psi, 1, length(phi(i,:)));
      // Envelope for current frequency replicated 
      // once per time instant and, if required,
      // enhanced through the application of 
      // a resonance enhancement exponent e
      envs = repmat(env(i,:).^e, length(t), 1);
      // Traveling waves for current frequency
      // Create a matrix with as many rows as
      // time instants and as many columns as 
      // positions along the cochlea
      y(:,:,i) = y(:,:,i) + envs .* sin(phis + psis);
   end

endfunction


// internal function
function [H, f] = cochlea_model(x, f)
   // Transfer function of a circuit model of the cochlea
   // 
   // Usage:
   //       H = cochlea_model(x, f)
   // where
   //       x: Scalar or vector of distances in m 
   //          from the oval window
   //       f: Scalar or vector of frequencies in Hz
   //       H: Transfer function matrix. There is 
   //          one row per frequency and one  column
   //          per distance from the oval window
   //
   // The model consists in a tansmission line with
   // variable parameters along the basilar membrane. The
   // transmission line can be viewed as a mechanical one
   // or an electrical one. The model has been proposed 
   // by Deutsch & M-Tzanakou (see the Reference section).
   // The model allows to compute the transfer function 
   // between the displacement of the oval window at
   // the base of the cochlea (which is roughly 
   // proportional to the sound pressure at the eardrum)
   // and a point at a distance x from the oval window
   // and a frequency x. Both x and f can be vectors.
   //
   // MOTES: If f is not provided, the frequencies of
   //        the centre of the one third octave bands
   //        are used. If x is not provided x is a 
   //        71-column vector of distances separated
   //        by 0.5 mm 
   //
   // References
   // Deutsch, Sid; Micheli-Tzanakou, Evangelina (1987)
   //    "Neuroelectric Systems". New york University
   //    Press. New York 
   // Deutsch, Sid; Deutsch, Alice (1993) "Understanding 
   //    the nervous system. An engineeering perspective".
   //    IEEE Press.New York
   // Miyara, Federico (2024) "Análisis de un modelo
   //    circuital de la cóclea" 
   //
   // ----------------------------------------
   // Federico Miyara
   // 2024-09-26

   // Length of basilar membrane from the oval window
   // to the helicotrema, in m 
   L = 0.035;
   
   // Characterístic resistance of the cochlea in ohm
   // or N*s/m
   Ro = 0.02;
   
   // Mean width of basilar membrane in m
   aBM = 0.0003;
   
   // Mean cross sectional area of vestibular scala in m2
   A = 1.5e-6;
   
   if argn(2)<2   
      // Vector of distances from the oval window in m
      x = 0:0.0005:L;
   else
      // Ensure x is a row vector
      x = x(:)'; 
   end
   if argn(2)<1
      // Frequency vector in Hz
      f = 20*10.^(0:0.1:3)';
   else
      // Ensure f is a column vector
      f = f(:); 
   end
   
   // Model current or perilymph velocity amplitude
   // in A o m/s
   Io = 1

   // Auxiliary variable (nonlinear distance 
   // transformation)
      w = sqrt(L - x);

   // Impedance per unit length (Miyara, ec. 12,
   // Deutsch et al., ec 9.9)
   a1 = -100;
   a2 = -20;
   A1 = 0.39;
   A2 = 0.1325;
   Zx = A1*f.^2*exp(a1*w) + %i* A2*f*exp(a2*w);

   // Integrals of exp(a1*w) and exp(a2*w) respect to x
   int_1 = 2/a1.^2 * exp(a1*w).*(1 - a1*w);
   int_2 = 2/a2.^2 * exp(a2*w).*(1 - a2*w);

   // Integral of impedance per unit length respect to x
   int_Zx = -( A1*f.^2 * int_1 + %i* A2*f * int_2);

   // Complex integration constant to adjust the value
   // of current to Io at x = 0  (Miyara, ecs. 26 a 28)
   K1 =  A1*2/a1^2/Ro * exp(a1*sqrt(L))*(1 - a1*sqrt(L))
   K2 =  A2*2/a2^2/Ro * exp(a2*sqrt(L))*(1 - a2*sqrt(L))
   K = K1*f.^2 + %i*K2*f;

   // Current (or velocity) as a function of x
   I = Io*exp(int_Zx/Ro + K*ones(x));
   
   // Transfer function as a function of x
   H = -A/aBM * Zx/Ro .* I/Io;

endfunction


