/*
   m__fp
   Any point-->FixedPointCurve
*/

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

#include "stagen.h"

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

#include "m__fp.h"
#include "_conti.h"	/* continuator'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_TEST         4
#define CL_UTEST	5
#define CL_STDFUN	6
#define CL_MULT		7

Local(StagenDataPtr) stagenData=NULL;	/* global copy of Starter's parameter passed by content */
Local(m__fp_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(Int2Ptr) active;	        	/* index of active parameter(s) */
Local(Int2) nappend;			/* number of appended user 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(Vector) x0=NULL;			/* approximate first point */
Local(Vector) v0=NULL;			/* approximate tangent vector at first point */
Local(FloatHiPtr) X=NULL;		/* current set of phases used by Decoder */
Local(FloatHiPtr) P=NULL;		/* current set of params used by Decoder */
Local(FloatHiPtr) M=NULL;		/* current set of multipliers */
Local(VAR) zero=0;
Local(Vector) StdFun=NULL;
Local(int) itmap;
Local(FloatHiPtr) F,DF,DFF,D2,DD2,D3F,DD3F;		/* work space for supdiff */

#define _itmap       staData->itmap   	/* number of superpositions */
#define _dx          staData->dx        /* increment */

#include "_supdiff.c"			    /* symbolic derivatives for superposition */

/* Working space for test functions  */

Local(Matrix) A0=NULL,A=NULL,B=NULL,C=NULL,AA=NULL;
Local(Vector) Re=NULL,Im=NULL,Mod=NULL,Arg=NULL;
Local(VAR) test[4],dx,ddx,dddx;
Local(Int2) i;               

#define BifVector 	    stagenData->BifDataOutPtr
#define ActiveUTestFuns     staData->ActiveUserTest
#define nUTest      	    stagenData->udFuncNum

/* Initialize function common for Starter and Decoder */
Local(Int2) Init(StagenDataPtr stagenptr) {
  /* Set local pointers to data area and functions array */
  stagenData=stagenptr;
  staData=(m__fp_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];
  itmap=_itmap;
  if (itmap<=0) {
    staUtil->ErrMessage("Superposition must be positive.");
    return 6;
  }
  /* Call udf SetInitPoint function, if any */
  if (!stagenData->ContDataGenLen && FuncParActive(staData->SetInitPoint)) {	/* there is UDF to set InitPoint */
    Int2 Index_X,Index_P;
    Index_X=stagenptr->ParIndex(staData,X);
    Index_P=stagenptr->ParIndex(staData,P);
    FuncParCall(staData->SetInitPoint)(staData);
    stagenptr->UpdatePar(Index_X);
    stagenptr->UpdatePar(Index_P);
  }
  /* Do all the checks related to active parameters and find them */
  for (i=0,nappend=0; i<nUTest; i++) if (ActiveUTestFuns[i]==UDF_APPEND) nappend++;
  active=MemNew(sizeof(Int2)*(nappend+1));    
  if (staUtil->CheckParameters(staData->Active,npar,active,1+nappend)) {
    active=MemFree(active);
    return 1;
  }
  /* Checks related to analytical/numerical differentiation of RHS */
  if (!stagenData->Der1 && staData->dx <= 0) {
    staUtil->ErrMessage("Analytical Jacoby matrix is undefined.\nRecompile system or set positive 'Increment'");
    active=MemFree(active);
    return 2;
  }
  if ((staData->ActiveSing[0] || staData->ActiveSing[1] || staData->ActiveSing[2]) && 
      !stagenData->Der2 && staData->dx <= 0) {
    staUtil->ErrMessage("Analytical Hessian matrix is undefined.\nRecompile system or set positive 'Increment'");
    active=MemFree(active);
    return 3;
  }
  if (staData->ActiveSing[2] && !stagenData->Der3 && staData->dx <= 0) {
    staUtil->ErrMessage("Analytical matrix of third derivatives is undefined.\nRecompile system or set positive 'Increment'");
    active=MemFree(active);
    return 4;
  }
  if (staData->dx <= 0 && nappend > 0) {
    staUtil->ErrMessage("Set positive 'Increment'");
    active=MemFree(active);
    return 5;
  }
  /* Create a copy parameters */
  Der2num=ind2_(nphase+npar-1,nphase+npar-1)+1;
  Der3num=ind3_(nphase+npar-1,nphase+npar-1,nphase+npar-1)+1;
  if (!P) {
    P=MemNew(sizeof(FloatHi)*npar);
    MemCopy(P,staData->P,sizeof(FloatHi)*npar);
    /* Allocate memory for copy of phase vars used by Decoder */
    X=MemNew(sizeof(FloatHi)*nphase);
    /* Allocate memory for values of rhs used by Decoder/Define/Tests */
    /* Allocate memory for eigenvalues */
    M=MemNew(sizeof(FloatHi)*4*nphase);
    x0=staFunc->CreateVector(nphase+1+nappend);
    x1=staFunc->CreateVector(nphase+1+nappend);
    v0=staFunc->CreateVector(nphase+1+nappend);
    Mod=staFunc->CreateVector(nphase);
    Arg=staFunc->CreateVector(nphase);
    Re=staFunc->CreateVector(nphase);
    Im=staFunc->CreateVector(nphase);
    StdFun=staFunc->CreateVector(1+2*nphase); /* L2_NORM, MIN&MAX */
    /* Allocate memory for matrices used in test functions and in the default processor */
    A0=staFunc->CreateMatrix(nphase,nphase+1+nappend);
    A=staFunc->CreateMatrix(nphase+nappend,nphase+1+nappend);
    AA=staFunc->CreateMatrix(nphase,nphase);
    B=staFunc->CreateMatrix(nphase+1+nappend,nphase+1+nappend);
    if (nphase >= 2 && staData->ActiveSing[3]) 
       C=staFunc->CreateMatrix(nphase*(nphase-1)/2,nphase*(nphase-1)/2);
    BifVector=MemNew(sizeof(FloatHi)*(4*nphase+3+3*nappend+1)); 
    F=MemNew(sizeof(FloatHi)*nphase);  
    DF=MemNew(sizeof(FloatHi)*nphase*(nphase+npar));  
    DFF=MemNew(sizeof(FloatHi)*nphase*(nphase+npar));  
    if (staData->ActiveSing[0]+staData->ActiveSing[1]+staData->ActiveSing[2]+staData->ActiveSing[3]) {
      D2=MemNew(sizeof(FloatHi)*nphase*Der2num);
      DD2=MemNew(sizeof(FloatHi)*nphase*Der2num);
      if (staData->ActiveSing[2]+staData->ActiveSing[3]) {
	D3F=MemNew(sizeof(FloatHi)*nphase*Der3num);
	DD3F=MemNew(sizeof(FloatHi)*nphase*Der3num);
      }
    }
  }
  dx=_dx;
  ddx=pow(dx,3.0/4.0);
  dddx=pow(dx,3.0/6.0);

  ((ContDataPtr)stagenptr->GenPtr)->ActiveSing=staData->ActiveSing;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveUserTest=staData->ActiveUserTest;
  ((ContDataPtr)stagenptr->GenPtr)->TypeTestFuns=staData->TypeTestFunc;
  ((ContDataPtr)stagenptr->GenPtr)->ZeroTestFuns=staData->ZeroTestFunc;

  return 0;
}
Local(void) Term (Boolean OK) {
    if (OK&&P&&(stagenData->Reason==RF_CLEAN)) {
      stagenData->ProcessPoint(stagenData->ProcFuncNum, x1, "Last point");
      stagenData->BifDataOutLen=0;
    }
    active=MemFree(active);
    if (P) {
      P=MemFree(P);
      X=MemFree(X);
      M=MemFree(M);
      x0=staFunc->DestroyVector(x0);
      x1=staFunc->DestroyVector(x1);
      v0=staFunc->DestroyVector(v0);
      Re=staFunc->DestroyVector(Re);
      Im=staFunc->DestroyVector(Im);
      Mod=staFunc->DestroyVector(Mod);
      Arg=staFunc->DestroyVector(Arg);
      StdFun=staFunc->DestroyVector(StdFun);
      A=staFunc->DestroyMatrix(A);
      AA=staFunc->DestroyMatrix(AA);
      A0=staFunc->DestroyMatrix(A0);
      B=staFunc->DestroyMatrix(B);
      C=staFunc->DestroyMatrix(C);
      BifVector=MemFree(BifVector);
      F=MemFree(F);
      DF=MemFree(DF);
      DFF=MemFree(DFF);
      if (staData->ActiveSing[0]+staData->ActiveSing[1]+staData->ActiveSing[2]+staData->ActiveSing[3]) {
        D2=MemFree(D2);
        DD2=MemFree(DD2);
	if (staData->ActiveSing[2]+staData->ActiveSing[3]) {
	  D3F=MemFree(D3F);
	  DD3F=MemFree(DD3F);
	}
      }
    }
}

/**********************************************/
/* Starters to continue the fixed point curve */

Local(Int2) SimpleInit(Vector X1, Vector X0, Vector V0);

/* o->fp Starter */
Entry(Int2) o_fp_Starter(StagenDataPtr stagenptr) {
Int2 i;
  if (!stagenptr) {	/* this is request to terminate */
    Term(TRUE);
    return 0;
  }
  /* initialize local variables */
  if (Init(stagenptr)) return 1;

  /* Create the first point for generator */

  if (!stagenptr->ContDataGenLen) {   /* "Continuation" */
    MemCopy(x1,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<=nappend; i++) x1[nphase+i]=staData->P[active[i]];
    /* Initialize x1 */
    if (SimpleInit(x1,x0,v0)) {
      staUtil->ErrMessage("No convergence to a fixed point from the initial point"); 
      Term(FALSE);
      return 1;  
    }
  }

  ((ContDataPtr)stagenptr->GenPtr)->X0=x0;
  ((ContDataPtr)stagenptr->GenPtr)->V0=v0;

  /* That's all */
  return 0;
}
Entry(Int2) fp_JacDefFixedPoint (Vector x, Matrix Jac);

/* Any regular point on the fp curve -> fp  Starter */

Entry(Int2) fp_Starter(StagenDataPtr stagenptr) {
  Int2 i,Ret;
  Vector v1;
  Matrix D,DE;
  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 */

  if (stagenData->BifDataInOk) _itmap=stagenData->BifDataInPtr[0];
  if (!stagenptr->ContDataGenLen) {   /* "Continuation" */
    v1=staFunc->CreateVector(nphase+1+nappend);
    D=staFunc->CreateMatrix(nphase+nappend,nphase+1+nappend);
    DE=staFunc->CreateMatrix(nphase+1+nappend,nphase+1+nappend);
    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<=nappend; i++) x0[nphase+i]=staData->P[active[i]];

    Ret=fp_JacDefFixedPoint (x0,D);
    for (i=0; i<nphase+1+nappend; i++) {
      v1[i]=1.0;
      staFunc->AppendMatrixVector(D,v1,DE);
      Ret=staFunc->Decomp(DE);
      if (Ret == 0) {
        v0[nphase+nappend]=stagenData->Forward ? 1 : -1;
        staFunc->Solve(DE,v0);
        break;
      };
      v1[i]=0.0;
    };
    if (Ret) {
       staUtil->ErrMessage(
       "Unable to find a tangent vector to the fixed point curve at the initial point"); 
       v1=staFunc->DestroyVector(v1);
       D=staFunc->DestroyMatrix(D);
       DE=staFunc->DestroyMatrix(DE);
       Term(FALSE);
       return 1;
    };
    staFunc->NormalizeVector(v0);

    v1=staFunc->DestroyVector(v1);
    D=staFunc->DestroyMatrix(D);
    DE=staFunc->DestroyMatrix(DE);
  }
  /* Initialize continuation */
  ((ContDataPtr)stagenptr->GenPtr)->X0=x0;
  ((ContDataPtr)stagenptr->GenPtr)->V0=v0;
  /* That's all */
  return 0;
}

/* bp-fp Starter  */

#define _branch   staData->branch     /* switch data  */

Entry(Int2) bp_fp_Starter(StagenDataPtr stagenptr) {
  Int2 i,n,dir,actold;
  Vector V1,V2;
  Boolean coin;
  
  if (!stagenptr) {	/* this is request to terminate */
    Term(TRUE);
    return 0;
  }
  /* initialize local variables */
  if (Init(stagenptr)) return 1;
  
  /* Create the first point for generator */
  if (!stagenptr->ContDataGenLen) {   /* "Continuation" */
    V1=staFunc->CreateVector(nphase+1+nappend);
    V2=staFunc->CreateVector(nphase+1+nappend);

    if (stagenData->BifDataInOk) {
  /* tangent vectors exist */
      _itmap=stagenData->BifDataInPtr[0];
      MemCopy(V1,stagenData->BifDataInPtr+1,sizeof(FloatHi)*(nphase+1+nappend));
      MemCopy(V2,stagenData->BifDataInPtr+nphase+1+nappend+1,sizeof(FloatHi)*(nphase+1+nappend));
      for (i=0,coin=TRUE; i<=nappend; i++) {
        actold=(Int2)stagenData->BifDataInPtr[2*(nphase+1+nappend)+i+1];
	if (active[i]!=actold) coin=FALSE;
      }   
      if (!coin) {
        staUtil->ErrMessage
        ("Restore active parameters to start from the branching point");  
        V1=staFunc->DestroyVector(V1);
        V2=staFunc->DestroyVector(V2);
        Term(FALSE); 
        return 1;
      }
    } else {
  /* no tangent vectors */
      staUtil->ErrMessage("No way to start from a user-supplied branching point");  
      V1=staFunc->DestroyVector(V1);
      V2=staFunc->DestroyVector(V2);
      Term(FALSE); 
      return 1;
    };

    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<=nappend; i++) x0[nphase+i]=staData->P[active[i]];
    dir=stagenData->Forward ? 1 : -1;
    if (_branch == 0) {
      /* primary branch */
      for (n=0; n<nphase+1+nappend; n++) v0[n]=dir*V1[n];
    } else {
      /* secondary branch */
      for (n=0; n<nphase+1+nappend; n++) v0[n]=dir*V2[n];
    }
    V1=staFunc->DestroyVector(V1);
    V2=staFunc->DestroyVector(V2);
  }

  /* Initialize the first point */

  ((ContDataPtr)stagenptr->GenPtr)->X0=x0;
  ((ContDataPtr)stagenptr->GenPtr)->V0=v0;

  /* That's all */
  return 0;
}

/* pd-fp Starter  */

Entry(Int2) pd_fp_Starter(StagenDataPtr stagenptr) {
  Int2 i,k,n,dir,actold;
  Vector V,V1,V2;
  Boolean coin;
  
  if (!stagenptr) {	/* this is request to terminate */
    Term(TRUE);
    return 0;
  }
  /* initialize local variables */
  if (Init(stagenptr)) return 1;
  
  /* Create the first point for generator */
  if (!stagenptr->ContDataGenLen) {   /* "Continuation" */
    V1=staFunc->CreateVector(nphase+1+nappend);
    V2=staFunc->CreateVector(nphase+1+nappend);
    V=staFunc->CreateVector(nphase);

    if (stagenData->BifDataInOk) {
  /* tangent vectors exist */
      _itmap=stagenData->BifDataInPtr[0];
      MemCopy(V,stagenData->BifDataInPtr+1,sizeof(FloatHi)*nphase);
      MemCopy(V1,stagenData->BifDataInPtr+nphase+1,sizeof(FloatHi)*(nphase+1+nappend));
      for (i=0,coin=TRUE; i<=nappend; i++) {
        actold=(Int2)stagenData->BifDataInPtr[2*nphase+1+nappend+i+1];
	if (active[i]!=actold) coin=FALSE;
      }   
      if (!coin) {
        staUtil->ErrMessage
        ("Restore active parameters to start from the period-doubling point");  
        V=staFunc->DestroyVector(V);
        V1=staFunc->DestroyVector(V1);
        V2=staFunc->DestroyVector(V2);
        Term(FALSE); 
        return 1;
      }
    } else {
  /* no tangent vectors */
      staUtil->ErrMessage("Starting from a user-supplied PD-point is NOT IMPLEMENTED");  
      V=staFunc->DestroyVector(V);
      V1=staFunc->DestroyVector(V1);
      V2=staFunc->DestroyVector(V2);
      Term(FALSE); 
      return 1;
    };

    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<=nappend; i++) x0[nphase+i]=staData->P[active[i]];
    dir=stagenData->Forward ? 1 : -1;
    if (_branch == 0) {
      /* primary branch */
      _itmap=stagenData->BifDataInPtr[0];
      for (n=0; n<nphase+1+nappend; n++) v0[n]=dir*V1[n];
    } else {
      /* double-period branch */
      _itmap=2*(stagenData->BifDataInPtr[0]);
      staFunc->CopyVectorElements(V,0,V2,0,nphase);
      staFunc->NormalizeVector(V2);
      for (n=0; n<nphase+1+nappend; n++) v0[n]=dir*V2[n];
    }
    V=staFunc->DestroyVector(V);
    V1=staFunc->DestroyVector(V1);
    V2=staFunc->DestroyVector(V2);
    itmap=_itmap;
  }

  /* Initialize the first point */

  ((ContDataPtr)stagenptr->GenPtr)->X0=x0;
  ((ContDataPtr)stagenptr->GenPtr)->V0=v0;

  /* That's all */
  return 0;
}

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

#define _maxit   staData->maxit      /* corrector data  */
#define _maxnewt staData->maxnewt
#define _epsx    staData->epsx
#define _epsf    staData->epsf


/***********************/
/*  Newton-chord - Moore-Penrose Corrector   */

static int newtmp (Int2 Ind,
	   Vector X0,             /* Predicted point  */ 
           Vector V0,             /* Predicted normalized tangent vector */          
	   Vector X,              /* Corrected point */ 
           Vector V,              /* Corrected normalized tanget vector */
	   Int2 *it)
{
    int ret=2;
    Int2 n=staFunc->GetVectorDim(X0);	/* space dimension */

    Vector F,W=NULL,Z=NULL;
    Matrix A, B0, B=NULL;
    VAR errx, errf;

    F = staFunc->CreateVector(n-1);
    A = staFunc->CreateMatrix (n-1,n);
    B0 = staFunc->CreateMatrix (n,n);

    staFunc->CopyVector(X0,X);
    staFunc->CopyVector(V0,V);

/*  Main loop  */
    for (*it=1; *it<=_maxit; (*it)++) {

	if (stagenData->DefFunc[0](X, F)) {ret=1; goto freemem;}
	errf = staFunc->NormOfVector(F);   /* f-error */
/*  compute or restore decomposed Jacobian  */
	if (*it <= _maxnewt) {
           Z = staFunc->CopyVector(X,Z);
	   if (stagenData->DefFunc[1](Z, A)) {ret=1; goto freemem;}
           B = staFunc->AppendMatrixVector(A,V,B);
           if (staFunc->Decomp(B)) {ret=3; goto freemem;}    /* singularity check */
           staFunc->CopyMatrix(B,B0); 

/*  compute and normalize the new kernel vector  */
           if (Ind == 0) {
              W=staFunc->MultiplyMatrixVector(A,V,W);
              staFunc->CopyVectorElements(W,0,Z,0,n-1);
	      Z[n-1]=0.0;
              staFunc->Solve(B,Z);
              staFunc->AddScaledVector(V,Z,-1,V);
              staFunc->NormalizeVector(V);
           }
          }
	else {
	   staFunc->CopyMatrix(B0,B);
	}
/*  solve linear problem for displacement  */
        staFunc->CopyVectorElements(F,0,Z,0,n-1);
        Z[n-1]=0.0;
	staFunc->Solve(B,Z);

	errx = staFunc->NormOfVector(Z);   /* x-error */      
        staFunc->AddScaledVector(X,Z,-1,X);    /* correction  */
	if (errx < _epsx && errf < _epsf) break;
    }
    if (*it <= _maxit)  ret=0;
freemem:;
    staFunc->DestroyVector(F);
    staFunc->DestroyVector(W);
    staFunc->DestroyVector(Z);
    staFunc->DestroyMatrix(A);
    staFunc->DestroyMatrix(B);
    staFunc->DestroyMatrix(B0);

    return ret;
}


/* Initializer for p->fp starter */
Local(Int2) SimpleInit(Vector X1, Vector X0, Vector V0)
{
  Int2 i,j,k=1,it,n;
  Vector X=NULL,V=NULL;
  n=staFunc->GetVectorDim(X1);
  X=staFunc->CreateVector(n);
  staFunc->CopyVector(X1,X);
  V=staFunc->CreateVector(n);
  for (i=n-1; i>=0; i--) {
    V[i]=stagenData->Forward ? 1 : -1;
    j=newtmp(0, X, V, X0, V0, &it);
    if (j==0) {k=0; break;}
    V[i]=0.0;
  }
  staFunc->DestroyVector(X);
  staFunc->DestroyVector(V);
  return (k);
}

/************************/
/* Some common routines */
Local(void) CopyToX(Vector vp) {
  MemCopy(X,vp,sizeof(FloatHi)*nphase);
}

Local(void) CopyToP(Vector vp) {
Int2 i;
  for (i=0; i<1+nappend; i++) P[active[i]]=vp[nphase+i];
}

Local(void) CopyToM(void) {
  MemCopy(M,Mod,sizeof(VAR)*nphase);
  MemCopy(M+nphase,Arg,sizeof(VAR)*nphase);
  MemCopy(M+2*nphase,Re,sizeof(VAR)*nphase);
  MemCopy(M+3*nphase,Im,sizeof(VAR)*nphase);
}

/***********/
/* Decoder */

Entry(Int2) fp_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]=&zero;
      indirect[CL_PHASE]=X;
      CopyToX(vpoint);
      indirect[CL_PAR]=P;
      CopyToP(vpoint);
      indirect[CL_MULT]=M;
      CopyToM();
      indirect[CL_TEST]=test;
      indirect[CL_STDFUN]=StdFun;
  }
  return 0;
}

/***************************/
/* User-defined functions  */
Entry(Int2) fp_DefUser (Vector x, Vector y) {
Int2 i,j;
  for (i=j=0; i<nUTest; i++) 
    if (ActiveUTestFuns[i]==UDF_APPEND) {
      (stagenData->udFunc)(x,i,&y[j]);
      j++;
    }
  return 0;
}

/*******/
/* RHS */
Entry(Int2) fp_DefRHS (Vector x, Vector y) {
Int2 i;
  CopyToP(x);
  staFunc->CopyVectorElements(x,0,y,0,nphase);
  for (i=1; i<=_itmap; i++) {
    CopyToX(y);
    stagenData->Rhs(X,P,&zero,y);
  }
  return 0;
}

/********************/
/* Jacobian of RHS  */
Local(Int2) JacRHS (Vector x, Matrix Jac) {
  if (stagenData->Der1) {
    SymJac(x,Jac);
  } else {
    if (staFunc->Der(fp_DefRHS, nphase, x, dx, Jac)) return 1;
  }
  return 0;
}

/**************************************/
/* Defining functions for fixed points */
Entry(Int2) fp_DefFixedPoint (Vector x, Vector y) {
Int2 i;
  fp_DefRHS(x,y);
  for (i=0; i<nphase; i++) y[i]-=x[i];
  if (nappend) fp_DefUser(x,y+nphase);
  return 0;
}

/*****************************************************/
/* Jacobian of defining functions for fixed points  */
Entry(Int2) fp_JacDefFixedPoint (Vector x, Matrix Jac) {
Int2 i;
  JacRHS(x, Jac);
  for (i=0; i<nphase; i++) Jac[i][i]-=1.0;  
  if (nappend) staFunc->Der(fp_DefUser, nappend, x, dx, Jac+nphase);
  return 0;
}
#define HessParElem(i,j,k) HessElem(i,nphase+active[j],nphase+active[k])

/*********************/
/* Default processor */
#define PI	4.0*atan(1.0)
#define PID2	2.0*atan(1.0)
#define SCALE	45.0/atan(1.0)

#ifdef WIN32
static int _USERENTRY fp_SortMult(const void *f, const void *s) {
#else
Entry(int) fp_SortMult(const void *f, const void *s) {
#endif
VAR Mod1,Mod2;
  Mod1=((FloatHiPtr)f)[0]*((FloatHiPtr)f)[0]+((FloatHiPtr)f)[1]*((FloatHiPtr)f)[1];
  Mod2=((FloatHiPtr)s)[0]*((FloatHiPtr)s)[0]+((FloatHiPtr)s)[1]*((FloatHiPtr)s)[1];
  if (Mod1==Mod2) {
    return 0;
  }
  return (Mod1<Mod2) ? 1 : -1;
}

Entry (Int2) fp_ProcDefault(Int2 Ind, Vector x, Vector v1, CharPtr *msg) {
VAR s;
Int2 i;
  if (Ind == -1) {  /* zero of user-defined function */
    BifVector[0]=(VAR)_itmap;
    stagenData->BifDataOutLen=sizeof(VAR);     
    return 0;
  }
  /* Assert: v1==NULL, msg==NULL */
  MemCopy(x1,x,sizeof(VAR)*(nphase+1+nappend));     /* for the last point */
  /* L2_NORM and Min&Max */
  for (i=0,s=0; i<nphase; i++) {
    s+=x[i]*x[i];
    StdFun[1+i]=x[i];
    StdFun[1+nphase+i]=x[i];
  }
  StdFun[0]=sqrt(s);
  /* Multipliers */
  if (staData->eigcomp) {
    if (JacRHS(x, A0)) return (1);     
    staFunc->CopyMatrixBlock(A0,0,0,AA,0,0,nphase,nphase); 
    staFunc->EigenValues(AA,Re,Im);
    for (i=0; i<nphase; i++) {
      M[2*i]=Re[i];
      M[2*i+1]=Im[i];
    }
    qsort(M,nphase,2*sizeof(VAR),fp_SortMult);
    for (i=0; i<nphase; i++) {
      Re[i]=M[2*i];
      Im[i]=M[2*i+1];
    }
    for (i=0; i<nphase; i++) {
      Mod[i]=sqrt(Re[i]*Re[i]+Im[i]*Im[i]);
      if (Mod[i]==0.0) {
        Arg[i]=0.0;
      } else {
        if (fabs(Re[i]/Mod[i]) > 0.5) {
          Arg[i]=atan(fabs(Im[i]/Re[i]));
        } else {
          if (Re[i]!=0.0) {
            Arg[i]=PID2-atan(fabs(Re[i]/Im[i]));
	  } else {
	    Arg[i]=PID2;
	  }
        }
      }
      if (Re[i] < 0) Arg[i]=PI-Arg[i];
      if (Im[i] < 0) Arg[i]=-Arg[i];
      Arg[i]*=SCALE;
    }
  }
  return 0;
}

/*********************************/
/* Test and processing functions */

Local(Char) buf[124];

/* Branching point (BP): Determinant of appended Jacobian matrix   */

Entry(Int2) fp_TestBranching(Vector x, Vector v,  FloatHiPtr res) {

  if (fp_JacDefFixedPoint(x, A)) return (1);
  B = staFunc->AppendMatrixVector(A, v, B);
  if (staFunc->Decomp(B)) {
    *res=0.0;
  } else {
    staFunc->GetMatrixData(B,res,NULL);  
  }
  test[0]=*res;
  return 0;
}

Entry(Int2) fp_ProcBranching(Int2 Ind, Vector x, Vector v1, CharPtr *msg) {

                        /*  A - global Jacobian matrix of the defining functions */
  Matrix D,DT;		/*  D - appended Jacobian matrix; DT=transposed(D); */
  Vector q,phi,psi;     /*  q - null-eigenvector of D, <q,q>=1, <v1,q>=0; 
                          phi - first (nphase+nappend) components of null-eigenvector of DT, <phi,phi>=1;  */
  Vector v2;            /* v2 - normalized tandent vector to the secondary branch */
  Vector p,v,w,x1,F,F1;
  int i,j,k,l;
  VAR koeff,f0,f1,f2,f3,b12,b22,eps=1.0e-2;

  stagenData->BifDataOutLen=0;

  if (Ind != 0) {
    sprintf(buf,"Branching (?)");
  } else {
    q=staFunc->CreateVector(nphase+1+nappend);
    p=staFunc->CreateVector(nphase+1+nappend);
    phi=staFunc->CreateVector(nphase+nappend);
    if (nappend) psi=staFunc->CreateVector(nappend);
    D=staFunc->CreateMatrix(nphase+1+nappend,nphase+1+nappend);
    DT=staFunc->CreateMatrix(nphase+1+nappend,nphase+1+nappend);
    v2=staFunc->CreateVector(nphase+1+nappend);
    v=staFunc->CreateVector(nphase+1+nappend);
    w=staFunc->CreateVector(nphase+1+nappend);
    x1=staFunc->CreateVector(nphase+1+nappend);
    F=staFunc->CreateVector(nphase+nappend);
    F1=staFunc->CreateVector(nappend);

/*  compute Jacobian and appended matrices */
    fp_JacDefFixedPoint(x,A);
    staFunc->AppendMatrixVector(A, v1, D);
    staFunc->TransposeMatrix(D,DT);

/*  compute null-vector and adjoint null-vector of D */ 
    if (staFunc->SngSolve(D,eps,q) != nphase+nappend) {
      sprintf(buf,"Branching: k=? (q not found)");
      goto out;
    } else {
      staFunc->NormalizeVector(q);
    }
    if (staFunc->SngSolve(DT,eps,p) != nphase+nappend) {
      sprintf(buf,"Branching: k=? (p not found)");
      goto out;
    } else {
      staFunc->NormalizeVector(p);
      staFunc->CopyVectorElements(p,0,phi,0,nphase+nappend);
      if (nappend) staFunc->CopyVectorElements(p,nphase,psi,0,nappend);
    }

/*  compute coefficients of the branching equation */
    if (stagenData->Der2)  {
        H=SymHess(x);	/* analytical second-order derivatives of RHS^itmap */
        for (b22=0.0,i=0; i<nphase; i++) {
           for (j=0; j<nphase; j++) b22+=phi[i]*HessElem(i,j,j)*q[j]*q[j];
           for (k=0; k<=nappend; k++) b22+=phi[i]*HessParElem(i,k,k)*q[nphase+k]*q[nphase+k];
        }
        for (i=0; i<nphase; i++) {
	   for (j=0; j<nphase; j++) {
	      /* Windows+Borland C: The standard loop
		 for (k=j+1; k < nphase; k++) b22+=2*phi[i]*HessElem(i,j,k)*q[j]*q[k];
		 would not work correctly with induction variable optimisation:
		 HessElem(0,0,1) does not even exists for nphase==1
              */
	      if (j+1 < nphase)
		for (k=j+1; k < nphase; k++) b22+=2*phi[i]*HessElem(i,j,k)*q[j]*q[k];
	      for (l=0; l<=nappend; l++) b22+=2*phi[i]*MixedHessElem(i,j,l)*q[j]*q[nphase+l];
           }
	   for (j=0; j<=nappend; j++) {
	       if(j+1 <= nappend)
		 for (k=j+1; k <= nappend; k++) b22+=2*phi[i]*HessParElem(i,j,k)*q[nphase+j]*q[nphase+k];
	   }
        }
	if (nappend) {       /* numerical second-order derivatives of user functions */
          fp_DefUser(x,F1);
          f0=staFunc->ScalProd(psi,F1);
          staFunc->AddScaledVector(x,q,ddx,x1);
          fp_DefUser(x1,F1);
          f1=staFunc->ScalProd(psi,F1);
          staFunc->AddScaledVector(x,q,-ddx,x1);
          fp_DefUser(x1,F1);
          f2=staFunc->ScalProd(psi,F1);
          b22+=(f1-2*f0+f2)/ddx/ddx;	
	}
/* printf("Analytical b22=%lg\n",b22);	 */     
	for (b12=0.0,i=0; i<nphase; i++) {
           for (j=0; j<nphase; j++) b12+=phi[i]*HessElem(i,j,j)*v1[j]*q[j];
           for (k=0; k<=nappend; k++) b12+=phi[i]*HessParElem(i,k,k)*v1[nphase+k]*q[nphase+k];
        }
        for (i=0; i<nphase; i++) {
	   for (j=0; j<nphase; j++) {
	      if (j+1<nphase)
		for (k=j+1; k<nphase; k++) b12+=phi[i]*HessElem(i,j,k)*(v1[j]*q[k]+v1[k]*q[j]);
	      for (l=0; l<=nappend; l++) 
	         b12+=phi[i]*MixedHessElem(i,j,l)*(v1[j]*q[nphase+l]+v1[nphase+l]*q[j]);
           }
	   for (j=0; j<=nappend; j++) {
	     if (j+1 <= nappend)
	      for (k=j+1; k <= nappend; k++)
	         b12+=phi[i]*HessParElem(i,j,k)*(v1[nphase+j]*q[nphase+k]+v1[nphase+k]*q[nphase+j]);
	   }
        }
	if (nappend) {
	  staFunc->AddScaledVector(v1,q,1.0,v);
          staFunc->AddScaledVector(v1,q,-1.0,w);
          staFunc->AddScaledVector(x,v,ddx/2,x1);
          fp_DefUser(x1,F1);
          f0=staFunc->ScalProd(psi,F1);
          staFunc->AddScaledVector(x,w,ddx/2,x1);
          fp_DefUser(x1,F1);
          f1=staFunc->ScalProd(psi,F1);
          staFunc->AddScaledVector(x,w,-ddx/2,x1);
          fp_DefUser(x1,F1);
          f2=staFunc->ScalProd(psi,F1);
          staFunc->AddScaledVector(x,v,-ddx/2,x1);
          fp_DefUser(x1,F1);
          f3=staFunc->ScalProd(psi,F1);
          b12+=(f0-f1-f2+f3)/ddx/ddx; 
	}	
/* printf("                   Analytical b12=%lg\n",b12);	*/      	
    } else {
      fp_DefFixedPoint(x,F);
      f0=staFunc->ScalProd(phi,F);
      staFunc->AddScaledVector(x,q,ddx,x1);
      fp_DefFixedPoint(x1,F);
      f1=staFunc->ScalProd(phi,F);
      staFunc->AddScaledVector(x,q,-ddx,x1);
      fp_DefFixedPoint(x1,F);
      f2=staFunc->ScalProd(phi,F);
      b22=(f1-2*f0+f2)/ddx/ddx;
/* printf("Numerical b22=%lg\n",b22); */	      
      staFunc->AddScaledVector(v1,q,1.0,v);
      staFunc->AddScaledVector(v1,q,-1.0,w);
      staFunc->AddScaledVector(x,v,ddx/2,x1);
      fp_DefFixedPoint(x1,F);
      f0=staFunc->ScalProd(phi,F);
      staFunc->AddScaledVector(x,w,ddx/2,x1);
      fp_DefFixedPoint(x1,F);
      f1=staFunc->ScalProd(phi,F);
      staFunc->AddScaledVector(x,w,-ddx/2,x1);
      fp_DefFixedPoint(x1,F);
      f2=staFunc->ScalProd(phi,F);
      staFunc->AddScaledVector(x,v,-ddx/2,x1);
      fp_DefFixedPoint(x1,F);
      f3=staFunc->ScalProd(phi,F);
      b12=(f0-f1-f2+f3)/ddx/ddx;
/* printf("                   Numerical b12=%lg\n",b12);*/	      
    }
    if (b12 == 0.0) {
      sprintf(buf,"Non-simple branching");
      goto out;
    }
    koeff=-b22/b12/2.0;
    staFunc->AddScaledVector(q,v1,koeff,v2);
    staFunc->NormalizeVector(v2);

/* provide bifurcation data for bp->fp */
    BifVector[0]=(VAR)_itmap;
    MemCopy(BifVector+1,v1,sizeof(FloatHi)*(nphase+1+nappend));
    MemCopy(BifVector+nphase+1+nappend+1,v2,sizeof(FloatHi)*(nphase+1+nappend));
    for (i=0; i<=nappend; i++) {
      BifVector[2*(nphase+1+nappend)+i+1]=(VAR)active[i];
    }
    stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+3*nappend+3+1);

    sprintf(buf,"Branching: k=%g",koeff);
out:;
    staFunc->DestroyMatrix(D);
    staFunc->DestroyMatrix(DT);
    staFunc->DestroyVector(q);
    staFunc->DestroyVector(p);
    staFunc->DestroyVector(phi);
    if (nappend) staFunc->DestroyVector(psi);
    staFunc->DestroyVector(v2);
    staFunc->DestroyVector(v);
    staFunc->DestroyVector(w);
    staFunc->DestroyVector(x1);
    staFunc->DestroyVector(F);
    staFunc->DestroyVector(F1);
  }
  *msg=buf;
  return 0;
}

/* Limit point (LP): p[active[0]]-component of the tangent vector to the curve or det */
Entry(Int2) fp_TestFold(Vector x, Vector v, FloatHiPtr res) {
int i;

  *res=1.0;
  if (nappend==0) {
    *res=v[nphase];
  } else {
    if (JacRHS(x, A0)) return 1;
    staFunc->CopyMatrixBlock(A0,0,0,AA,0,0,nphase,nphase);     
    for (i=0; i<nphase; i++) AA[i][i]-=1.0;
    if (staFunc->Decomp(AA)) {
      *res=0.0;
    } else {
      staFunc->GetMatrixData(AA,res,NULL);
    }
  return 0;
  }
  test[1]=*res;
  return 0;
}

Entry(Int2) fp_ProcFold(Int2 Ind, Vector x, Vector v, CharPtr *msg) {

  Matrix A,B;		/* A - RHS Jacobian matrix; B=transposed(A); */
  Vector V,W;		/* V - eigenvector of A, AV=V, <V,V>=1; 
                           W - eigenvector of B, BW=W, <W,V>=1;  */
  Vector x1,F;
  int i,j,k;
  VAR f0,f1,f2,eps=1.0e-2;

  stagenData->BifDataOutLen=0;

  if (Ind == 0) {
    V=staFunc->CreateVector(nphase);
    W=staFunc->CreateVector(nphase);
    A=staFunc->CreateMatrix(nphase,nphase);
    B=staFunc->CreateMatrix(nphase,nphase);
    x1=staFunc->CreateVector(nphase+1+nappend);
    F=staFunc->CreateVector(nphase);

/*  compute Jacobian matrix */
    JacRHS(x,A0);
    staFunc->CopyMatrixBlock(A0,0,0,A,0,0,nphase,nphase);

/*  provide bifurcation data for lp->lp starter */
/*  BifVector[0]:	                _itmap */
/*  BifVector[1] to BifVector[nphase]:  fold vector */
    staFunc->CopyVectorElements(v,0,V,0,nphase);
    staFunc->NormalizeVector(V);
    BifVector[0]=(VAR)_itmap;
    MemCopy(BifVector+1,V,sizeof(FloatHi)*nphase);
    stagenData->BifDataOutLen=sizeof(FloatHi)*(nphase+1);

/*  compute adjoint eigenvector */
    staFunc->TransposeMatrix(A,B);
    for (i=0; i<nphase; i++) B[i][i]-=1.0;
    
    if (staFunc->SngSolve(B,eps,W) != nphase-1) {
      sprintf(buf,"Limit point: a=?");
    } else {
      staFunc->NormalizeVectorRelative(W,V);

/*  compute second-order normal form coefficient  */
      if (stagenData->Der2)  {
        H=SymHess(x);
        for (f0=0.0,i=0; i<nphase; i++) {
           for (j=0; j<nphase; j++) {
              for (k=0; k<nphase; k++) 
                 f0+=W[i]*HessianElem(i,j,k)*V[j]*V[k];
           }
        }
        sprintf(buf,"Limit point: a=%g",f0);

      } else {
        for (i=0; i<=nappend; i++) v[nphase+i]=0.0;
        fp_DefRHS(x,F);
        f0=staFunc->ScalProd(W,F);
        staFunc->AddScaledVector(x,v,ddx,x1);
        fp_DefRHS(x1,F);
        f1=staFunc->ScalProd(W,F);
        staFunc->AddScaledVector(x,v,-ddx,x1);
        fp_DefRHS(x1,F);
        f2=staFunc->ScalProd(W,F);
        sprintf(buf,"Limit point: a=%g",(f1-2*f0+f2)/ddx/ddx);
      }
    }
    staFunc->DestroyMatrix(B);
    staFunc->DestroyMatrix(A);
    staFunc->DestroyVector(W);
    staFunc->DestroyVector(V);
    staFunc->DestroyVector(x1);
    staFunc->DestroyVector(F);
  } else {
    sprintf(buf,"Limit point (?)");
  }
  *msg=buf;
  return 0;
}
/* Period doubling (PD):  det(A+I) */
Entry(Int2) fp_TestPD(Vector x, Vector v, FloatHiPtr res) {
int i;
  *res=1.0;
  if (JacRHS(x, A0)) return 1;
  staFunc->CopyMatrixBlock(A0,0,0,AA,0,0,nphase,nphase);     
  for (i=0; i<nphase; i++) AA[i][i]+=1.0;
  if (staFunc->Decomp(AA)) {
      *res=0.0;
  } else {
      staFunc->GetMatrixData(AA,res,NULL);
  }
  test[2]=*res;
  return 0;
}
Entry(Int2) fp_ProcPD(Int2 Ind, Vector x, Vector v, CharPtr *msg) {
Matrix A,B;
Vector x0,x1,v0,V,W;
VAR eps=1.0e-2;
Int2 i;

  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Period-doubling (?)");
  else {
    x0=staFunc->CreateVector(nphase+1+nappend);
    x1=staFunc->CreateVector(nphase+1+nappend);
    v0=staFunc->CreateVector(nphase+1+nappend);
    V=staFunc->CreateVector(nphase);
    W=staFunc->CreateVector(nphase);
    A=staFunc->CreateMatrix(nphase,nphase);
    B=staFunc->CreateMatrix(nphase,nphase);

/*  compose Jacobian matrix */
    JacRHS(x,A0);
    staFunc->CopyMatrixBlock(A0,0,0,A,0,0,nphase,nphase);

/* compute eigenvector */
    staFunc->CopyMatrix(A,B);
    for (i=0; i<nphase; i++) B[i][i]+=1.0;
    if (staFunc->SngSolve(B,eps,V) != nphase-1) {
      sprintf(buf,"Period-doubling: critical eigenvector (?)");
      goto final;
    };
    staFunc->NormalizeVector(V);

/*  provide bifurcation data for pd->* starters */
/*  BifVector[0]:	                _itmap */
/*  BifVector[1] to BifVector[nphase]:  flip vector AV=-V */
/*  BifVector[nphase+1] to BifVector[2*nphase+1+nappend]: tangent vector to the fp-curve */
/*  BifVector[2*nphase+2+nappend] to BifVector[2*nphase+2+2*nappend]:  active[] */
    BifVector[0]=(VAR)_itmap;
    MemCopy(BifVector+1,V,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+nphase+1,v,sizeof(FloatHi)*(nphase+1+nappend));
    for (i=0; i<=nappend; i++) {
      BifVector[2*nphase+1+nappend+i+1]=(VAR)active[i];
    }
    stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+1+2*nappend+1+1);

/*  compute adjoint eigenvector */
    staFunc->TransposeMatrix(A,B);
    for (i=0; i<nphase; i++) B[i][i]+=1.0;
    if (staFunc->SngSolve(B,eps,W) != nphase-1) {
      sprintf(buf,"Period-doubling: adjoint critical eigenvector (?)");
      goto final;
    };
    staFunc->NormalizeVectorRelative(W,V);

/* compute the coefficient in the normal form x'=-x+cx^3 */
 
    staFunc->ClearVector(v0);
    staFunc->CopyVectorElements(V,0,v0,0,nphase);
    staFunc->CopyVector(x,x0);

{
Vector a,F,F1,F2;
VAR b,c;
Int2 i,j,k,l;
    a=staFunc->CreateVector(nphase);
    F=staFunc->CreateVector(nphase);
    F1=staFunc->CreateVector(nphase);
    F2=staFunc->CreateVector(nphase);

/* compute vector a=B(V,V) */
    if (stagenData->Der2) {
        H=SymHess(x);
        for (i=0; i<nphase; i++) {
           for (a[i]=0.0,j=0; j<nphase; j++) {
              for (k=0; k<nphase; k++) a[i]+=HessianElem(i,j,k)*V[j]*V[k];    /* a=B(V,V) */
           }
        }
    } else {
        staFunc->AddScaledVector(x0,v0,ddx,x1);
        fp_DefRHS(x1,a);
        fp_DefRHS(x0,F);
	staFunc->AddScaledVector(a,F,-2.0,a);
        staFunc->AddScaledVector(x0,v0,-ddx,x1);
        fp_DefRHS(x1,F);
	staFunc->AddScaledVector(a,F,1.0,a);
	staFunc->MultiplyVectorScalar(a,1.0/ddx/ddx,a);
    };

/* solve the system for a */
    staFunc->CopyMatrix(A,B);
    for (i=0; i<nphase; i++) B[i][i]-=1.0;
    staFunc->Decomp(B);
    staFunc->Solve(B,a);

/* compute F=B(V,a) */
    if (stagenData->Der2) {
      for (i=0; i<nphase; i++) {
        for (F[i]=0.0,j=0; j<nphase; j++) {
          for (k=0; k<nphase; k++) F[i]+=HessianElem(i,j,k)*V[j]*a[k];    
        }
      }
    } else {
      staFunc->AddScaledVector(V,a,1.0,F1);
      staFunc->CopyVectorElements(F1,0,v0,0,nphase);
      staFunc->AddScaledVector(x0,v0,0.5*ddx,x1);
      fp_DefRHS(x1,F1);
      staFunc->AddScaledVector(x0,v0,-0.5*ddx,x1);
      fp_DefRHS(x1,F);
      staFunc->AddScaledVector(F1,F,1.0,F2);

      staFunc->AddScaledVector(V,a,-1.0,F1);
      staFunc->CopyVectorElements(F1,0,v0,0,nphase);
      staFunc->AddScaledVector(x0,v0,0.5*ddx,x1);
      fp_DefRHS(x1,F1);
      staFunc->AddScaledVector(x0,v0,-0.5*ddx,x1);
      fp_DefRHS(x1,F);
      staFunc->AddScaledVector(F1,F,1.0,F1);
      staFunc->AddScaledVector(F2,F1,-1.0,F);
      staFunc->MultiplyVectorScalar(F,1.0/ddx/ddx,F);   
    }
/* compute c */
    if (stagenData->Der3) {
      D3=SymD3(x);
      for (c=0.0,i=0; i<nphase; i++) {
         for (j=0; j<nphase; j++) {
            for (k=0; k<nphase; k++) {
	       for (l=0; l<nphase; l++) {
	         c+=W[i]*D3Elem(i,j,k,l)*V[j]*V[k]*V[l];
	       }
	    }
         }
      }
    } else {
	staFunc->CopyVectorElements(V,0,v0,0,nphase);
        staFunc->AddScaledVector(x0,v0,3*dddx/2,x1);
        fp_DefRHS(x1,F1);
	c=staFunc->ScalProd(W,F1);
        staFunc->AddScaledVector(x0,v0,dddx/2,x1);
        fp_DefRHS(x1,F1);
	c-=3.0*staFunc->ScalProd(W,F1);
        staFunc->AddScaledVector(x0,v0,-dddx/2,x1);
        fp_DefRHS(x1,F1);
	c+=3.0*staFunc->ScalProd(W,F1);
        staFunc->AddScaledVector(x0,v0,-3*dddx/2,x1);
        fp_DefRHS(x1,F1);
	c-=staFunc->ScalProd(W,F1);
        c/=(dddx*dddx*dddx);   
     }
     b=staFunc->ScalProd(W,F);
 
    sprintf(buf,"Period-doubling: c=%g",-b/2+c/6);

    a=staFunc->DestroyVector(a);
    F=staFunc->DestroyVector(F);
    F1=staFunc->DestroyVector(F1);
    F2=staFunc->DestroyVector(F2);
}
final:
  x0=staFunc->DestroyVector(x0);
  x1=staFunc->DestroyVector(x1);
  v0=staFunc->DestroyVector(v0);
  V=staFunc->DestroyVector(V);
  W=staFunc->DestroyVector(W);
  A=staFunc->DestroyMatrix(A);
  B=staFunc->DestroyMatrix(B);
  }
  *msg=buf;
  return 0;
}

/* Neimark-Sacker (NS): det(A*A-I) */
Entry(Int2) fp_TestNS(Vector x, Vector v, FloatHiPtr res) {
  int i,j,p,q,r,s;

  *res=1.0;
  if (nphase == 1) goto end;

  if (JacRHS(x, A0)) return 1;
  staFunc->ClearMatrix(C);
  i=-1;
  for (p=1; p<nphase; p++) {
    for (q=0; q<p; q++) {
      i++;
      j=-1;
      for (r=1; r<nphase; r++) {
        for (s=0; s<r; s++) {
          j++;
          C[i][j]=A0[p][r]*A0[q][s]-A0[q][r]*A0[p][s];
	  if (i==j) C[i][j]-=1.0;
        }
      }
    }
  }
  if (staFunc->Decomp(C)) {
    *res=0.0;
  } else {
    staFunc->GetMatrixData(C,res,NULL);
  }
end:;
  test[3]=*res;
  return 0;
}

Entry(Int2) fp_ProcNS(Int2 Ind, Vector x, Vector v, CharPtr *msg) {

Int2 i,j,k,l,imin,jmin,Ret;
VAR d,s,s1,r,r1,h,beta,gamma,phi,lambda,lambda1,omega,lambda2,omega2,theta,eps=0.01;
Matrix AT,BB;
Vector V,Qr,Qi,Pr,Pi;
Vector aa,bb,cc,dd,rr;
   
  AT=staFunc->CreateMatrix(nphase,nphase);
  BB=staFunc->CreateMatrix(2*nphase,2*nphase);
  V=staFunc->CreateVector(2*nphase);
  Qr=staFunc->CreateVector(nphase);
  Qi=staFunc->CreateVector(nphase);
  Pr=staFunc->CreateVector(nphase);
  Pi=staFunc->CreateVector(nphase);
  aa=staFunc->CreateVector(nphase);
  bb=staFunc->CreateVector(nphase);
  cc=staFunc->CreateVector(nphase);
  dd=staFunc->CreateVector(nphase);
  rr=staFunc->CreateVector(nphase);

  BifVector[0]=(VAR)_itmap;
  stagenData->BifDataOutLen=sizeof(VAR);
  if (Ind == 0) {
    /*  compose Jacobian matrix */
    JacRHS(x,A0);
    staFunc->CopyMatrixBlock(A0,0,0,AA,0,0,nphase,nphase);
    staFunc->EigenValues(AA,Re,Im);
    for (i=0; i<nphase; i++) {
      M[2*i]=Re[i];
      M[2*i+1]=Im[i];
    }
    qsort(M,nphase,2*sizeof(VAR),fp_SortMult);
    for (i=0; i<nphase; i++) {
      Re[i]=M[2*i];
      Im[i]=M[2*i+1];
    }

    /* find the pair of multipliers on the unit circle */
    s=fabs(Re[0]*Re[1]-Im[0]*Im[1]-1.0)+fabs(Re[0]*Im[1]+Re[1]*Im[0]);
    imin=0;
    jmin=1;
    for (i=0; i<nphase; i++) {
      for (j=i+1; j<nphase; j++) {
        d=fabs(Re[i]*Re[j]-Im[i]*Im[j]-1.0)+fabs(Re[i]*Im[j]+Re[j]*Im[i]);
        if (d<s) {
          imin=i;
          jmin=j;
          s=d;
        }
      }
    }
    if (fabs(Im[imin])+fabs(Im[jmin])==0) {
/* Neutral saddle */
      lambda=Re[imin];
      lambda1=Re[jmin];
      /* compute real eigenvectors associated with lambda and lambda1=1/lambda */
      staFunc->CopyMatrixBlock(A0,0,0,AA,0,0,nphase,nphase);
      for (i=0; i<nphase; i++) AA[i][i]-=lambda;
      if (staFunc->SngSolve(AA,eps,Qr) != nphase-1) {
        sprintf(buf,"Neutral saddle (?) (cannot find eigenvector)");
        goto final;
      } 
      staFunc->CopyMatrixBlock(A0,0,0,AA,0,0,nphase,nphase);
      for (i=0; i<nphase; i++) AA[i][i]-=lambda1;
      if (staFunc->SngSolve(AA,eps,Qi) != nphase-1) {
        sprintf(buf,"Neutral saddle (?) (cannot find eigenvector)");
        goto final;
      } 
      staFunc->NormalizeVector(Qr);   
      staFunc->NormalizeVectorRelative(Qi,Qr);    /* ensures smallest possible angle */
      staFunc->NormalizeVector(Qi);

      /* provide bifurcation data for ns->ns starter */
      BifVector[0]=(VAR)_itmap;
      staFunc->AddScaledVector(Qr,Qi,1.0,Pr);
      staFunc->NormalizeVector(Pr);
      MemCopy(BifVector+1,Pr,sizeof(FloatHi)*nphase);
      staFunc->AddScaledVector(Qr,Qi,-1.0,Pr);
      staFunc->NormalizeVector(Pr);
      MemCopy(BifVector+nphase+1,Pr,sizeof(FloatHi)*nphase);
      BifVector[2*nphase+1]=(lambda+lambda1)/2;    /* kappa */
      stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+1+1);      
      sprintf(buf,"Neutral saddle: lambda=%g",lambda);

    } else {
/* Neimark-Sacker */
      if (fabs(Re[imin]) > 0.5) {
        theta=atan(fabs(Im[imin]/Re[imin]));
      } else {
        theta=PID2-atan(fabs(Re[imin]/Im[imin]));
      }
      if (Re[imin] < 0) theta=PI-theta;
      if (Im[imin] < 0) theta=-theta;
      lambda=cos(theta);
      omega=sin(theta);
      lambda2=cos(2*theta);
      omega2=sin(2*theta);

      /* compute real and imaginary part of the eigenvector */
      staFunc->CopyMatrixBlock(A0,0,0,BB,0,0,nphase,nphase);
      staFunc->CopyMatrixBlock(A0,0,0,BB,nphase,nphase,nphase,nphase);
      for (i=0; i<nphase; i++) {
        BB[i][i+nphase]=omega;
        BB[i][i]-=lambda;
        BB[i+nphase][i]=-omega;
        BB[i+nphase][i+nphase]-=lambda;
      }
      if (staFunc->SngSolve(BB,eps,V) != 2*nphase-2) {
        sprintf(buf,"Neimark-Sacker (?) (singular BB)");
        goto final;
      }
      staFunc->CopyVectorElements(V,0,Qr,0,nphase);
      staFunc->CopyVectorElements(V,nphase,Qi,0,nphase);
      /* normalize real and imaginary parts <Qr,Qr>=1, <Qr,Qi>=0 */ 
        d=staFunc->ScalProd(Qr,Qr);
        s=staFunc->ScalProd(Qi,Qi);
        r=staFunc->ScalProd(Qr,Qi);
/*
          beta=0.5*(d-s)/r;
          phi=0.5*atan(-1.0/beta);
*/
        if (r!=0) {
          beta=sqrt(fabs(4.0*r*r+(s-d)*(s-d)));
          if (2.0*fabs(2.0*r) > beta) {
            phi=PID2-atan(fabs((s-d)/(2.0*r)));
          } else {
            phi=atan(fabs(2.0*r/(s-d)));
          }
          if (r < 0) phi=PI-phi;
          if ((s-d) < 0) phi=-phi;
          phi/=2.0;
          for (i=0; i<nphase; i++) {
            V[i]=Qr[i]*cos(phi)-Qi[i]*sin(phi);
            V[i+nphase]=Qr[i]*sin(phi)+Qi[i]*cos(phi);
          }
        }
        staFunc->CopyVectorElements(V,0,Qr,0,nphase);
        d=staFunc->ScalProd(Qr,Qr);
        gamma=1/sqrt(d);
        staFunc->MultiplyVectorScalar(V,gamma,V);
        staFunc->CopyVectorElements(V,0,Qr,0,nphase);
        staFunc->CopyVectorElements(V,nphase,Qi,0,nphase);
/*
staFunc->PrintVector(Qr,"Qr=");
staFunc->PrintVector(Qi,"Qi=");
*/
      /* provide bifurcation data for ns->* starters */
       BifVector[0]=(VAR)_itmap;
       MemCopy(BifVector+1,V,sizeof(FloatHi)*2*nphase);
       BifVector[2*nphase+1]=cos(theta);  /* kappa */
       stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+1+2);

      /* compute real and imaginary part of the adjoint eigenvector */
       staFunc->CopyMatrixBlock(A0,0,0,AA,0,0,nphase,nphase);
       staFunc->TransposeMatrix(AA,AT);
       staFunc->ClearMatrix(BB);
       staFunc->CopyMatrixBlock(AT,0,0,BB,0,0,nphase,nphase);
       staFunc->CopyMatrixBlock(AT,0,0,BB,nphase,nphase,nphase,nphase);
       for (i=0; i<nphase; i++) {
         BB[i][i+nphase]=-omega;
         BB[i][i]-=lambda;
         BB[i+nphase][i]=omega;
         BB[i+nphase][i+nphase]-=lambda;
       }
       if (staFunc->SngSolve(BB,eps,V) != 2*nphase-2) {
         sprintf(buf,"Neimark-Sacker: theta=%g, a=?",theta);
         goto final;
       }
       staFunc->CopyVectorElements(V,0,Pr,0,nphase);
       staFunc->CopyVectorElements(V,nphase,Pi,0,nphase);

      /* normalize the adjoint vector <Pr,Qr>+<Pi,Qi>=1, <Pr,Qi>-<Pi,Qr>=0 */
       d=staFunc->ScalProd(Pr,Qr);
       s=staFunc->ScalProd(Pi,Qi);
       r=staFunc->ScalProd(Pr,Qi);
       h=staFunc->ScalProd(Pi,Qr);
/*       phi=atan((r-h)/(d+s)); */
       beta=sqrt(fabs((r-h)*(r-h)+(d+s)*(d+s)));
       if (2.0*fabs(d+s) > beta) {
         phi=atan(fabs((r-h)/(d+s)));
       } else {
         phi=PID2-atan(fabs((d+s)/(r-h)));
       }
       if ((r-h) < 0) phi=PI-phi;
       if ((d+s) < 0) phi=-phi;

       for (i=0; i<nphase; i++) {
         V[i]=Pr[i]*cos(phi)-Pi[i]*sin(phi);
         V[i+nphase]=Pr[i]*sin(phi)+Pi[i]*cos(phi);
       }
       staFunc->CopyVectorElements(V,0,Pr,0,nphase);
       staFunc->CopyVectorElements(V,nphase,Pi,0,nphase);
       r=staFunc->ScalProd(Pr,Qr);
       h=staFunc->ScalProd(Pi,Qi);
       beta=1.0/(r+h);
       staFunc->MultiplyVectorScalar(V,beta,V);
       staFunc->MultiplyVectorScalar(Pr,beta,Pr);
       staFunc->MultiplyVectorScalar(Pi,beta,Pi);

/*
staFunc->PrintVector(Pr,"Pr=");
staFunc->PrintVector(Pi,"Pi=");
*/
      /* compute the cubic normal form coefficient a(0) */
      /* quadratic terms */
      if (stagenData->Der2)  {
        H=SymHess(x);
	for (i=0; i<nphase; i++) {
           for (s=0.0,r=0.0,h=0.0,j=0; j<nphase; j++) {
              for (k=0; k<nphase; k++) {
                 s+=HessianElem(i,j,k)*Qr[j]*Qr[k];
                 r+=HessianElem(i,j,k)*Qi[j]*Qi[k];
                 h+=HessianElem(i,j,k)*Qr[j]*Qi[k];
	      }
           }
	   aa[i]=s;	
	   bb[i]=r;	
	   cc[i]=-2.0*h;	
	}
      } else {
        staFunc->ClearVector(v0);
        staFunc->CopyVectorElements(Qr,0,v0,0,nphase);
        fp_DefRHS(x,aa);
        staFunc->AddScaledVector(x,v0,ddx,x1);
        fp_DefRHS(x1,bb);
        staFunc->AddScaledVector(x,v0,-ddx,x1);
        fp_DefRHS(x1,cc);
	staFunc->AddScaledVector(bb,aa,-2.0,rr);
	staFunc->AddScaledVector(rr,cc,1.0,aa);
	staFunc->MultiplyVectorScalar(aa,1.0/ddx/ddx,aa);

        staFunc->CopyVectorElements(Qi,0,v0,0,nphase);
        fp_DefRHS(x,rr);
        staFunc->AddScaledVector(x,v0,ddx,x1);
        fp_DefRHS(x1,bb);
        staFunc->AddScaledVector(x,v0,-ddx,x1);
        fp_DefRHS(x1,cc);
	staFunc->AddScaledVector(bb,rr,-2.0,rr);
	staFunc->AddScaledVector(rr,cc,1.0,bb);
	staFunc->MultiplyVectorScalar(bb,1.0/ddx/ddx,bb);

        staFunc->AddScaledVector(Qr,Qi,1.0,rr);
        staFunc->CopyVectorElements(rr,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,rr);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,cc);
	staFunc->AddScaledVector(rr,cc,1.0,dd);

        staFunc->AddScaledVector(Qr,Qi,-1.0,rr);
        staFunc->CopyVectorElements(rr,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,rr);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,cc);
	staFunc->AddScaledVector(rr,cc,1.0,rr);
        staFunc->AddScaledVector(dd,rr,-1.0,cc);
	staFunc->MultiplyVectorScalar(cc,-2.0/ddx/ddx,cc);   /* cc=-2*c */
      }
            
      staFunc->CopyMatrixBlock(A0,0,0,AA,0,0,nphase,nphase);
      for (i=0; i<nphase; i++) AA[i][i]-=1.0;
      staFunc->AddScaledVector(aa,bb,1.0,rr);
      Ret=staFunc->Decomp(AA);
      if (Ret != 0) {
        sprintf(buf,"Neimark-Sacker: theta=%g, a=? (singular AA)",theta);
	goto final;
      }
      staFunc->Solve(AA,rr);
      staFunc->MultiplyVectorScalar(rr,-2.0,rr);		/* rr=-2*r */

      staFunc->AddScaledVector(bb,aa,-1.0,dd);
      staFunc->CopyVectorElements(dd,0,V,0,nphase);
      staFunc->CopyVectorElements(cc,0,V,nphase,nphase);
      staFunc->ClearMatrix(BB);
      staFunc->CopyMatrixBlock(A0,0,0,BB,0,0,nphase,nphase);
      staFunc->CopyMatrixBlock(A0,0,0,BB,nphase,nphase,nphase,nphase);
      for (i=0; i<nphase; i++) {
        BB[i][i+nphase]=omega2;
        BB[i][i]-=lambda2;
        BB[i+nphase][i]=-omega2;
        BB[i+nphase][i+nphase]-=lambda2;
      }
      Ret=staFunc->Decomp(BB);
      if (Ret != 0) {
        sprintf(buf,"Neimark-Sacker: theta=%g, a=? (singular BB)",theta);
	goto final;
      }
      staFunc->Solve(BB,V);
      staFunc->CopyVectorElements(V,0,aa,0,nphase);          	/* aa=Sr */
      staFunc->CopyVectorElements(V,nphase,bb,0,nphase);	/* bb=Si */

      if (stagenData->Der2)  {
        for (s=0.0,i=0; i<nphase; i++) {
           for (j=0; j<nphase; j++) {
              for (k=0; k<nphase; k++) 
	         s+=HessianElem(i,j,k)*(
		   Pr[i]*(Qr[j]*rr[k]+Qr[j]*aa[k]+Qi[j]*bb[k])+
		   Pi[i]*(Qi[j]*rr[k]+Qr[j]*bb[k]-Qi[j]*aa[k]));
           }
        }
        for (s1=0.0,i=0; i<nphase; i++) {
           for (j=0; j<nphase; j++) {
              for (k=0; k<nphase; k++) 
	         s1+=HessianElem(i,j,k)*(
		   Pr[i]*(Qi[j]*rr[k]+Qr[j]*bb[k]-Qi[j]*aa[k])-
		   Pi[i]*(Qr[j]*rr[k]+Qr[j]*aa[k]+Qi[j]*bb[k]));
           }
        }
      } else {
        staFunc->ClearVector(v0);

/*  -2*<Pr,B(Qr,r)>  */
        staFunc->AddScaledVector(Qr,rr,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pr,dd);

        staFunc->AddScaledVector(Qr,rr,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
/*  -2*<Pr,B(Qi,r)>  */
        staFunc->AddScaledVector(Qi,rr,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pr,dd);

        staFunc->AddScaledVector(Qi,rr,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pr,dd);
/*  -2*<Pi,B(Qi,r)>  */
        staFunc->AddScaledVector(Qi,rr,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);

        staFunc->AddScaledVector(Qi,rr,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);
/*  2*<Pi,B(Qr,r)>  */
        staFunc->AddScaledVector(Qr,rr,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pi,dd);

        staFunc->AddScaledVector(Qr,rr,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pi,dd);
/*  <Pr,B(Qr,Sr)>  */
        staFunc->AddScaledVector(Qr,aa,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pr,dd);

        staFunc->AddScaledVector(Qr,aa,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
/*  <Pr,B(Qr,Si)>  */
        staFunc->AddScaledVector(Qr,bb,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pr,dd);

        staFunc->AddScaledVector(Qr,bb,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pr,dd);
/*  <Pr,B(Qi,Si)>  */
        staFunc->AddScaledVector(Qi,bb,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pr,dd);

        staFunc->AddScaledVector(Qi,bb,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
/* -<Pr,B(Qi,Sr)>  */
        staFunc->AddScaledVector(Qi,aa,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pr,dd);

        staFunc->AddScaledVector(Qi,aa,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pr,dd);
/*  <Pi,B(Qr,Si)>  */
        staFunc->AddScaledVector(Qr,bb,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);

        staFunc->AddScaledVector(Qr,bb,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);
/*  -<Pi,B(Qr,Sr)>  */
        staFunc->AddScaledVector(Qr,aa,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pi,dd);

        staFunc->AddScaledVector(Qr,aa,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pi,dd);
/*  -<Pi,B(Qi,Sr)>  */
        staFunc->AddScaledVector(Qi,aa,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);

        staFunc->AddScaledVector(Qi,aa,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);
	s/=(ddx*ddx); 
/*  -<Pi,B(Qi,Si)>  */
        staFunc->AddScaledVector(Qi,bb,1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pi,dd);

        staFunc->AddScaledVector(Qi,bb,-1.0,dd);
        staFunc->CopyVectorElements(dd,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        fp_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pi,dd);
	s1/=(ddx*ddx); 
      }
     /* cubic terms */
     if (stagenData->Der3) {
        D3=SymD3(x);
        for (r=0.0,i=0; i<nphase; i++) {
           for (j=0; j<nphase; j++) {
              for (k=0; k<nphase; k++) {
	         for (l=0; l<nphase; l++) {
	           r+=D3Elem(i,j,k,l)*(Pr[i]*Qr[j]*(Qr[k]*Qr[l]+Qi[k]*Qi[l])+
		                       Pi[i]*Qi[l]*(Qi[j]*Qi[k]+Qr[j]*Qr[k]));
	         }
	      }
           }
        }	
        for (r1=0.0,i=0; i<nphase; i++) {
           for (j=0; j<nphase; j++) {
              for (k=0; k<nphase; k++) {
	         for (l=0; l<nphase; l++) {
	           r1+=D3Elem(i,j,k,l)*(Pr[i]*Qi[j]*(Qr[k]*Qr[l]+Qi[k]*Qi[l])-
		                        Pi[i]*Qr[l]*(Qi[j]*Qi[k]+Qr[j]*Qr[k]));
	         }
	      }
           }
        }	
     } else {
        staFunc->ClearVector(v0);

/* (2/3)<Pr,C(Qr,Qr,Qr)> */	
	staFunc->CopyVectorElements(Qr,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r=(2.0/3.0)*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        fp_DefRHS(x1,dd);
	r-=2.0*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        fp_DefRHS(x1,dd);
	r+=2.0*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r-=(2.0/3.0)*staFunc->ScalProd(Pr,dd);
/* (2/3)<Pr,C(Qi,Qi,Qi)> */	
	staFunc->CopyVectorElements(Qi,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1=(2.0/3.0)*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1-=2.0*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1+=2.0*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1-=(2.0/3.0)*staFunc->ScalProd(Pr,dd);
	
/* (2/3)<Pi,C(Qi,Qi,Qi)> */
	staFunc->CopyVectorElements(Qi,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r+=(2.0/3.0)*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        fp_DefRHS(x1,dd);
	r-=2.0*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        fp_DefRHS(x1,dd);
	r+=2.0*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r-=(2.0/3.0)*staFunc->ScalProd(Pi,dd);
/* -(2/3)<Pi,C(Qr,Qr,Qr)> */
	staFunc->CopyVectorElements(Qr,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1-=(2.0/3.0)*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1+=2.0*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1-=2.0*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1+=(2.0/3.0)*staFunc->ScalProd(Pi,dd);
	
/* (1/6) <Pr+Pi,C(Qr+Qi,Qr+Qi,Qr+Qi)> */
        staFunc->AddScaledVector(Qr,Qi,1.0,aa);
        staFunc->CopyVectorElements(aa,0,v0,0,nphase);
        staFunc->AddScaledVector(Pr,Pi,1.0,bb);
        staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r+=staFunc->ScalProd(bb,dd)/6;
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        fp_DefRHS(x1,dd);
	r-=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        fp_DefRHS(x1,dd);
	r+=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r-=staFunc->ScalProd(bb,dd)/6;
/* (1/6) <Pr-Pi,C(Qr+Qi,Qr+Qi,Qr+Qi)> */
        staFunc->AddScaledVector(Qr,Qi,1.0,aa);
        staFunc->CopyVectorElements(aa,0,v0,0,nphase);
        staFunc->AddScaledVector(Pr,Pi,-1.0,bb);
        staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1+=staFunc->ScalProd(bb,dd)/6;
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1-=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1+=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1-=staFunc->ScalProd(bb,dd)/6;

/* (1/6) <Pr-Pi,C(Qr-Qi,Qr-Qi,Qr-Qi)> */
        staFunc->AddScaledVector(Qr,Qi,-1.0,aa);
        staFunc->CopyVectorElements(aa,0,v0,0,nphase);
        staFunc->AddScaledVector(Pr,Pi,-1.0,bb);
        staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r+=staFunc->ScalProd(bb,dd)/6;
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        fp_DefRHS(x1,dd);
	r-=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        fp_DefRHS(x1,dd);
	r+=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r-=staFunc->ScalProd(bb,dd)/6;
        r/=(dddx*dddx*dddx);                                         
/* -(1/6) <Pr+Pi,C(Qr-Qi,Qr-Qi,Qr-Qi)> */
        staFunc->AddScaledVector(Qr,Qi,-1.0,aa);
        staFunc->CopyVectorElements(aa,0,v0,0,nphase);
        staFunc->AddScaledVector(Pr,Pi,1.0,bb);
        staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1-=staFunc->ScalProd(bb,dd)/6;
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1+=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1-=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        fp_DefRHS(x1,dd);
	r1+=staFunc->ScalProd(bb,dd)/6;
        r1/=(dddx*dddx*dddx);                                         
     }	       
     sprintf(buf,"Neimark-Sacker: theta=%g, a=%g",theta,(lambda*(s+r)+omega*(s1+r1))/2);
    }
  } else {
    sprintf(buf,"Neimark-Sacker or neutral saddle (?)");
  }
final:;
  AT=staFunc->DestroyMatrix(AT);
  BB=staFunc->DestroyMatrix(BB);
  V=staFunc->DestroyVector(V);
  Qr=staFunc->DestroyVector(Qr);
  Qi=staFunc->DestroyVector(Qi);
  Pr=staFunc->DestroyVector(Pr);
  Pi=staFunc->DestroyVector(Pi);
  aa=staFunc->DestroyVector(aa);
  bb=staFunc->DestroyVector(bb);
  cc=staFunc->DestroyVector(cc);
  dd=staFunc->DestroyVector(dd);
  rr=staFunc->DestroyVector(rr);
  *msg=buf;
  return 0;
}


/**********************/
/* Trivial Adaptation */
Entry(Int2) fp_Adapter(Vector X, Vector V) {
  return (0);
}