/********************************************************************
*  HyperRossler.h
*  說明：超Rossler吸引子
*  2004.11.12 zya
*********************************************************************/

#include <stdlib.h>    
#include "matlab.h"

//超Rossler吸引子參數
#define rossler_a 0.25
#define rossler_b 3
#define rossler_c 0.5
#define rossler_d 0.05

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

    mlfEnterNewContext(0,2,tm,ym);

    mlfAssign(&ypm, mlfDoubleMatrix(4,1,NULL,NULL));
    y = mxGetPr(ym);
    yp = mxGetPr(ypm);
	//定義超Rossler映射
	yp[0] = -y[1]-y[2];
	yp[1] = y[0]+rossler_a*y[1]+y[3];
	yp[2] = rossler_b+y[0]*y[2];
	yp[3] = -rossler_c*y[2]+rossler_d*y[3];

    mlfRestorePreviousContext(0,2,tm,ym);
    return mlfReturnValue(ypm);

}

static int _hyperrossler_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 MFuncTab[] = 
{
	{ "hyperrossler", (mlfFuncp) hyperrossler, _hyperrossler_thunk_fcn_ },
	{ 0, 0, 0 }
};

int hyperrossler(double *X, double *Y, double *Z,double *W,double *yinit, int seqlen)
{
    mxArray *tm = NULL, *ym = NULL, *tsm = NULL, *ysm = NULL;
    double *tspan= new double[seqlen];
    double *hyperrossler_X;
    double *hyperrossler_Y;
	double *hyperrossler_Z;
	double *hyperrossler_W;
    //double *t;
    int n,i;
    //設定時間間隔
	for (i=0; i<seqlen; i++)
	{
		tspan[i] = i*0.1;
	}

    mlfEnterNewContext(0, 0);

    mlfFevalTableSetup ( MFuncTab );

    mlfAssign(&tsm, mlfDoubleMatrix(seqlen, 1, tspan, NULL));
	//超Rossler映射初始值
    mlfAssign(&ysm, mlfDoubleMatrix(1, 4, yinit, NULL));

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

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

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