%			
%	PHYSICAL PARAMETERS ASSUMED FOR THE HUMAN COCHLEAR PARTITION	
%
%	The following parameters are plotted:
%
%	1) mass per unit length [Kg/m]
%	2) basilar membrane stiffness per unit length [Kg/(m·sec²)]
%	3) absolute viscosity per unit length [Kg/(m·sec)]
%	4) shear viscosity coefficient per unit length [Kg·m/sec]
%	

N=300;	% Number of points on the BM

x=(1:N)/N;

BETA=0.0335;  % BM length in meters, 1 meter = 29.85 BETA

%--------------------------- MASS ------------------------------

m = Corti_mass(x);	% organ of Corti mass in Kg/BETA
m=m/BETA;	% organ of Corti mass in Kg/m

%----------------------- STIFFNESS -------------------------

ko = 1.675e+3;		% BM stiffness at base in Kg/(sec^2*BETA)
expk = -3.6*log(10); 	% BM stiffness exponent (3.6 orders of magnitude
			% decrease from base to apex)
k=ko*exp(expk*x);
     k=k/BETA;			% BM stiffness in Kg/(sec^2*m)

%------------------------ ABSOLUTE VISCOSITY ----------------------
% 'absolute' means with respect to motions orthogonal to the BM plane
% in order to be distinguished from 'shear' viscosity, which depends on
% relative motions between adjacent portions of BM 

ho = 10e-3; 	% Abs. BM viscosity at base in Kg/(BETA sec)
		% (about 30 times water viscosity)
 
exph = -log(4);	% BM decreases 4 times

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% abs. viscosity is corrected by terms accounting for irregularities
% at the BM ends:  
 
Ce=5;		% relative damping at the helicotrema
Ce=Ce*exp(-1/0.03+exph);	
Cb=5;		% relative damping at the base
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h =ho*( exp(exph*x) + Ce*exp(x/0.03) + Cb*exp(-x/0.005) );

% Absolute viscosity vector in  Kg/(BETA sec)

h=h/BETA; % Absolute viscosity vector in  Kg/(m sec)

% -------------------- SHEAR VISCOSITY ---------------------

micron = 29.85e-6; % 1 micron in BETA
h0 = 40*micron;	   % 40 micron,  
		   % putative height of the organ of Corti at base
		   % (effective for shear viscosity) 

h = h0*exp(log(4)*x);	% effective organ of Corti thickness
			% (grows 4 times from base to apex)

w0 = 50*micron;	  % 50 micron, mean radial width of the organ of Corti at base
		  % (effective for shear viscosity) 

w = w0*exp(log(4)*x); % effective organ of Corti radial width
		      % (grows 4 times from base to apex)

eta = 1e-3/29.85; %= 3.35e-5 Kg/(BETA*sec) 
		  % 10^{-3} Kg/(m*sec) = viscosity coefficient of water
		  % at 20 Celsius degrees (very close to 1 centipoise)	

s = eta*w.*h; 	  % Sectional shear viscosity constant, assuming that
		  % adjacent sectional segments of the organ of Corti
		  % are viscous like water. Unit dimension of s is Kg*BETA/sec
  	
s=35*s;	% Shear viscosity is multiplied by the phenomenologic factor 35
	% this value is consistent with observed values of travelling
	% wave width in various conditions
 
s=s*BETA; % Now unit dimension of s is Kg*m/sec

% ------------------- PLOTS ----------------

figure(1)
%axis([0 1 0 6e-6]); 
subplot(221),
plot(x,m);
title('Mass');
ylabel('m(x)  [Kg/m]');
% text(.2, .875, 'a)','Units', 'normalized');
%-------------------------
%axis([0 1 0 0.2]);
subplot(223),
plot(x,h)
title('Abs. viscosity');
ylabel('h(x)  [Kg/(m·sec)]');
xlabel('Fractional distance from stapes');
% text(.7, .875, 'b)','Units', 'normalized');
%--------------------------
%axis([0 1 0 1e-10]);
subplot(224),
plot(x,s)
title('Shear viscosity');
ylabel('s(x)  [Kg·m/sec]');
xlabel('Fractional distance from stapes');
% text(.2, .375, 'c)','Units', 'normalized');
%-------------------------------
%axis([0 1 1 7]);
subplot(222),
semilogy(x,k)
title('Stiffness');
ylabel('k(x)  [Kg/(m·sec²)]');
% text(.7, .375, 'd)','Units', 'normalized');
axis;
