function [tw, TMsite, TMindex] = twave(x,G,Gs,Sh,lambda,Hz,flag,Csh)
%
%	Computes Travelling Waves with TM resonating near CF
%
%	[TW, TM_SITE, TM_INDEX] = TWAVE(X,G,Gs,Sh,Lambda,Hz,Flag,Csh)
%
%	X	= BM point vector in [0, 1] (can be non uniformely spaced);
%	G, Gs	= stapes propagator and BM Green's function computed by GREENF.M
%	Sh	= Sectional shear viscosity operator computed by GREENF.M
%	Lambda	= outer hair cell undamping control parameter
%	Hz	= stimulus frequency in Hertz
% 	Flag	= 0 for const. displ. , 1 for const. velocity at the stapes 
%	Csh	= adimensional viscosity factor (relative to water)
% 	TM_SITE = BM site at which TM resonates
% 	TM_INDEX = index at which TM resonates
%	Values for mass is provided by function Corti_mass(x), all terms are  
%	computed as averages over i-th interval dx(i)=x(i)-x(i-1), with xo=0.
%        
	
			

t0=clock;
disp('	Function TWAVE is at work ...');
omega = 2*pi*Hz;
flag=-flag;
omp=omega^flag;
om2=omega^2;
iom= sqrt(-1)*omega;

N = length(x);

x=x(:)'; 				% make sure x is a row.

dx=[x(1),x(2:N)-x(1:N-1)];		% interpoint distances


disp('	Computing mass ...');

m = Corti_mass(x);	% Unit dimension is Kg/BETA
u = ones(size(x));

disp('	Computing stiffness ...')

%---------------------------------------------------------------------------
%  Default values of parameters
     expk = -3.6*log(10); 	% Stiffness exponent  (Base E)
     ko  = 2e+3;		% Kg/(BETA*sec^2)
% ------------------------------------------------------------------

k = ko*exp(expk*x);	 % Stiffness vector in Kg/(sec^2*BETA)

stif=[k(1)-ko, k(2:N)-k(1:N-1)]./(expk*dx); % Averaging stiffness locally
					    % by integration over dx:

stif=stif/om2; % Stiffness divided by square frequency (unit dimension = Kg/BETA)

% ------------------------------------------------------------------
disp('	Computing damping ...')

[damp, damp0] = damping(x);
% Here damp and damp0 have unit dimension Kg/(BETA*sec)

damp  = -damp/iom;
damp0 = -damp0/iom; % Now their dimension is Kg/BETA
	
OCSh = (Csh/iom)*Sh;	% Organ of Corti shear viscosity term
			% Csh is adimensional.
			% OCSh has unit dimension Kg/BETA

if lambda==0,
	TMsite=-1;	% to signal that TM is uneffective
    TMindex=-1;
	disp('	Computing kernel ...')
	Kernel = G + OCSh + diag(m-stif-damp);
else,
	disp('	Computing OHC undamping term')	
	[res, TMsite, TMindex]=tmres(x, omega);
	lam=exp(+0.025*x + 0.15*(x-0.35).^2);
	% lambda values are corrected to uniformise output profiles 
	undamp=(lambda*damp0).*lam.*res;
	% undamp=(lambda*damp0).*res;
	disp('	Computing kernel ...')
	Kernel = G + OCSh + diag(m-stif-damp-undamp);
end


disp('	Computing solution ...');
Ampl = 1e-3; % Stapes oscillation amplitude (arbitrary length unit)

tw=((-omp*Gs*Ampl)/Kernel);
% tw amplitude has the same unit dimension as Ampl

tw=tw(:);		% To make sure the output is a column 

disp(['	CPU time taken by TWAVE: ' num2str(etime(clock,t0)) ' sec.'])