
function [P,dx] = peaks_hc(x, plotflag);
%
% 	PEAKS_hc.M  computes the Green's function peak matrix of the human 
%       cochlea. Data is loaded from MAN_DATA.MAT.
%
%		[P, dx] = PEAKS_hc(x);
%
% 	x  =  BM point vector.
%	plotflag = If set to 1, plots are showed (def. 0)
% 
%	dx = Interpoint distances
%	P  = Peak matrix
%	
%	Matrix P is not normalised (see PSIGMA.M) and must be multiplied 
%	by b(x)*b(x'), where b(x) is the BM width. The dimension of P
%	is  1/BETA  in BETA units (1 BETA = 33.5 mm).
%
%	Files MAN_DATA.MAT, BM_Width.mat, Cochlear_radius.m  and Psigma.m are called.
	


if nargin <2,
	plotflag=0;
end

xp=0;	% initialisation of arrays to be loaded from MAN_DATA.MAT 
yp=0;
xm=0;
ym=0;

disp('	Function PEAKS is at work');
BM_Width=0;	% initialisation of basilar membrane width array
N=length(x);
x=x(:);		% To make sure that x is a column
c=[x(1);x(2:N)+x(1:N-1)]/2;	% averaged coordinates

load BM_Width.mat -ascii	% Loading BM width data. This yields 2-column matrix 'bmw'

b=cochlear_radius(BM_Width(:,1),BM_Width(:,2),x); 	% Interpolating BM-width data	

b=b/2; 		% The effective BM width is 1/2 the geometrical width		

load man_data 	% Loading cochlear radius. This yields xp, yp, xm, ym  

rp=cochlear_radius(xp, yp, c); % Interpolating upper cochlear radius

rm=cochlear_radius(xm, ym, c); % Interpolating lower cochlear radius


r = (rp+rm)/2;   % Computing averaged coclear radius at point c(i) 

% ---------------------------------------------------------


xo=[0;x(1:N-1)];
dx=x-xo;	% interval lengths

% The following routine computes peaks without reflection effect at the stapes

P=zeros(N,N);

for j=1:N,
	P(:,j)=psigma(x-c(j),b(j),r(j))-psigma(xo-c(j),b(j),r(j));
end

%	P(:,j)  is column 'j'; 

if plotflag==1,
	figure(1)
	clf
	plot(x, P(:,1:20:N) ),
	title('BM-BM Green''s function peaks'),
	pause(3)
end	
	
