 
 function g=gamma(x,m,s,c,a)

g1=s/sqrt(2*pi)*exp(-(x-m-c/a)^2/(2*s^2))+(m+c/a-x)*(1-normcdf(x,m+c/a,s));
g=a*x-g1;
