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" 
   //
   // See also: cochlea_travelling_wave()
   // 
   // ----------------------------------------
   // Federico Miyara
   // 2024-09-26
   
   // Missing arguments
   rhs = argn(2);
   if rhs<2 
      f = 1000;
   end
   if rhs<1
      x = 0:0.0005:0.035;
   end

   // 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



