/*********************************************************************
 * This file is part of the PRAPI library.
 *
 * Copyright (C) 2001-2002 Topi Mäenpää
 * All rights reserved.
 *
 * This program is free software. You can redistribute and/or modify
 * it under the terms of the free software licence found in the
 * accompanying file "COPYING". The licence terms must always be
 * redistributed with this source file. The above copyright notice
 * must be reproduced in all modified and unmodified copies of this
 * source file.
 *
 * $Revision: 1.5 $
 *********************************************************************/

#include "Gaussian.h"

using namespace util;

namespace prapi
{
	Matrix<double> Gaussian::create2DGaussian(double radius, int size, double width)
	{
		Matrix<double> result(size, false);
		double step = radius*2/(size-1);

		for (int i=0;i<size;i++)
			for (int j=0;j<size;j++)
				{
					double x=-radius+i*step;
					double y=-radius+j*step;
					result(i,j) = exp(-x*x-y*y)/width;
				}
		return result;
	}

	Matrix<double> Gaussian::create2DNormal(double radius, int size, double deviation)
	{
		Matrix<double> result(size, false);
		double step = radius*2/(size-1);
		double twoD2 = deviation*deviation*2;
		double normalizer = step*step/(twoD2*M_PI);

		for (int i=0;i<size;i++)
			for (int j=0;j<size;j++)
				{
					double x=-radius+i*step;
					double y=-radius+j*step;
					result(i,j) = normalizer * exp(-(x*x+y*y)/twoD2);
				}
		return result;
	}

	List<double> Gaussian::create1DGaussian(double range, int size, double width)
	{
		List<double> result(size);
		double step = range*2/(size-1);

		for (int i=0;i<size;i++)
			{
				double x=-range+i*step;
				result += exp(-x*x)/width;
			}
		return result;
	}

	double Gaussian::erf(double x)
	{
		int xSign = 1;
		if (x < 0)
			{
				x = -x;
				xSign = -1;
			}

		if (x > 5) return xSign;
		int sign = 1;
		unsigned int n = 0;
		double factorial = 1, power = x, xSquared = x*x;
		double sum = 0, increment = 0;
		do
			{
				if (n)
					factorial *= n;
				increment = power / (factorial*((n<<1)+1));
				sum += increment*sign;
				sign = -sign;
				n++;
				power *= xSquared;
			} while (n <= 170 && increment > 1e-12);

		return xSign * M_2_SQRTPI * sum;
	}
}