/* Gaussian.m  	by Paul Kunz	January 1993
 * Subclass of PFunction to manage a Gaussian function
 *
 * Copyright (C)  1993  The Board of Trustees of
 * The Leland Stanford Junior University. All Rights Reserved.
 */

#import "Gaussian.h"

const char Gaussian_h_rcsid[] = GAUSSIAN_H_ID;
const char Gaussian_m_rcsid[] = "$Id: Gaussian.m,v 1.2 1993/01/31 20:17:38 pfkeb Exp $";

#import <objc/Storage.h>
#include <math.h>

#define SQR2P	2.50662828


@implementation Gaussian 

static double gaussian( double x, double binW, double *par )
{
    double 	norm 	= par[0];
    double 	xzero	= par[1];
    double 	sigma	= par[2];
    double 	value;
    double	t;
    
    if ( sigma != 0.0 ) {
        t = (x - xzero) / sigma ;
	t = 0.5 * t*t;
	if ( fabs(t) > 50.0 ) {
	    value = 0.0;
	} else {
	    value = exp( -t ) / (SQR2P * sigma );
	}
    } else {
        if ( x != xzero ) {
	    value = 0.0;
	} else {
	    value = 1.0;
	}
    }

    return binW * norm * value;
}

- init
{
    [super init];
    [[[self setTitle:"Gaussian"] setFunctionPtr:gaussian] setNumberArgs:3];
    [self registerFunc];	/* must be called after the above */
    [self setArgName:"Norm" at:0];
    [self setArgName:"Mean" at:1];
    [self setArgName:"Sigma" at:2];
    return self;
}

- setInitialValues
{
    fitParm             *parm;
    graphtype_t		type;
    
    [super setInitialValues];
    if ( !disp ) return self;
    
    type = h_getDispType(disp);
    
    parm = [variedParms elementAt:0];
    if ( type == HISTOGRAM ) {
	parm->value       = disp->bins.totals[1][0];
	parm->lower_limit = 0.0;
	parm->upper_limit = 2.0 * parm->value; 
	parm->step_size   = 0.01 * parm->upper_limit;
    } else {
        parm->lower_limit = disp->yAxis.low;
        parm->upper_limit = 4.0 * disp->yAxis.high;
        parm->value       = 0.5 * parm->upper_limit;
	parm->step_size   = 0.01 * (disp->yAxis.high - disp->yAxis.low);
    }
    
    parm = [variedParms elementAt:1];
    parm->lower_limit = disp->xAxis.low;
    parm->upper_limit = disp->xAxis.high;
    parm->step_size   = 0.01 * (disp->xAxis.high - disp->xAxis.low);
    parm->value       = parm->lower_limit 
                        + 0.50 * (disp->xAxis.high - disp->xAxis.low);
    
    parm = [variedParms elementAt:2];
    parm->upper_limit = 0.2 * (disp->xAxis.high - disp->xAxis.low);
    parm->lower_limit = 0.02 * parm->upper_limit;
    parm->step_size   = 0.01 * (parm->upper_limit - parm->lower_limit);
    parm->value = 0.5 * (parm->upper_limit - parm->lower_limit );
    return self;
}

@end