function [lat,coeff] = latencyfit_jp(y,t,fig)

%TIME IN MILLISECONDS!
%Must be col vectors

% ----------------------------------------------
% clear all;
% 
% generate some example data
% coeff = [0.004    0.11    0.07    80      14]; % typical values of V1 modulation
% coeff = [0.04    0.2     1.3     50     7]; % typical values of V1 PSTH
% 
% a = coeff(1);
% c = coeff(2);
% d = coeff(3);
% m = coeff(4);
% s = coeff(5); 
% 
% t = [-300:1:600]';
% y = d * exp(m*a+0.5*s^2*a^2-a.*t).*normcdf(t,m+s^2*a,s)+c.*normcdf(t,m,s);
% 
% y = y+0.1*randn(size(y)); % add some noise
% 
% figure;plot(t,y,'b');legend({'data'});

% make educated guess for starting values (instead of these you can also
% use a list of starting values, for example based on previous data,
% and loop over these and pick the one with the best fit)
[Y,I] = max(y);
m     = t(I);
st_=[0.004 Y/2 Y/2 m 200];
% use curve tracing toolbox
ft_ = fittype('d * exp(m*a+0.5*s^2*a^2-a.*t).*normcdf(t,m+s^2*a,s)+c.*normcdf(t,m,s);',...
    'dependent',{'r'},'independent',{'t'},...
    'coefficients',{'a', 'c', 'd', 'm', 's'});
fo_ = fitoptions('method','NonlinearLeastSquares','MaxFunEvals',1000,'Lower',[0,0,0,0,0]);
set(fo_,'Startpoint',st_);

try
    [cfun,gof,output] = fit(t,y,ft_,fo_);
catch
    %Fit went wrong
    lat = NaN;
    coeff = NaN(1,5);
    return
end
[coeff]=coeffvalues(cfun);

if 0 % use optimization toolbox instead
    % r = d * exp(m*a+0.5*s^2*a^2-a.*t).*normcdf(t,m+s^2*a,s)+c.*normcdf(t,m,s);
    predicted = @(a,tdata) a(3) * exp(a(4)*a(1)+0.5*a(5)^2*a(1)^2-a(1).*t).*normcdf(t,a(4)+a(5)^2*a(1),a(5))+a(2).*normcdf(t,a(4),a(5));
    options = optimset('MaxFunEvals',2000);    
    x0=st_;
    lb=[0,0,0,0,0];
    ub=[Inf,Inf,Inf,Inf,Inf];
    [coeff]=lsqcurvefit(predicted,x0,t,y,lb,ub,options);
end

% coefficients of fit
a = coeff(1);
c = coeff(2);
d = coeff(3);
m = coeff(4);
s = coeff(5);
yf = d * exp(m*a+0.5*s^2*a^2-a.*t).*normcdf(t,m+s^2*a,s)+c.*normcdf(t,m,s);

[mf,mfix] = max(yf);
ixs = [1:mfix]; % only search before max peak
[l33,l33ix]=min(abs(yf(ixs)-mf*0.33));

lat = t(l33ix);

if fig
    figure; hold on;
    plot(t,y,'b',t,yf,'r');
    legend({'data','fit'});
    plot([t(l33ix),t(l33ix)],get(gca,'YLim')); % point where 33% of max response reached
end