% moh83 - Finite-difference solution of 2-D cochlear model from MOH83

function fdm83()
global r1 r2 r4 srd y m;
%
% prepare for plotting
%
L=2.55;
%
figure(1);clf;
subplot(211);hold on;axis([0 L -110 -50]);
title('admittance');
ylabel('magnitude (dB)');
subplot(212);hold on;axis([0 L -0.5 0.5]);
xlabel('distance from stapes, cm');
ylabel('phase (cyc)');
%
figure(2);clf;
subplot(211);hold on;axis([0 L -3000 1000]);
title('impedance');
ylabel('real');
subplot(212);hold on;axis([0 L -3000 1000]);
xlabel('distance from stapes, cm');
ylabel('imaginary');
%
figure(3);clf;
subplot(211);hold on;axis([0 L 120 240]);
title('presure');
ylabel('magnitude (dB)');
subplot(212);hold on;axis([0 L -16 1]);
xlabel('distance from stapes, cm');
ylabel('phase (cyc)');
%
figure(4);clf;
subplot(211);hold on;axis([0 L -40 80]);
title('displacement');
ylabel('magnitude (dB)');
subplot(212);hold on;axis([0 L -16 1]);
xlabel('distance from stapes, cm');
ylabel('phase (cyc)');
drawnow;

% set parameter values
m=6; n=409; xl=2.55; yh=0.1; rho=1;

% declare arrays
a=zeros(m,m,n);
p=zeros(m,n);

% useful constructs
dx=xl/(n-1);
dy=yh/(m-1);
r1=(dx/dy)^2;
r2=r1+r1;
r4=2+r2;

% loop over frequencies, f=1000*2^fi
for fi=0:5,
   
   % set up for solution
   w=2*pi*1000*(2^fi);
   s=i*w;
   srd=i*4*rho*w*dy*r1;
   p0p=2*rho*w^2;
   x=transpose(linspace(0,xl,n));
   
   % compute admittance for all x
   z=Zbm(x,s,1);
   y=z.^(-1);
   
   % solve block matrix equation
   %  work down from top
  	p(:,1)=-2*dx*p0p;   
   a(:,:,1)=mxadd(1);
	b=inv(reshape(a(:,:,1),m,m));
	a(:,:,1)=reshape(-2*b,1,m,m);
  	p(:,1)=b*p(:,1);   
   for k=2:n-1,
      a(:,:,k)=a(:,:,k-1)+mxadd(k);
		b=inv(reshape(a(:,:,k),m,m));
		a(:,:,k)=reshape(-b,1,m,m);
  		p(:,k)=b*p(:,k-1);   
   end
   %  work up from bottom
   for k=n-2:-1:1,
      b=reshape(a(:,:,k),m,m);
		p(:,k)=p(:,k)-b*p(:,k+1);
   end
   d=y.*p(1,:).'./s;
   
   % plot data
   figure(1);
   subplot(211);plot(x,dbmag(y));
   subplot(212);plot(x,phase(y));
   figure(2);
   subplot(211);plot(x,real(z));
   subplot(212);plot(x,imag(z));
   figure(3);
   subplot(211);plot(x,dbmag(p(1,:)));
   subplot(212);plot(x,phase(p(1,:)));
   figure(4);
   subplot(211);plot(x,dbmag(d));
   subplot(212);plot(x,phase(d));
   drawnow;
end

function aa=mxadd(k)
global r1 r2 r4 srd y m;
aa=zeros(m,m,1);
aa(1,1,1)=r4+y(k)*srd;
aa(1,2,1)=-r2;
for j=2:m-1,
   aa(j,j-1,1)=-r1;
   aa(j,j,1)=r4;
   aa(j,j+1,1)=-r1;
end
aa(m,m-1,1)=-r2;
aa(m,m,1)=r4;

function y=dbmag(x)
y=20*real(log10(max(1e-9,x)));

function y=phase(x)
y=unwrap(imag(log(max(1e-9,x))))/(2*pi);

function z=Zbm(x,s,g)
K1=1e9*exp(-2.4*x);
R1=400+1e3*exp(-1.2);
M1=1e-3*exp(1.2*x);
K2=250*exp(-2.2*x);
R2=1e-4*exp(1.1*x)+8e-4*exp(-2.2*x);
M2=3e-9*exp(2.2*x);
R3=1;
Zb=K1/s + R1 + M1*s;
Za=-(R3)./(K2/s + R2 + M2*s);
z=Zb+g*Za;