double dU = 500.0; // the upper endpoint of the interval
	double dL = 0.0;	// the lower endpoint of the interval
	double z;
	int nSteps = 100; // the integrating interval will be broken into nsteps for calculation. this will affect the fitting speed and precision.
	double dz = (dU-dL)/(double)nSteps;
	double dFz;	  // F(z);
	double dFzz;  // F(z+dz);
	double dInteg = 0.0;
	for (int ii=0; ii<nSteps; ii++)
	{
		if(ii==0) 
			dFz = exp(-z)/((P1+P2*x)*(P1+P2*x)+(z-P3*x)*(z-P3*x));
		else
			dFz = dFzz;
		z += dz;
		dFzz = exp(-z)/((P1+P2*x)*(P1+P2*x)+(z-P3*x)*(z-P3*x));
		dInteg += 0.5*(dFz+dFzz)*dz;
	}
	y = P4*P4*dInteg/(P4*P4+x*x);