/***************************************************************
 * jHepWork.
 * Original Developer
 *     Dale Visser (dale@visser.name)
 *     S.Chekanov
 *     Based on nuclear reactor simulator software (open source)
 **************************************************************/

package jhpro.tseries;

import org.apache.commons.math3.util.FastMath;

/**
 * Gaussian filter for time series analysis.
 * 
 * @author S.Chekanov.
 * 
 */

public class GaussianFilter {

	/**
	 * magic number for calculating
	 */
	static private final double a = 0.93911;

	/**
	 * magic number for calculating
	 */
	static private final double b = 2.77066;

	/**
	 * magic number for calculating
	 */
	static final double SIGMA_TO_FWHM = 2.354;

	private double[] inSpectrum, outSpectrum;
	private double width;

	/**
	 * Initialialize the fitter
	 * @param inSpectrum input spectrum.
	 * @param width width
	 */
	public GaussianFilter(double[] inSpectrum, double width) {
		this.inSpectrum = inSpectrum;
		this.width = width;
		processSpectrum();
	}

	/**
	 * Process the spectrum.
	 */
	private void processSpectrum() {
		int length = inSpectrum.length;
		outSpectrum = new double[length];
		for (int newChannel = 0; newChannel < length; newChannel++) {
			double newCounts = 0.0;
			double scaleFactor = a * inSpectrum[newChannel] / width;
			int minCh = (int) FastMath.round(FastMath.max(0, newChannel - 3 * width));
			int maxCh = (int) FastMath.round(FastMath.min(length - 1, newChannel + 3
					* width));
			for (int oldChannel = minCh; oldChannel <= maxCh; oldChannel++) {
				newCounts += FastMath.exp(-b
						* FastMath.pow((oldChannel - newChannel) / width, 2));
			}
			newCounts *= scaleFactor;
			outSpectrum[newChannel] = newCounts;
		}
	}

	/**
	 * Get filtered result.
	 * @return
	 */
	public double[] getFilteredSpectrum() {
		return outSpectrum;
	}

}