a mex file for converting p coordinate into isentropic coordinate. improve the performance by about 500 times

p2isen.c

#include "mex.h"
#include "matrix.h"

void   mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
	double *vcomp3dInterp, *bin3d; /* pointers to the input matrixces */
	double *flux2dIsen;            /* pointer to output matrix */
        double *numIndex;              /* num of index array is third input */
 
	mwSignedIndex lon,lat,p,numIsen;   /* matrix dimensions */
        mwSize  i,j,k;                     /* use in for loop   */

	vcomp3dInterp = mxGetPr(prhs[0]);
	bin3d         = mxGetPr(prhs[1]);
	numIndex      = mxGetPr(prhs[2]);

        lon     = numIndex[0];
	lat     = numIndex[1];
	p       = numIndex[2];
	numIsen = numIndex[3];

	plhs[0]       = mxCreateDoubleMatrix(lat, numIsen, mxREAL);
	flux2dIsen    = mxGetPr(plhs[0]);

        for (i = 1; i <= lon; i++)
	{ for (j = 1; j <= lat; j++)
	  {  for (k=1; k <= p;  k++)
	    {
		if (bin3d[i+lon*(j-1)+lon*lat*(k-1)-1] > 0)
		{
		   flux2dIsen[j+lat*((int)bin3d[i+lon*(j-1)+lon*lat*(k-1)-1]-1)-1] = 
                      flux2dIsen[j+lat*((int)bin3d[i+lon*(j-1)+lon*lat*(k-1)-1]-1)-1] + 
                      vcomp3dInterp[i+lon*(j-1)+lon*lat*(k-1)-1];
                }}}}

}