function  [cochlea]=coch_input_V4(Nfreq, azi, nstimuli, Is, Freq_flag);

% nonlinear, saturating, azimuth dependent activity for both cochleae:
% the threshold is put at 90 deg left and right 
% 
% This function calculates the activity values of the cochlea layers (left and right) which are presented in a forward
% direction to the audmap layer
% function that generates Gaussian receptive field profile:
% F = exp( -[(x-x0)^2]/[2sigma^2] )

if nargin < 5
   Freq_flag=0;
end
if nargin < 4
  Is=1.0;
end  
if nargin < 3
  nstimuli=1;;
end
if nargin < 2
  azi=40;
end
if nargin < 1
  Nfreq = 30;
end

global Gain_freq; 
global Exp_freq; 
global Sig_freq;  
global Freq_min;
global Freq_max;

Gain_freq = 1.0;
Exp_freq  = 40.0;
Sig_freq  = 4;
Freq_min=100;
Freq_max=10000;

 cochlea=zeros(2,Nfreq);
 for k=1:1:nstimuli
	 
	 %/* determine cochlear location in frequency map */
	 io = randi(Nfreq,1,1);
	 m = (Nfreq - 1) / (log(Freq_max) - log(Freq_min));
	 logf =  ((io - 1.0)/m) + log(Freq_min);
	 cf =  exp(logf);
	 
	 
	 for i=1:1:Nfreq
		 % changed dependency blocking is freqency dependent        %
		 % High freq much diminution                                %
		 % Low frequency less diminution    ((float)Exp_freq+i*10)  %
		 IL=Is-exp( +(1/(Exp_freq + i*0))*(azi - 90));
		 if(IL<0) IL=0.0; end;
		 IR=Is-exp( -(1/(Exp_freq + i*0))*(azi + 90));
		 if(IR<0) IR=0.0; end;
		
		 rate_r = Gain_freq  * IR * exp( -((i-io)*(i-io)) / (2*(Sig_freq*Sig_freq)) ); 
		 if (rate_r<IR) rate_L = IR; end;
		 cochlea(1,i) =  cochlea(1,i)+ rate_r;
		 
		 rate_l = Gain_freq  * IL * exp( -((i-io)*(i-io)) / (2*(Sig_freq*Sig_freq)) );
		 if (rate_l<IL) rate_L = IL; end;
		 cochlea(2,i) = cochlea(2,i) + rate_l;
		 	 
	 end
	 
 end
% figure(1)
% subplot(221)
%      hold off
%      plot((cochlea(1,:)),'k');
% 	 hold on
% 	 plot((cochlea(2,:)),'r');
end %% CLOSE FUNCTION