function [res, xmin, imax]  = tmres(x, om, plotflag);

%	FUNCTION [RES, IMIN, IMAX] = TMRES(X, OM, PLOTFLAG);
% 
%	X = Points on the BM
%	OM = Input radian frequency
%	PLOTFLAG = Plot flag (if 1, plot tectorial membrane resonance profile)
%
%	RES = TM resonance transfer function (=1 at OM) 
%	XMIN = X value at which abs(RES)=1; 
%	IMAX = index at which abs(RES) = 1; 



if nargin<3,
	plotflag=0;
end	

	micron = 29.85e-6;	% 1 micron in BETA
	tmh0 = 15*micron;	% 15 micron,  
		   		% putative height of TM at base
		   		% (effective for shear viscosity) 

	tmh = tmh0*exp(log(4)*x);	% TM mean height 
					% (grows 4 times from base to apex)

	tmw0 = 70*micron;	  % 80 micron,  
				  % putative width of the TM at base
				  % (effective for shear viscosity) 

	tmw = tmw0*exp(log(4)*x);	% TM mean width 
					% (grows 4 times from base to apex)

	tmvsc = 4e-3; % Kg/(BETA sec) % viscosity of TM at base

	tmvsc=tmvsc*exp(log(4)*x);       %  grows 4 times

	rho = 0.0376;	% = 1000/(29.85)^3 Kg/BETA^3, density of water  

	tm_mass=rho*(tmh.*tmw);     % TM mass (constant) per unit length

	omb=2*pi*18000; 	    % TM resonance frequency at the base
	oma=2*pi*60;		    % TM resonance frequency at the apex	
	omexp = log(oma/omb);
	wm=omb*exp(omexp*x);
	

	bigamma=tmvsc./tm_mass;   % the TM-RL and TM-LIMBUS damping constant

	% figure(1)
	% plot(x,bigamma);
	% pause

	wm2=wm.*wm;
	
	om2=om*om;
	iom=om*sqrt(-1);

	res = -iom*bigamma./((wm2 - om2) +iom*bigamma); 
		
	[val, imax]=max(abs(res));
	
	
	xmin=x(imax);

if plotflag==1,
	% damp=damping(x);
	% undamp=damp.*res;
	figure(1)
	subplot(211), plot(x, abs(res), xmin, 0,'*c', 'linewidth', 5);
	title('TM transfer function modulus');
	ylabel('Amplitude'); 
	subplot(212), plot(x, angle(iom*res), xmin, -pi/2, '*c', 'linewidth', 5);
	title('TM transfer function phase');
	ylabel('Phase [rad]')
	xlabel('Fractional distance from stapes');
end	