function [TW, X] = tw(lambda, Hz, N, Csh, vflag, spacexp)
%
%	Function TW computes Travelling Waves for a set of frequencies. 
%
%       	[TW, X] = TW(Lambda, Hz, N, Csh, Vflag, SpaceExp)
%
%	Lambda = OHC undamping control parameter (0 < Lambda < 1.12) [def. 0]
%	Hz = input frequency vector [def. 100*[1, 2, 4, 8, 16, 32, 64, 128, 200]] 
%	N  = number of oscillating points on the basilar membrane (BM) [def. 300]
%	Csh = shear-damping control coefficient [def.25]
%	Vflag = 0 for constant stapes displacement, 1 for constant stapes velocity [def. 0]
%	SpaceExp = BM point spacing increase exponentially with this exponent [def. 1.5]	
%
%	MARKS:	(*) = location of tectorial membrane resonance
%		(o) = location of travelling wave peak  
% 
%	The function returns the result in matrix TW; the i-th COLUMN of
%	which is the travelling wave elicited by a tone of frequency Hz(i)
%	input to stapes.
%	  The columns of matrix X contain the corresponding x-axis points used 
%	for computing solutions.
%	  The following auxiliary routines are needed:
%	TWAVE.M, GREENF.M, PSIGMA.M, GPEAKS.M, G0.M, INTERPOL.M,
%	MAN_DATA.MAT, BMW_MAN.MAT, VARGRID.M, FR2X.M, X2FR.M,  CF.MAT 
%
%				MATHEMATICAL NOTE
%
%	Function TW calls function TWAVE to solve for tw the linear equation
%
%		[-w²(G+M) + iw(H + S) + K + U(y)] tw = -w² Gs,
%	
%	where: w = radian frequency of stapes input; G = Green's function
%	matrix; M, H and K = organ of Corti mass, absolute viscosity and
%	stiffness resp. (diagonale matrices); S = shear viscosity operator  
%	(non diagonal matrix); U(y) = OHC undamping force term, depending
%	on the tectorial membrane displacement y; Gs = stapes force times
%	stapes  Green's function. At each BM point x, U(y) depends on w as:
%
%			-iw lambda h(x)/[1 + f(x,w)],
%
%	where h(x) is a diagonal element of H and f(x,w) vanishes if x is
%	close to the BM location having w as characteristic radian frequency. 
%


% DEFAULT SETTING 
if nargin < 6
	spacexp=1.5;
end
if nargin < 5,
	vflag=0;
end
if nargin < 4,
	Csh=25;
end
if nargin < 3,
	N=300;
end
if nargin < 2,
	Hz = 100*[1, 2, 4, 8, 16, 32, 64, 128, 200]; 
end

if nargin < 1,
	lambda=0;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if(length(lambda) > 1), 
	disp('parameter lambda in TW must be a scalar');
	return;
end

if(min(size(Hz)) > 1),
	disp('parameter Hz in TW cannot be a matrix');
	return;
end

if N<30,
	disp('parameter N in TW is too small: irregular behaviour expected')
end

disp(' ')
disp(' ')
disp('******************* TRAVELLING WAVES AT THEIR BEST! ******************');
disp(' ')


NC=length(Hz);	% Number of Travelling Waves to be computed.

x=vargrid(N, spacexp);
x=x(:);

N=length(x);

disp(['   Number of points: ' int2str(N)]);

TW =zeros(N,NC);	% Allocate space for solutions
X  =zeros(N,NC);

%%%%%%%%%%%%%%%%%%%% Initialisations for progressive settings %%%%%%%%%%%
TMx	 = [];	% tectorial membrane resonance sites
TWx	 = [];	% tectorial membrane resonance levels
CFsites	 = [];	% Characteristic frequency sites
CFphases = [];	% Phase values at characteristic frequency
TMphases = [];
TMsites = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a = zeros(N,NC);
phas = zeros(N,NC);

u=ones(size(x));



%%%%%%%%%%%%%%%  Solve problem NC times %%%%%%%%%%%%%%%%%%%%%%

for i=1:NC,
	disp(' ');
	disp(['	Stimulus frequency: ' num2str(Hz(i)) ' Hz']);
	% compute propagators, but only if necessary
	
	X(:,i)=x;		
	if i==1,
		G=0;
		if exist('PROP.mat')==2,
			disp('	Loading from PROP.mat');
			load PROP.mat,
		end
		if exist('Xg') == 1,
			test = ceil(1/Xg(1))-ceil(1/x(1));
		else
			test=1;
            Xg = 1;   %*something*. So it can be saved.
		end
		if (size(G) ~= [N,N]) & test ~=0,
			[Gs, G, Sh]=greenf(x);
			save PROP.mat Gs G Sh Xg
		end
	end
	
	%%%%%%%%%%%%%%%%%%%%%%% SOLVE EQUATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%
	[TW(:,i), TMmax, TMindex] = twave(x,G,Gs,Sh,lambda,Hz(i),vflag,Csh);
	a(:,i)=abs(TW(:,i));
	mxa=0.9*max(max(a));
	if TMmax >=0,
		TMx=[TMx,TMmax];
		TWx=[TWx, mxa];
	end

	phas(:,i)=unwrap(angle(TW(:,i)))/(2*pi);	

	%%%%%%%%%%%%%%%%%%%  PLOTS  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
	if i==1,
		figure(1); 
		clf;
	end;
	subplot(211)
	plot(X, a, '-r', TMx, TWx, '*b'),

	if lambda==0,
		if i==1, 
	    		title(['Passive Travelling Wave at ', num2str(Hz(1)), 'Hz']);
	  	else,
			if vflag==0,
		  		title('Passive Travelling Waves (constant stapes displacement)');
			else,
			  	title('Passive Travelling Waves (constant stapes velocity)');
			end
		end
	else,
		if i==1, 
	    		title(['Active Travelling Wave at ', num2str(Hz(1)), ' Hz - Undamping = ' num2str(lambda)]);
		else,
			title(['Active Travelling Waves - Undamping = ', num2str(lambda)]);
 		end
	end

	if vflag==0,	
		ylabel('TW-Ampl./Stapes-displ.'),
	else
		ylabel('TW-Ampl./Stapes-vel.'),
	end

	[ma,ja]=max(a(:,1:i));	% values and indices of maximum TW aplitudes up to now

	for j=1:i,
		mj=ma(j);
		xj=X(ja(j),j);
		text(xj,mj,[num2str(Hz(j))]);
	end

	CFphases = [CFphases, phas(ja(i),i)];
	CFsites = [CFsites, xj];

	if TMmax > 0,
		TMphases = [TMphases, phas(TMindex,i)];
		TMsites = [TMsites, TMmax];
	end
	subplot(212), 
	if vflag==1,
		if TMmax <=0, 
			plot(X, phas+pi/2, '-w', CFsites, CFphases, 'ob');
		else
  			plot(X, phas+pi/2, '-w', CFsites, CFphases, 'ob',...
			 	TMsites, TMphases, '*r');
		end
	else
		if TMmax <=0, 
			plot(X, phas, '-w', CFsites, CFphases, 'ob');
		else
  			plot(X, phas, '-w', CFsites, CFphases, 'ob', ...
				TMsites, TMphases, '*r');
		end
	end
	 
	xlabel('Fractional distance from stapes'), 
	ylabel('Phase [cycles]'),
	pause(1)
end % end cycle


