1.12.1 Mathematics

Mathematics

Normalize

The following example shows how to pick a point in a data plot (curve) and then normalize all curves in the layer to the same value at that point. This snippet of code assumes a graph layer with multiple curves is active and all curves share the same X values. This assumption is typical in spectroscopy.

GraphLayer gl = Project.ActiveLayer();
if( !gl )
    return;

// Allow user to click and select one particular point of one particular curve
GetGraphPoints mypts;
mypts.SetFollowData(true);
mypts.GetPoints(1, gl);

vector vx, vy;
vector<int> vn;
if(mypts.GetData(vx, vy, vn) == 1)
{
    // Save index and y value of picked point
    int nxpicked = vn[0] - 1;
    double dypicked = vy[0];
    
    // Loop over all data plots in layer
    foreach( DataPlot dp in gl.DataPlots )
    {
        // Get the data range and then the y column for current plot
        XYRange xy;
        Column cy;
        if(dp.GetDataRange(xy) && xy.GetYColumn(cy))
        {
            // Get a vector reference to y values from the y column                
            vectorbase &vycurrent = cy.GetDataObject();
            
            // Scale vector so y value matches user-picked point
            vycurrent *= dypicked/vycurrent[nxpicked];
        }
    }
}

InterpolationExtrapolationInterpolation/Extrapolation

The ocmath_interpolate function is used to do interpolation/extrapolation with modes of Linear, Spline and B-Spline.

// Make sure there are 4 columns in active worksheet
// The first two columns are source xy data, 
// 3rd column has input x data, 4th column to put output y.
Worksheet    wks = Project.ActiveLayer();    
wks.SetSize(-1, 4); 

DataRange drSource;
drSource.Add(wks, 0, "X"); // 1st column - source x data
drSource.Add(wks, 1, "Y"); // 2nd column - source y data

vector vSrcx, vSrcy;
drSource.GetData(&vSrcx, 0);
drSource.GetData(&vSrcy, 1);

DataRange drOut;
drOut.Add(wks, 2, "X"); // 3rd column - input x data
drOut.Add(wks, 3, "Y"); // 4th column - interpolated y data

vector vOutx, vOuty;
drOut.GetData(&vOutx, 0);

int    nSrcSize = vSrcx.GetSize();
int    nOutSize = vOutx.GetSize();
vOuty.SetSize(nOutSize);

int nMode = INTERP_TYPE_BSPLINE;
double dSmoothingFactor = 1;
int iRet = ocmath_interpolate(vOutx, vOuty, nOutSize, vSrcx, vSrcy, nSrcSize, 
nMode, dSmoothingFactor);

drOut.SetData(&vOuty, &vOutx);

Integration

Origin C provides access to NAG's integral routines to perform integration. With Origin C and NAG you can do integration on a normal integrand, an integrand with parameters, an integrand with oscillation, an infinite integral, higher dimension integration, and more. The following examples show how to do integration with NAG.

Your Origin C code will need to include the NAG header file at least once before your code calls any NAG functions.

#include <OC_nag.h> // NAG declarations

Simple Integral Function

Integral, with One Integration Variable The first example shows how to do a basic integration on a simple integrand with only one integration variable.

// NAG_CALL denotes proper calling convention. You may treat it 
// like a function pointer and define your own integrand
double NAG_CALL func(double x,Nag_User *comm)
{
        int *use_comm = (int *)comm->p;
        return (x*sin(x*30.0)/sqrt(1.0-x*x/(PI*PI*4.0)));
}

void nag_d01sjc_ex()
{
	double a = 0.0;
	double b = PI * 2.0;  // integration interval
	
	double epsabs, abserr, epsrel, result;
	// you may use epsabs and epsrel and this quantity to enhance your desired 
	// precision when not enough precision encountered
	epsabs = 0.0;
	epsrel = 0.0001;
	
	// The max number of sub-intervals needed to evaluate the function in the 
	// integral. For most cases 200 to 500 is adequate and recommmended.
	int max_num_subint = 200;
	
	Nag_QuadProgress qp;
	NagError fail;   
	
        Nag_User comm;
	static int use_comm[1] = {1};
	comm.p = (Pointer)&use_comm;
        d01sjc(func, a, b, epsabs, epsrel, max_num_subint, &result, &abserr, 
	&qp, &comm, &fail);
	
	
	// For the error other than the following three errors which are due to 
	// bad input parameters or allocation failure. You will need to free
	// the memory allocation before calling the integration routine again
	// to avoid memory leakage
	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);
	}
	
	printf("%g\n", result); 
}

Integral Function with Parameters

Integral, with Multiple Parameters The next example shows how to define and perform integration on an integrand with parameters. Notice that the parameters are passed to the integrator by a user-defined structure. This avoids having to use static variables as parameters of the integrand, and makes it thread-safe.

This example can also be adapted to use NAG's infinite integrator. For instance, by enabling the line calling the infinite integrator d01smc function, the example can be used to perform infinite integration.

struct user // integrand parameters
{
    double A;
    double Xc;
    double W; 
};

// Function supplied by user, return the value of the integrand at a given x.
static double NAG_CALL f_callback(double x, Nag_User *comm) 
{
    struct user *param = (struct user *)(comm->p);
	
    return param->A * exp(-2 * (x - param->Xc) * (x - param->Xc)
        / param->W / param->W) / (param->W * sqrt(PI / 2));
}

Now, we set parameter values for the function and define the additional parameters necessary to perform the integration. The integration is then performed by a single function call, passing the parameters as arguments.

void nag_d01sjc_ex()
{
    double a = 0.0;
    double b = 2.0; // integration interval
    
    // The following variables are used to control
    // the accuracy and precision of the integration.
    double epsabs = 0.0;      // absolute accuracy, set negative to use relative
    double epsrel = 0.0001;   // relative accuracy, set negative to use absolute
    int max_num_subint = 200; // max sub-intervals, 200 to 500 is recommended
    
    // Result keeps the approximate integral value returned by the algorithm
    // abserr is an estimate of the error which should be an upper bound 
    // for |I - result| where I is the integral value
    double result, abserr;
    
    // The structure of type Nag_QuadProgress, it contains pointers 
    // allocated memory internally with max_num_subint elements
    Nag_QuadProgress qp;
    
    // The NAG error parameter (structure)
    NagError fail;
    
    // Parameters passed to integrand by NAG user communication struct
    struct user param;
    param.A  = 1.0;
    param.Xc = 0.0;
    param.W  = 1.0;

    Nag_User comm;
    comm.p = (Pointer)&param;
    
    // Perform integration
    // There are 3 kinds of infinite boundary types you can use in Nag infinite 
    // integrator Nag_LowerSemiInfinite, Nag_UpperSemiInfinite, Nag_Infinite
    /*
    d01smc(f_callback, Nag_LowerSemiInfinite, b, epsabs, epsrel, max_num_subint, 
    &result, &abserr, &qp, &comm, &fail);
    */
    d01sjc(f_callback, a, b, epsabs, epsrel, max_num_subint, 
    &result, &abserr, &qp, &comm, &fail);
        
    // check the error by printing out error message
    if (fail.code != NE_NOERROR)
        printf("%s\n", fail.message);
    
    // For errors other than the following three errors which are due to 
    // bad input parameters, or allocation failure,
    // you will need to free the memory allocation before calling the 
    // integration routine again to avoid memory leakage.
    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);
    }
    
    printf("%g\n", result);
}

Multi-dimension Integral Function

Integral, Multi-dimension For integrals of dimension higher than 2, you can call the NAG integrator function d01wcc to perform the integration.

Our user defined call back function will be passed to the NAG d01wcc function.

double NAG_CALL f_callback(int n, double* z, Nag_User *comm)
{ 
	double tmp_pwr;
	tmp_pwr = z[1]+1.0+z[3];
	return z[0]*4.0*z[2]*z[2]*exp(z[0]*2.0*z[2])/(tmp_pwr*tmp_pwr);
}

Main function:

void nag_d01wcc_ex()
{	
	// Input variables
	int ndim = NDIM;  // the integral dimension
	double a[4], b[4];
	for(int ii=0; ii < 4; ++ii)  // integration interval
	{
		a[ii] = 0.0;
		b[ii] = 1.0;
	}
	int minpts = 0;	
	int maxpts = MAXPTS;  // maximum number of function evaluation 
	double eps = 0.0001; // set the precision
	
	// Output variable
	double finval, acc;
	Nag_User comm;	
	NagError fail;	
	
	d01wcc(ndim, f_callback, a, b, &minpts, maxpts, eps, &finval, &acc, 
	&comm, &fail);
	
	if (fail.code != NE_NOERROR)
		printf("%s\n", fail.message);
	
	if (fail.code == NE_NOERROR || fail.code == NE_QUAD_MAX_INTEGRAND_EVAL)
	{
		printf("Requested accuracy =%12.2e\n", eps);
		printf("Estimated value    =%12.4f\n", finval);  
		printf("Estimated accuracy =%12.2e\n", acc);
	}
}

Differentiation

The ocmath_derivative function is used to do simple derivative calculations without smoothing. The function is declared in ocmath.h as shown below.

int ocmath_derivative(
    const double* pXData, double* pYData, uint nSize, DWORD dwCntrl = 0);

The function ignores all missing values and computes the derivative by taking the average of the two slopes between the point and each of its neighboring data points. If the dwCntrl argument uses the default value of 0, the function fills in the average when data changes direction.

if( OE_NOERROR == ocmath_derivative(vx, vy, vx.GetSize()) )
    out_str("successfully");

If dwCntrl is set to DERV_PEAK_AS_ZERO, the function fills in zero if data changes direction.

if( OE_NOERROR == ocmath_derivative(vx, vy, vx.GetSize(), DERV_PEAK_AS_ZERO) )
    out_str("successfully");