trapz_dke_yp

C/C++ source

0001 /*=================================================================
0002  *
0003  * trapz_dke_yp.C  .MEX file corresponding to trapz_dke_yp.M
0004  *
0005  * Numerical integration by the trapeze's method of a function y = f(x).
0006  *Steps may be irregular. (see the .m equivalent trapz_dke_yp.m)
0007  *
0008  *    Input:
0009  *
0010  *        - f: [x,f1,f2,...,fm], x is the abciss vector, fi are m functions [n,m+1]
0011  *
0012  *    Output:
0013  *
0014  *        - sf: [sf1,sf2,...,sfm], sum of fi's [n,m] 
0015  *
0016  *
0017  *by Y.PEYSSON CEA-DRFC (yves.peysson@cea.fr) and Joan Decker CEA-DRFC (joan.decker@cea.fr)
0018  *
0019  * This is a MEX-file for DKE.  
0020  *
0021  *=================================================================
0022  */
0023 #include <math.h>
0024 #include <stdlib.h>
0025 #include "mex.h"
0026 
0027 /* Input Arguments */
0028 
0029 #define    F_IN    prhs[0]
0030 
0031 /* Output Arguments */
0032 
0033 #define    SF_OUT    plhs[0]
0034 
0035 #if !defined(MAX)
0036 #define    MAX(A, B)    ((A) > (B) ? (A) : (B))
0037 #endif
0038 
0039 #if !defined(MIN)
0040 #define    MIN(A, B)    ((A) < (B) ? (A) : (B))
0041 #endif
0042 
0043 void trapz(sf,f,m,n)
0044 double    sf[];
0045 double    f[];
0046 unsigned int m;
0047 unsigned int n;        
0048 {
0049 unsigned int i;
0050 unsigned int j;
0051 unsigned int k;
0052 double    *delta;
0053 double    dd;
0054 
0055     delta = mxMalloc(m*sizeof(double)+1);
0056 
0057     if (m > 2) {
0058           for (i=0;ifor (j=1;jfor (i=0;ielse if (m == 2) {
0068        for (j=1;jelse if (m == 1) {
0072        for (j=1;jvoid mexFunction( int nlhs, mxArray *plhs[], 
0085           int nrhs, const mxArray*prhs[] )
0086      
0087 { 
0088     double *sf; 
0089     double *f; 
0090     unsigned int m,n; 
0091     
0092     /* Check for proper number of arguments */
0093     
0094     if (nrhs != 1) { 
0095     mexErrMsgTxt("One input argument required."); 
0096     } else if (nlhs > 1) {
0097     mexErrMsgTxt("Too many output arguments."); 
0098     } 
0099     
0100     /* Check the dimensions of F */
0101     
0102     m = mxGetM(F_IN);/*number of rows*/ 
0103     n = mxGetN(F_IN);/*number of columns*/
0104     if (!mxIsDouble(F_IN)) {
0105     mexErrMsgTxt("trapz_dke_yp requires double type input data"); 
0106     } 
0107     
0108     if (n==0 || m==0) {
0109     mexErrMsgTxt("Empty matrix in trapz_dke_yp"); 
0110     } 
0111     
0112     /* Create a matrix for the return argument */ 
0113     SF_OUT = mxCreateDoubleMatrix(1, n-1, mxREAL); 
0114     
0115     /* Assign pointers to the various parameters */ 
0116     sf = mxGetPr(SF_OUT);
0117     f = mxGetPr(F_IN); 
0118            
0119     /* Do the actual computations in a subroutine */
0120     
0121     trapz(sf,f,m,n); 
0122     return;
0123     
0124 }
0125 
0126