/************************************************************
*  Lorenz.h
*  說明：lorenz吸引子
*  2004.11.12 zya
************************************************************/

#include <stdlib.h>    
#include "matlab.h"
//lorenz吸引子參數
#define SIGMA 10
#define RHO 28
#define BETA 8/3

mxArray *lorenz(mxArray *tm, mxArray *ym)
{
    mxArray *ypm = NULL;
    double *y, *yp;

    mlfEnterNewContext(0,2,tm,ym);

    mlfAssign(&ypm, mlfDoubleMatrix(3,1,NULL,NULL));
    y = mxGetPr(ym);
    yp = mxGetPr(ypm);
    //定義lorenz映射
    yp[0] = -BETA*y[0] + y[1]*y[2];
    yp[1] = -SIGMA*y[1] + SIGMA*y[2];
    yp[2] = -y[0]*y[1] + RHO*y[1] - y[2];
    
    mlfRestorePreviousContext(0,2,tm,ym);
    return mlfReturnValue(ypm);
}

static int _lorenz_thunk_fcn_( mlfFuncp pFunc, int nlhs, mxArray **lhs, int nrhs, mxArray **rhs )
{
    typedef mxArray * (*PFCN_1_2)( mxArray * , mxArray * );
    mxArray *Out;

    if ( nlhs > 1 || nrhs > 2 ) {
        return(0);
    }
    
    Out = (*((PFCN_1_2)pFunc))(
                               nrhs > 0 ? rhs[0] : 0,
                               nrhs > 1 ? rhs[1] : 0
                               );
    
    if (nlhs > 0) {
        lhs[0] = Out;
    }
    
    return(1);
}

static mlfFuncTabEnt MFuncTab1[] = 
{
	{ "lorenz", (mlfFuncp) lorenz, _lorenz_thunk_fcn_ },
	{ 0, 0, 0 }
};

int lorenz(double *X, double *Y, double *Z,double *yinit, int seqlen)
{
    mxArray *tm = NULL, *ym = NULL, *tsm = NULL, *ysm = NULL;
    double *tspan= new double[seqlen];
    double *lorenz_X;
    double *lorenz_Y;
	double *lorenz_Z;
    //double *t;
    int n,i;

	for (i=0; i<seqlen; i++)
	{
		tspan[i] = i*0.01;
	}

    mlfEnterNewContext(0, 0);

    mlfFevalTableSetup ( MFuncTab1 );

    mlfAssign(&tsm, mlfDoubleMatrix(seqlen, 1, tspan, NULL));
    mlfAssign(&ysm, mlfDoubleMatrix(1, 3, yinit, NULL));

    mlfOde45(mlfVarargout(&tm, &ym, NULL),
             mxCreateString("lorenz"), tsm, ysm, NULL, NULL);

    n = mxGetM(tm);
    //t = mxGetPr(tm);
	lorenz_X = mxGetPr(ym);
    lorenz_Y = lorenz_X+n;
	lorenz_Z = lorenz_Y+n;
	for (i=0; i<seqlen; i++)  //mat數據以行為先的形式存儲，與c語言等數據存儲方式不同
	{
		X[i] = lorenz_X[i];
		Y[i] = lorenz_Y[i];
		Z[i] = lorenz_Z[i];
		//cout<<X[i];
	}
    /* Free the matrices. */
    mxDestroyArray(tsm);
    mxDestroyArray(ysm);
    mxDestroyArray(tm);
    mxDestroyArray(ym);
	delete [] tspan;

    mlfRestorePreviousContext(0,0);
    return(EXIT_SUCCESS);
}