void _nlsfIntegral(
// Fit Parameter(s):
double mu, double sigma,
// Independent Variable(s):
double x,
// Dependent Variable(s):
double& y)
{
	// Beginning of editable part
		double dIntegral = 0.0;
		double dInt = 0.0;
		double t = 0.0; // the integration interval ranges from 0.0 to infinity
		double dPrecision = 1e-6; // the precision to control the integraton to infinity 
		double dt = 1.0; // this is the step size for the integration. 
		//the following lines actually perform a integrate to the function, the trapezoidal rule is used.
		double F1, F2, B1, B2;
		if (mu<0) 
		{
			y = NANUM;
			return;
		}
		do
		{
			if(t == 0.0) 
			{
				B1 = 0.0; // lim   B(t,mu,sigma) = 0;
				          // t->0
			}
			else B1 = exp( -0.5*(ln(t/mu))^2/sigma/sigma ) / sqrt(2.0*PI)/t;
			
			B2 = exp( -0.5*(ln((t+dt)/mu))^2/sigma/sigma ) / sqrt(2.0*PI)/(t+dt);
			F1 = sin(x+t); // F(x) = sin(x+t);
			F2 = sin(x+t+dt);
						
		    dInt = 0.5 * (B1*F1 + B2*F2) * dt ; 		
			dIntegral += dInt;
			t +=dt;
		}while( dInt/dIntegral > dPrecision );
		y=dIntegral;
	// End of editable part
}

#include <oc_nag8.h>

struct user
{
double b, tau, fitX;  // fitX the independent variable of fitting function
};

static double NAG_CALL f_callback(double x, Nag_User *comm)  // x is the independent variable of the integrand
{
 
struct user *sp = (struct user *)(comm->p);
 
double b, tau, fitX; // temp variable to accept the parameters in the Nag_User communication struct
b = sp->b;
tau = sp->tau;
fitX = sp->fitX;
 
return 1/sqrt(x)*exp(-b/x-x/tau)*cos(2.21e5*fitX*x);
}

void _nlsffitint(
// Fit Parameter(s):
double P, double tau, double rhoN, double wN, double tN, double L,
// Independent Variable(s):
double x,
// Dependent Variable(s):
double& y)
{
	// Beginning of editable part
	double a = 5.2532E-6*sqrt(rhoN)/(wN*tN);
	double b = 7.209E8*L^2*rhoN;
	
	double epsabs = 0.00001, epsrel = 1e-3, result, abserr;
	Integer max_num_subint = 100; 
	
	Nag_QuadProgress qp;
	static NagError fail;
	
	Nag_User comm;	
	struct user s;
	s.b = b;
	s.tau = tau;
	s.fitX = x;
	comm.p = (Pointer)&s;
	    
	d01smc(f_callback, Nag_UpperSemiInfinite, 0, epsabs, epsrel, max_num_subint, &result, &abserr, &qp, &comm, &fail);
	
	if (fail.code != NE_INT_ARG_LT && fail.code != NE_BAD_PARAM && fail.code != NE_ALLOC_FAIL)
	{
		NAG_FREE(qp.sub_int_beg_pts);
		NAG_FREE(qp.sub_int_end_pts);
		NAG_FREE(qp.sub_int_result);
		NAG_FREE(qp.sub_int_error);
	}
	
	y = P^2*a*result;
	// End of editable part
}