function t_inh = timewarp(t_hom)
% TIMEWARP generates a set of spike times of an inhomogeneous
%          Poisson process with rate given by the function lambda_t.m
  
% Input: 
% t_hom: vector of firing times (spike train) of a homogeneous Poisson process
%  
% Output
% t_inh: vector of firing times (spike train) of an inhomogeneous Poisson
%        process
%  
% Parameters:
lambda_0 = 100; %rate of the homogeneous Poisson process in t_hom
r = 100;  % mean firing rate of the inh. Poisson process (units of Hz)
f = 3.3; % oscillation frequency of the inh. Poisson process (units of Hz)
DT =  0.002; % integration time step


% body of function  starts here

N = length(t_hom); % number of spikes

time = 0; %  initialization of time, ...
n=1;      %  spike index, and ...
I = 0;    %  integral

for n=1:N     % go through all spikes
  while I < t_hom(n)   % integrate until a spike has been reached
    delta_I = DT *lambda_t(r,f,time)/lambda_0;
    I = I + delta_I; % update integral (of course there are much
		     % better integration routines ...)
    time = time + DT;% update time 
  end  
  t_inh(n) = time;                      % zero order approximation

  delta = DT * (t_hom(n) - I)/delta_I;  % first order correction    
  t_inh(n) = t_inh(n) + delta ;                    
end % for n=1:N