%A simple model of the glacial cycles:
%
%   function [To V]=Hmodel(t,threshold,Vo,mu,sig)
%
% Input:
%  t, time in ky
%  threshold, the threshold value across which terminations initiate
%  Vo, initial amount of icevolume
%  sig, the standard deviation of the random component of ice
%       volume chnange.
%
% Output:
%  To, the time of terminations
%  V,  the time hisotry of ice volume
%
%Peter Huybers, 2011
%phuybers@fas.harvard.edu

function [V T]=Hmodel(t,threshold,Vo,mu,sig);

if nargin<4,
  mu=1;
  sig=0;
end;

a=mu;

%---------
%run model
V=Vo;
Tb=[];
Te=[];
for ct=2:length(t); 
  V(ct)=V(ct-1)+a+sig*randn; 
  if V(ct)>threshold(ct),
    a=-V(ct)/10;
  end;
  if V(ct)<=0,
    a=mu;
  end;
end;

%-----------------------------------------------------------
%determine the timing of termination from rates of change in
%ice-volume, as done witht the actual delta18O record.
dV=diff(smoothPH(V,11)); 
pl=find(dV(2:end-1)<-5 & dV(2:end-1)<dV(1:end-2) & dV(2:end-1)<=dV(3:end));
pl=pl+1;
T=t(pl);