/**
* \file return_mex_arguments.cpp
* \author basile.graf@epfl.ch
* 
* This file contains functions for returning parameters to MATLAB
*/


#include <vector>

#include "Datas.h"
#include "return_mex_arguments.h"

#include "../CVODEincludes.h"  // ODE solver includes

using namespace std;


#ifndef PI_DEF
#define pi 3.141592653589793
#define PI 3.141592653589793
#endif

#define Ith(v,i)    NV_Ith_S(v,i)       /* Ith numbers components 1..NEQ */




 vector <double> sout00;
 vector <double> sout01;
 vector <double> sout02;
 vector <double> sout03;
 vector <double> sout04;
 vector <double> sout05;
 vector <double> sout06;
 vector <double> sout07;
 vector <double> sout08;
 vector <double> sout09;
 vector <double> sout10;
 vector <double> sout11;
 vector <double> timev;

void store_state(double time, N_Vector state_at_t)
{
	sout00.push_back(Ith(state_at_t,0));
	sout01.push_back(Ith(state_at_t,1));
	sout02.push_back(Ith(state_at_t,2));
	sout03.push_back(Ith(state_at_t,3));
	sout04.push_back(Ith(state_at_t,4));
	sout05.push_back(Ith(state_at_t,5));
	sout06.push_back(Ith(state_at_t,6));
	sout07.push_back(Ith(state_at_t,7));
	sout08.push_back(Ith(state_at_t,8));
	sout09.push_back(Ith(state_at_t,9));
	sout10.push_back(Ith(state_at_t,10));
	sout11.push_back(Ith(state_at_t,11));
	timev.push_back(time);
}

void return_args(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[], ODE_data *f_dat)
{
	double *outArray;
	double *outTime;
	double *numStep;
	long n, m, N, M;

	N=12;
	M=sout00.size();

	if (nlhs>0)
	{
		plhs[0] = mxCreateDoubleMatrix(N, M, mxREAL);
		outArray = mxGetPr(plhs[0]);
		for (m=M-1; m>=0; m--)
		{
			outArray[m*N+0] = sout00.back(); sout00.pop_back();
			outArray[m*N+1] = sout01.back(); sout01.pop_back();
			outArray[m*N+2] = sout02.back(); sout02.pop_back();
			outArray[m*N+3] = sout03.back(); sout03.pop_back();
			outArray[m*N+4] = sout04.back(); sout04.pop_back();
			outArray[m*N+5] = sout05.back(); sout05.pop_back();
			outArray[m*N+6] = sout06.back(); sout06.pop_back();
			outArray[m*N+7] = sout07.back(); sout07.pop_back();
			outArray[m*N+8] = sout08.back(); sout08.pop_back();
			outArray[m*N+9] = sout09.back(); sout09.pop_back();
			outArray[m*N+10]= sout10.back(); sout10.pop_back();
			outArray[m*N+11]= sout11.back(); sout11.pop_back();
		}
	}
	if (nlhs>1)
	{
		plhs[1] = mxCreateDoubleMatrix(1, M, mxREAL);
		outTime = mxGetPr(plhs[1]);
		for (m=M-1; m>=0; m--)
		{
			outTime[m] = timev.back(); timev.pop_back();
		}
	}
	if (nlhs>2)
	{
		plhs[2] = mxCreateDoubleMatrix(1, 1, mxREAL);
		numStep = mxGetPr(plhs[2]);
		numStep[0]=f_dat->CumulNumOfSteps;
	}
}

