function f=phi(sigma);
sigma=sigma.*(sigma<0);
f=0.5*(-sigma+sinh(sigma));
