%
%         a Hopf bifurcation for driven outer hair cells mechanics
%
%     in Hudspeth et al. 2010 formulated as:   
%       dz/dt = -(C-Cc-iwC)*z(t) - (a+ib)|z|quit^2*z + F(t)/G
%          critical when C=Cc. Oscillates when C-Cc = -10
%           near criticality: C-Cc = 0.1 Then: compressing nonlinearity....
%
%    Note that (see Izhkevich book pp 210):
%
%        the system 
%          dz/dt = (c(b) + iw(b))z + (a+id)z |z|^2    (w(b) = frequency)
%       is equivalent to the polar coordinate formulation:
%
%           dr/dt = c(b)r + ar^3
%           dphi/dt = w(b) + dr^2
%
%      and the cartesian system
%
%           du/dt = c(b)u - w(b)v + (au-dv)(u^2+v^2)
%           dv/dt = w(b)u + c(b)v + (du+av)(u^2+v^2)
%
function hopf_ohc2(dC, a, b, w0, F0)

if nargin<5 A = 0.2; end
if nargin<4 w0 = 80*2*pi; end
if nargin<3 b = 1; end
if nargin<2 a = 1; end
if nargin<1 dC = 0.1; end     % near criticality

dt=0.0001;
t = (0:dt:5.0);    % time in 0.1 ms intervals, i.e. up to 5 kHz
x=zeros(1,length(t));
y=x;
x(1) = 0.02;   % initial condition
y(1) = 0.02;
z=complex(x,y);
p=complex(a,b);

F0 = zeros(1,length(t));
F0(1000 : 2500) = A;    % 150 ms of sound stimulus
F=F0.*cos(w0*t);   % the driving frequency

figure(1)
clf
plot(x(1), y(1),'ko','markersize',8);
hold on
plot(0, 0,'r+','markersize',8);

for n=2:length(t)-1
    z2 = (abs(z(n-1)))^2;                      % |z|^2
    z(n) = z(n-1) - (dC - 1i*w0)*z(n-1)*dt - p * z2 *z(n-1)*dt + F(n-1)*dt;   % complex plane
end
z(end)=z(end-1);
zamp = abs(z(end));

plot(z, 'r-','linewidth', 1.5);     % plot in the complex plane
axis([-1.2*zamp 1.2*zamp -1.2*zamp 1.2*zamp]);
axis square

figure(2); clf;
plot(t, real(z), 'r-', 'linewidth',2)
hold on
axis([0 1 -1.2*zamp 1.2*zamp]);
%plot(t, imag(z), 'k--', 'linewidth',1)

x=real(z);
figure(3); clf    % spectrum of real(z)
ftx=fft(x(40000:end),4*4096);   ftx=ftx(1:4*2048);   % steady state ...
ftxa=abs(ftx); 
df = 5000/(4*2048);
freq=[0:df:5000-df];
loglog(freq, ftxa,'r-', 'linewidth',2);

hold on
axis([0.1 2000 1 1.2*max(ftxa)]);
ftx=fft(x,4*4096);   ftx=ftx(1:4*2048);   % the whole signal
ftxa=abs(ftx); 
loglog(freq, ftxa,'b--')

end