/* $Revision: 1.2 $ */
/*
 *
 * YPRIME.C Sample .MEX file corresponding to YPRIME.M
 *					Solves simple 3 body orbit problem 
 *
 * The calling syntax is:
 *
 *		[yp] = yprime(t, y)
 *
 * Copyright (c) 1984-1998 by The MathWorks, Inc.
 * All Rights Reserved.
 */

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

/* Input Arguments */

#define X_IN	prhs[0]
#define C_IN	prhs[1]
#define V_IN	prhs[2]
#define L_IN	prhs[3]
#define NORM 	prhs[4]
#define MB_TYP_IN prhs[5]
#define M_PAR_IN prhs[6]
/* Output Arguments */

#define Y_OUT plhs[0]

static void calcOut(
			 double *f,
			 double *x,
			 double *cp,
			 double *vp,
			 double *lp,
			 int		nPoints,
			 int		nIn,
			 int		nOut,
			 int		nRules,
			 int		norm,
			 int		memb_type,
			 double		m)
			 
{
	double *xMinC=NULL, *memb=NULL, *y=NULL;
	double sumMemb;
	int thePoint, xOff, fOff;
	
	/* Allocate the space for the temporary variables */
	xMinC=(double*)mxCalloc(nRules*nIn,sizeof(double));
	y=(double*)mxCalloc(nRules*nOut,sizeof(double));
	memb=(double*)mxCalloc(nRules,sizeof(double));
	/* mxCalloc throws an exception to Matlab if there is a problem */
	
	
	for (thePoint=0,xOff=0,fOff=0;thePoint<nPoints;thePoint++,xOff+=nIn,fOff+=nOut)
	{
		process(f + fOff, xMinC, memb, y, &sumMemb, x + xOff, cp, vp, lp, nIn, nOut, nRules, norm, memb_type, m);
	}
	
	mxFree(xMinC);
	mxFree(memb);
	mxFree(y);
	
 }

void mexFunction(
								 int nlhs,	mxArray *plhs[],
								 int nrhs, const mxArray *prhs[])
								 
{
	double		*out;
	double		*x,*c,*v,*l;
	double		m;
	const int		*dim_array_C;
	int			nRules,nIn,nOut,nPoints,memb_type, norm;
 
					 
	/* Check for proper number of arguments */
	
	if (nrhs < 4) {
		mexErrMsgTxt("outMex requires four input arguments.");
	} else	if (nrhs > 7) {
		mexErrMsgTxt("outMex don't handle more than seven input arguments.");
	} else if (nlhs > 1) {
		mexErrMsgTxt("outMex requires one output argument.");
	}
	
	
	dim_array_C = mxGetDimensions(C_IN);
	nIn = dim_array_C[0];
	nRules = dim_array_C[1];
	nOut = mxGetM(L_IN);
	nPoints = mxGetN(X_IN);
	
 
	/* Create a matrix for the return argument */
	
	Y_OUT = mxCreateDoubleMatrix(nOut, nPoints, mxREAL);
	
	
	/* Assign pointers to the various parameters */
	
	out = mxGetPr(Y_OUT);
	
	x = mxGetPr(X_IN);
	c = mxGetPr(C_IN);
	v = mxGetPr(V_IN);
	l = mxGetPr(L_IN);
	
	if (nrhs > 4)
		norm = *mxGetPr(NORM);
	else
		norm = 1;
	
	if (nrhs > 5)
		memb_type = *mxGetPr(MB_TYP_IN);
	else
		memb_type = 1;
	
	if (nrhs > 6)
		m = *mxGetPr(M_PAR_IN);
 	else
		m = 2;
	
	/* Do the actual computations in a subroutine */
	
	calcOut(out,x,c,v,l,nPoints,nIn,nOut,nRules,norm,memb_type,m);
	return;
}


