/*
   da__o
   Any point-->OrbitCurve
*/

#include "common.h"
#include "corefunc.h"

#include "stagen.h"

#include "_util.h"
#include "_linalg.h"

#include "da__o.h"

/*
   IMPORTANT: all headers of integrators must have identical two first fields
   in their data structures: X0 and ActiveUserTest.
*/
#include "_integr.h"	/* integrator's (i.e. generator's) data */


/* Numbers of visible name classes;
   the order must be the same as in /names section */
#define CL_TIME		1
#define CL_PHASE	2
#define CL_PAR		3
#define CL_UTEST	4


Local(StagenDataPtr) stagenData=NULL;	/* global copy of Starter's parameter passed by content */
Local(da__o_DataPtr) staData=NULL;	/* global copy of our own data area */
Local(StdLinAlgDataPtr) staFunc=NULL;	/* global copy of pointer to standard linea algebra functions */
Local(UtilitiesDataPtr) staUtil=NULL;	/* global copy of pointer to utility functions */
Local(Int2) npar,nphase;		/* dims of par and phase spaces */
Local(Vector) x1=NULL;			/* ptr to vector - the first point produced by starter */
Local(VAR) t;				/* time */
Local(FloatHiPtr) P=NULL;		/* current set of params used by Decoder */
Local(FloatHiPtr) JJJ;                  /* current complete Jacobian */
Local(FloatHiPtr) MMM;                  /* current complete Mass Matrix */

/* Initialize function common for Starter and Decoder */
Local(Int2) Init(StagenDataPtr stagenptr) {
  /* Set local pointers to data area and functions array */
  stagenData=stagenptr;
  staData=(da__o_DataPtr)stagenptr->StaPtr;
  /* Use standard linea algebra */
  staFunc=(StdLinAlgDataPtr)stagenptr->StdLinAlg;
  staUtil=(UtilitiesDataPtr)stagenptr->Utilities;
  /* Get dimensions of Phase and Par spaces */
  npar=*stagenptr->ClassDim[CL_PAR];
  nphase=*stagenptr->ClassDim[CL_PHASE];

  
  /* Call udf SetInitPoint function, if any */
  if (!stagenData->ContDataGenLen && FuncParActive(staData->SetInitPoint)) {	/* there is UDF to set InitPoint */
    Int2 Index_X,Index_P,Index_T;
    Index_X=stagenptr->ParIndex(staData,X);
    Index_P=stagenptr->ParIndex(staData,P);
    Index_T=stagenptr->ParIndex(staData,T);
    FuncParCall(staData->SetInitPoint)(staData);
    stagenptr->UpdatePar(Index_X);
    stagenptr->UpdatePar(Index_P);
    stagenptr->UpdatePar(Index_T);
  }
  /* Create a copy parameters */
  if (!P) {
    P=MemNew(sizeof(FloatHi)*npar);
    MemCopy(P,staData->P,sizeof(FloatHi)*npar);
    /* Allocate memory for values of rhs used by Decoder/Define/Tests */
    x1=staFunc->CreateVector(nphase+1);
  }
  stagenData->BifDataOutLen=0;
  ((IntegrDataPtr)stagenptr->GenPtr)->ActiveUserTest=staData->ActiveUserTest;
  return 0;
}

Local(void) Term (Boolean OK) {
    if (OK&&P&&(stagenData->Reason==RF_CLEAN)) stagenData->ProcessPoint(stagenData->ProcFuncNum, x1, "Last point");
    if (P) {
      P=MemFree(P);
      x1=staFunc->DestroyVector(x1);
    }
}


/* Any point  -> Orbit Starter */

Entry(Int2) da_o_Starter(StagenDataPtr stagenptr) {
  Int2 i,Ret;
  if (!stagenptr) {	/* this is request to terminate */
    Term(TRUE);
    return 0;
  }
  if (Init(stagenptr)) return 1;	/* initialize local variables */
  
  /* Create the first point for generator */
  x1[0]=*stagenptr->ClassDim[CL_TIME] ? *staData->T : *staData->Tdef;
  MemCopy(x1+1,staData->X,sizeof(FloatHi)*nphase);  
  ((IntegrDataPtr)stagenptr->GenPtr)->X0=x1;
  /* That's all */
  return 0;
}

/***********/
/* Decoder */
Entry(Int2) da_o_Decoder(VoidPtr vpoint, FloatHiPtr PNTR indirect, Int2 oper) {
  switch (oper) {
    case DECODER_INIT:	/* point is a pointer to StagenData structure */
    case DECODER_TERM:
      break;
    case DECODER_DIM:	/* returns number of M-points in given G-point */
      return 1;
    default:		/* extract next M-point from G-point pointed by vpoint */
      /* vpoint is a Vector */
      indirect[CL_TIME]=(VAR *)vpoint;
      indirect[CL_PHASE]=(VAR *)vpoint+1;
      indirect[CL_PAR]=P;
  }
  return 0;
}

/**************/
/* System RHS */
Entry(Int2) da_o_DefRHS (Vector x, Vector y) {
  stagenData->Rhs(x+1,P,&x[0],y+1);    /* 0-component is for time */
  y[0]=0;                           	
  return 0;
}

/**************/
/* System Jacobian */
Entry(Int2) da_o_DefJAC (Vector x, Matrix Jac) {
  int i;
  JJJ=stagenData->Der1(x+1,P,&x[0]);    /* 0-component is for time */
  for(i=0;i<nphase;i++) MemCopy(Jac[i],JJJ+(nphase+npar)*i,sizeof(FloatHi)*nphase);
  return 0;
}

/**************/
/* Mass Matrix */
Entry(Int2) da_o_DefMAS (Matrix Mas) {
  stagenData->AuxFunc[5](P,Mas);    /* 0-component is for time */

/**{{
typedef double mt[2][2][2];
double *m;
int i,j,k;
typedef EntryPtr(double *,LhsDerPtr,(double *p));
m=((LhsDerPtr)(stagenData->AuxFunc[6]))(P);
  printf("\n");
for (i=0; i<2; i++) {
  for (j=0; j<2; j++) {
    for (k=0; k<2; k++)
      printf ("%i %i %i: %g\n",i,j,k,m[2*2*i+2*j+k]);
  }
  printf("\n");
}
}}**/

  return 0;
}

/*********************/
/* Default processor */
Entry(Int2) da_o_ProcDefault(Int2 Ind, Vector x, Vector v1, CharPtr *msg) {
  /* Assert: v1==NULL, msg==NULL, Ind==0 */
  MemCopy(x1,x,sizeof(VAR)*(nphase+1));     /* for the last point */
  return 0;
}