%PUPILC  integrates the delay-differential equation 
%for the pupil light reflex with Ornstein-Uhlenbeck noise on 
%the parameter c. Program averages results over navg realizations. 
%The program uses the integral Euler-Maruyama algorithm proposed in
%R.F. Fox et al., Physical Review A 38, 5938 (1988). 
%Solution can be seen using plot(t,u);
%Histogram of solution can be seen using plot(rhox,rhoy);
%
%   Copyright 2003
%   Centre for Nonlinear Dynamics in Physiology and Medicine

tau=0.3; 
%delay
tot=20000; 
%total number of integration time steps per realizations
navgs=10;  
%total number of realizations
N=100; 
%subdivide each delay interval by this number 
dt=tau/N; 
%actual integration time step
trans=5; 
%number of delays to discard as transients
n=10.0;    
%bifurcation parameter
alf=3.21; 
% alpha
the=50.0; 
%theta
k=0.0; 
% parameter k
c=200; 
% parameter c 
dnz=1.0; 
%intensity of the Gaussian white noise process
tcor=0.001; 
%correlation time of the Ornstein-Uhlenbeck process
xinit=0.0; 
%initial condition
t=ones(1,tot+N); 
%time vector
bins=200;  
%number of bins to use for histogram of solution 
rhoy=zeros(1,bins); 
%histogram of the solution
transtime=trans*N+1; 
%number of integration time steps discarded as transients
efac=exp(-dt/tcor);
gausfac=sqrt(dnz*(1-efac*efac)/tcor);
for realiznumber=1:navgs,
    realiznumber
    randn('state',sum(100*clock));
    eta=sqrt(dnz/tcor)*randn; %initial condition for Ornstein-Uhlenbeck process
for i=1:N+1,
    t(1,i)=-tau+(i-1)*dt;
end;    
u=xinit*ones(1,tot+N);
for i=N+2:tot+N,
    t(1,i)=-tau+(i-1)*dt;
    gausx=randn*gausfac;
    etah=eta*efac+gausx;
    fb=(c+etah)/(1+(u(1,i-1-N)/the)^n)+k;
    u(1,i)=u(1,i-1)+dt*(fb-alf*u(1,i-1));
    eta=etah;
end;
rhoy=rhoy+hist(u(transtime:tot+N),bins);
end;
[rhodum,rhox]=hist(u(transtime:tot+N),bins);

