yprime.c

MATLAB扩展编程by 孙永康.ZIP.有需要的朋友请来下载.

#include <math.h>
#include "mex.h"

/* 输入变量 */

#define	T_IN	prhs[0]
#define	Y_IN	prhs[1]


/* 输出变量 */

#define	YP_OUT	plhs[0]

#if !defined(MAX)
#define	MAX(A, B)	((A) > (B) ? (A) : (B))
#endif

#if !defined(MIN)
#define	MIN(A, B)	((A) < (B) ? (A) : (B))
#endif

#define PI 3.14159265

static	double	mu = 1/82.45;
static	double	mus = 1 - 1/82.45;

/* 计算核心 */
static void yprime(
		   double	yp[],
		   double	*t,
 		   double	y[]
		   )
{
	double	r1,r2;
	
	r1 = sqrt((y[0]+mu)*(y[0]+mu) + y[2]*y[2]);
	r2 = sqrt((y[0]-mus)*(y[0]-mus) + y[2]*y[2]);
	
	/* 如果有除0操作，显示警告信息 */
	if (r1 == 0.0 || r2 == 0.0 )
		mexWarnMsgTxt("Division by zero!\n");
	
	yp[0] = y[1];
	yp[1] = 2*y[3]+y[0]-mus*(y[0]+mu)/(r1*r1*r1)-mu*(y[0]-mus)/(r2*r2*r2);
	yp[2] = y[3];
	yp[3] = -2*y[1] + y[2] - mus*y[2]/(r1*r1*r1) - mu*y[2]/(r2*r2*r2);
}

/* mex程序接口函数 */
void mexFunction( int nlhs, mxArray *plhs[],
				  int nrhs, const mxArray*prhs[] )

{
    double *yp;
    double *t,*y;
    unsigned int m,n;

    /* 检查输入变量个数 */

    if (nrhs != 2)
		mexErrMsgTxt("Two input arguments required.");
    else if (nlhs > 1)
		mexErrMsgTxt("Too many output arguments.");

    /* 检查Y的维数：4 X 1 或 1 X 4 */
    m = mxGetM(Y_IN);
    n = mxGetN(Y_IN);
    if (!mxIsDouble(Y_IN) || mxIsComplex(Y_IN) ||
		(MAX(m,n) != 4)   || (MIN(m,n) != 1))
		mexErrMsgTxt("YPRIME requires that Y be a 4 x 1 vector.");

    /* 为输出变量创建mxArray对象 */
    YP_OUT = mxCreateDoubleMatrix(m, n, mxREAL);

    /* 获取输入、输出变量的数据指针 */
    yp = mxGetPr(YP_OUT);
    t  = mxGetPr(T_IN);
    y  = mxGetPr(Y_IN);

    /* 调用计算程序 */
    yprime(yp,t,y);
}