% Pbmodel1.m
%
%
%
%   		1D Linear Cochlear Model for the gerbil
%
%					R. Naidu and D. Mountain
%				Cochlear Biophysics Laboratory
%			Boston University Hearing Research Center
%
% The model computes Pressure difference and the BM volume
% velocity by picking a number for the fluid volume velocity
% through the helicotrema and then working back towards
% the base.  At each section the fluid volume velocity is computed
% by adding the fluid volume velocity from the previous section to
% the basilar membrane volume velocity.  The pressure drop across the
% fluid due to this velocity is then computed and added to the pressure
% difference across the cochlear partition.  Since the model is linear,
% the results above can then be scaled to the pressure difference at the
% stapes.

%The following are for gerbil, other frequency-place maps are available at:
%   http://earlab.bu.edu/anatomy/cochlea/freq_place_maps/Greenwood.htm

L   = 1.21;		% length of cochlea (cm)
a   = 398.;		%Greenwood coefficient
b   = 2.2;		%Greenwood coefficient
k   = 0.631;	%Greenwood coefficient
 
N = 512; 		% No. of sections
T = L/N; 		% spatial step size

rho = 1.02;		% Fluid density (gm/cm^3)

P  = zeros(1,N);	% Pressure difference array (dyne/cm)
Zp = zeros(1,N);	% Cochlear Partition Impedance array (cgs Ohms)
Vbm= zeros(1,N);	% Basilar membrane volume velocity array (cm^3/s)
						% To get to mechanical velocity (cm/s), this needs to be mulitiplied
						% by 1/(T*Wbm) where Wbm is the effective width of the BM

f = input('Frequency of tone ==>');
q = input('Quality factor (try 10) ==>');

w = 2*pi*f;

for i = 1:N					% Compute model parameters
   
	Wbm(i) = 0.0175;		%Width of the BM (cm)
	y(i)   = 0.01;			%Thickness of the organ of Corti (cm)
	As(i)  = .0002;		%Cross-sectional area of scalae (cm^2)
	Q(i)   = q;				%Quality factor of the RLC circuit
	x(i)   = i*T;			%Distance from base (cm)

	Ls(i)  = rho*T/As(i);						%Mass of scalae
	CF(i)  = a*(10^(((L-x(i))/L)*b) - k);	%Frequency-place map
	Cbm(i) =  exp(-30.9 + 4.04*x(i));		%Volume compliance of cochlear partition
	Lbm(i) = 1/(Cbm(i)*(2*pi*CF(i))^2);		%Effective Mass of cochlear partition
	Rbm(i) = (1/Q(i))*(Lbm(i)/Cbm(i))^0.5;	%viscous resistance of cochlear partition
       
%Cochlear partition impedance
	Zp(i)  = (-w^2*Lbm(i)*Cbm(i) + 1 + Rbm(i)*j*w*Cbm(i))/(j*w*Cbm(i));   %j = sqrt(-1)
	end
      
% define the conditions at the helicotrema modeling it as a resistance:
%
Rh  = (Ls(N)/Cbm(N))^0.5;		%Match the cochlear impedance	
Vh  = 1;								%Pick a velocity at the helicotrema
Ph = Rh*Vh;
            
for i = N:-1:1						% Start from the apex and solve for P and vBM
    if i == N
       P(i) = Ph;    			% Boundary conditions at helicotrema
       V(i) = Vh + P(i)/Zp(i);
      
    else 
       P(i) = P(i+1) + V(i+1)*j*w*Ls(i+1);
       V(i) = V(i+1) + P(i)/Zp(i);
    end

    Vbm(i)= P(i)/Zp(i);
   
end
norm_Vbm = abs(Vbm)/abs(P(1));	% Scale for input pressure of 1 dyne/cm^2 inside stapes
norm_P  = abs(P)/abs(P(1));		% which corresponds to approximately 40 dB SPL
norm_Zp = abs(Zp);


theta_P=unwrap(angle(P))-angle(P(1));
theta_Vbm=unwrap(angle(Vbm)-angle(P(1)));
theta_Zp=unwrap(angle(Zp));

MagFig=figure;				%Start generating the magnitude figure (assume 1024x768 display
set(MagFig, 'Position',[5  200   500   500])	%Position it on the left side of screen
shg;

subplot(3,1,1)
semilogy(x, norm_Vbm);
axis([0 1.4 1e-10 1e-6]);
grid;
ylabel('BM Volume Velocity');
title('Magnitude (cgs units)');

subplot(3,1,2)
semilogy(x, norm_P);
axis([0 1.4 0.01 1]);
grid;
ylabel('Pressure Difference');

subplot(3,1,3)
semilogy(x, norm_Zp);
axis([0 1.4 1e6 1e10]);
grid;
xlabel('Distance from Base (cm)');
ylabel('CP Impedance');

PhaseFig=figure;				%Start generating the magnitude figure (assume 1024x768 display
set(PhaseFig, 'Position',[512  200   500   500])	%Position it on the left side of screen

subplot(3,1,1)
plot(x, theta_Vbm);
%axis([0 1.4 -5 5]);
grid;
ylabel('BM Velocity');
title('Phase (radians)');

subplot(3,1,2)
plot(x, theta_P);
%axis([0 1.4 -4 0 ]);
grid;
ylabel('Pressure Difference');

subplot(3,1,3)
plot(x, theta_Zp);
axis([0 1.4 -2 2]);
grid;
xlabel('Distance from Base (cm)');
ylabel('CP Impedance');
shg;