What this does is the following.

In matlab you write a function that you want to integrate. It can be anything. Because it is a matlab function, any parameters in it must be either hard coded or global variables set elsewhere.

You call the mexfunction (NRIntergrate.cpp) to integrate the function between the limits (You also need a copy of Function.h)

I call the function like this

integral = NRIntegrate('Gaussian','Adaptive',lower,upper,tol) ;

Where Gaussian is the name of my function (Gaussian.m)

The gaussians parameters,its mean, sigma and size are stored internally in that function.

If you look in the code you can see how I use a functor for the function which internally accesses the matlab function.

Hope this is useful

Will Handler

NRIntegrate.cpp

#include <nr3matlab.h>
#include <quadrature.h>
#include <interp_1d.h>
#include <romberg.h>
#include <adapt.h>
#include "Function.h"

/*
Input to use various numerical recipes techniques from Matlab, we pass
in a function to integrate

INPUTS: Function, 1 dimensional
Method, which can be one of the following

Trapz
Simpson
Romberg
Adaptive
OpenRomberg


Lowerlimit
UpperLImit
Tolerance
We also
AUTHOR W. Handler August 2008
*/


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
//////////////////////////////////////////////////////////////////////
//Test the inputs and make sure the all make sense, maximum number of
//inputs has to be equal to 2
//
//This is most of the work, testing the inputs and putting them into
//arrays, the actual call to NR code is simple.
//////////////////////////////////////////////////////////////////////



if ( nrhs != 5 ) {
cerr << "The inputs are not correct!\n The correct format to use this function is\n... NRIntegrate( function, Method, lower limit, upper limit, tolerance ) \n";
return;
}

//now process the func name
if ( !mxIsChar( prhs[0] ) ){
mexErrMsgTxt("The function name needs to be given as a string");
}
const char *funcname = (char*) mxArrayToString(prhs[0]);

//mow process the method
if ( !mxIsChar( prhs[1] ) ){
mexErrMsgTxt("The method needs to be given as a string");
}
const char *method = (char*) mxArrayToString(prhs[1]);


//now process the limits
//LowerLimit
if ( !mxIsNumeric(prhs[2]) || mxIsComplex(prhs[2]) ) {
mexErrMsgTxt("The lower limit has to be a real numeric scalar");
}
int nd = mxGetNumberOfDimensions(prhs[2]);
int *size =(int*) mxGetDimensions(prhs[2]);
if ( nd != 2 )
mexErrMsgTxt("The lower limit has too many or too few dimensions, it should be a vector in nmatlab");

if ( !( size[0] == 1 && size[1] == 1) )
mexErrMsgTxt("The lower limit has to be a scalar");
Doub LowerLimit = *((double*)mxGetData(prhs[2]));

//UpperLImit
if ( !mxIsNumeric(prhs[3]) || mxIsComplex(prhs[3]) ) {
mexErrMsgTxt("The lower limit has to be a real numeric scalar");
}
nd = mxGetNumberOfDimensions(prhs[3]);
size =(int*) mxGetDimensions(prhs[3]);
if ( nd != 2 )
mexErrMsgTxt("The UpperLimit has too many or too few dimensions, it should be a vector in nmatlab");

if ( !( size[0] == 1 && size[1] == 1) )
mexErrMsgTxt("The Upper LImit has to be a scalar");
Doub UpperLimit = *((double*)mxGetData(prhs[3]));


//Tolerance
if ( !mxIsNumeric(prhs[4]) || mxIsComplex(prhs[4]) ) {
mexErrMsgTxt("The tolerance has to be a real numeric scalar");
}
nd = mxGetNumberOfDimensions(prhs[4]);
size =(int*) mxGetDimensions(prhs[4]);
if ( nd != 2 )
mexErrMsgTxt("The tolerance has too many or too few dimensions, it should be a scalar in nmatlab");

if ( !( size[0] == 1 && size[1] == 1) )
mexErrMsgTxt("The tolerance has to be a scalar");
Doub TOL = *((double*)mxGetData(prhs[4]));



if ( LowerLimit > UpperLimit) mexErrMsgTxt("Limits do not make sense\n");

//mexPrintf("Doing an integration of %s between %lf and %lf using %s with a tolerance of %e\n", funcname, LowerLimit, UpperLimit, method, TOL);

//Perform the integration
Function func((const char*) funcname);

Doub integral, oldintegral = 0.0;

if ( !strcmp("Simpson", method) ){
integral = qsimp(func, LowerLimit, UpperLimit, TOL);
} else if ( !strcmp("Trapz", method) ){
integral = qtrap(func,LowerLimit, UpperLimit,TOL);
} else if ( !strcmp("Romberg", method) ) {
integral = qromb(func,LowerLimit, UpperLimit,TOL);
} else if ( !strcmp("Adaptive", method) ) {
Adapt s(TOL);
integral = s.integrate(func,LowerLimit, UpperLimit);
} else if ( !strcmp("OpenRomberg", method) ) {
Midpnt<Function> q(func,LowerLimit,UpperLimit);
integral = qromo(q,TOL);
} else {
mexErrMsgTxt("Method for integration not recognised");
}



VecDoub output(1, plhs[0]);
output.resize(1);
output[0] = integral;

}




Function.h
#ifndef Function
#define Function__h


/*
Functor to a matlab function for use with numrec code
Will Handler: July 2008

*/
class Function
{
public:
//FitFunction::~FitFunction(){
// fSinglePoint->mxDestroyArray();
//} //
Function(const char* funcname); //constructor
// ~FitFunction(); //destructor
Doub operator()(Doub Point);
//{
//*input = Point;
//mexCallMATLAB(1,&fOutput,1,&fInput,fName);
//output = ( (double*) mxGetData(fOutput));
//return *output;
//}
//this overloaded operator is used to make the class a functor for use by
//some of the fitting routines


void SetFuncName(const char *aFuncName);

private:
char * fName; //name of function to call in Matlab space

mxArray* fInput; //array of mxArrays to be used with mxCallMATLAB
mxArray* fOutput;
double *fFuncVals;//pointer to calculated function values
double *input;
double *output;

};


Function::Function(const char *funcname)
{
//CONSTRUCTOR FOR FUNCTION CLASS
//creat the arrays to use, need to be destroyed in hte destructor
fInput = mxCreateDoubleMatrix(1,1, mxREAL);
fOutput = mxCreateDoubleMatrix(1,1, mxREAL);
this->SetFuncName(funcname);
input = (double*) mxGetPr(fInput);
output= (double*) mxGetData(fOutput);
}

Doub Function::operator()(Doub Point)
{
//call the matlab fnction by name and returns its value, it should hgave only one input.

*input = Point;
mexCallMATLAB(1,&fOutput,1,&fInput,fName);
output = (double*) mxGetData(fOutput);
return *output;
}

void Function::SetFuncName(const char *aName){
fName =(char*) aName;
}

#endif