% Finite-difference solution of the 2-D cochlear model
% Neely (1981, JASA 69, 1386-1393)

function fdm81()
global r1 r2 r4 srd y m;
figure(1);clf;L=3.5;
subplot(211);hold on;axis([0 L -130 -30]);
title('admittance');
ylabel('magnitude');
subplot(212);hold on;axis([0 L -1 1]);
xlabel('distance from stapes, cm');
ylabel('phase/pi');
figure(2);clf;
subplot(211);hold on;axis([0 L 70 170]);
title('presure');
ylabel('magnitude');
subplot(212);hold on;axis([0 L -9 1]);
xlabel('distance from stapes, cm');
ylabel('phase/pi');
figure(3);clf;
subplot(211);hold on;axis([0 L -70 30]);
title('displacement');
ylabel('magnitude');
subplot(212);hold on;axis([0 L -9 1]);
xlabel('distance from stapes, cm');
ylabel('phase/pi');
drawnow;

% set parameter values
m=8; n=246; xl=3.5; yh=0.1; rho=1;
k0=1e9; r0=200; m0=0.15; k1=-2;

% 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 10 frequencies
for fi=0:9,
   w=2*pi*400*(1.414^fi);
   s=i*w;
   srd=i*4*rho*w*dy*r1;
   p0p=2*rho*w^2;
   
   % compute admittance for all x
   x=transpose(linspace(0,xl,n));
   y=(k0*exp(k1*x)/s + r0 + m0*s).^(-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,dbmag(p(1,:)));
   subplot(212);plot(x,phase(p(1,:)));
   figure(3);
   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))))/pi;