function z = igrnd(N, theta,c)
%Inverse Gaussian generator where mean of the IG is theta>0 and the diffusion term is c>0. 
% First passage time to a barrier at theta for process dX(t)=dt+cdW(t)
% N = number of simulated values
one1 = ones(1,N);
v = chi2rnd(1,1, N);
z1 = theta*one1+ c^2*v/2  - 0.5*c^2*sqrt(4*theta*v/c^2+ v.^2);
U = rand(1,N);
H = (U <= theta./(theta+z1));
z = z1.*H + (theta^2./z1).*(one1 - H);  % z  has the inverse gaussian distribution