function ff = stnormfx(xx)
% stnormfx takes a value from a std. normal distr., xx, 
% and returns its cumulative distribution function value
% (Abramowitz and Stegun 1964).
   ci = 0.196854; % these are approximation constants
   cii = 0.115194;
   ciii = 0.000344;
   civ = 0.019527;

   if xx >= 0; z= xx; else z = -xx; end;
   a = 1 + (ci*z) + (cii*z*z);
   b = (ciii*z*z*z) + (civ*z*z*z*z);
   w = a+b;
   if xx >= 0   
       ff= 1 - 1/(2*w*w*w*w);
   else   
		 ff = 1 - (1 - 1/(2*w*w*w*w)); 
   end;