//
// gaussian.cc
//
// C++ program for normal computations
// based on the algorithm 26.2.17 in
// Abromowitz and Stegun, Handbook of Mathematical
// Functions, Dover, NY
//
#include <math.h>
#include "gaussian.h"
distdens Gaussian(double quantity,double mean,double variance) {
  distdens toreturn;
  double q,sigma,t;
  sigma = sqrt(variance);
  q = (quantity - mean)/sigma;
  toreturn.dens = exp(-q*q/2.0)*.3989422804014327/sigma;
  if(q>0.0) {
    t = 1.0/(1.0+0.2316419*q);
    toreturn.dist = 1.0-exp(-q*q/2.0)*.3989422804014327*
      t * (0.319381530 +
      t * (-0.356563782 +
      t * (1.781477937 +
      t * (-1.821255978 +
      t * 1.330274429))));}
  else {
    t = 1.0/(1.0-0.2316419*q);
    toreturn.dist = exp(-q*q/2.0)*.3989422804014327*
      t * (0.319381530 +
      t * (-0.356563782 +
      t * (1.781477937 +
      t * (-1.821255978 +
      t * 1.330274429))));}
  return(toreturn);}