const double sqrtln2 = 0.832554611;
	double	T[] = {0.0,	.314240376,	.947788391,	1.59768264,	

2.27950708,	3.02063703,	3.8897249};
	double	C[] = {0.0,	1.01172805,	-.75197147,	1.2557727e-2,	

1.00220082e-2,	-2.42068135e-4,	5.00848061e-7};
	double	S[] = {0.0,	1.39323,	.231152406,	-.155351466,	

6.21836624e-3,	9.19082986e-5,	-6.27525958e-7};
	double	wr = 0.0;
	double	wi = 0.0;
	double	X, Y, y1, y2, y3, d, d1, d2, d3, d4, r, r2;
	int	i;

	double	temp;

	X = 2*sqrtln2*(x - x0)/wG;
	Y = sqrtln2*wL/wG;

	y1 = Y + 1.5;
	y2 = y1 * y1;

	if(Y > 0.85 || fabs(X) < 18.1*Y+1.65)
	{
		for(i=1; i<=6; i++)
		{
			r = X - T[i];
			d = 1.0/(r*r + y2);
        	d1 = y1*d;
        	d2 = r*d;
			r = X + T[i];
        	d = 1.0/(r*r + y2);
        	d3 = y1*d;
        	d4 = r*d;
        	wr = wr + C[i]*(d1 + d3) - S[i]*(d2 - d4);
        	wi = wi + C[i]*(d2 + d4) + S[i]*(d1 - d3);
    	}
    }
    else
    {
		if(fabs(X) < 12.0)
		{
			wr = exp(-X*X);
		}
		y3 = Y + 3.0;
		for(i=1; i<=6; i++)
		{
			r = X - T[i];
			r2 = r * r;
			d = 1.0/(r2 + y2);
			d1 = y1 * d;
			d2 = r * d;
			wr = wr + Y*(C[i]*(r*d2 - 1.5*d1) + S[i]*y3*d2)/(r2 + 2.25);
			r = X + T[i];
			r2 = r * r;
			d = 1.0/(r2 + y2);
			d3 = y1 * d;
			d4 = r * d;
			wr = wr + Y*(C[i]*(r*d4 - 1.5*d3) - S[i]*y3*d4)/(r2 + 2.25);
			wi = wi + C[i]*(d2 + d4) + S[i]*(d1 - d3);
		}
    }
    temp = A*wr*0.939437278/wG;

y = y0 + temp * temp;

sort(x_y_curve);
smooth(x_y_curve);
y0 = yatxmin( x_y_curve );
wG = 0.85 * fwhm( x_y_curve, y0 );
wL = wG;
A = area( x_y_curve, y0 );
if ( A > 0 ) x0 = xatymax( x_y_curve );
else x0 = xatymin( x_y_curve );