function [Gs, G, Sh] = greenf(x, plotflag);
%
%	FUNCTION [Gs, G, Sh] = GREENF(X, PlotFlag);    
%	
%	GREENF.m  calls  Greenwood.m,  GPEAKS.M and BM_Width.MAT to produce hydrodynamic
%	propagators (Greenwood's functions) Gs, G, the organ of Corti shear operator Sh.
%	
%	X  = BM point vector covering interval 0-1. 
%	PlotFlag = if set to 1, Green's functions plots are showed [def=0]
% 
%	Gs = Stapes-to-BM  Green's function (force generated by stapes acceleration
%	     and transmitted by fluid pressure field on a unit length BM segment).
%	G  = BM-to-BM Green's functions (force generated by the acceleration of
%	     a BM segment on the other BM segments, left-multiplied by integration
%	     differentials dx).
%	Sh = sectional shear viscosity operator [= d../dx s(x) d../dx, where
%	     s(x) is the sectional shear viscosity of the organ of Corti assumed
%	     as viscous as water (needs multiplication by relative viscosity)]. 	
%
%	Unequally spaced X elements will be accepted also. 	
%	The unit dimension of Gs and G is Kg/BETA, that of Sh is 
%	Kg/(BETA*sec). BETA = 33.5 mm for the human cochlea.
%					


% The density of water is 10^3 Kg/m^3. But lengths are expressed as
% fractions of the cochlear duct length (33.5 mm). Let's call BETA
% the unit length in our system. As 1 BETA = 33.5 x 10^{-3} meters
% then 1 m = 29.85 BETA. So we can express the density of water 
% in our system as:
 		
if nargin < 2
	plotflag=0;
end

rho = 0.0376;	% = 1000/(29.85)^3 Kg/BETA^3 

disp('	Building Green''s function peaks ...')

time=clock;

[P, dx] = peaks_hc(x);	% dx is a column representing the differentials of
			% the fractional distance from stapes x, for general
			% spacing.  
			% P = unshaped singular part of Green's function
			% times integration differential factors dx. 
			% P is computed as the differential after x of a 
			% symmetric sigmoidal function (see peaks_hc.m)

disp(['	CPU time taken by PEAKS: ', num2str(etime(clock,time)),' sec']);

disp('	Computing coarse Green''s functions ...');

[Gs, G]=Greenwood(x);		% unshaped coarse profile of the Green's function 

BM_Width=0;		% array initialisation: bmw_man(:,1)=abscissae,
			% bmw_man(:,2)=ordinates
 
load BM_Width.mat -ascii	% Loading BM-width data from ascii *.mat file

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

b=b/2;	% the effective BM width is assumed to be 1/2 of the geometric width

b=b(:);

u=ones(size(dx));

disp('	Computing Green''s functions ...');

G=G.*(dx*u')+ P;    % matrix  dx*u' is the Kronecker product of unit
		    % vector u by the integration differential array dx
		    % providing G with integration differential factors
		    % (already accounted for in P). Here G is adimensional 


G=G.*(b*(rho*b'));  % Matrix form factor on account of the BM width and 
		    % water density. Unit dimension of G is Kg/BETA
		     
S=pi*(2/33.5)^2;    % Stapes surface in BETA^2

Gs=(rho*S*b').*Gs;		

disp('	Computing sectional shear viscosity operator ...');

N=length(dx);
invdx=u./dx;		% inverse differentials

Der=diag(invdx,0);
Derl = diag(invdx(1:N-1),1) - Der; % First derivative at left.
				   % Unit dimension = 1/BETA
Derr = diag(invdx(2:N),-1) - Der;  % First derivative at right.
				   % Unit dimension = 1/BETA

x=x(:);

% --------------- Shear viscosity computations -----------------------------

micron = 29.85e-6; % 1 micron in BETA
h0 = 40*micron;	   % 40 micron,  
		   % putative height of the organ of Corti at base
		   % (effective for shear viscosity) 

h = h0*exp(log(4)*x);	% effective organ of Corti thickness
			% (grows 4 times from base to apex)

w0 = 50*micron;	  % 50 micron, mean radial width of the organ of Corti at base
		  % (effective for shear viscosity) 

w = w0*exp(log(4)*x); % effective organ of Corti radial width
		      % (grows 4 times from base to apex)

eta = 1e-3/29.85; %= 3.35e-5 Kg/(BETA*sec) 
		  % 10^{-3} Kg/(m*sec) = viscosity coefficient of water
		  % at 20 Celsius degrees (very close to 1 centipoise)	

s = eta*w.*h; 	  % Sectional shear viscosity constant, assuming that
		  % adjacent sectional segments of the organ of Corti
		  % are viscous like water. Unit dimension of s is Kg*BETA/sec
% ----------------------------------------------------------------------------	

Sh=Derl*diag(s)*Derr; % sectional shear viscosity operator
		      % Unit dimension of Sh is Kg/(BETA*sec) 

%%%%%%%%%%%%%  PLOTS %%%%%%%%%%%%%%%%%%%
if plotflag==1,
	figure(1)
	subplot(211), 
	plot(x, G(1:20:N, :),'-r', 'linewidth', 1),
	title('BM-BM Greenwood''s functions (sampled)'),
	subplot(212),
	plot(x, Gs,'-r', 'linewidth', 1),
	title('Stapes-BM Greenwood''s function')
	xlabel('Fractional distance from stapes');
end