function [G0S, G0M] = Greenwood(x, plotflag)

%	Computes the coarse propagators over BM site array X for a 
%	realistic profile of the cochlear duct.

%		[G0S, G0M] = G0(X, PlotFlag)
%
%	X   = BM point vector (points may be irregularly spaced!)
%	PlotFlag = if set to 1, plots are showed [def.0]
%
%	G0S = stapes propagator; 
%   G0M = basilar membrane propagator.
%	
%	G0S and G0M are the result of the integration over
%	the inverse of the coclear duct section surfaces and are
%	NOT multiplied by the BM form factors. This latter step
%	is accomplished by function GREENF.M. The dimension of 
%	G0S and G0M is  1/BETA  in BETA units (1 BETA = 33.5 mm).



if nargin < 2,
	plotflag=0;
end

% initialising the coordinates of cochlear semiduct radii: 
xp=0;
xm=0;
yp=0;
ym=0;

disp('	Reading human cochlea data...');

x=x(:)'; 	% To make sure that x is a row
N=length(x);

load man_data	% File contains cochlear duct profile data 
		% All lengths are scaled. 1 unit == 33.5 mm
		% Data from the upper side:  xp, yp 
		% Data from the lower side:  xm, ym
		% ap, am = cross--sectional areas
		% "p" is for upper profile, "m" is for lower.


Np=max(size(xp))-1;
Nm=max(size(xm))-1;     


ypm= (yp(1:Np)+yp(2:Np+1))/2;   % middle point values (number = Np)	
Sp = (pi/2)*ypm.^2;		% Find cross section areas at middle points
				% from human cochlea data, assuming cochlear
				% duct cross sections are two half circles
ymm= (ym(1:Nm)+ym(2:Nm+1))/2;
Sm = (pi/2)*ymm.^2;	 

%%%%%%%%%%%%%%%%%%%% COMMENT THIS TO AVOID PLOTTING %%%%%%%%%%%%%%%%
% if plotflag==1,
%	figure(1)
%	plot(xp,yp,xm,-ym)
%	xlabel('Fractional distance from stapes');
%	ylabel('Semiduct radii in BM length units');
%	title('Human cochlea profile'),
%	pause(5)
%end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% disp('	Computing G0 ...');

dxp = xp(2:Np+1)-xp(1:Np);	% Find interval lengths
dxm = xm(2:Nm+1)-xm(1:Nm);

Intp=zeros(size(1:Np+1));		% Initialise
Intm=zeros(size(1:Nm+1));

dx_over_Sp = dxp./Sp;		% Find integrand
dx_over_Sm = dxm./Sm;

Intp(Np+1) = (0.5*pi)*yp(Np+1)/Sp(Np);	% Contribution from upper side of helicotrema;
Intp(Nm+1) = (0.5*pi)*ym(Nm+1)/Sm(Nm); 	% Contribution from lower side of helicotrema;

for i=Np:-1:1,
	Intp(i) = Intp(i+1)+dx_over_Sp(i);
end

for i=Nm:-1:1,
	Intm(i) = Intm(i+1)+dx_over_Sm(i);
end

% Interpolation of Intp and Intm: 


G0S=cochlear_radius(xp(1:Np+1),Intp,x) + cochlear_radius(xm(1:Nm+1),Intm,x);  

G0M=G0S'*ones(size(G0S));
G0M=tril(G0M) + tril(G0M)'- diag(G0S);

%%%%%%%%%%%%%%%%%%%%%% COMMENT THIS TO AVOID PLOTTING %%%%%%%%%%%%%%%
if plotflag==1,
	figure(1)
	plot(x', G0M(:,(1:10:N)),'.')
	title('Unshaped coarse BM-BM Greenwood''s functions'), 
	pause(5)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


disp('	Loading BM width data ...')
load BM_Width.mat -ascii

%%%%%%%%%%%%%%%%%%%%% COMMENT THIS TO AVOID PLOTTING %%%%%%%%%%%%%%
if plotflag==1,
	figure(2)
	clf,
	axis([0 18.5 0 .2]/33.5);
	plot(BM_Width(:,1), BM_Width(:,2),BM_Width(:,1),BM_Width(:,2),'.c', 'linewidth', 5),
	title('Basilar membrane width');
	xlabel('Fractional distance from stapes');
	ylabel('BM width in BM length units');
	pause(5)
	axis;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b=cochlear_radius(BM_Width(:,1),BM_Width(:,2),x);	% b is a column
G0=G0M.*(b'*b);
if plotflag==1,
	figure(3)
	plot(x',G0(:,1:10:N)), title('Coarse BM-BM Greenwood''s functions')
end