#include <values.h>

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

#include "stagen.h"

#include "_linalg.h"

#include "_util.h"

#include "_iter.h"

/* Ptr to Data area passed to generator by content.
   Its structure is defined in stagen.h */
Local(StagenDataPtr) genData;

#define NumSing (stagen->ProcFuncNum-1)	/* 0-th is a default processor */

/* The pointer to Generator's functional parameters
   points to the following structure
   (names of functions must be listed in appropriate
    generator's description files under /genfunc section
    in proper order). */
typedef struct {
  PrintVectorPtr	PrintVector;
  CreateVectorPtr       CreateVector;
  DestroyVectorPtr      DestroyVector;
  GetVectorDimPtr       GetVectorDim;
  CopyVectorPtr         CopyVector;
  ScalProdPtr           ScalProd;
  NormalizeVectorPtr    NormalizeVector;
  AddScaledVectorPtr    AddScaledVector;
  CopyVectorElementsPtr CopyVectorElements;
} ContFuncPars, PNTR ContFuncParsPtr;

/* Pointer to a structure containing pointers to functional parameters */
Local(ContFuncParsPtr) contFunc;

/* Pointer to a continuator's own data area */
Local(ContDataPtr) contData;

/***************************************/

#define RHS        genData->DefFunc[0]

/***************************/
#define Nmax    contData->Nmax	        /* Number of iterations */
#define ActiveUTestFuns     contData->ActiveUserTest
					/* continuation data */
#define nUTest        genData->udFuncNum

#define DefaultProcessor(p) ((ProcFuncPtr)genData->ProcFunc[0])(0,p,NULL,NULL);

Local(VAR) CurTint;		/* current maximum of time */

/***************************************/

typedef EntryPtr(Int2,ProcFuncPtr,(Int2 Ind, Vector xx,  Vector vv, CharPtr *msg));

/****************************/


/* Message Processing */

Local(Int2) MsgProc(Pointtype type, Vector p, char *text) {
  /* PostProcessing */
  return genData->ProcessPoint(type, p, text);
}

/***************/
/* StopCheck   */
Local(int) np;                 /* number of points along orbit segment */
Local(int) npc;                /* total number of points along orbit */

Local(int) Checker(Vector X) {
  npc++;
  if (++np>Nmax) {
    CurTint+=Nmax;
    return 1;
  }	
  return 0;
}

Local(Vector) X,F;

Local(void) Term (void) {
  contFunc->DestroyVector(X);
  contFunc->DestroyVector(F);
}

/* Entry point to the generator */
Entry(Int2) OrbitMap(StagenDataPtr stagen) {


  Int2 n;		        /* space dimension */
  VAR TV;
  int i,ii;
  int Cont=1,Ret=0;

  if (!stagen) {
    Term();
    return Ret;
  }
  genData=stagen;
  contData=(ContDataPtr)stagen->GenPtr;
  contFunc=(ContFuncParsPtr)stagen->GenFunc;

  n=contFunc->GetVectorDim(contData->X0);

  X = contFunc->CreateVector(n);
  F = contFunc->CreateVector(n);

/* Orbit "continuation" */
  if (genData->ContDataGenLen) { 
    MemCopy(X,genData->ContDataGenPtr,n*sizeof(VAR));
    CurTint=genData->ContDataGenPtr[n+1];
    npc=(int)genData->ContDataGenPtr[n];
    MemFree(genData->ContDataGenPtr);  
  } else {
    CurTint=Nmax;
    contFunc->CopyVector(contData->X0,X);
    npc=X[0]+1;
  }
  np=1;

  /* User test function values at the first point  */
  for (i=0; i<nUTest; i++) 
    if (ActiveUTestFuns[i]) {
      ii=(genData->udFunc)(X,i,&TV);
      if (ii) { 
        MsgProc(0,X,"Undefined user test function in the first point");
        Ret=2;
        Cont=0;
	break;
      }
    }
  /* Regular step */
  while (Cont) {

    /* Issue the message */
    if (MsgProc(0,X,NULL)) {
       Cont=0;
       continue;
    }

    /* STEP */
    RHS(X,F);
    contFunc->CopyVector(F,X);
    X[0]=npc;
        
    /* Compute values of all test function at the next regular point */
    for (i=0; i<nUTest; i++) 
      if (ActiveUTestFuns[i]) genData->udFunc(X,i,&TV);

    /* Stop check  */
    switch (Checker(X)) {
      case 1:		/* number of points */
        Cont=0;
	MsgProc(0,X,NULL);        
	continue;
    }
  } /* while */
  
  DefaultProcessor(X);
  /* assert: V is the tangent vector at the last found point X */
  genData->ContDataGenLen=(n+2)*sizeof(VAR);
  genData->ContDataGenPtr=MemNew(genData->ContDataGenLen);
  MemCopy(genData->ContDataGenPtr,X,n*sizeof(VAR));
  genData->ContDataGenPtr[n]=npc;
  genData->ContDataGenPtr[n+1]=CurTint;
  
  Term();

  return (Ret);
}