%
%  the Zwislocki/Von B?k?sy model
%
x=[0:0.1:3.5]; % cm
f=[25, 50, 100, 200, 400, 800, 1600, 3200, 6400, 12800]; % Hz
arg1 = exp(2.5*x)-1;
arg2 = exp(1.25*x)-1;
BM=zeros(length(f), length(x));
figure(1)
clf

for n=1:10
    c1=1.5*10^(-10)*f(n)^2;
    % c2=5*10^(-4)*sqrt(f(n));
    c2=1.84*10^(-3)*sqrt(f(n));
    arg3 = x - c1*arg1-c2*arg2;
    BM(n,:) = exp(arg3);
    maxima(n) = find(BM(n,:)>=max(BM(n,:)));
    mx=max(BM(n,:));
    BM(n,:)=BM(n,:)/mx;
end

for n=1:2:9
    subplot(5,1,floor(n/2)+1);
    plot(x, 0.9*BM(n,:), 'k-','linewidth',2)
    s=sprintf('%4d Hz',f(n));
    text(0.1, 0.85*max(BM(n,:)), s,'fontsize',11,'fontweight','bold')
    set(gca,'ytick',[]); set(gca,'xticklabel',[]);
    set(gca,'fontsize',13, 'fontweight','bold') 
end
xlabel('Distance from stapes (mm)','fontsize',14, 'fontweight','bold')
set(gca, 'xtick',[0 1.0 2.0 3.0]);
set(gca,'xticklabel',[' 0';'10';'20';'30']);
print -djpeg Fig5-13a

figure(2)
clf
semilogx(f,-maxima, 'ko-','markerfacecolor','k', 'linewidth',2)
axis([0 13000 -40 0]); set(gca, 'yticklabel',[]);
set(gca,'ytick',[-40 -30 -20 -10 0]);
set(gca,'yticklabel',['40';'30';'20';'10';' 0']);

set(gca,'fontsize',13, 'fontweight','bold')
xlabel('Frequency (Hz)','fontsize',14, 'fontweight','bold')
ylabel('Distance from stapes (mm)','fontsize',14, 'fontweight','bold')
print -djpeg Fig5-13b

%
%   the pressure wave
%
     
  x=[0:0.1:3.5]; % cm
  RelP=zeros(length(f), length(x)); 
  arg1=exp(2.5*x)-1;   arg2 = exp(1.25*x)-1;
  
  for n=1:10
    c1=1.5*10^(-10)*f(n)^2;
    %c2=5*10^(-4)*sqrt(f(n));
    c2=1.84*10^(-3)*sqrt(f(n));
    arg4 = -0.5*x - c1*arg1-c2*arg2;
    RelP(n,:) = exp(arg4);
    maxima(n) = find(RelP(n,:)>=max(RelP(n,:)));
    mx=max(RelP(n,:));
    RelP(n,:)=RelP(n,:)/mx;
    end
    
figure(3); clf;
for n=1:2:9
    subplot(5,1,floor(n/2)+1);
    plot(x, 0.9*RelP(n,:), 'k-','linewidth',2)
    s=sprintf('%4d Hz',f(n));
    text(3, 0.85*max(RelP(n,:)), s,'fontsize',11,'fontweight','bold')
    set(gca,'ytick',[]); set(gca,'xticklabel',[]);
    set(gca, 'xtick', [0 1 2 3]);
    set(gca,'fontsize',13, 'fontweight','bold') 
end
xlabel('Distance from stapes (mm)','fontsize',14, 'fontweight','bold')
text(-0.2, 2, 'Relative pressure [a.u.]','fontsize',14, 'fontweight','bold','rotation',90)
set(gca,'xticklabel',[' 0';'10';'20';'30']);
print -djpeg Fig5-13c

figure(4); clf;
Phi = zeros(length(f),length(x));
arg5 = exp(x)-1;

for n=1:10
    Phi(n,:)=-1.12*0.001*f(n)*arg5;
    end

for n=2:2:10
    plot(x, Phi(n,:), 'k-','linewidth',1)
    hold on
    s=sprintf('%4d Hz',f(n));
    for m=[5 10 15 20 25 30]
        plot(x(m), Phi(n,m), 'ko','markersize',6)      
    end
    axis([0 3.5 -4*pi 0.2])
        if n==2 text(2.8,Phi(n,28)+0.3,s,'fontsize',11,'fontweight','bold'); end
        if n==4 text(2.8,Phi(n,28), s,'fontsize',11,'fontweight','bold'); end
        if n==6 text(2.0,Phi(n,20), s,'fontsize',11,'fontweight','bold'); end
        if n==8 text(1.0,Phi(n,10), s,'fontsize',11,'fontweight','bold'); end
        if n==10 text(0.5,Phi(n,5), s,'fontsize',11,'fontweight','bold'); end
    set(gca,'xticklabel',[]);
    set(gca, 'xtick', [0 1 2 3]);
    set(gca,'fontsize',13, 'fontweight','bold') 
end
xlabel('Distance from stapes (mm)','fontsize',14, 'fontweight','bold')
ylabel('Relative phase (rad)','fontsize',14, 'fontweight','bold','rotation',90)
set(gca,'xticklabel',[' 0';'10';'20';'30']);
print -djpeg Fig5-13d



