/*
   m__ns2.c
   Neimarck Sacker(neutral saddle) point--> Neimarck Sacker(neutral saddle) Curve 
   using double-bordering (A^2 - 2kA + I) method
   December 17, 1998
*/

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

#include "stagen.h"

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

#include "m__ns2.h"
#include "_conti.h"	/* continuer'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

#define _eps     staData->eps
#define _dx      staData->dx
#define _branch  staData->branch   

#define PI	4.0*atan(1.0)
#define PID2	2.0*atan(1.0)   

Local(StagenDataPtr) stagenData=NULL;	/* global copy of Starter's parameter passed by content */
Local(m__ns_DataPtr) staData=NULL;	/* global copy of our own data area */
Local(StdLinAlgDataPtr) staFunc=NULL;	/* global copy of pointer standard linear algebra struct */
Local(UtilitiesDataPtr) staUtil=NULL;	/* global copy of pointer to utility functions struct */
Local(Int2Ptr) active=NULL;        	/* 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) x0=NULL;			/* approximate first point */
Local(Vector) x1=NULL;			/* current 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 */

#include "_supdiff.c"

/* Working space for defining and test functions  */

Local(Matrix) A0,A,B,C,C1,C2,AA,AT,AH,BB,BP,BT,D,DE,G,GM,BP1,BT1,JM1,JM2,TVM,TVMT,HM,GM1,GM2,Jac1,Jac2;
Local(Vector) vn,xp,xp1,V,W,F,F1,Vp,L,VV,VB,V2,W2,Q1,Q2,Qr,Qi,PP1,PP2,Pr,Pi,aa,bb,cc,dd,rr;
Local(Vector) G0,G1,G2,G3,TV1,TV2,HV,HV1,HV2,HW1,HW2,VB1,VB2,V1,W1,RL,K1,K2,K3,OK1,OK2,PIVOT;
Local(Vector) Re=NULL,Im=NULL,Mod=NULL,Arg=NULL;
Local(VAR) test[8],dx,ddx,dddx;
Local(Int2) m1,n1,m2,n2;

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

/* Initialize function common for Starter and Decoder */
Local(Int2) Init(StagenDataPtr stagenptr) {
Int2 i,m;
  /* Set local pointers to data area and functions array */
  stagenData=stagenptr;
  staData=(m__ns_DataPtr)stagenptr->StaPtr;
  /* Use standard linear 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+2));    
  if (staUtil->CheckParameters(staData->Active,npar,active,2+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;
  }
  /* Create a copy parameters */
  m=nphase*(nphase-1)/2;
  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 eigenvalues */
    M=MemNew(sizeof(FloatHi)*4*nphase);
    Re=staFunc->CreateVector(nphase);
    Im=staFunc->CreateVector(nphase);
    Mod=staFunc->CreateVector(nphase);
    Arg=staFunc->CreateVector(nphase); 
    StdFun=staFunc->CreateVector(1+2*nphase); /* L2_NORM, MIN&MAX */       
    /* Allocate memory for working space used in defining and test functions */
    A0=staFunc->CreateMatrix(nphase,nphase+2+nappend);
    A=staFunc->CreateMatrix(nphase,nphase);
    B=staFunc->CreateMatrix(nphase,nphase);
    BP=staFunc->CreateMatrix(nphase+2,nphase+2);
    BT=staFunc->CreateMatrix(nphase+2,nphase+2);
    BP1=staFunc->CreateMatrix(nphase+1,nphase+1);
    BT1=staFunc->CreateMatrix(nphase+1,nphase+1);
    
    V1=staFunc->CreateVector(nphase+2);
    W1=staFunc->CreateVector(nphase+2);    
    V2=staFunc->CreateVector(nphase+2);
    W2=staFunc->CreateVector(nphase+2);
    VB=staFunc->CreateVector(nphase);
    VB1=staFunc->CreateVector(nphase+2);
    VB2=staFunc->CreateVector(nphase+2);    
    VV=staFunc->CreateVector(2*nphase);
    Q1=staFunc->CreateVector(nphase+2);
    Q2=staFunc->CreateVector(nphase+2);
    PP1=staFunc->CreateVector(nphase+2);
    PP2=staFunc->CreateVector(nphase+2);
    RL=staFunc->CreateVector(m+1);
    K1=staFunc->CreateVector(nphase);
    K2=staFunc->CreateVector(nphase);
    K3=staFunc->CreateVector(nphase);
    OK1=staFunc->CreateVector(m+2);
    OK2=staFunc->CreateVector(m+2);
    
    G0=staFunc->CreateVector(nphase+3+nappend);
    G1=staFunc->CreateVector(nphase+3+nappend);
    G2=staFunc->CreateVector(nphase+3+nappend);
    G3=staFunc->CreateVector(nphase+3+nappend);
    JM1=staFunc->CreateMatrix(nphase+2,nphase+2);
    JM2=staFunc->CreateMatrix(nphase+3,nphase+3);  
    TV1=staFunc->CreateVector(nphase+2);   
    TV2=staFunc->CreateVector(nphase+3);
    TVM=staFunc->CreateMatrix(nphase+3,3);
    TVMT=staFunc->CreateMatrix(3,nphase+3);
    HM=staFunc->CreateMatrix(3,3);
    HV=staFunc->CreateVector(3);
    HV1=staFunc->CreateVector(nphase);
    HV2=staFunc->CreateVector(nphase);
    HW1=staFunc->CreateVector(nphase);
    HW2=staFunc->CreateVector(nphase);
    GM1=staFunc->CreateMatrix(nphase+3+nappend,4);
    GM2=staFunc->CreateMatrix(nphase+3+nappend,4);
    G=staFunc->CreateMatrix(2,2);
    AA=staFunc->CreateMatrix(nphase,nphase);
    BB=staFunc->CreateMatrix(2*nphase,2*nphase);
    C=staFunc->CreateMatrix(nphase+nappend,nphase+3+nappend);
    C1=staFunc->CreateMatrix(nphase,nphase+3+nappend);
    C2=staFunc->CreateMatrix(2,nphase+3+nappend);
    x0=staFunc->CreateVector(nphase+3+nappend);
    x1=staFunc->CreateVector(nphase+3+nappend);
    v0=staFunc->CreateVector(nphase+3+nappend);
    vn=staFunc->CreateVector(nphase);
    xp=staFunc->CreateVector(nphase+2+nappend);
    xp1=staFunc->CreateVector(nphase+2+nappend);
    V=staFunc->CreateVector(nphase);
    W=staFunc->CreateVector(nphase);
    L=staFunc->CreateVector(nphase);
    Qr=staFunc->CreateVector(nphase);
    Qi=staFunc->CreateVector(nphase);
    Vp=staFunc->CreateVector(nphase+2+nappend);
    BifVector=MemNew(sizeof(FloatHi)*(4*nphase+3));
    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);
      }
    }
    if (nphase==1) staData->ActiveSing[0]=FALSE;
    if ((nphase>1) && staData->ActiveSing[0]) {  /* Degenerate Neimark-Sacker detection ON*/
      AT=staFunc->CreateMatrix(nphase,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);
    }
    if (nphase<=3) staData->ActiveSing[4]=FALSE;
    if ((nphase>3) && staData->ActiveSing[4]) {  /* Double NS detection ON */
      AT=staFunc->CreateMatrix(nphase,nphase);
      AH=staFunc->CreateMatrix(nphase-2,nphase-2);
      DE=staFunc->CreateMatrix((nphase-2)*(nphase-3)/2,(nphase-2)*(nphase-3)/2);
      F1=staFunc->CreateVector(nphase);
    }
  }
  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(Int2) JacRHS (Vector x, Matrix J);

Local(void) Term (Boolean OK) {
  stagenData->BifDataOutLen=0;
  if (OK&&P&&(stagenData->Reason==RF_CLEAN)) {
/* saving global vectors V1,V2,W1 and W2 used in defining functions */
    stagenData->ContDataStaLen=4*(nphase+2)*sizeof(FloatHi);
    stagenData->ContDataStaPtr=MemNew(stagenData->ContDataStaLen);
    MemCopy(stagenData->ContDataStaPtr,V1,(nphase+2)*sizeof(FloatHi));
    MemCopy(stagenData->ContDataStaPtr+(nphase+2),W1,(nphase+2)*sizeof(FloatHi));
    MemCopy(stagenData->ContDataStaPtr+2*(nphase+2),V2,(nphase+2)*sizeof(FloatHi));
    MemCopy(stagenData->ContDataStaPtr+3*(nphase+2),W2,(nphase+2)*sizeof(FloatHi));    

    stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+2);
    BifVector[0]=itmap;
    MemCopy(BifVector+1,VB1,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+nphase+1,VB2,sizeof(FloatHi)*nphase);
    BifVector[2*nphase+1]=x1[nphase+2+nappend];
    stagenData->ProcessPoint(stagenData->ProcFuncNum, x1, "Last point");
  } else {
    stagenData->ContDataStaLen=0;
    stagenData->BifDataOutLen=0;
  } 
  x0=staFunc->DestroyVector(x0);
  x1=staFunc->DestroyVector(x1);
  v0=staFunc->DestroyVector(v0);
  active=MemFree(active); 
  if (P) {
    P=MemFree(P);
    X=MemFree(X);
    M=MemFree(M);
    Re=staFunc->DestroyVector(Re);
    Im=staFunc->DestroyVector(Im);
    Mod=staFunc->DestroyVector(Mod);
    Arg=staFunc->DestroyVector(Arg); 
    StdFun=staFunc->DestroyVector(StdFun);              
    A=staFunc->DestroyMatrix(A);
    A0=staFunc->DestroyMatrix(A0);
    B=staFunc->DestroyMatrix(B);
    BP=staFunc->DestroyMatrix(BP);
    BT=staFunc->DestroyMatrix(BT);
    BP1=staFunc->DestroyMatrix(BP1);
    BT1=staFunc->DestroyMatrix(BT1);
    V1=staFunc->DestroyVector(V1);
    W1=staFunc->DestroyVector(W1);    
    V2=staFunc->DestroyVector(V2);
    W2=staFunc->DestroyVector(W2);
    Q1=staFunc->DestroyVector(Q1);
    Q2=staFunc->DestroyVector(Q2);
    PP1=staFunc->DestroyVector(PP1);
    PP2=staFunc->DestroyVector(PP2);
    AA=staFunc->DestroyMatrix(AA);
    BB=staFunc->DestroyMatrix(BB);
    VV=staFunc->DestroyVector(VV);
    VB=staFunc->DestroyVector(VB);
    VB1=staFunc->DestroyVector(VB1);
    VB2=staFunc->DestroyVector(VB2);    
    RL=staFunc->DestroyVector(RL);    
    K1=staFunc->DestroyVector(K1);
    K2=staFunc->DestroyVector(K2);
    K3=staFunc->DestroyVector(K3);
    OK1=staFunc->DestroyVector(OK1);
    OK2=staFunc->DestroyVector(OK2); 
    G0=staFunc->DestroyVector(G0);
    G1=staFunc->DestroyVector(G1);      
    G2=staFunc->DestroyVector(G2);      
    G3=staFunc->DestroyVector(G3);      
    JM1=staFunc->DestroyMatrix(JM1);
    JM2=staFunc->DestroyMatrix(JM2);
    TV1=staFunc->DestroyVector(TV1);
    TV2=staFunc->DestroyVector(TV2);
    TVM=staFunc->DestroyMatrix(TVM);
    TVMT=staFunc->DestroyMatrix(TVMT);
    HM=staFunc->DestroyMatrix(HM);
    HV=staFunc->DestroyVector(HV);
    HV1=staFunc->DestroyVector(HV1);
    HV2=staFunc->DestroyVector(HV2);
    HW1=staFunc->DestroyVector(HW1);
    HW2=staFunc->DestroyVector(HW2);
    GM1=staFunc->DestroyMatrix(GM1);
    GM2=staFunc->DestroyMatrix(GM2);
    G=staFunc->DestroyMatrix(G);
  
    C=staFunc->DestroyMatrix(C);
    C1=staFunc->DestroyMatrix(C1);
    C2=staFunc->DestroyMatrix(C2);
    vn=staFunc->DestroyVector(vn);
    xp=staFunc->DestroyVector(xp);
    xp1=staFunc->DestroyVector(xp1);
    V=staFunc->DestroyVector(V);
    W=staFunc->DestroyVector(W);
    Vp=staFunc->DestroyVector(Vp);
    L=staFunc->DestroyVector(L);
    Qr=staFunc->DestroyVector(Qr);
    Qi=staFunc->DestroyVector(Qi);
    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);
      }
    }
    if ((nphase>1) && staData->ActiveSing[0]) {  /* Degenerate NS detection ON*/
      AT=staFunc->DestroyMatrix(AT);
      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);
    }
    if ((nphase>3) && staData->ActiveSing[4]) {  /* Double NS detection ON*/
      AT=staFunc->DestroyMatrix(AT);
      AH=staFunc->DestroyMatrix(AH);
      DE=staFunc->DestroyMatrix(DE);
      F1=staFunc->DestroyVector(F1);
    }
  }
} 

Local(Int2) ns2_JacDefFixedPoint (Vector x, Matrix J);
Local(Int2) JacDefEig(Vector x, Matrix J);
Entry(Int2) ns2_JacDefNS (Vector x, Matrix J);
Entry(Int2) ns2_DefNS (Vector x, Vector J);

#ifdef WIN32
static int _USERENTRY ns2_SortMult(const void *f, const void *s) {
#else
Entry(int) ns2_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;
}              

/************/
/* Starters */
/* ns -> ns2 Starter */
Entry(Int2) ns_ns2_Starter(StagenDataPtr stagenptr) {
  Int2 i,j,k,kp1,l,Ret,p,q;
  Vector v1,v2;
  Matrix D,DE,S;
  VAR t,scal;
  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) {       /* extension */
    MemCopy(V1,stagenptr->ContDataStaPtr,sizeof(FloatHi)*(nphase+2));
    MemCopy(W1,stagenptr->ContDataStaPtr+(nphase+2),sizeof(FloatHi)*(nphase+2));
    MemCopy(V2,stagenptr->ContDataStaPtr+2*(nphase+2),sizeof(FloatHi)*(nphase+2));    
    MemCopy(W2,stagenptr->ContDataStaPtr+3*(nphase+2),sizeof(FloatHi)*(nphase+2));    
    MemFree(stagenptr->ContDataStaPtr);
    stagenptr->ContDataStaLen=0;
    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) x0[nphase+i]=staData->P[active[i]];                             
    x0[nphase+2+nappend]=BifVector[2*nphase+1]; 
  } else { /* forward or backward */
    v1=staFunc->CreateVector(nphase+3+nappend);
    v2=staFunc->CreateVector(nphase+3+nappend);
    D=staFunc->CreateMatrix(nphase+2+nappend,nphase+3+nappend);
    DE=staFunc->CreateMatrix(nphase+3+nappend,nphase+3+nappend);

    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) x0[nphase+i]=staData->P[active[i]];
    JacRHS(x0,C1);

    if (stagenData->BifDataInOk) {
/* Message(MSG_OK,"Eigenvectors exist"); */
/* bifurcation data exist */
/*  BifVector[0]:                      _itmap */
/*  BifVector[1] to [2*nphase]:        NS-vectors */
/*  BifVector[2*nphase+1]:	       kappa */
      MemCopy(VB1,stagenData->BifDataInPtr+1,sizeof(FloatHi)*nphase); 
      MemCopy(VB2,stagenData->BifDataInPtr+nphase+1,sizeof(FloatHi)*nphase); 
      x0[nphase+2+nappend]=stagenData->BifDataInPtr[2*nphase+1]; 
    } else { 
     /* bifurcation data to be computed */
     /* Message(MSG_OK,"Eigenvectors have to be computed");  */
      Ret=0;
{int i,j,imin,jmin;
 double r,s,d,phi,beta,gamma,lambda,lambda1,omega,theta;
    JacRHS(x0,C1);
    staFunc->CopyMatrixBlock(C1,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),ns2_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(C1,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) {
        staUtil->ErrMessage("Unable to find eigenvector at the initial point");
        Ret=1;
        goto polufinal;
      } 
      staFunc->CopyMatrixBlock(C1,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) {
        staUtil->ErrMessage("Unable to find eigenvector at the initial point");
        Ret=1;
        goto polufinal;
      } 
      staFunc->NormalizeVector(Qr);   
      staFunc->NormalizeVectorRelative(Qi,Qr);
      staFunc->NormalizeVector(Qi);
      
      staFunc->AddScaledVector(Qr,Qi,1.0,V);
      staFunc->NormalizeVector(V);
      staFunc->AddScaledVector(Qr,Qi,-1.0,L);             
      staFunc->NormalizeVector(L); 
      x0[nphase+2+nappend]=-lambda*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);
         
      /* compute real and imaginary part of the eigenvector */
      staFunc->ClearMatrix(BB);
      staFunc->CopyMatrixBlock(C1,0,0,BB,0,0,nphase,nphase);
      staFunc->CopyMatrixBlock(C1,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,VV) != 2*nphase-2) {
        staUtil->ErrMessage("Unable to find eigenvectors at the initial point");
        Ret=1;
        goto polufinal;
      }
      staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
      staFunc->CopyVectorElements(VV,nphase,Qi,0,nphase);

      /* normalize real and imaginary parts */
      d=staFunc->ScalProd(Qr,Qr);
      s=staFunc->ScalProd(Qi,Qi);
      r=staFunc->ScalProd(Qr,Qi);
      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++) {
          VV[i]=Qr[i]*cos(phi)-Qi[i]*sin(phi);
          VV[i+nphase]=Qr[i]*sin(phi)+Qi[i]*cos(phi);
        }
      }

      staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
      d=staFunc->ScalProd(Qr,Qr);
      gamma=1/sqrt(d);
      staFunc->MultiplyVectorScalar(VV,gamma,VV);
      staFunc->CopyVectorElements(VV,0,V,0,nphase);
      staFunc->CopyVectorElements(VV,nphase,L,0,nphase);
      x0[nphase+2+nappend]=cos(theta);
      staFunc->NormalizeVector(L);
    }
polufinal:
      if (Ret) {
        v1=staFunc->DestroyVector(v1);
        v2=staFunc->DestroyVector(v2);
        D=staFunc->DestroyMatrix(D);
        DE=staFunc->DestroyMatrix(DE);
        Term(FALSE);
        return 1;
      }
    }
}
  if (JacRHS(x0,C1)) return 1;  /* for safety */
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  staFunc->MultiplyMatrixMatrix(A,A,BP); 
  staFunc->AddScaledMatrix(BP,A,-2.0*x0[nphase+2+nappend],BP);
  for (i=0;i<nphase;i++) {
   BP[i][i]+=1.0;
  }
/* searching 2 rows with maximal pivots */  
  if (nphase >= 2) {
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) {
      BP1[i][nphase]=i;
      BP1[nphase][i]=i;
    } 
 /* Calculation of PIVOTS - complete pivoting */
    for (k=0;k<nphase-2;k++) {
      t=0.0;
      for (i=k;i<nphase;i++) {
	for (j=k;j<nphase;j++) {
    	  if (fabs(BP1[i][j]) > t) {
            p = i;
            q = j;
            t = fabs(BP1[i][j]);
          }           
	}
      }      
      if (k!=q) { /* switch colom k with colom  q */        
        for (i=0;i<=nphase;i++) {
          t=BP1[i][k];
          BP1[i][k]=BP1[i][q];
          BP1[i][q]=t;
	}
      }
      if (k!=p) { /* switch row k with row p */       
        for (i=0;i<=nphase;i++) {
          t=BP1[k][i];         
          BP1[k][i]=BP1[p][i];
          BP1[p][i]=t;
	}
      }        
      if (fabs(BP1[k][k]) == 0.0) {
        staUtil->ErrMessage("(A^2 + k I) has rank defect at least three");
        Term(FALSE);
        return 1;
      }
      kp1=k+1;
      for (i=kp1;i<nphase;i++) { /* perform elimination */
        t=-BP1[i][k]/BP1[k][k];
        for (j=kp1;j<nphase;j++) {
          BP1[i][j]+=t*BP1[k][j];
	}
      }
    } 
    p=(Int2)BP1[nphase-2][nphase];
    q=(Int2)BP1[nphase-1][nphase];   
    BP[p][nphase]=1.0;
    BP[q][nphase+1]=1.0;
    p=(Int2)BP1[nphase][nphase-2];
    q=(Int2)BP1[nphase][nphase-1];   
    BP[nphase][p]=1.0;
    BP[nphase+1][q]=1.0;    

/* solve systems to obtain new values for V1 and V2 */
    staFunc->ClearVector(V1);
    staFunc->ClearVector(V2);
    staFunc->ClearVector(W1);
    staFunc->ClearVector(W2);
    staFunc->CopyMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(VB1);
    VB1[nphase]=1.0;
    staFunc->Solve(BT,VB1);
    staFunc->CopyVector(VB1,V1);
    staFunc->ClearVector(VB2);
    VB2[nphase+1]=1.0;
    staFunc->Solve(BT,VB2);
    staFunc->CopyVector(VB2,V2);

/* solve transposed systems to obtain new values for W1 and W2*/
    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(PP1);
    PP1[nphase]=1.0;
    staFunc->Solve(BT,PP1);
    staFunc->CopyVector(PP1,W1);
    staFunc->ClearVector(PP2);
    PP2[nphase+1]=1.0;
    staFunc->Solve(BT,PP2);
    staFunc->CopyVector(PP2,W2);
                
    for (i=0;i<nphase;i++) {
      BP[nphase][i]=V1[i];
      BP[i][nphase]=W1[i];
    }
    for (i=0;i<nphase+2;i++) {
      BP[nphase+1][i]=V2[i];
      BP[i][nphase+1]=W2[i];
    }  
  } else {
     BP[0][nphase] = 1.0;
     BP[1][0] = 1.0;
     BP[nphase+1][nphase+1] = 1.0;
     staFunc->ClearVector(V1);
     staFunc->ClearVector(V2);
     staFunc->ClearVector(W1);
     staFunc->ClearVector(W2);
     V1[0]=1.0;     
     W1[0]=1.0;
     W2[nphase+1]=1.0;
    } 

    staFunc->CopyMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(VB1);
    VB1[nphase]=1.0;
    staFunc->Solve(BT,VB1);  
    staFunc->ClearVector(VB2);
    VB2[nphase+1]=1.0;
    staFunc->Solve(BT,VB2);

    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(PP1);
    PP1[nphase]=1.0;
    staFunc->Solve(BT,PP1);
    staFunc->ClearVector(PP2);
    PP2[nphase+1]=1.0;
    staFunc->Solve(BT,PP2);
/* compute derivatives of g11,g12,g21,g22 */
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);   
    if (stagenData->Der2) {
      H=SymHess(x0);
      staFunc->ClearVector(HV1);
      staFunc->ClearVector(HV2);
      staFunc->ClearVector(HW1);
      staFunc->ClearVector(HW2);
      staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
      staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
      staFunc->MultiplyMatrixVector(A,HV1,HV1);
      staFunc->MultiplyMatrixVector(A,HV2,HV2);
      staFunc->TransposeMatrix(A,B);
      staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
      staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
      staFunc->MultiplyMatrixVector(B,HW1,HW1);
      staFunc->MultiplyMatrixVector(B,HW2,HW2);  
      S=staFunc->CreateMatrix(2,2); 
      for (l=0; l<nphase+2+nappend; l++) {
        for (i=0,S[0][0]=0.0,S[0][1]=0.0,S[1][0]=0.0,S[1][1]=0.0;i<nphase;i++) {
          for (j=0; j<nphase; j++) {
    	    S[0][0]-=HW1[i]*VB1[j]*HessianElem(i,j,l)+PP1[i]*HV1[j]*HessianElem(i,j,l);
            S[0][0]-=2.0*x0[nphase+2+nappend]*PP1[i]*VB1[j]*HessianElem(i,j,l);
    	    S[0][1]-=HW1[i]*VB2[j]*HessianElem(i,j,l)+PP1[i]*HV2[j]*HessianElem(i,j,l);
            S[0][1]-=2.0*x0[nphase+2+nappend]*PP1[i]*VB2[j]*HessianElem(i,j,l);
    	    S[1][0]-=HW2[i]*VB1[j]*HessianElem(i,j,l)+PP2[i]*HV1[j]*HessianElem(i,j,l);
            S[1][0]-=2.0*x0[nphase+2+nappend]*PP2[i]*VB1[j]*HessianElem(i,j,l);
    	    S[1][1]-=HW2[i]*VB2[j]*HessianElem(i,j,l)+PP2[i]*HV2[j]*HessianElem(i,j,l);            
            S[1][1]-=2.0*x0[nphase+2+nappend]*PP2[i]*VB2[j]*HessianElem(i,j,l);
          }
        }
        G0[l]=S[0][0];
        G1[l]=S[0][1];
        G2[l]=S[1][0];
        G3[l]=S[1][1];  
      }
      S=staFunc->DestroyMatrix(S);
      G0[nphase+2+nappend]=2.0*staFunc->ScalProd(PP1,HV1);
      G1[nphase+2+nappend]=2.0*staFunc->ScalProd(PP1,HV2);
      G2[nphase+2+nappend]=2.0*staFunc->ScalProd(PP2,HV1);
      G3[nphase+2+nappend]=2.0*staFunc->ScalProd(PP2,HV2);
    } else {
        staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase); 
        staFunc->ClearVector(HV1);
        staFunc->ClearVector(HV2);
        staFunc->ClearVector(HW1);
        staFunc->ClearVector(HW2);
        staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
        staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
        staFunc->MultiplyMatrixVector(A,HV1,HV1);
        staFunc->MultiplyMatrixVector(A,HV2,HV2);
        staFunc->TransposeMatrix(A,B);
        staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
        staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
        staFunc->MultiplyMatrixVector(B,HW1,HW1);
        staFunc->MultiplyMatrixVector(B,HW2,HW2);  
        Jac1=staFunc->CreateMatrix(nphase,nphase);
        Jac2=staFunc->CreateMatrix(nphase,nphase);
        for (l=0; l<nphase+2+nappend; l++) {
          staFunc->ClearVector(xp);
          staFunc->CopyVectorElements(x0,0,xp,0,nphase+2);
          xp[l]+=_dx;
          staFunc->ClearVector(xp1);
          staFunc->CopyVectorElements(x0,0,xp1,0,nphase+2);
          xp1[l]-=_dx;
          JacRHS(xp,C1);
          staFunc->CopyMatrixBlock(C1,0,0,Jac1,0,0,nphase,nphase);
          JacRHS(xp1,C1);
          staFunc->CopyMatrixBlock(C1,0,0,Jac2,0,0,nphase,nphase);
       /* g11 */
          staFunc->CopyVectorElements(VB1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2); 
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G0[l]=-staFunc->ScalProd(HW1,K3);
          G0[l]-=2.0*x0[nphase+2+nappend]*staFunc->ScalProd(PP1,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G0[l]-=staFunc->ScalProd(K3,HV1);                                   
       /* g12 */
          staFunc->CopyVectorElements(VB2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1);  
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G1[l]=-staFunc->ScalProd(HW1,K3);
          G1[l]-=2.0*x0[nphase+2+nappend]*staFunc->ScalProd(PP1,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G1[l]-=staFunc->ScalProd(K3,HV2);                                   
       /* g21 */
          staFunc->CopyVectorElements(VB1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G2[l]=-staFunc->ScalProd(HW2,K3);
          G2[l]-=2.0*x0[nphase+2+nappend]*staFunc->ScalProd(PP2,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);   
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);   
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G2[l]-=staFunc->ScalProd(K3,HV1);                                   
       /* g22 */
          staFunc->CopyVectorElements(VB2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G3[l]=-staFunc->ScalProd(HW2,K3);
          G3[l]-=2.0*x0[nphase+2+nappend]*staFunc->ScalProd(PP2,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);   
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G3[l]-=staFunc->ScalProd(K3,HV2);                                             
        }
        Jac1=staFunc->DestroyMatrix(Jac1);
        Jac2=staFunc->DestroyMatrix(Jac2);
        G0[nphase+2+nappend]=2.0*staFunc->ScalProd(PP1,HV1);
        G1[nphase+2+nappend]=2.0*staFunc->ScalProd(PP1,HV2);
        G2[nphase+2+nappend]=2.0*staFunc->ScalProd(PP2,HV1);
        G3[nphase+2+nappend]=2.0*staFunc->ScalProd(PP2,HV2);        
      }
    staFunc->ClearMatrix(BT);
    staFunc->ClearMatrix(JM1);
    staFunc->CopyMatrixBlock(C1,0,0,BT,0,0,nphase,nphase+2);
    for (i=0; i<nphase; i++) BT[i][i]-=1.0;
/* placement of ones determined by partial pivoting */
    PIVOT=staFunc->CreateVector(nphase+2);
    staFunc->ClearVector(PIVOT);
    staFunc->TransposeMatrix(BT,JM1);
    PIVOT[nphase+1]=1.0;
    for (k=0;k<=nphase;k++) {
     kp1=k+1;
     p = k;
     for (i=kp1;i<=nphase+1;i++) {
      if (fabs(JM1[i][k]) > fabs(JM1[p][k]))
        p = i;
     }
     PIVOT[k]=p;
     if (p!=k) PIVOT[nphase+1]=-PIVOT[nphase+1];
     t=JM1[p][k];
     JM1[p][k]=JM1[k][k];
     JM1[k][k]=t;
     if (t==0.0) continue;
     for (i=kp1;i<=nphase+1;i++) JM1[i][k]=-JM1[i][k]/t;
     for (j=kp1;j<=nphase+1;j++) {
       t=JM1[p][j];
       JM1[p][j]=JM1[k][j];
       JM1[k][j]=t;
       if (t!=0.0) for (i=kp1;i<=nphase+1;i++) JM1[i][j]+=JM1[i][k]*t;
     }
    }
    staFunc->ClearVector(TV1);
    for (i=0;i<nphase+2;i++) TV1[i]=i;
    for (k=0;k<nphase+1;k++) {
      i=(Int2)PIVOT[k];
      t=TV1[i];
      TV1[i]=TV1[k];
      TV1[k]=t;
    }
    p=(Int2)TV1[nphase];
    q=(Int2)TV1[nphase+1];  
    PIVOT=staFunc->DestroyVector(PIVOT);
/* end of pivot information */

    staFunc->ClearMatrix(JM1);
    staFunc->CopyMatrixBlock(C1,0,0,JM1,0,0,nphase,nphase+2);
    for (i=0; i<nphase; i++) JM1[i][i]-=1.0;
    JM1[nphase][p]=1.0;
    JM1[nphase+1][q]=1.0;
    staFunc->ClearVector(TV1);
    TV1[nphase]=1.0;
    staFunc->Decomp(JM1);
    staFunc->Solve(JM1,TV1);
    staFunc->ClearMatrix(TVM);
    for (i=0;i<nphase+2;i++) TVM[i][0] = TV1[i];
    staFunc->ClearVector(TV1);
    TV1[nphase+1]=1.0;  
    staFunc->Solve(JM1,TV1);
    for (i=0;i<nphase+2;i++) TVM[i][1] = TV1[i];    
    TVM[nphase+2][2]=1.0;

    staFunc->TransposeMatrix(TVM,TVMT);
    staFunc->MultiplyMatrixMatrix(TVMT,TVM,HM); 
    staFunc->Decomp(HM);
    staFunc->MultiplyMatrixVector(TVMT,G0,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G0);  
    staFunc->MultiplyMatrixVector(TVMT,G1,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G1);      
    staFunc->MultiplyMatrixVector(TVMT,G2,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G2);
    staFunc->MultiplyMatrixVector(TVMT,G3,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G3);       
  
/* Determine component with maximum norm => gives m1,n1 */
    G[0][0] = staFunc->NormOfVector(G0);
    G[0][1] = staFunc->NormOfVector(G1);
    G[1][0] = staFunc->NormOfVector(G2);
    G[1][1] = staFunc->NormOfVector(G3);    
    m1=0;
    n1=0;      
    if (G[0][1] > G[m1][n1]) n1=1;  
    if (G[1][0] > G[m1][n1]) { m1=1; n1=0; }
    if (G[1][1] > G[m1][n1]) { m1=1; n1=1; }
    for (i=0;i<nphase+3;i++) {
      GM1[i][0] = G0[i]; 
      GM1[i][1] = G1[i];
      GM1[i][2] = G2[i];
      GM1[i][3] = G3[i];
    }
/* G0 has biggest norm ; necessary for projection */ 
    for (i=0;i<nphase+3;i++) {
      G0[i] = GM1[i][2*m1+n1];
    }      
    staFunc->NormalizeVector(G0);
/* G1,G2,G3 are projected on G0 */ 
    for (j=0;j<4;j++) {
      if (j != (2*m1+n1)) {
       for (i=0;i<nphase+3;i++) { /* neem een vector, niet de grootste */
          G1[i] = GM1[i][j];
       }
       scal = staFunc->ScalProd(G1,G0);
       staFunc->AddScaledVector(G1,G0,-scal,G1);
       for (i=0;i<nphase+3;i++) { /* zet nieuwe vector terug */
         GM1[i][j] = G1[i];
       }
      }
    }      
/* compute biggest norm between G1,G2 and G3  => gives m2,n2 */
    for (i=0;i<nphase+3;i++) {
      G0[i] = GM1[i][0];
      G1[i] = GM1[i][1];
      G2[i] = GM1[i][2];
      G3[i] = GM1[i][3];
    }
    G[0][0] = staFunc->NormOfVector(G0);
    G[0][1] = staFunc->NormOfVector(G1);
    G[1][0] = staFunc->NormOfVector(G2);
    G[1][1] = staFunc->NormOfVector(G3);    
    G[m1][n1]=-1.0;
    m2=0;
    n2=0;
    if (G[0][1] > G[m2][n2]) n2=1;  
    if (G[1][0] > G[m2][n2]) { m2=1; n2=0; }
    if (G[1][1] > G[m2][n2]) { m2=1; n2=1; }    

/* Compute a tangent vector to the NS curve */
    ns2_DefNS(x0,v0);
    Ret=ns2_JacDefNS(x0,D);
    for (i=0; i<nphase+3+nappend; i++) {
      v1[i]=1.0;
      staFunc->AppendMatrixVector(D,v1,DE);
      Ret=staFunc->Decomp(DE);
      if (Ret == 0) {
        v2[nphase+2+nappend]=stagenData->Forward ? 1 : -1;
        staFunc->Solve(DE,v2);
        staFunc->CopyVector(v2,v0);
        break;
      };
      v1[i]=0.0;
    };
    if (Ret) {
       staUtil->ErrMessage
         ("Unable to find a tangent vector to the Neimark-Sacker curve");   
       v1=staFunc->DestroyVector(v1);
       v2=staFunc->DestroyVector(v2);
       D=staFunc->DestroyMatrix(D);
       DE=staFunc->DestroyMatrix(DE);
       Term(FALSE);
       return 1;
    }  
    staFunc->NormalizeVector(v0);
    v1=staFunc->DestroyVector(v1);
    v2=staFunc->DestroyVector(v2);
    D=staFunc->DestroyMatrix(D);
    DE=staFunc->DestroyMatrix(DE);
  }
  /* Initialize continuation */
  ((ContDataPtr)stagenptr->GenPtr)->X0=x0;
  ((ContDataPtr)stagenptr->GenPtr)->V0=v0;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveSing=staData->ActiveSing;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveUserTest=staData->ActiveUserTest;
  ((ContDataPtr)stagenptr->GenPtr)->TypeTestFuns=staData->TypeTestFunc;
  ((ContDataPtr)stagenptr->GenPtr)->ZeroTestFuns=staData->ZeroTestFunc;
  /* That's all */
  
  return 0;
}

/* Fold (flip) -> ns2 Starter */
Entry(Int2) f_ns2_Starter(StagenDataPtr stagenptr) {
  Int2 i,j,k,kp1,l,Ret,p,q;
  Vector v1,v2;
  Matrix D,DE,S;
  VAR t,scal;
  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) {       /* "Extension" */
    MemCopy(V1,stagenptr->ContDataStaPtr,sizeof(FloatHi)*(nphase+2));
    MemCopy(W1,stagenptr->ContDataStaPtr+(nphase+2),sizeof(FloatHi)*(nphase+2));
    MemCopy(V2,stagenptr->ContDataStaPtr+2*(nphase+2),sizeof(FloatHi)*(nphase+2));    
    MemCopy(W2,stagenptr->ContDataStaPtr+3*(nphase+2),sizeof(FloatHi)*(nphase+2));    
    MemFree(stagenptr->ContDataStaPtr);
    stagenptr->ContDataStaLen=0;
    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) x0[nphase+i]=staData->P[active[i]];                             
    x0[nphase+2+nappend]=cos(BifVector[2*nphase+1]); 
  } else { 
    v1=staFunc->CreateVector(nphase+3+nappend);
    v2=staFunc->CreateVector(nphase+3+nappend);
    D=staFunc->CreateMatrix(nphase+2+nappend,nphase+3+nappend);
    DE=staFunc->CreateMatrix(nphase+3+nappend,nphase+3+nappend);

    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) x0[nphase+i]=staData->P[active[i]];
    JacRHS(x0,C1);

    if (stagenData->BifDataInOk) {
/* Message(MSG_OK,"Eigenvectors exist"); 			*/
/* bifurcation data exist 					*/
/*  BifVector[0]:                      _itmap                   */
/*  BifVector[nphase+1] to [3*nphase]: NS vectors 		*/
/*  BifVector[3*nphase+1]:	       kappa                    */
      MemCopy(V,stagenData->BifDataInPtr+nphase+1,sizeof(FloatHi)*nphase);
      MemCopy(L,stagenData->BifDataInPtr+2*nphase+1,sizeof(FloatHi)*nphase);
      staFunc->NormalizeVector(L);
      x0[nphase+2+nappend]=stagenData->BifDataInPtr[3*nphase+1];      
    } else { 
/* bifurcation data to be computed */
/* Message(MSG_OK,"Eigenvectors have to be computed");  */
      Ret=0;
{int i,j,imin,jmin;
 double r,s,d,phi,beta,gamma,lambda,lambda1,lambda2,omega,omega2,theta;
    JacRHS(x0,C1);
    staFunc->CopyMatrixBlock(C1,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),ns2_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(C1,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) {
        staUtil->ErrMessage("Unable to find eigenvector at the initial point");
        Ret=1;
        goto polufinal;
      } 
      staFunc->CopyMatrixBlock(C1,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) {
        staUtil->ErrMessage("Unable to find eigenvector at the initial point");
        Ret=1;
        goto polufinal;
      } 
      staFunc->NormalizeVector(Qr);   
      staFunc->NormalizeVectorRelative(Qi,Qr);
      staFunc->NormalizeVector(Qi);
      
      staFunc->AddScaledVector(Qr,Qi,1.0,V);
      staFunc->NormalizeVector(V);
      staFunc->AddScaledVector(Qr,Qi,-1.0,L);             
      staFunc->NormalizeVector(L); 
      x0[nphase+2+nappend]=-lambda*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(C1,0,0,BB,0,0,nphase,nphase);
      staFunc->CopyMatrixBlock(C1,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,VV) != 2*nphase-2) {
        staUtil->ErrMessage("Unable to find eigenvectors at the initial point");
        Ret=1;
        goto polufinal;
      }
      staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
      staFunc->CopyVectorElements(VV,nphase,Qi,0,nphase);

      /* normalize real and imaginary parts */
      d=staFunc->ScalProd(Qr,Qr);
      s=staFunc->ScalProd(Qi,Qi);
      r=staFunc->ScalProd(Qr,Qi);
      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++) {
          VV[i]=Qr[i]*cos(phi)-Qi[i]*sin(phi);
          VV[i+nphase]=Qr[i]*sin(phi)+Qi[i]*cos(phi);
        }
      }

      staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
      d=staFunc->ScalProd(Qr,Qr);
      gamma=1/sqrt(d);
      staFunc->MultiplyVectorScalar(VV,gamma,VV);
      staFunc->CopyVectorElements(VV,0,V,0,nphase);
      staFunc->CopyVectorElements(VV,nphase,L,0,nphase);
      x0[nphase+2+nappend]=cos(theta);
      staFunc->NormalizeVector(L);
    }
polufinal:
      if (Ret) {
        v1=staFunc->DestroyVector(v1);
        v2=staFunc->DestroyVector(v2);
        D=staFunc->DestroyMatrix(D);
        DE=staFunc->DestroyMatrix(DE);
        Term(FALSE);
        return 1;
      }
    }
}
  if (JacRHS(x0,C1)) return 1;  /* for safety */
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  staFunc->MultiplyMatrixMatrix(A,A,BP); 
  staFunc->AddScaledMatrix(BP,A,-2.0*x0[nphase+2+nappend],BP);
  for (i=0;i<nphase;i++) {
   BP[i][i]+=1.0;
  }
  
/* searching 2 rows with maximal pivots */  
  if (nphase >= 2) {
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) {
      BP1[i][nphase]=i;
      BP1[nphase][i]=i;
    } 
    /* Calculation of PIVOTS - complete pivoting */
    for (k=0;k<nphase-2;k++) {
      t=0.0;
      for (i=k;i<nphase;i++) {
	for (j=k;j<nphase;j++) {
    	  if (fabs(BP1[i][j]) > t) {
            p = i;
            q = j;
            t = fabs(BP1[i][j]);
          }           
	}
      }      
      if (k!=q) { /* switch colom k with colom  q */        
        for (i=0;i<=nphase;i++) {
          t=BP1[i][k];
          BP1[i][k]=BP1[i][q];
          BP1[i][q]=t;
	}
      }
      if (k!=p) { /* switch row k with row p */       
        for (i=0;i<=nphase;i++) {
          t=BP1[k][i];         
          BP1[k][i]=BP1[p][i];
          BP1[p][i]=t;
	}
      }        
      if (fabs(BP1[k][k]) == 0.0) {
        staUtil->ErrMessage("(A^2 + k I) has rank defect at least three");
        Term(FALSE);
        return 1;
      }
      kp1=k+1;
      for (i=kp1;i<nphase;i++) { /* perform elimination */
        t=-BP1[i][k]/BP1[k][k];
        for (j=kp1;j<nphase;j++) {
          BP1[i][j]+=t*BP1[k][j];
	}
      }
    } 
    p=(Int2)BP1[nphase-2][nphase];
    q=(Int2)BP1[nphase-1][nphase];   
    BP[p][nphase]=1.0;
    BP[q][nphase+1]=1.0;
    p=(Int2)BP1[nphase][nphase-2];
    q=(Int2)BP1[nphase][nphase-1];   
    BP[nphase][p]=1.0;
    BP[nphase+1][q]=1.0;    
       
/* solve systems to obtain new values for V1 and V2 */
    staFunc->ClearVector(V1);
    staFunc->ClearVector(V2);
    staFunc->ClearVector(W1);
    staFunc->ClearVector(W2);
    staFunc->CopyMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(VB1);
    VB1[nphase]=1.0;
    staFunc->Solve(BT,VB1);
    staFunc->CopyVector(VB1,V1);
    staFunc->ClearVector(VB2);
    VB2[nphase+1]=1.0;
    staFunc->Solve(BT,VB2);
    staFunc->CopyVector(VB2,V2);

/* solve transposed systems to obtain new values for W1 and W2 */
    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(PP1);
    PP1[nphase]=1.0;
    staFunc->Solve(BT,PP1);
    staFunc->CopyVector(PP1,W1);
    staFunc->ClearVector(PP2);
    PP2[nphase+1]=1.0;
    staFunc->Solve(BT,PP2);
    staFunc->CopyVector(PP2,W2);
                
    for (i=0;i<nphase;i++) {
      BP[nphase][i]=V1[i];
      BP[i][nphase]=W1[i];
    }
    for (i=0;i<nphase+2;i++) {
      BP[nphase+1][i]=V2[i];
      BP[i][nphase+1]=W2[i];
    }  
  } else {
     BP[0][nphase] = 1.0;
     BP[1][0] = 1.0;
     BP[nphase+1][nphase+1] = 1.0;
     staFunc->ClearVector(V1);
     staFunc->ClearVector(V2);
     staFunc->ClearVector(W1);
     staFunc->ClearVector(W2);
     V1[0]=1.0;     
     W1[0]=1.0;
     W2[nphase+1]=1.0;
    } 

    staFunc->CopyMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(VB1);
    VB1[nphase]=1.0;
    staFunc->Solve(BT,VB1);  
    staFunc->ClearVector(VB2);
    VB2[nphase+1]=1.0;
    staFunc->Solve(BT,VB2);

    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(PP1);
    PP1[nphase]=1.0;
    staFunc->Solve(BT,PP1);
    staFunc->ClearVector(PP2);
    PP2[nphase+1]=1.0;
    staFunc->Solve(BT,PP2);
/* compute derivatives of g11,g12,g21,g22 */
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);   
    if (stagenData->Der2) {
      H=SymHess(x0);
      staFunc->ClearVector(HV1);
      staFunc->ClearVector(HV2);
      staFunc->ClearVector(HW1);
      staFunc->ClearVector(HW2);
      staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
      staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
      staFunc->MultiplyMatrixVector(A,HV1,HV1);
      staFunc->MultiplyMatrixVector(A,HV2,HV2);
      staFunc->TransposeMatrix(A,B);
      staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
      staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
      staFunc->MultiplyMatrixVector(B,HW1,HW1);
      staFunc->MultiplyMatrixVector(B,HW2,HW2);  
      S=staFunc->CreateMatrix(2,2); 
      for (l=0; l<nphase+2+nappend; l++) {
        for (i=0,S[0][0]=0.0,S[0][1]=0.0,S[1][0]=0.0,S[1][1]=0.0;i<nphase;i++) {
          for (j=0; j<nphase; j++) {
    	    S[0][0]-=HW1[i]*VB1[j]*HessianElem(i,j,l)+PP1[i]*HV1[j]*HessianElem(i,j,l);
    	    S[0][1]-=HW1[i]*VB2[j]*HessianElem(i,j,l)+PP1[i]*HV2[j]*HessianElem(i,j,l);
    	    S[1][0]-=HW2[i]*VB1[j]*HessianElem(i,j,l)+PP2[i]*HV1[j]*HessianElem(i,j,l);
    	    S[1][1]-=HW2[i]*VB2[j]*HessianElem(i,j,l)+PP2[i]*HV2[j]*HessianElem(i,j,l);            
          }
        }
        G0[l]=S[0][0];
        G1[l]=S[0][1];
        G2[l]=S[1][0];
        G3[l]=S[1][1];  
      }
      S=staFunc->DestroyMatrix(S);
      staFunc->ClearVector(HV1);
      staFunc->ClearVector(HV2);
      staFunc->ClearVector(HW1);
      staFunc->ClearVector(HW2);
      staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
      staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
      staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
      staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
      G0[nphase+2+nappend]=-staFunc->ScalProd(HW1,HV1);
      G1[nphase+2+nappend]=-staFunc->ScalProd(HW1,HV2);
      G2[nphase+2+nappend]=-staFunc->ScalProd(HW2,HV1);
      G3[nphase+2+nappend]=-staFunc->ScalProd(HW2,HV2);
    } else {
        staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase); 
        staFunc->ClearVector(HV1);
        staFunc->ClearVector(HV2);
        staFunc->ClearVector(HW1);
        staFunc->ClearVector(HW2);
        staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
        staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
        staFunc->MultiplyMatrixVector(A,HV1,HV1);
        staFunc->MultiplyMatrixVector(A,HV2,HV2);
        staFunc->TransposeMatrix(A,B);
        staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
        staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
        staFunc->MultiplyMatrixVector(B,HW1,HW1);
        staFunc->MultiplyMatrixVector(B,HW2,HW2);  
        Jac1=staFunc->CreateMatrix(nphase,nphase);
        Jac2=staFunc->CreateMatrix(nphase,nphase);
        for (l=0; l<nphase+2+nappend; l++) {
          staFunc->ClearVector(xp);
          staFunc->CopyVectorElements(x0,0,xp,0,nphase+2);
          xp[l]+=_dx;
          staFunc->ClearVector(xp1);
          staFunc->CopyVectorElements(x0,0,xp1,0,nphase+2);
          xp1[l]-=_dx;
          JacRHS(xp,C1);
          staFunc->CopyMatrixBlock(C1,0,0,Jac1,0,0,nphase,nphase);
          JacRHS(xp1,C1);
          staFunc->CopyMatrixBlock(C1,0,0,Jac2,0,0,nphase,nphase);
       /* g11 */
          staFunc->CopyVectorElements(VB1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2); 
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G0[l]=-staFunc->ScalProd(HW1,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G0[l]-=staFunc->ScalProd(K3,HV1);                                   
       /* g12 */
          staFunc->CopyVectorElements(VB2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1);  
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G1[l]=-staFunc->ScalProd(HW1,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G1[l]-=staFunc->ScalProd(K3,HV2);                                   
       /* g21 */
          staFunc->CopyVectorElements(VB1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G2[l]=-staFunc->ScalProd(HW2,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);   
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);   
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G2[l]-=staFunc->ScalProd(K3,HV1);                                   
       /* g22 */
          staFunc->CopyVectorElements(VB2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G3[l]=-staFunc->ScalProd(HW2,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);   
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G3[l]-=staFunc->ScalProd(K3,HV2);                                             
        }
        Jac1=staFunc->DestroyMatrix(Jac1);
        Jac2=staFunc->DestroyMatrix(Jac2);
        G0[nphase+2+nappend]=-staFunc->ScalProd(HW1,HV1);
        G1[nphase+2+nappend]=-staFunc->ScalProd(HW1,HV2);
        G2[nphase+2+nappend]=-staFunc->ScalProd(HW2,HV1);
        G3[nphase+2+nappend]=-staFunc->ScalProd(HW2,HV2);        
      }
    staFunc->ClearMatrix(BT);
    staFunc->ClearMatrix(JM1);
    staFunc->CopyMatrixBlock(C1,0,0,BT,0,0,nphase,nphase+2);
/* placement of ones determined by partial pivoting */
    PIVOT=staFunc->CreateVector(nphase+2);
    staFunc->ClearVector(PIVOT);
    staFunc->TransposeMatrix(BT,JM1);
    PIVOT[nphase+1]=1.0;
    for (k=0;k<=nphase;k++) {
     kp1=k+1;
     p = k;
     for (i=kp1;i<=nphase+1;i++) {
      if (fabs(JM1[i][k]) > fabs(JM1[p][k]))
        p = i;
     }
     PIVOT[k]=p;
     if (p!=k) PIVOT[nphase+1]=-PIVOT[nphase+1];
     t=JM1[p][k];
     JM1[p][k]=JM1[k][k];
     JM1[k][k]=t;
     if (t==0.0) continue;
     for (i=kp1;i<=nphase+1;i++) JM1[i][k]=-JM1[i][k]/t;
     for (j=kp1;j<=nphase+1;j++) {
       t=JM1[p][j];
       JM1[p][j]=JM1[k][j];
       JM1[k][j]=t;
       if (t!=0.0) for (i=kp1;i<=nphase+1;i++) JM1[i][j]+=JM1[i][k]*t;
     }
    }
    staFunc->ClearVector(TV1);
    for (i=0;i<nphase+2;i++) TV1[i]=i;
    for (k=0;k<nphase+1;k++) {
      i=(Int2)PIVOT[k];
      t=TV1[i];
      TV1[i]=TV1[k];
      TV1[k]=t;
    }
    p=(Int2)TV1[nphase];
    q=(Int2)TV1[nphase+1];  
    PIVOT=staFunc->DestroyVector(PIVOT);
/* end of pivot information */

    staFunc->ClearMatrix(JM1);
    staFunc->CopyMatrixBlock(C1,0,0,JM1,0,0,nphase,nphase+2);
    JM1[nphase][p]=1.0;
    JM1[nphase+1][q]=1.0;
    staFunc->ClearVector(TV1);
    TV1[nphase]=1.0;
    staFunc->Decomp(JM1);
    staFunc->Solve(JM1,TV1);
    staFunc->ClearMatrix(TVM);
    for (i=0;i<nphase+2;i++) TVM[i][0] = TV1[i];
    staFunc->ClearVector(TV1);
    TV1[nphase+1]=1.0;  
    staFunc->Solve(JM1,TV1);
    for (i=0;i<nphase+2;i++) TVM[i][1] = TV1[i];    
    TVM[nphase+2][2]=1.0;

    staFunc->TransposeMatrix(TVM,TVMT);
    staFunc->MultiplyMatrixMatrix(TVMT,TVM,HM); 
    staFunc->Decomp(HM);
    staFunc->MultiplyMatrixVector(TVMT,G0,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G0);  
    staFunc->MultiplyMatrixVector(TVMT,G1,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G1);      
    staFunc->MultiplyMatrixVector(TVMT,G2,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G2);
    staFunc->MultiplyMatrixVector(TVMT,G3,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G3);       
  
/* Determine component with maximum norm => gives m1,n1 */
    G[0][0] = staFunc->NormOfVector(G0);
    G[0][1] = staFunc->NormOfVector(G1);
    G[1][0] = staFunc->NormOfVector(G2);
    G[1][1] = staFunc->NormOfVector(G3);    
    m1=0;
    n1=0;      
    if (G[0][1] > G[m1][n1]) n1=1;  
    if (G[1][0] > G[m1][n1]) { m1=1; n1=0; }
    if (G[1][1] > G[m1][n1]) { m1=1; n1=1; }
    for (i=0;i<nphase+3;i++) {
      GM1[i][0] = G0[i]; 
      GM1[i][1] = G1[i];
      GM1[i][2] = G2[i];
      GM1[i][3] = G3[i];
    }
/* G0 has biggest norm ; necessary for projection */ 
    for (i=0;i<nphase+3;i++) {
      G0[i] = GM1[i][2*m1+n1];
    }      
    staFunc->NormalizeVector(G0);
/* G1,G2,G3 are projected on G0 */ 
    for (j=0;j<4;j++) {
      if (j != (2*m1+n1)) {
       for (i=0;i<nphase+3;i++) { /* neem een vector, niet de grootste */
          G1[i] = GM1[i][j];
       }
       scal = staFunc->ScalProd(G1,G0);
       staFunc->AddScaledVector(G1,G0,-scal,G1);
       for (i=0;i<nphase+3;i++) { /* zet nieuwe vector terug */
         GM1[i][j] = G1[i];
       }
      }
    }      
/* compute biggest norm between G1,G2 and G3  => gives m2,n2 */
    for (i=0;i<nphase+3;i++) {
      G0[i] = GM1[i][0];
      G1[i] = GM1[i][1];
      G2[i] = GM1[i][2];
      G3[i] = GM1[i][3];
    }
    G[0][0] = staFunc->NormOfVector(G0);
    G[0][1] = staFunc->NormOfVector(G1);
    G[1][0] = staFunc->NormOfVector(G2);
    G[1][1] = staFunc->NormOfVector(G3);    
    G[m1][n1]=-1.0;
    m2=0;
    n2=0;
    if (G[0][1] > G[m2][n2]) n2=1;  
    if (G[1][0] > G[m2][n2]) { m2=1; n2=0; }
    if (G[1][1] > G[m2][n2]) { m2=1; n2=1; }    
   
/* Compute a tangent vector to the NS curve */
    ns2_DefNS(x0,v0);
    Ret=ns2_JacDefNS(x0,D);
    for (i=0; i<nphase+3+nappend; i++) {
      v1[i]=1.0;
      staFunc->AppendMatrixVector(D,v1,DE);
      Ret=staFunc->Decomp(DE);
      if (Ret == 0) {
        v2[nphase+2+nappend]=stagenData->Forward ? 1 : -1;
        staFunc->Solve(DE,v2);
        staFunc->CopyVector(v2,v0);
        break;
      };
      v1[i]=0.0;
    };
    if (Ret) {
       staUtil->ErrMessage
         ("Unable to find a tangent vector to the Neimark-Sacker curve");   
       v1=staFunc->DestroyVector(v1);
       v2=staFunc->DestroyVector(v2);
       D=staFunc->DestroyMatrix(D);
       DE=staFunc->DestroyMatrix(DE);
       Term(FALSE);
       return 1;
    }  
    staFunc->NormalizeVector(v0);
    v1=staFunc->DestroyVector(v1);
    v2=staFunc->DestroyVector(v2);
    D=staFunc->DestroyMatrix(D);
    DE=staFunc->DestroyMatrix(DE);
  }
  /* Initialize continuation */
  ((ContDataPtr)stagenptr->GenPtr)->X0=x0;
  ((ContDataPtr)stagenptr->GenPtr)->V0=v0;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveSing=staData->ActiveSing;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveUserTest=staData->ActiveUserTest;
  ((ContDataPtr)stagenptr->GenPtr)->TypeTestFuns=staData->TypeTestFunc;
  ((ContDataPtr)stagenptr->GenPtr)->ZeroTestFuns=staData->ZeroTestFunc;
  /* That's all */
  return 0;
}

/* R1 -> ns2 Starter */
Entry(Int2) r1_ns2_Starter(StagenDataPtr stagenptr) {
  Int2 i,j,k,kp1,l,Ret,p,q;
  Vector v1,v2;
  Matrix D,DE,S;
  VAR t,scal;
  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) {       /* extension */
    MemCopy(V1,stagenptr->ContDataStaPtr,sizeof(FloatHi)*(nphase+2));
    MemCopy(W1,stagenptr->ContDataStaPtr+(nphase+2),sizeof(FloatHi)*(nphase+2));
    MemCopy(V2,stagenptr->ContDataStaPtr+2*(nphase+2),sizeof(FloatHi)*(nphase+2));    
    MemCopy(W2,stagenptr->ContDataStaPtr+3*(nphase+2),sizeof(FloatHi)*(nphase+2));    
    MemFree(stagenptr->ContDataStaPtr);
    stagenptr->ContDataStaLen=0;
    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) x0[nphase+i]=staData->P[active[i]];                             
    x0[nphase+2+nappend]=BifVector[2*nphase+1];
  } else { 
    v1=staFunc->CreateVector(nphase+3+nappend);
    v2=staFunc->CreateVector(nphase+3+nappend);
    D=staFunc->CreateMatrix(nphase+2+nappend,nphase+3+nappend);
    DE=staFunc->CreateMatrix(nphase+3+nappend,nphase+3+nappend);

    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) x0[nphase+i]=staData->P[active[i]];
    JacRHS(x0,C1);

    if (stagenData->BifDataInOk) {
/* Message(MSG_OK,"Eigenvectors exist"); */
/* bifurcation data exist */
/*  BifVector[0]:                      _itmap */
/*  BifVector[1] to [nphase]:          fold vector V (A*V=V, <V,V>=1) */
/*  BifVector[nphase+1]:               kappa=1.0 */
      MemCopy(V,stagenData->BifDataInPtr+1,sizeof(FloatHi)*nphase);
      x0[nphase+2+nappend]=1.0;      
    } else { 
/* bifurcation data to be computed */
/* Message(MSG_OK,"Eigenvectors have to be computed");  */
      Ret=0;
      JacRHS(x0,C1);
      staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
      for (i=0; i<nphase; i++) AA[i][i]-=1.0;
      staFunc->CopyMatrixBlock(AA,0,0,BB,0,0,nphase,nphase);
      staFunc->CopyMatrixBlock(AA,0,0,BB,nphase,nphase,nphase,nphase);
      for (i=0; i<nphase; i++) BB[i+nphase][i]=-1;
      i=staFunc->SngSolve(BB,_eps,VV);
      if (i != 2*nphase-2) {
        staUtil->ErrMessage("Unable to find (generalized)eigenvectors at the initial point");
        Ret=1;
        goto final;
      }
      staFunc->CopyVectorElements(VV,0,V,0,nphase);       /* eigenvector */
      staFunc->NormalizeVector(V);
      staFunc->CopyVectorElements(VV,nphase,W,0,nphase);  /* generalized eigenvector */
      staFunc->AddScaledVector(W,V,-staFunc->ScalProd(V,W),W);
      staFunc->NormalizeVector(W);

      staFunc->CopyVectorElements(V,0,L,0,nphase);
      staFunc->CopyVectorElements(W,0,V,0,nphase);
      x0[nphase+2+nappend]=1.0;
final:
      if (Ret) {
        v1=staFunc->DestroyVector(v1);
        v2=staFunc->DestroyVector(v2);
        D=staFunc->DestroyMatrix(D);
        DE=staFunc->DestroyMatrix(DE);
        Term(FALSE);
        return 1;
      }
}

  if (JacRHS(x0,C1)) return 1;  /* for safety */
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  staFunc->MultiplyMatrixMatrix(A,A,BP); 
  staFunc->AddScaledMatrix(BP,A,-2.0*x0[nphase+2+nappend],BP);
  for (i=0;i<nphase;i++) {
   BP[i][i]+=1.0;
  }
  
/* searching 2 rows with maximal pivots */  
  if (nphase >= 2) {
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) {
      BP1[i][nphase]=i;
      BP1[nphase][i]=i;
    } 
    /* Calculation of PIVOTS - complete pivoting */
    for (k=0;k<nphase-2;k++) {
      t=0.0;
      for (i=k;i<nphase;i++) {
	for (j=k;j<nphase;j++) {
    	  if (fabs(BP1[i][j]) > t) {
            p = i;
            q = j;
            t = fabs(BP1[i][j]);
          }           
	}
      }      
      if (k!=q) { /* switch colom k with colom  q */        
        for (i=0;i<=nphase;i++) {
          t=BP1[i][k];
          BP1[i][k]=BP1[i][q];
          BP1[i][q]=t;
	}
      }
      if (k!=p) { /* switch row k with row p */       
        for (i=0;i<=nphase;i++) {
          t=BP1[k][i];         
          BP1[k][i]=BP1[p][i];
          BP1[p][i]=t;
	}
      }        
      if (fabs(BP1[k][k]) == 0.0) {
        staUtil->ErrMessage("(A^2 + k I) has rank defect at least three");
        Term(FALSE);
        return 1;
      }
      kp1=k+1;
      for (i=kp1;i<nphase;i++) { /* perform elimination */
        t=-BP1[i][k]/BP1[k][k];
        for (j=kp1;j<nphase;j++) {
          BP1[i][j]+=t*BP1[k][j];
	}
      }
    } 
    p=(Int2)BP1[nphase-2][nphase];
    q=(Int2)BP1[nphase-1][nphase];   
    BP[p][nphase]=1.0;
    BP[q][nphase+1]=1.0;
    p=(Int2)BP1[nphase][nphase-2];
    q=(Int2)BP1[nphase][nphase-1];   
    BP[nphase][p]=1.0;
    BP[nphase+1][q]=1.0;    
       
/* solve systems to obtain new values for V1 and V2 */
    staFunc->ClearVector(V1);
    staFunc->ClearVector(V2);
    staFunc->ClearVector(W1);
    staFunc->ClearVector(W2);
    staFunc->CopyMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(VB1);
    VB1[nphase]=1.0;
    staFunc->Solve(BT,VB1);
    staFunc->CopyVector(VB1,V1);
    staFunc->ClearVector(VB2);
    VB2[nphase+1]=1.0;
    staFunc->Solve(BT,VB2);
    staFunc->CopyVector(VB2,V2);

/* solve transposed systems to obtain new values for W1 and W2 */
    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(PP1);
    PP1[nphase]=1.0;
    staFunc->Solve(BT,PP1);
    staFunc->CopyVector(PP1,W1);
    staFunc->ClearVector(PP2);
    PP2[nphase+1]=1.0;
    staFunc->Solve(BT,PP2);
    staFunc->CopyVector(PP2,W2);
                
    for (i=0;i<nphase;i++) {
      BP[nphase][i]=V1[i];
      BP[i][nphase]=W1[i];
    }
    for (i=0;i<nphase+2;i++) {
      BP[nphase+1][i]=V2[i];
      BP[i][nphase+1]=W2[i];
    }  
  } else {
     BP[0][nphase] = 1.0;
     BP[1][0] = 1.0;
     BP[nphase+1][nphase+1] = 1.0;
     staFunc->ClearVector(V1);
     staFunc->ClearVector(V2);
     staFunc->ClearVector(W1);
     staFunc->ClearVector(W2);
     V1[0]=1.0;     
     W1[0]=1.0;
     W2[nphase+1]=1.0;
    } 

    staFunc->CopyMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(VB1);
    VB1[nphase]=1.0;
    staFunc->Solve(BT,VB1);  
    staFunc->ClearVector(VB2);
    VB2[nphase+1]=1.0;
    staFunc->Solve(BT,VB2);

    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(PP1);
    PP1[nphase]=1.0;
    staFunc->Solve(BT,PP1);
    staFunc->ClearVector(PP2);
    PP2[nphase+1]=1.0;
    staFunc->Solve(BT,PP2);
/* compute derivatives of g11,g12,g21,g22 */
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);   
    if (stagenData->Der2) {
      H=SymHess(x0);
      staFunc->ClearVector(HV1);
      staFunc->ClearVector(HV2);
      staFunc->ClearVector(HW1);
      staFunc->ClearVector(HW2);
      staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
      staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
      staFunc->MultiplyMatrixVector(A,HV1,HV1);
      staFunc->MultiplyMatrixVector(A,HV2,HV2);
      staFunc->TransposeMatrix(A,B);
      staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
      staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
      staFunc->MultiplyMatrixVector(B,HW1,HW1);
      staFunc->MultiplyMatrixVector(B,HW2,HW2);  
      S=staFunc->CreateMatrix(2,2); 
      for (l=0; l<nphase+2+nappend; l++) {
        for (i=0,S[0][0]=0.0,S[0][1]=0.0,S[1][0]=0.0,S[1][1]=0.0;i<nphase;i++) {
          for (j=0; j<nphase; j++) {
    	    S[0][0]-=HW1[i]*VB1[j]*HessianElem(i,j,l)+PP1[i]*HV1[j]*HessianElem(i,j,l);
    	    S[0][1]-=HW1[i]*VB2[j]*HessianElem(i,j,l)+PP1[i]*HV2[j]*HessianElem(i,j,l);
    	    S[1][0]-=HW2[i]*VB1[j]*HessianElem(i,j,l)+PP2[i]*HV1[j]*HessianElem(i,j,l);
    	    S[1][1]-=HW2[i]*VB2[j]*HessianElem(i,j,l)+PP2[i]*HV2[j]*HessianElem(i,j,l);            
          }
        }
        G0[l]=S[0][0];
        G1[l]=S[0][1];
        G2[l]=S[1][0];
        G3[l]=S[1][1];  
      }
      S=staFunc->DestroyMatrix(S);
      staFunc->ClearVector(HV1);
      staFunc->ClearVector(HV2);
      staFunc->ClearVector(HW1);
      staFunc->ClearVector(HW2);
      staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
      staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
      staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
      staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
      G0[nphase+2+nappend]=-staFunc->ScalProd(HW1,HV1);
      G1[nphase+2+nappend]=-staFunc->ScalProd(HW1,HV2);
      G2[nphase+2+nappend]=-staFunc->ScalProd(HW2,HV1);
      G3[nphase+2+nappend]=-staFunc->ScalProd(HW2,HV2);
    } else {
        staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase); 
        staFunc->ClearVector(HV1);
        staFunc->ClearVector(HV2);
        staFunc->ClearVector(HW1);
        staFunc->ClearVector(HW2);
        staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
        staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
        staFunc->MultiplyMatrixVector(A,HV1,HV1);
        staFunc->MultiplyMatrixVector(A,HV2,HV2);
        staFunc->TransposeMatrix(A,B);
        staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
        staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
        staFunc->MultiplyMatrixVector(B,HW1,HW1);
        staFunc->MultiplyMatrixVector(B,HW2,HW2);  
        Jac1=staFunc->CreateMatrix(nphase,nphase);
        Jac2=staFunc->CreateMatrix(nphase,nphase);
        for (l=0; l<nphase+2+nappend; l++) {
          staFunc->ClearVector(xp);
          staFunc->CopyVectorElements(x0,0,xp,0,nphase+2);
          xp[l]+=_dx;
          staFunc->ClearVector(xp1);
          staFunc->CopyVectorElements(x0,0,xp1,0,nphase+2);
          xp1[l]-=_dx;
          JacRHS(xp,C1);
          staFunc->CopyMatrixBlock(C1,0,0,Jac1,0,0,nphase,nphase);
          JacRHS(xp1,C1);
          staFunc->CopyMatrixBlock(C1,0,0,Jac2,0,0,nphase,nphase);
       /* g11 */
          staFunc->CopyVectorElements(VB1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2); 
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G0[l]=-staFunc->ScalProd(HW1,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G0[l]-=staFunc->ScalProd(K3,HV1);                                   
       /* g12 */
          staFunc->CopyVectorElements(VB2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1);  
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G1[l]=-staFunc->ScalProd(HW1,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G1[l]-=staFunc->ScalProd(K3,HV2);                                   
       /* g21 */
          staFunc->CopyVectorElements(VB1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G2[l]=-staFunc->ScalProd(HW2,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);   
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);   
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G2[l]-=staFunc->ScalProd(K3,HV1);                                   
       /* g22 */
          staFunc->CopyVectorElements(VB2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G3[l]=-staFunc->ScalProd(HW2,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);   
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G3[l]-=staFunc->ScalProd(K3,HV2);                                             
        }
        Jac1=staFunc->DestroyMatrix(Jac1);
        Jac2=staFunc->DestroyMatrix(Jac2);
        G0[nphase+2+nappend]=-staFunc->ScalProd(HW1,HV1);
        G1[nphase+2+nappend]=-staFunc->ScalProd(HW1,HV2);
        G2[nphase+2+nappend]=-staFunc->ScalProd(HW2,HV1);
        G3[nphase+2+nappend]=-staFunc->ScalProd(HW2,HV2);        
      }
    staFunc->ClearMatrix(BT);
    staFunc->ClearMatrix(JM1);
    staFunc->CopyMatrixBlock(C1,0,0,BT,0,0,nphase,nphase+2);
/* placement of ones determined by partial pivoting */
    PIVOT=staFunc->CreateVector(nphase+2);
    staFunc->ClearVector(PIVOT);
    staFunc->TransposeMatrix(BT,JM1);
    PIVOT[nphase+1]=1.0;
    for (k=0;k<=nphase;k++) {
     kp1=k+1;
     p = k;
     for (i=kp1;i<=nphase+1;i++) {
      if (fabs(JM1[i][k]) > fabs(JM1[p][k]))
        p = i;
     }
     PIVOT[k]=p;
     if (p!=k) PIVOT[nphase+1]=-PIVOT[nphase+1];
     t=JM1[p][k];
     JM1[p][k]=JM1[k][k];
     JM1[k][k]=t;
     if (t==0.0) continue;
     for (i=kp1;i<=nphase+1;i++) JM1[i][k]=-JM1[i][k]/t;
     for (j=kp1;j<=nphase+1;j++) {
       t=JM1[p][j];
       JM1[p][j]=JM1[k][j];
       JM1[k][j]=t;
       if (t!=0.0) for (i=kp1;i<=nphase+1;i++) JM1[i][j]+=JM1[i][k]*t;
     }
    }
    staFunc->ClearVector(TV1);
    for (i=0;i<nphase+2;i++) TV1[i]=i;
    for (k=0;k<nphase+1;k++) {
      i=(Int2)PIVOT[k];
      t=TV1[i];
      TV1[i]=TV1[k];
      TV1[k]=t;
    }
    p=(Int2)TV1[nphase];
    q=(Int2)TV1[nphase+1];  
    PIVOT=staFunc->DestroyVector(PIVOT);
/* end of pivot information */

    staFunc->ClearMatrix(JM1);
    staFunc->CopyMatrixBlock(C1,0,0,JM1,0,0,nphase,nphase+2);
    JM1[nphase][p]=1.0;
    JM1[nphase+1][q]=1.0;
    staFunc->ClearVector(TV1);
    TV1[nphase]=1.0;
    staFunc->Decomp(JM1);
    staFunc->Solve(JM1,TV1);
    staFunc->ClearMatrix(TVM);
    for (i=0;i<nphase+2;i++) TVM[i][0] = TV1[i];
    staFunc->ClearVector(TV1);
    TV1[nphase+1]=1.0;  
    staFunc->Solve(JM1,TV1);
    for (i=0;i<nphase+2;i++) TVM[i][1] = TV1[i];    
    TVM[nphase+2][2]=1.0;

    staFunc->TransposeMatrix(TVM,TVMT);
    staFunc->MultiplyMatrixMatrix(TVMT,TVM,HM); 
    staFunc->Decomp(HM);
    staFunc->MultiplyMatrixVector(TVMT,G0,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G0);  
    staFunc->MultiplyMatrixVector(TVMT,G1,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G1);      
    staFunc->MultiplyMatrixVector(TVMT,G2,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G2);
    staFunc->MultiplyMatrixVector(TVMT,G3,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G3);       
  
/* Determine component with maximum norm => gives m1,n1 */
    G[0][0] = staFunc->NormOfVector(G0);
    G[0][1] = staFunc->NormOfVector(G1);
    G[1][0] = staFunc->NormOfVector(G2);
    G[1][1] = staFunc->NormOfVector(G3);    
    m1=0;
    n1=0;      
    if (G[0][1] > G[m1][n1]) n1=1;  
    if (G[1][0] > G[m1][n1]) { m1=1; n1=0; }
    if (G[1][1] > G[m1][n1]) { m1=1; n1=1; }
    for (i=0;i<nphase+3;i++) {
      GM1[i][0] = G0[i]; 
      GM1[i][1] = G1[i];
      GM1[i][2] = G2[i];
      GM1[i][3] = G3[i];
    }
/* G0 has biggest norm ; necessary for projection */ 
    for (i=0;i<nphase+3;i++) {
      G0[i] = GM1[i][2*m1+n1];
    }      
    staFunc->NormalizeVector(G0);
/* G1,G2,G3 are projected on G0 */ 
    for (j=0;j<4;j++) {
      if (j != (2*m1+n1)) {
       for (i=0;i<nphase+3;i++) { /* neem een vector, niet de grootste */
          G1[i] = GM1[i][j];
       }
       scal = staFunc->ScalProd(G1,G0);
       staFunc->AddScaledVector(G1,G0,-scal,G1);
       for (i=0;i<nphase+3;i++) { /* zet nieuwe vector terug */
         GM1[i][j] = G1[i];
       }
      }
    }      
/* compute biggest norm between G1,G2 and G3  => gives m2,n2 */
    for (i=0;i<nphase+3;i++) {
      G0[i] = GM1[i][0];
      G1[i] = GM1[i][1];
      G2[i] = GM1[i][2];
      G3[i] = GM1[i][3];
    }
    G[0][0] = staFunc->NormOfVector(G0);
    G[0][1] = staFunc->NormOfVector(G1);
    G[1][0] = staFunc->NormOfVector(G2);
    G[1][1] = staFunc->NormOfVector(G3);    
    G[m1][n1]=-1.0;
    m2=0;
    n2=0;
    if (G[0][1] > G[m2][n2]) n2=1;  
    if (G[1][0] > G[m2][n2]) { m2=1; n2=0; }
    if (G[1][1] > G[m2][n2]) { m2=1; n2=1; }    
   
/* Compute a tangent vector to the NS curve */
    ns2_DefNS(x0,v0);
    Ret=ns2_JacDefNS(x0,D);
    for (i=0; i<nphase+3+nappend; i++) {
      v1[i]=1.0;
      staFunc->AppendMatrixVector(D,v1,DE);
      Ret=staFunc->Decomp(DE);
      if (Ret == 0) {
        v2[nphase+2+nappend]=stagenData->Forward ? 1 : -1;
        staFunc->Solve(DE,v2);
        staFunc->CopyVector(v2,v0);
        break;
      };
      v1[i]=0.0;
    };
    if (Ret) {
       staUtil->ErrMessage
         ("Unable to find a tangent vector to the Neimark-Sacker curve");   
       v1=staFunc->DestroyVector(v1);
       v2=staFunc->DestroyVector(v2);
       D=staFunc->DestroyMatrix(D);
       DE=staFunc->DestroyMatrix(DE);
       Term(FALSE);
       return 1;
    }  
    staFunc->NormalizeVector(v0);
    v1=staFunc->DestroyVector(v1);
    v2=staFunc->DestroyVector(v2);
    D=staFunc->DestroyMatrix(D);
    DE=staFunc->DestroyMatrix(DE);
  }
  /* Initialize continuation */
  ((ContDataPtr)stagenptr->GenPtr)->X0=x0;
  ((ContDataPtr)stagenptr->GenPtr)->V0=v0;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveSing=staData->ActiveSing;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveUserTest=staData->ActiveUserTest;
  ((ContDataPtr)stagenptr->GenPtr)->TypeTestFuns=staData->TypeTestFunc;
  ((ContDataPtr)stagenptr->GenPtr)->ZeroTestFuns=staData->ZeroTestFunc;
  /* That's all */
  return 0;
}

/* R2 -> ns2 Starter */
Entry(Int2) r2_ns2_Starter(StagenDataPtr stagenptr) {

  Int2 i,j,k,kp1,l,Ret,p,q;
  Vector v1,v2;
  Matrix D,DE,S;
  VAR t,scal;
  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) {       /* extension */
    MemCopy(V1,stagenptr->ContDataStaPtr,sizeof(FloatHi)*(nphase+2));
    MemCopy(W1,stagenptr->ContDataStaPtr+(nphase+2),sizeof(FloatHi)*(nphase+2));
    MemCopy(V2,stagenptr->ContDataStaPtr+2*(nphase+2),sizeof(FloatHi)*(nphase+2));    
    MemCopy(W2,stagenptr->ContDataStaPtr+3*(nphase+2),sizeof(FloatHi)*(nphase+2));    
    MemFree(stagenptr->ContDataStaPtr);
    stagenptr->ContDataStaLen=0;
    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) x0[nphase+i]=staData->P[active[i]];                             
    x0[nphase+2+nappend]=BifVector[2*nphase+1];
  } else { 
    v1=staFunc->CreateVector(nphase+3+nappend);
    v2=staFunc->CreateVector(nphase+3+nappend);
    D=staFunc->CreateMatrix(nphase+2+nappend,nphase+3+nappend);
    DE=staFunc->CreateMatrix(nphase+3+nappend,nphase+3+nappend);

    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) x0[nphase+i]=staData->P[active[i]];
    JacRHS(x0,C1);

    if (stagenData->BifDataInOk) {
/* Message(MSG_OK,"Eigenvectors exist"); */
/* bifurcation data exist */
/*  BifVector[0]:                      _itmap */
/*  BifVector[1] to [nphase]:          fold vector V (A*V=-V, <V,V>=1) */
/*  BifVector[nphase+1]:               kappa */
      MemCopy(V,stagenData->BifDataInPtr+1,sizeof(FloatHi)*nphase);
      x0[nphase+2+nappend]=-1.0;      
    } else { 
/* bifurcation data to be computed */
/* Message(MSG_OK,"Eigenvectors have to be computed");  */
      Ret=0;
      JacRHS(x0,C1);
      staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
      for (i=0; i<nphase; i++) AA[i][i]+=1.0;
      staFunc->CopyMatrixBlock(AA,0,0,BB,0,0,nphase,nphase);
      staFunc->CopyMatrixBlock(AA,0,0,BB,nphase,nphase,nphase,nphase);
      for (i=0; i<nphase; i++) BB[i+nphase][i]=-1;
      i=staFunc->SngSolve(BB,_eps,VV);
      if (i != 2*nphase-2) {
        staUtil->ErrMessage("Unable to find (generalized)eigenvectors at the initial point");
        Ret=1;
        goto final;
      }
      staFunc->CopyVectorElements(VV,0,V,0,nphase);       /* eigenvector */
      staFunc->NormalizeVector(V);
      staFunc->CopyVectorElements(VV,nphase,W,0,nphase);  /* generalized eigenvector */
      staFunc->AddScaledVector(W,V,-staFunc->ScalProd(V,W),W);
      staFunc->NormalizeVector(W);

      staFunc->CopyVectorElements(V,0,L,0,nphase);
      staFunc->CopyVectorElements(W,0,V,0,nphase);
      x0[nphase+2+nappend]=-1.0;
final:
      if (Ret) {
        v1=staFunc->DestroyVector(v1);
        v2=staFunc->DestroyVector(v2);
        D=staFunc->DestroyMatrix(D);
        DE=staFunc->DestroyMatrix(DE);
        Term(FALSE);
        return 1;
      }
}

  if (JacRHS(x0,C1)) return 1;  /* for safety */
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  staFunc->MultiplyMatrixMatrix(A,A,BP); 
  staFunc->AddScaledMatrix(BP,A,-2.0*x0[nphase+2+nappend],BP);
  for (i=0;i<nphase;i++) {
   BP[i][i]+=1.0;
  }
  
/* searching 2 rows with maximal pivots */  
  if (nphase >= 2) {
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) {
      BP1[i][nphase]=i;
      BP1[nphase][i]=i;
    } 
    /* Calculation of PIVOTS - complete pivoting */
    for (k=0;k<nphase-2;k++) {
      t=0.0;
      for (i=k;i<nphase;i++) {
	for (j=k;j<nphase;j++) {
    	  if (fabs(BP1[i][j]) > t) {
            p = i;
            q = j;
            t = fabs(BP1[i][j]);
          }           
	}
      }      
      if (k!=q) { /* switch colom k with colom  q */        
        for (i=0;i<=nphase;i++) {
          t=BP1[i][k];
          BP1[i][k]=BP1[i][q];
          BP1[i][q]=t;
	}
      }
      if (k!=p) { /* switch row k with row p */       
        for (i=0;i<=nphase;i++) {
          t=BP1[k][i];         
          BP1[k][i]=BP1[p][i];
          BP1[p][i]=t;
	}
      }        
      if (fabs(BP1[k][k]) == 0.0) {
        staUtil->ErrMessage("(A^2 + k I) has rank defect at least three");
        Term(FALSE);
        return 1;
      }
      kp1=k+1;
      for (i=kp1;i<nphase;i++) { /* perform elimination */
        t=-BP1[i][k]/BP1[k][k];
        for (j=kp1;j<nphase;j++) {
          BP1[i][j]+=t*BP1[k][j];
	}
      }
    } 
    p=(Int2)BP1[nphase-2][nphase];
    q=(Int2)BP1[nphase-1][nphase];   
    BP[p][nphase]=1.0;
    BP[q][nphase+1]=1.0;
    p=(Int2)BP1[nphase][nphase-2];
    q=(Int2)BP1[nphase][nphase-1];   
    BP[nphase][p]=1.0;
    BP[nphase+1][q]=1.0;    
       
/* solve systems to obtain new values for V1 and V2 */
    staFunc->ClearVector(V1);
    staFunc->ClearVector(V2);
    staFunc->ClearVector(W1);
    staFunc->ClearVector(W2);
    staFunc->CopyMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(VB1);
    VB1[nphase]=1.0;
    staFunc->Solve(BT,VB1);
    staFunc->CopyVector(VB1,V1);
    staFunc->ClearVector(VB2);
    VB2[nphase+1]=1.0;
    staFunc->Solve(BT,VB2);
    staFunc->CopyVector(VB2,V2);

/* solve transposed systems to obtain new values for W1 and W2 */
    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(PP1);
    PP1[nphase]=1.0;
    staFunc->Solve(BT,PP1);
    staFunc->CopyVector(PP1,W1);
    staFunc->ClearVector(PP2);
    PP2[nphase+1]=1.0;
    staFunc->Solve(BT,PP2);
    staFunc->CopyVector(PP2,W2);
                
    for (i=0;i<nphase;i++) {
      BP[nphase][i]=V1[i];
      BP[i][nphase]=W1[i];
    }
    for (i=0;i<nphase+2;i++) {
      BP[nphase+1][i]=V2[i];
      BP[i][nphase+1]=W2[i];
    }  
  } else {
     BP[0][nphase] = 1.0;
     BP[1][0] = 1.0;
     BP[nphase+1][nphase+1] = 1.0;
     staFunc->ClearVector(V1);
     staFunc->ClearVector(V2);
     staFunc->ClearVector(W1);
     staFunc->ClearVector(W2);
     V1[0]=1.0;     
     W1[0]=1.0;
     W2[nphase+1]=1.0;
    } 

    staFunc->CopyMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(VB1);
    VB1[nphase]=1.0;
    staFunc->Solve(BT,VB1);  
    staFunc->ClearVector(VB2);
    VB2[nphase+1]=1.0;
    staFunc->Solve(BT,VB2);

    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(PP1);
    PP1[nphase]=1.0;
    staFunc->Solve(BT,PP1);
    staFunc->ClearVector(PP2);
    PP2[nphase+1]=1.0;
    staFunc->Solve(BT,PP2);
/* compute derivatives of g11,g12,g21,g22 */
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);   
    if (stagenData->Der2) {
      H=SymHess(x0);
      staFunc->ClearVector(HV1);
      staFunc->ClearVector(HV2);
      staFunc->ClearVector(HW1);
      staFunc->ClearVector(HW2);
      staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
      staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
      staFunc->MultiplyMatrixVector(A,HV1,HV1);
      staFunc->MultiplyMatrixVector(A,HV2,HV2);
      staFunc->TransposeMatrix(A,B);
      staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
      staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
      staFunc->MultiplyMatrixVector(B,HW1,HW1);
      staFunc->MultiplyMatrixVector(B,HW2,HW2);  
      S=staFunc->CreateMatrix(2,2); 
      for (l=0; l<nphase+2+nappend; l++) {
        for (i=0,S[0][0]=0.0,S[0][1]=0.0,S[1][0]=0.0,S[1][1]=0.0;i<nphase;i++) {
          for (j=0; j<nphase; j++) {
    	    S[0][0]-=HW1[i]*VB1[j]*HessianElem(i,j,l)+PP1[i]*HV1[j]*HessianElem(i,j,l);
    	    S[0][1]-=HW1[i]*VB2[j]*HessianElem(i,j,l)+PP1[i]*HV2[j]*HessianElem(i,j,l);
    	    S[1][0]-=HW2[i]*VB1[j]*HessianElem(i,j,l)+PP2[i]*HV1[j]*HessianElem(i,j,l);
    	    S[1][1]-=HW2[i]*VB2[j]*HessianElem(i,j,l)+PP2[i]*HV2[j]*HessianElem(i,j,l);            
          }
        }
        G0[l]=S[0][0];
        G1[l]=S[0][1];
        G2[l]=S[1][0];
        G3[l]=S[1][1];  
      }
      S=staFunc->DestroyMatrix(S);
      staFunc->ClearVector(HV1);
      staFunc->ClearVector(HV2);
      staFunc->ClearVector(HW1);
      staFunc->ClearVector(HW2);
      staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
      staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
      staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
      staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
      G0[nphase+2+nappend]=-staFunc->ScalProd(HW1,HV1);
      G1[nphase+2+nappend]=-staFunc->ScalProd(HW1,HV2);
      G2[nphase+2+nappend]=-staFunc->ScalProd(HW2,HV1);
      G3[nphase+2+nappend]=-staFunc->ScalProd(HW2,HV2);
    } else {
        staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase); 
        staFunc->ClearVector(HV1);
        staFunc->ClearVector(HV2);
        staFunc->ClearVector(HW1);
        staFunc->ClearVector(HW2);
        staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
        staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
        staFunc->MultiplyMatrixVector(A,HV1,HV1);
        staFunc->MultiplyMatrixVector(A,HV2,HV2);
        staFunc->TransposeMatrix(A,B);
        staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
        staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
        staFunc->MultiplyMatrixVector(B,HW1,HW1);
        staFunc->MultiplyMatrixVector(B,HW2,HW2);  
        Jac1=staFunc->CreateMatrix(nphase,nphase);
        Jac2=staFunc->CreateMatrix(nphase,nphase);
        for (l=0; l<nphase+2+nappend; l++) {
          staFunc->ClearVector(xp);
          staFunc->CopyVectorElements(x0,0,xp,0,nphase+2);
          xp[l]+=_dx;
          staFunc->ClearVector(xp1);
          staFunc->CopyVectorElements(x0,0,xp1,0,nphase+2);
          xp1[l]-=_dx;
          JacRHS(xp,C1);
          staFunc->CopyMatrixBlock(C1,0,0,Jac1,0,0,nphase,nphase);
          JacRHS(xp1,C1);
          staFunc->CopyMatrixBlock(C1,0,0,Jac2,0,0,nphase,nphase);
       /* g11 */
          staFunc->CopyVectorElements(VB1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2); 
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G0[l]=-staFunc->ScalProd(HW1,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G0[l]-=staFunc->ScalProd(K3,HV1);                                   
       /* g12 */
          staFunc->CopyVectorElements(VB2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1);  
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G1[l]=-staFunc->ScalProd(HW1,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G1[l]-=staFunc->ScalProd(K3,HV2);                                   
       /* g21 */
          staFunc->CopyVectorElements(VB1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G2[l]=-staFunc->ScalProd(HW2,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);   
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);   
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G2[l]-=staFunc->ScalProd(K3,HV1);                                   
       /* g22 */
          staFunc->CopyVectorElements(VB2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G3[l]=-staFunc->ScalProd(HW2,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);   
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G3[l]-=staFunc->ScalProd(K3,HV2);                                             
        }
        Jac1=staFunc->DestroyMatrix(Jac1);
        Jac2=staFunc->DestroyMatrix(Jac2);
        G0[nphase+2+nappend]=-staFunc->ScalProd(HW1,HV1);
        G1[nphase+2+nappend]=-staFunc->ScalProd(HW1,HV2);
        G2[nphase+2+nappend]=-staFunc->ScalProd(HW2,HV1);
        G3[nphase+2+nappend]=-staFunc->ScalProd(HW2,HV2);        
      }
    staFunc->ClearMatrix(BT);
    staFunc->ClearMatrix(JM1);
    staFunc->CopyMatrixBlock(C1,0,0,BT,0,0,nphase,nphase+2);
/* placement of ones determined by partial pivoting */
    PIVOT=staFunc->CreateVector(nphase+2);
    staFunc->ClearVector(PIVOT);
    staFunc->TransposeMatrix(BT,JM1);
    PIVOT[nphase+1]=1.0;
    for (k=0;k<=nphase;k++) {
     kp1=k+1;
     p = k;
     for (i=kp1;i<=nphase+1;i++) {
      if (fabs(JM1[i][k]) > fabs(JM1[p][k]))
        p = i;
     }
     PIVOT[k]=p;
     if (p!=k) PIVOT[nphase+1]=-PIVOT[nphase+1];
     t=JM1[p][k];
     JM1[p][k]=JM1[k][k];
     JM1[k][k]=t;
     if (t==0.0) continue;
     for (i=kp1;i<=nphase+1;i++) JM1[i][k]=-JM1[i][k]/t;
     for (j=kp1;j<=nphase+1;j++) {
       t=JM1[p][j];
       JM1[p][j]=JM1[k][j];
       JM1[k][j]=t;
       if (t!=0.0) for (i=kp1;i<=nphase+1;i++) JM1[i][j]+=JM1[i][k]*t;
     }
    }
    staFunc->ClearVector(TV1);
    for (i=0;i<nphase+2;i++) TV1[i]=i;
    for (k=0;k<nphase+1;k++) {
      i=(Int2)PIVOT[k];
      t=TV1[i];
      TV1[i]=TV1[k];
      TV1[k]=t;
    }
    p=(Int2)TV1[nphase];
    q=(Int2)TV1[nphase+1];  
    PIVOT=staFunc->DestroyVector(PIVOT);
/* end of pivot information */

    staFunc->ClearMatrix(JM1);
    staFunc->CopyMatrixBlock(C1,0,0,JM1,0,0,nphase,nphase+2);
    JM1[nphase][p]=1.0;
    JM1[nphase+1][q]=1.0;
    staFunc->ClearVector(TV1);
    TV1[nphase]=1.0;
    staFunc->Decomp(JM1);
    staFunc->Solve(JM1,TV1);
    staFunc->ClearMatrix(TVM);
    for (i=0;i<nphase+2;i++) TVM[i][0] = TV1[i];
    staFunc->ClearVector(TV1);
    TV1[nphase+1]=1.0;  
    staFunc->Solve(JM1,TV1);
    for (i=0;i<nphase+2;i++) TVM[i][1] = TV1[i];    
    TVM[nphase+2][2]=1.0;

    staFunc->TransposeMatrix(TVM,TVMT);
    staFunc->MultiplyMatrixMatrix(TVMT,TVM,HM); 
    staFunc->Decomp(HM);
    staFunc->MultiplyMatrixVector(TVMT,G0,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G0);  
    staFunc->MultiplyMatrixVector(TVMT,G1,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G1);      
    staFunc->MultiplyMatrixVector(TVMT,G2,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G2);
    staFunc->MultiplyMatrixVector(TVMT,G3,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G3);       
  
/* Determine component with maximum norm => gives m1,n1 */
    G[0][0] = staFunc->NormOfVector(G0);
    G[0][1] = staFunc->NormOfVector(G1);
    G[1][0] = staFunc->NormOfVector(G2);
    G[1][1] = staFunc->NormOfVector(G3);    
    m1=0;
    n1=0;      
    if (G[0][1] > G[m1][n1]) n1=1;  
    if (G[1][0] > G[m1][n1]) { m1=1; n1=0; }
    if (G[1][1] > G[m1][n1]) { m1=1; n1=1; }
    for (i=0;i<nphase+3;i++) {
      GM1[i][0] = G0[i]; 
      GM1[i][1] = G1[i];
      GM1[i][2] = G2[i];
      GM1[i][3] = G3[i];
    }
/* G0 has biggest norm ; necessary for projection */ 
    for (i=0;i<nphase+3;i++) {
      G0[i] = GM1[i][2*m1+n1];
    }      
    staFunc->NormalizeVector(G0);
/* G1,G2,G3 are projected on G0 */ 
    for (j=0;j<4;j++) {
      if (j != (2*m1+n1)) {
       for (i=0;i<nphase+3;i++) { /* neem een vector, niet de grootste */
          G1[i] = GM1[i][j];
       }
       scal = staFunc->ScalProd(G1,G0);
       staFunc->AddScaledVector(G1,G0,-scal,G1);
       for (i=0;i<nphase+3;i++) { /* zet nieuwe vector terug */
         GM1[i][j] = G1[i];
       }
      }
    }      
/* compute biggest norm between G1,G2 and G3  => gives m2,n2 */
    for (i=0;i<nphase+3;i++) {
      G0[i] = GM1[i][0];
      G1[i] = GM1[i][1];
      G2[i] = GM1[i][2];
      G3[i] = GM1[i][3];
    }
    G[0][0] = staFunc->NormOfVector(G0);
    G[0][1] = staFunc->NormOfVector(G1);
    G[1][0] = staFunc->NormOfVector(G2);
    G[1][1] = staFunc->NormOfVector(G3);    
    G[m1][n1]=-1.0;
    m2=0;
    n2=0;
    if (G[0][1] > G[m2][n2]) n2=1;  
    if (G[1][0] > G[m2][n2]) { m2=1; n2=0; }
    if (G[1][1] > G[m2][n2]) { m2=1; n2=1; }    
   
/* Compute a tangent vector to the NS curve */
    ns2_DefNS(x0,v0);
    Ret=ns2_JacDefNS(x0,D);
    for (i=0; i<nphase+3+nappend; i++) {
      v1[i]=1.0;
      staFunc->AppendMatrixVector(D,v1,DE);
      Ret=staFunc->Decomp(DE);
      if (Ret == 0) {
        v2[nphase+2+nappend]=stagenData->Forward ? 1 : -1;
        staFunc->Solve(DE,v2);
        staFunc->CopyVector(v2,v0);
        break;
      };
      v1[i]=0.0;
    };
    if (Ret) {
       staUtil->ErrMessage
         ("Unable to find a tangent vector to the Neimark-Sacker curve");   
       v1=staFunc->DestroyVector(v1);
       v2=staFunc->DestroyVector(v2);
       D=staFunc->DestroyMatrix(D);
       DE=staFunc->DestroyMatrix(DE);
       Term(FALSE);
       return 1;
    }  
    staFunc->NormalizeVector(v0);
    v1=staFunc->DestroyVector(v1);
    v2=staFunc->DestroyVector(v2);
    D=staFunc->DestroyMatrix(D);
    DE=staFunc->DestroyMatrix(DE);
  }
  /* Initialize continuation */
  ((ContDataPtr)stagenptr->GenPtr)->X0=x0;
  ((ContDataPtr)stagenptr->GenPtr)->V0=v0;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveSing=staData->ActiveSing;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveUserTest=staData->ActiveUserTest;
  ((ContDataPtr)stagenptr->GenPtr)->TypeTestFuns=staData->TypeTestFunc;
  ((ContDataPtr)stagenptr->GenPtr)->ZeroTestFuns=staData->ZeroTestFunc;
  /* That's all */
  return 0;
}

/*************************************************************/
/* Defining conditions for the eigenvector: (A*A+kappa*E)V=0 */
/* [xxx] : x[2*nphase+2+nappend] nok */
Entry(Int2) DefEig (Vector x, Vector y) {

  if (JacRHS(x,C1)) return 1;	    /* A,C1,V,W - global */
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);

  staFunc->CopyVectorElements(x, nphase+2+nappend, V, 0, nphase);
  staFunc->MultiplyMatrixVector(A,V,W);
  staFunc->MultiplyMatrixVector(A,W,V);
  staFunc->CopyVectorElements(x, nphase+2+nappend, W, 0, nphase);
  staFunc->AddScaledVector(V,W,x[2*nphase+2+nappend],V);
  staFunc->CopyVectorElements(V, 0, y, 0, nphase);
  return 0;
}

/****************************************************/
/* Jacobian of defining conditions for eigenvectors */
Local(Int2) JacDefEig (Vector x, Matrix Jac) {
int i,j,k,l;
double s;
  staFunc->ClearMatrix(Jac);
  if (stagenData->Der2) {
    /* creation of DefEig Jacobian matrix Jac using symbolic second derivatives */
    H=SymHess(x);
    if (JacRHS(x,C1)) return 1;
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) {
      for (l=0; l<nphase+2+nappend; l++) {
        for (s=0.0,j=0; j<nphase; j++) {
          for (k=0; k<nphase; k++)
            s+=(HessianElem(i,l,j)*A[j][k]+HessianElem(j,l,k)*A[i][j])*x[k+nphase+2+nappend];       
        }
        Jac[i][l]=s;
      }
    }
    for (i=0; i<nphase; i++) {
      for (l=0; l<nphase; l++) {
        for (s=0.0,j=0; j<nphase; j++) s+=A[i][j]*A[j][l];
        if (i==l) s+=x[2*nphase+2+nappend];
        Jac[i][l+nphase+2+nappend]=s;
      }
      Jac[i][2*nphase+2+nappend]=x[i+nphase+2+nappend];
    }
  } else {
    if (staFunc->Der(DefEig, nphase, x, dx, Jac)) return 1;
  }
  return(0);
}

Local(Int2) ns2_JacDefFixedPoint (Vector x, Matrix Jac);

/* nn -> ns2 Starter */
Entry(Int2) nn_ns2_Starter(StagenDataPtr stagenptr) {
  Int2 i,j,k,kp1,l,Ret,p,q;
  Vector v1,v2;
  Matrix D,DE,S;
  VAR t,scal;
  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) {       /* extension */
    MemCopy(V1,stagenptr->ContDataStaPtr,sizeof(FloatHi)*(nphase+2));
    MemCopy(W1,stagenptr->ContDataStaPtr+(nphase+2),sizeof(FloatHi)*(nphase+2));
    MemCopy(V2,stagenptr->ContDataStaPtr+2*(nphase+2),sizeof(FloatHi)*(nphase+2));    
    MemCopy(W2,stagenptr->ContDataStaPtr+3*(nphase+2),sizeof(FloatHi)*(nphase+2));    
    MemFree(stagenptr->ContDataStaPtr);
    stagenptr->ContDataStaLen=0;
    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) x0[nphase+i]=staData->P[active[i]];                             
    x0[nphase+2+nappend]=BifVector[2*nphase+1]; 
  } else { /* forward or backward */
    v1=staFunc->CreateVector(nphase+3+nappend);
    v2=staFunc->CreateVector(nphase+3+nappend);
    D=staFunc->CreateMatrix(nphase+2+nappend,nphase+3+nappend);
    DE=staFunc->CreateMatrix(nphase+3+nappend,nphase+3+nappend);

    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<2+nappend; i++) x0[nphase+i]=staData->P[active[i]];
    JacRHS(x0,C1);

    if (stagenData->BifDataInOk) { /* bifurcation data exist */
/* Primary branch: */
/*  BifVector[0]:                          _itmap */
/*  BifVector[1] to [nphase]:              1-NS continuation vector */
/*  BifVector[nphase+1] to [2*nphase]:     1-global vector L */
/*  BifVector[2*nphase+1]:                 1-kappa */
/* Secondary branch: */
/*  BifVector[2*nphase+2] to [3*nphase+1]: 2-NS continuation vector */
/*  BifVector[3*nphase+2] to [4*nphase+1]: 2-global vector L */
/*  BifVector[4*nphase+2]:                 2-kappa*/

      if (_branch == 0) { /* primary branch */
        MemCopy(VB1,stagenData->BifDataInPtr+1,sizeof(FloatHi)*nphase); 
        MemCopy(VB2,stagenData->BifDataInPtr+nphase+1,sizeof(FloatHi)*nphase); 
        x0[nphase+2+nappend]=stagenData->BifDataInPtr[2*nphase+1]; 
      } else { /* secondary branch */
        MemCopy(VB1,stagenData->BifDataInPtr+2*nphase+2,sizeof(FloatHi)*nphase); 
        MemCopy(VB2,stagenData->BifDataInPtr+3*nphase+2,sizeof(FloatHi)*nphase); 
        x0[nphase+2+nappend]=stagenData->BifDataInPtr[4*nphase+2]; 
      }
    } else { /* bifurcation data to be computed */
      Ret=0;
{int i,j,imin,jmin,imin0,jmin0,imin1,jmin1;
 double r,s,d,phi,beta,gamma,lambda,lambda1,omega,theta;
    JacRHS(x0,C1);
    staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
    staFunc->EigenValues(AA,Re,Im);
    /* find two pairs of multipliers with unit product */
    s=fabs(Re[0]*Re[1]-Im[0]*Im[1]-1.0)+fabs(Re[0]*Im[1]+Re[1]*Im[0]);
    imin0=0;
    jmin0=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) {
          imin0=i;
          jmin0=j;
          s=d;
        }
      }
    }
    imin1=imin0+1;
    if (imin1==nphase) imin1=0; 
    jmin1=jmin0+1;
    if (jmin1==nphase) jmin1=0; 
    s=fabs(Re[imin1]*Re[jmin1]-Im[imin1]*Im[jmin1]-1.0)+fabs(Re[imin1]*Im[jmin1]+Re[imin1]*Im[jmin1]);
    for (i=0; i<nphase; i++) {
      if (i!=imin0) {
        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) {
            imin1=i;
            jmin1=j;
            s=d;
          }
        }
      }
    }
    if (_branch == 0) { /* primary branch */
      imin=imin0;
      jmin=jmin0;
    } else { /* secondary branch */
      imin=imin1;
      jmin=jmin1;
    }

    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(C1,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) {
        staUtil->ErrMessage("Unable to find eigenvector at the initial point");
        Ret=1;
        goto polufinal;
      } 
      staFunc->CopyMatrixBlock(C1,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) {
        staUtil->ErrMessage("Unable to find eigenvector at the initial point");
        Ret=1;
        goto polufinal;
      } 
      staFunc->NormalizeVector(Qr);   
      staFunc->NormalizeVectorRelative(Qi,Qr);
      staFunc->NormalizeVector(Qi);
      
      staFunc->AddScaledVector(Qr,Qi,1.0,V);
      staFunc->NormalizeVector(V);
      staFunc->AddScaledVector(Qr,Qi,-1.0,L);             
      staFunc->NormalizeVector(L); 
      x0[nphase+2+nappend]=(lambda+lambda1)/2.0;
    } 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);
         
      /* compute real and imaginary part of the eigenvector */
      staFunc->ClearMatrix(BB);
      staFunc->CopyMatrixBlock(C1,0,0,BB,0,0,nphase,nphase);
      staFunc->CopyMatrixBlock(C1,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,VV) != 2*nphase-2) {
        staUtil->ErrMessage("Unable to find eigenvectors at the initial point");
        Ret=1;
        goto polufinal;
      }
      staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
      staFunc->CopyVectorElements(VV,nphase,Qi,0,nphase);

      /* normalize real and imaginary parts */
      d=staFunc->ScalProd(Qr,Qr);
      s=staFunc->ScalProd(Qi,Qi);
      r=staFunc->ScalProd(Qr,Qi);
      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++) {
          VV[i]=Qr[i]*cos(phi)-Qi[i]*sin(phi);
          VV[i+nphase]=Qr[i]*sin(phi)+Qi[i]*cos(phi);
        }
      }

      staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
      d=staFunc->ScalProd(Qr,Qr);
      gamma=1/sqrt(d);
      staFunc->MultiplyVectorScalar(VV,gamma,VV);
      staFunc->CopyVectorElements(VV,0,V,0,nphase);
      staFunc->CopyVectorElements(VV,nphase,L,0,nphase);
      x0[nphase+2+nappend]=cos(theta);
      staFunc->NormalizeVector(L);
    }
polufinal:
      if (Ret) {
        v1=staFunc->DestroyVector(v1);
        v2=staFunc->DestroyVector(v2);
        D=staFunc->DestroyMatrix(D);
        DE=staFunc->DestroyMatrix(DE);
        Term(FALSE);
        return 1;
      }
    }
}

  if (JacRHS(x0,C1)) return 1;  /* for safety */
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  staFunc->MultiplyMatrixMatrix(A,A,BP); 
  staFunc->AddScaledMatrix(BP,A,-2.0*x0[nphase+2+nappend],BP);
  for (i=0;i<nphase;i++) {
   BP[i][i]+=1.0;
  }
/* searching 2 rows with maximal pivots */  
  if (nphase >= 2) {
    staFunc->ClearMatrix(BP1);
    staFunc->CopyMatrixBlock(BP,0,0,BP1,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) {
      BP1[i][nphase]=i;
      BP1[nphase][i]=i;
    } 
 /* Calculation of PIVOTS - complete pivoting */
    for (k=0;k<nphase-2;k++) {
      t=0.0;
      for (i=k;i<nphase;i++) {
	for (j=k;j<nphase;j++) {
    	  if (fabs(BP1[i][j]) > t) {
            p = i;
            q = j;
            t = fabs(BP1[i][j]);
          }           
	}
      }      
      if (k!=q) { /* switch colom k with colom  q */        
        for (i=0;i<=nphase;i++) {
          t=BP1[i][k];
          BP1[i][k]=BP1[i][q];
          BP1[i][q]=t;
	}
      }
      if (k!=p) { /* switch row k with row p */       
        for (i=0;i<=nphase;i++) {
          t=BP1[k][i];         
          BP1[k][i]=BP1[p][i];
          BP1[p][i]=t;
	}
      }        
      if (fabs(BP1[k][k]) < _eps) {
        staUtil->ErrMessage("(A^2 - 2*k*A + I) has rank defect at least three");
        Term(FALSE);
        return 1;
      }
      kp1=k+1;
      for (i=kp1;i<nphase;i++) { /* perform elimination */
        t=-BP1[i][k]/BP1[k][k];
        for (j=kp1;j<nphase;j++) {
          BP1[i][j]+=t*BP1[k][j];
	}
      }
    } 
    p=(Int2)BP1[nphase-2][nphase];
    q=(Int2)BP1[nphase-1][nphase];   
    BP[p][nphase]=1.0;
    BP[q][nphase+1]=1.0;
    p=(Int2)BP1[nphase][nphase-2];
    q=(Int2)BP1[nphase][nphase-1];   
    BP[nphase][p]=1.0;
    BP[nphase+1][q]=1.0;    

/* solve systems to obtain new values for V1 and V2 */
    staFunc->ClearVector(V1);
    staFunc->ClearVector(V2);
    staFunc->ClearVector(W1);
    staFunc->ClearVector(W2);
    staFunc->CopyMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(VB1);
    VB1[nphase]=1.0;
    staFunc->Solve(BT,VB1);
    staFunc->CopyVector(VB1,V1);
    staFunc->ClearVector(VB2);
    VB2[nphase+1]=1.0;
    staFunc->Solve(BT,VB2);
    staFunc->CopyVector(VB2,V2);

/* solve transposed systems to obtain new values for W1 and W2 */
    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(PP1);
    PP1[nphase]=1.0;
    staFunc->Solve(BT,PP1);
    staFunc->CopyVector(PP1,W1);
    staFunc->ClearVector(PP2);
    PP2[nphase+1]=1.0;
    staFunc->Solve(BT,PP2);
    staFunc->CopyVector(PP2,W2);
                
    for (i=0;i<nphase;i++) {
      BP[nphase][i]=V1[i];
      BP[i][nphase]=W1[i];
    }
    for (i=0;i<nphase+2;i++) {
      BP[nphase+1][i]=V2[i];
      BP[i][nphase+1]=W2[i];
    }  
  } else {
     BP[0][nphase] = 1.0;
     BP[1][0] = 1.0;
     BP[nphase+1][nphase+1] = 1.0;
     staFunc->ClearVector(V1);
     staFunc->ClearVector(V2);
     staFunc->ClearVector(W1);
     staFunc->ClearVector(W2);
     V1[0]=1.0;     
     W1[0]=1.0;
     W2[nphase+1]=1.0;
    } 

    staFunc->CopyMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(VB1);
    VB1[nphase]=1.0;
    staFunc->Solve(BT,VB1);  
    staFunc->ClearVector(VB2);
    VB2[nphase+1]=1.0;
    staFunc->Solve(BT,VB2);

    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(PP1);
    PP1[nphase]=1.0;
    staFunc->Solve(BT,PP1);
    staFunc->ClearVector(PP2);
    PP2[nphase+1]=1.0;
    staFunc->Solve(BT,PP2);

/* compute derivatives of g11,g12,g21,g22 */
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);   
    if (stagenData->Der2) {
      H=SymHess(x0);
      staFunc->ClearVector(HV1);
      staFunc->ClearVector(HV2);
      staFunc->ClearVector(HW1);
      staFunc->ClearVector(HW2);
      staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
      staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
      staFunc->MultiplyMatrixVector(A,HV1,HV1);
      staFunc->MultiplyMatrixVector(A,HV2,HV2);
      staFunc->TransposeMatrix(A,B);
      staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
      staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
      staFunc->MultiplyMatrixVector(B,HW1,HW1);
      staFunc->MultiplyMatrixVector(B,HW2,HW2);  
      S=staFunc->CreateMatrix(2,2); 
      for (l=0; l<nphase+2+nappend; l++) {
        for (i=0,S[0][0]=0.0,S[0][1]=0.0,S[1][0]=0.0,S[1][1]=0.0;i<nphase;i++) {
          for (j=0; j<nphase; j++) {
    	    S[0][0]-=HW1[i]*VB1[j]*HessianElem(i,j,l)+PP1[i]*HV1[j]*HessianElem(i,j,l);
            S[0][0]-=2.0*x0[nphase+2+nappend]*PP1[i]*VB1[j]*HessianElem(i,j,l);
    	    S[0][1]-=HW1[i]*VB2[j]*HessianElem(i,j,l)+PP1[i]*HV2[j]*HessianElem(i,j,l);
            S[0][1]-=2.0*x0[nphase+2+nappend]*PP1[i]*VB2[j]*HessianElem(i,j,l);
    	    S[1][0]-=HW2[i]*VB1[j]*HessianElem(i,j,l)+PP2[i]*HV1[j]*HessianElem(i,j,l);
            S[1][0]-=2.0*x0[nphase+2+nappend]*PP2[i]*VB1[j]*HessianElem(i,j,l);
    	    S[1][1]-=HW2[i]*VB2[j]*HessianElem(i,j,l)+PP2[i]*HV2[j]*HessianElem(i,j,l);            
            S[1][1]-=2.0*x0[nphase+2+nappend]*PP2[i]*VB2[j]*HessianElem(i,j,l);
          }
        }
        G0[l]=S[0][0];
        G1[l]=S[0][1];
        G2[l]=S[1][0];
        G3[l]=S[1][1];  
      }
      S=staFunc->DestroyMatrix(S);
      G0[nphase+2+nappend]=2.0*staFunc->ScalProd(PP1,HV1);
      G1[nphase+2+nappend]=2.0*staFunc->ScalProd(PP1,HV2);
      G2[nphase+2+nappend]=2.0*staFunc->ScalProd(PP2,HV1);
      G3[nphase+2+nappend]=2.0*staFunc->ScalProd(PP2,HV2);
    } else {
        staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase); 
        staFunc->ClearVector(HV1);
        staFunc->ClearVector(HV2);
        staFunc->ClearVector(HW1);
        staFunc->ClearVector(HW2);
        staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
        staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
        staFunc->MultiplyMatrixVector(A,HV1,HV1);
        staFunc->MultiplyMatrixVector(A,HV2,HV2);
        staFunc->TransposeMatrix(A,B);
        staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
        staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
        staFunc->MultiplyMatrixVector(B,HW1,HW1);
        staFunc->MultiplyMatrixVector(B,HW2,HW2);  
        Jac1=staFunc->CreateMatrix(nphase,nphase);
        Jac2=staFunc->CreateMatrix(nphase,nphase);
        for (l=0; l<nphase+2+nappend; l++) {
          staFunc->ClearVector(xp);
          staFunc->CopyVectorElements(x0,0,xp,0,nphase+2);
          xp[l]+=_dx;
          staFunc->ClearVector(xp1);
          staFunc->CopyVectorElements(x0,0,xp1,0,nphase+2);
          xp1[l]-=_dx;
          JacRHS(xp,C1);
          staFunc->CopyMatrixBlock(C1,0,0,Jac1,0,0,nphase,nphase);
          JacRHS(xp1,C1);
          staFunc->CopyMatrixBlock(C1,0,0,Jac2,0,0,nphase,nphase);
       /* g11 */
          staFunc->CopyVectorElements(VB1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2); 
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G0[l]=-staFunc->ScalProd(HW1,K3);
          G0[l]-=2.0*x0[nphase+2+nappend]*staFunc->ScalProd(PP1,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G0[l]-=staFunc->ScalProd(K3,HV1);                                   
       /* g12 */
          staFunc->CopyVectorElements(VB2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1);  
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G1[l]=-staFunc->ScalProd(HW1,K3);
          G1[l]-=2.0*x0[nphase+2+nappend]*staFunc->ScalProd(PP1,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G1[l]-=staFunc->ScalProd(K3,HV2);                                   
       /* g21 */
          staFunc->CopyVectorElements(VB1,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G2[l]=-staFunc->ScalProd(HW2,K3);
          G2[l]-=2.0*x0[nphase+2+nappend]*staFunc->ScalProd(PP2,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);   
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);   
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G2[l]-=staFunc->ScalProd(K3,HV1);                                   
       /* g22 */
          staFunc->CopyVectorElements(VB2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
          staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G3[l]=-staFunc->ScalProd(HW2,K3);
          G3[l]-=2.0*x0[nphase+2+nappend]*staFunc->ScalProd(PP2,K3);
          staFunc->TransposeMatrix(Jac1,B);
          staFunc->CopyVectorElements(PP2,0,K3,0,nphase);
          staFunc->MultiplyMatrixVector(B,K3,K1);  
          staFunc->TransposeMatrix(Jac2,B);
          staFunc->MultiplyMatrixVector(B,K3,K2);   
          staFunc->AddScaledVector(K1,K2,-1.0,K3);
          staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
          G3[l]-=staFunc->ScalProd(K3,HV2);                                             
        }
        Jac1=staFunc->DestroyMatrix(Jac1);
        Jac2=staFunc->DestroyMatrix(Jac2);
        G0[nphase+2+nappend]=2.0*staFunc->ScalProd(PP1,HV1);
        G1[nphase+2+nappend]=2.0*staFunc->ScalProd(PP1,HV2);
        G2[nphase+2+nappend]=2.0*staFunc->ScalProd(PP2,HV1);
        G3[nphase+2+nappend]=2.0*staFunc->ScalProd(PP2,HV2);        
      }
    staFunc->ClearMatrix(BT);
    staFunc->ClearMatrix(JM1);
    staFunc->CopyMatrixBlock(C1,0,0,BT,0,0,nphase,nphase+2);
    for (i=0; i<nphase; i++) BT[i][i]-=1.0;
/* placement of ones determined by partial pivoting */
    PIVOT=staFunc->CreateVector(nphase+2);
    staFunc->ClearVector(PIVOT);
    staFunc->TransposeMatrix(BT,JM1);
    PIVOT[nphase+1]=1.0;
    for (k=0;k<=nphase;k++) {
     kp1=k+1;
     p = k;
     for (i=kp1;i<=nphase+1;i++) {
      if (fabs(JM1[i][k]) > fabs(JM1[p][k]))
        p = i;
     }
     PIVOT[k]=p;
     if (p!=k) PIVOT[nphase+1]=-PIVOT[nphase+1];
     t=JM1[p][k];
     JM1[p][k]=JM1[k][k];
     JM1[k][k]=t;
     if (t==0.0) continue;
     for (i=kp1;i<=nphase+1;i++) JM1[i][k]=-JM1[i][k]/t;
     for (j=kp1;j<=nphase+1;j++) {
       t=JM1[p][j];
       JM1[p][j]=JM1[k][j];
       JM1[k][j]=t;
       if (t!=0.0) for (i=kp1;i<=nphase+1;i++) JM1[i][j]+=JM1[i][k]*t;
     }
    }
    staFunc->ClearVector(TV1);
    for (i=0;i<nphase+2;i++) TV1[i]=i;
    for (k=0;k<nphase+1;k++) {
      i=(Int2)PIVOT[k];
      t=TV1[i];
      TV1[i]=TV1[k];
      TV1[k]=t;
    }
    p=(Int2)TV1[nphase];
    q=(Int2)TV1[nphase+1];  
    PIVOT=staFunc->DestroyVector(PIVOT);
/* end of pivot information */

    staFunc->ClearMatrix(JM1);
    staFunc->CopyMatrixBlock(C1,0,0,JM1,0,0,nphase,nphase+2);
    for (i=0; i<nphase; i++) JM1[i][i]-=1.0;
    JM1[nphase][p]=1.0;
    JM1[nphase+1][q]=1.0;
    staFunc->ClearVector(TV1);
    TV1[nphase]=1.0;
    staFunc->Decomp(JM1);
    staFunc->Solve(JM1,TV1);
    staFunc->ClearMatrix(TVM);
    for (i=0;i<nphase+2;i++) TVM[i][0] = TV1[i];
    staFunc->ClearVector(TV1);
    TV1[nphase+1]=1.0;  
    staFunc->Solve(JM1,TV1);
    for (i=0;i<nphase+2;i++) TVM[i][1] = TV1[i];    
    TVM[nphase+2][2]=1.0;

    staFunc->TransposeMatrix(TVM,TVMT);
    staFunc->MultiplyMatrixMatrix(TVMT,TVM,HM); 
    staFunc->Decomp(HM);
    staFunc->MultiplyMatrixVector(TVMT,G0,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G0);  
    staFunc->MultiplyMatrixVector(TVMT,G1,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G1);      
    staFunc->MultiplyMatrixVector(TVMT,G2,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G2);
    staFunc->MultiplyMatrixVector(TVMT,G3,HV);  
    staFunc->Solve(HM,HV);
    staFunc->MultiplyMatrixVector(TVM,HV,G3);       
  
/* Determine component with maximum norm => gives m1,n1 */
    G[0][0] = staFunc->NormOfVector(G0);
    G[0][1] = staFunc->NormOfVector(G1);
    G[1][0] = staFunc->NormOfVector(G2);
    G[1][1] = staFunc->NormOfVector(G3);    
    m1=0;
    n1=0;      
    if (G[0][1] > G[m1][n1]) n1=1;  
    if (G[1][0] > G[m1][n1]) { m1=1; n1=0; }
    if (G[1][1] > G[m1][n1]) { m1=1; n1=1; }
    for (i=0;i<nphase+3;i++) {
      GM1[i][0] = G0[i]; 
      GM1[i][1] = G1[i];
      GM1[i][2] = G2[i];
      GM1[i][3] = G3[i];
    }
/* G0 has biggest norm ; necessary for projection */ 
    for (i=0;i<nphase+3;i++) {
      G0[i] = GM1[i][2*m1+n1];
    }      
    staFunc->NormalizeVector(G0);
/* G1,G2,G3 are projected on G0 */ 
    for (j=0;j<4;j++) {
      if (j != (2*m1+n1)) {
       for (i=0;i<nphase+3;i++) { /* neem een vector, niet de grootste */
          G1[i] = GM1[i][j];
       }
       scal = staFunc->ScalProd(G1,G0);
       staFunc->AddScaledVector(G1,G0,-scal,G1);
       for (i=0;i<nphase+3;i++) { /* zet nieuwe vector terug */
         GM1[i][j] = G1[i];
       }
      }
    }      
/* compute biggest norm between G1,G2 and G3  => gives m2,n2 */
    for (i=0;i<nphase+3;i++) {
      G0[i] = GM1[i][0];
      G1[i] = GM1[i][1];
      G2[i] = GM1[i][2];
      G3[i] = GM1[i][3];
    }
    G[0][0] = staFunc->NormOfVector(G0);
    G[0][1] = staFunc->NormOfVector(G1);
    G[1][0] = staFunc->NormOfVector(G2);
    G[1][1] = staFunc->NormOfVector(G3);    
    G[m1][n1]=-1.0;
    m2=0;
    n2=0;
    if (G[0][1] > G[m2][n2]) n2=1;  
    if (G[1][0] > G[m2][n2]) { m2=1; n2=0; }
    if (G[1][1] > G[m2][n2]) { m2=1; n2=1; }    

/* Compute a tangent vector to the NS curve */
    ns2_DefNS(x0,v0);
    Ret=ns2_JacDefNS(x0,D);
    for (i=0; i<nphase+3+nappend; i++) {
      v1[i]=1.0;
      staFunc->AppendMatrixVector(D,v1,DE);
      Ret=staFunc->Decomp(DE);
      if (Ret == 0) {
        v2[nphase+2+nappend]=stagenData->Forward ? 1 : -1;
        staFunc->Solve(DE,v2);
        staFunc->CopyVector(v2,v0);
        break;
      };
      v1[i]=0.0;
    };
    if (Ret) {
       staUtil->ErrMessage
         ("Unable to find a tangent vector to the Neimark-Sacker curve");   
       v1=staFunc->DestroyVector(v1);
       v2=staFunc->DestroyVector(v2);
       D=staFunc->DestroyMatrix(D);
       DE=staFunc->DestroyMatrix(DE);
       Term(FALSE);
       return 1;
    }  
    staFunc->NormalizeVector(v0);
    v1=staFunc->DestroyVector(v1);
    v2=staFunc->DestroyVector(v2);
    D=staFunc->DestroyMatrix(D);
    DE=staFunc->DestroyMatrix(DE);
  }
  /* Initialize continuation */
  ((ContDataPtr)stagenptr->GenPtr)->X0=x0;
  ((ContDataPtr)stagenptr->GenPtr)->V0=v0;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveSing=staData->ActiveSing;
  ((ContDataPtr)stagenptr->GenPtr)->ActiveUserTest=staData->ActiveUserTest;
  ((ContDataPtr)stagenptr->GenPtr)->TypeTestFuns=staData->TypeTestFunc;
  ((ContDataPtr)stagenptr->GenPtr)->ZeroTestFuns=staData->ZeroTestFunc;
  /* That's all */
  return 0;
}

/************************/
/* 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<2+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) ns2_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) ns2_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) ns2_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;
}

/*************************************/
/* Defining function for fixed point */
Local(Int2) ns2_DefFixedPoint (Vector x, Vector y) {
Int2 i;
  ns2_DefRHS(x,y);
  for (i=0; i<nphase; i++) y[i]-=x[i];
  if (nappend) ns2_DefUser(x,y+nphase);
  return 0;
}

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

/****************************************************/
/* Jacobian of defining conditions for fixed points */
Local(Int2) ns2_JacDefFixedPoint (Vector x, Matrix Jac) {
Int2 i;
  staFunc->ClearMatrix(Jac);
  JacRHS(x, Jac);
  for (i=0; i<nphase; i++) Jac[i][i]-=1.0;
  if (nappend) staFunc->Der(ns2_DefUser, nappend, x, dx, Jac+nphase);  
  return 0;
}

/*************************/
/* NS defining condition */
Entry(Int2) ns2_DefBord (Vector x, Vector y) {
Int2 i,Ret;

  if (JacRHS(x,C1)) return 1;	  /* A,C1,BP,V1,V2,W1,W2 - global */
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  staFunc->MultiplyMatrixMatrix(A,A,BP);
  staFunc->AddScaledMatrix(BP,A,-2.0*x[nphase+2+nappend],BP);
  for (i=0;i<nphase;i++) {
    BP[i][i]+=1.0;
  } 
/* Fill in the bordered matrix (A^2 - 2 cos(theta) A + I) */
  for (i=0;i<nphase;i++) {
    BP[nphase][i]=V1[i];
    BP[i][nphase]=W1[i];
  }
  for (i=0;i<nphase+2;i++) {
    BP[nphase+1][i]=V2[i];
    BP[i][nphase+1]=W2[i];
  }   
  staFunc->CopyMatrix(BP,BT);
  Ret=staFunc->Decomp(BT);
  if (Ret) return 1;
  staFunc->ClearVector(VB1);
  VB1[nphase]=1.0;
  staFunc->Solve(BT,VB1);
  staFunc->CopyVector(VB1,Q1);
  staFunc->ClearVector(VB2);
  VB2[nphase+1]=1.0;
  staFunc->Solve(BT,VB2);
  staFunc->CopyVector(VB2,Q2);

  G[0][0]=VB1[nphase];    /* g11 */
  G[1][0]=VB1[nphase+1];  /* g21 */
  G[0][1]=VB2[nphase];    /* g12 */
  G[1][1]=VB2[nphase+1];  /* g22 */
  
  y[0]=G[m1][n1];
  y[1]=G[m2][n2];
  return 0;
}

/*************************************/
/* Jacobian of NS defining condition */
Local(Int2) JacDefBord (Vector x, Matrix Jac) {
int i,j,l;
Int2 Ret;
Matrix S;

  staFunc->ClearMatrix(Jac);
  if (stagenData->Der2) {
/* creation of DefBord Jacobian matrix Jac using symbolic second derivatives */
    H=SymHess(x);
/* take global vectors VB1 and VB2 computed in DefBord */
/* compute PP1 cand PP2 using the tranposed squared Jacobian */
    staFunc->TransposeMatrix(BP,BT);
    Ret=staFunc->Decomp(BT);
    if (Ret) return 1;
    staFunc->ClearVector(PP1);
    PP1[nphase]=1.0;
    staFunc->Solve(BT,PP1);
    staFunc->ClearVector(PP2);
    PP2[nphase+1]=1.0;
    staFunc->Solve(BT,PP2);
 
    staFunc->ClearVector(HV1);
    staFunc->ClearVector(HV2);
    staFunc->ClearVector(HW1);
    staFunc->ClearVector(HW2);     
    staFunc->CopyVectorElements(VB1,0,VB,0,nphase);
    staFunc->MultiplyMatrixVector(A,VB,HV1);
    staFunc->CopyVectorElements(VB2,0,VB,0,nphase);
    staFunc->MultiplyMatrixVector(A,VB,HV2);
    staFunc->TransposeMatrix(A,B);
    staFunc->CopyVectorElements(PP1,0,VB,0,nphase);
    staFunc->MultiplyMatrixVector(B,VB,HW1);
    staFunc->CopyVectorElements(PP2,0,VB,0,nphase);
    staFunc->MultiplyMatrixVector(B,VB,HW2);  
    S=staFunc->CreateMatrix(2,2); 
    for (l=0; l<nphase+2+nappend; l++) {
      for (i=0,S[0][0]=0.0,S[0][1]=0.0,S[1][0]=0.0,S[1][1]=0.0;i<nphase;i++) {
        for (j=0; j<nphase; j++) {
    	  S[0][0]-=HW1[i]*VB1[j]*HessianElem(i,j,l)+PP1[i]*HV1[j]*HessianElem(i,j,l);
    	  S[0][1]-=HW1[i]*VB2[j]*HessianElem(i,j,l)+PP1[i]*HV2[j]*HessianElem(i,j,l);
    	  S[1][0]-=HW2[i]*VB1[j]*HessianElem(i,j,l)+PP2[i]*HV1[j]*HessianElem(i,j,l);
    	  S[1][1]-=HW2[i]*VB2[j]*HessianElem(i,j,l)+PP2[i]*HV2[j]*HessianElem(i,j,l);            
    	  S[0][0]+=2.0*x[nphase+2+nappend]*PP1[i]*VB1[j]*HessianElem(i,j,l);
    	  S[0][1]+=2.0*x[nphase+2+nappend]*PP1[i]*VB2[j]*HessianElem(i,j,l);    	  
    	  S[1][0]+=2.0*x[nphase+2+nappend]*PP2[i]*VB1[j]*HessianElem(i,j,l);    	  
    	  S[1][1]+=2.0*x[nphase+2+nappend]*PP2[i]*VB2[j]*HessianElem(i,j,l);    	  
        }

      }
      Jac[0][l]=S[m1][n1];
      Jac[1][l]=S[m2][n2];
    }
    S=staFunc->DestroyMatrix(S);
    G0[nphase+2+nappend]=2.0*staFunc->ScalProd(HW1,VB1);
    G1[nphase+2+nappend]=2.0*staFunc->ScalProd(HW1,VB2);
    G2[nphase+2+nappend]=2.0*staFunc->ScalProd(HW2,VB1);
    G3[nphase+2+nappend]=2.0*staFunc->ScalProd(HW2,VB2);
    GM1[nphase+2+nappend][0]=G0[nphase+2+nappend];
    GM1[nphase+2+nappend][1]=G1[nphase+2+nappend];
    GM1[nphase+2+nappend][2]=G2[nphase+2+nappend];
    GM1[nphase+2+nappend][3]=G3[nphase+2+nappend];
    Jac[0][nphase+2+nappend]=GM1[nphase+2+nappend][2*m1+n1];
    Jac[1][nphase+2+nappend]=GM1[nphase+2+nappend][2*m2+n2];
  } else {
    if (staFunc->Der(ns2_DefBord, 2, x, dx, Jac)) return 1;
  }
  return(0);
}

/***********************************/
/* Defining functions for NS curve */

Entry(Int2) ns2_DefNS (Vector x, Vector y) {
  if (ns2_DefFixedPoint(x,y)) return 1;
  if (ns2_DefBord(x,y+nphase+nappend)) return 2;
  return 0;
}

/***********************************************/
/* Jacobian of defining functions for NS curve */  
Entry(Int2) ns2_JacDefNS (Vector x, Matrix Jac) {
  staFunc->ClearMatrix(Jac);
  if (ns2_JacDefFixedPoint(x,C)) return 1;
  staFunc->CopyMatrixBlock(C,0,0,Jac,0,0,nphase+nappend,nphase+2+nappend);
  if (JacDefBord(x,C2)) return 2;
  staFunc->CopyMatrixBlock(C2,0,0,Jac,nphase+nappend,0,2,nphase+3+nappend);

  return 0;
}

/*********************/
/* Default processor */
#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) ns2_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, Ind==0 */
  ns2_DefNS (x, xp1);  /* compute global Q for test functions */
  MemCopy(x1,x,sizeof(FloatHi)*(nphase+3+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 */
typedef EntryPtr(Int2,TestFuncPtr,(Vector xx, Vector vv, FloatHiPtr TV));
typedef EntryPtr(Int2,ProcFuncPtr,(Int2 Ind, Vector xx, Vector vv, CharPtr *msg));

Local(Char) buf[80];

/* first Lyapunov coefficient */
Entry(Int2) ns2_TestDegenNS(Vector x, Vector v, FloatHiPtr res) {
VAR beta,gamma,phi,lambda,lambda2,omega,omega2,d,s,s1,r,r1,h;
Int2 i,j,k,l,Ret;  

  *res=1.0;
  if (nphase < 2) goto end;
  if (fabs(x[nphase+2+nappend]) >= 1.0) { 	/* neutral saddle */
    *res=777;
  } else {   				/* Neimark-Sacker */
    lambda=x[nphase+2+nappend];         /* cos(theta) = kappa */
    omega=sqrt(1.0-lambda*lambda);      /* sin(theta) */
    lambda2=2*lambda*lambda-1.0	;	/* cos(2*theta) */
    omega2=2*omega*lambda;		/* sin(2*theta) */
    
    if (JacRHS(x,C1)) return 1;	    /* C1,A,Qr,Qi - global */

    /* compute real(Qr) and imaginary(Qi) part of the eigenvector */
    staFunc->ClearMatrix(BB);
    staFunc->CopyMatrixBlock(C1,0,0,BB,0,0,nphase,nphase);
    staFunc->CopyMatrixBlock(C1,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,VV) != 2*nphase-2) {
      *res=1.0;
      goto end;
    }
    staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
    staFunc->CopyVectorElements(VV,nphase,Qi,0,nphase);

    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) A[i][i]-=lambda;

    staFunc->MultiplyMatrixVector(A,Qr,Qi);
    staFunc->MultiplyVectorScalar(Qi,-1.0/omega,Qi);
    staFunc->CopyVectorElements(Qr,0,VV,0,nphase);
    staFunc->CopyVectorElements(Qi,0,VV,nphase,nphase);
    
    /* normalize real and imaginary parts <Qr,Qr>=1, <Qi,Qi>=0 */
    d=staFunc->ScalProd(Qr,Qr);
    s=staFunc->ScalProd(Qi,Qi);
    r=staFunc->ScalProd(Qr,Qi);
    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++) {
          VV[i]=Qr[i]*cos(phi)-Qi[i]*sin(phi);
          VV[i+nphase]=Qr[i]*sin(phi)+Qi[i]*cos(phi);
        }
    }
    staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
    d=staFunc->ScalProd(Qr,Qr);
    gamma=1/sqrt(d);
    staFunc->MultiplyVectorScalar(VV,gamma,VV);
    staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
    staFunc->CopyVectorElements(VV,nphase,Qi,0,nphase);   

    /* compute real and imaginary part of the adjoint eigenvector */
    staFunc->CopyMatrixBlock(C1,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,VV) != 2*nphase-2) {
       *res=1.0;
       test[0]=*res;
       return 1;
    }
    staFunc->CopyVectorElements(VV,0,Pr,0,nphase);
    staFunc->CopyVectorElements(VV,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++) {
      VV[i]=Pr[i]*cos(phi)-Pi[i]*sin(phi);
      VV[i+nphase]=Pr[i]*sin(phi)+Pi[i]*cos(phi);
    }

    staFunc->CopyVectorElements(VV,0,Pr,0,nphase);
    staFunc->CopyVectorElements(VV,nphase,Pi,0,nphase);
    r=staFunc->ScalProd(Pr,Qr);
    h=staFunc->ScalProd(Pi,Qi);
    beta=1.0/(r+h);
    staFunc->MultiplyVectorScalar(VV,beta,VV);
    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);
        ns2_DefRHS(x,aa);
        staFunc->AddScaledVector(x,v0,ddx,x1);
        ns2_DefRHS(x1,bb);
        staFunc->AddScaledVector(x,v0,-ddx,x1);
        ns2_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);
        ns2_DefRHS(x,rr);
        staFunc->AddScaledVector(x,v0,ddx,x1);
        ns2_DefRHS(x1,bb);
        staFunc->AddScaledVector(x,v0,-ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,rr);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,rr);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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(C1,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) {
        *res=1.0;
        test[0]=*res;
        return 1;
      }
      staFunc->Solve(AA,rr);
      staFunc->MultiplyVectorScalar(rr,-2.0,rr);		/* rr=-2*r */

      staFunc->AddScaledVector(bb,aa,-1.0,dd);
      staFunc->CopyVectorElements(dd,0,VV,0,nphase);
      staFunc->CopyVectorElements(cc,0,VV,nphase,nphase);
      staFunc->ClearMatrix(BB);
      staFunc->CopyMatrixBlock(C1,0,0,BB,0,0,nphase,nphase);
      staFunc->CopyMatrixBlock(C1,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) {
        *res=1.0;
        test[0]=*res;
        return 1;
      }
      staFunc->Solve(BB,VV);
      staFunc->CopyVectorElements(VV,0,aa,0,nphase);          	/* aa=Sr */
      staFunc->CopyVectorElements(VV,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);
        ns2_DefRHS(x1,dd);
	s=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s1=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s1-=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	s1+=staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	r=(2.0/3.0)*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r-=2.0*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r+=2.0*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	r1=(2.0/3.0)*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r1-=2.0*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r1+=2.0*staFunc->ScalProd(Pr,dd);
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	r+=(2.0/3.0)*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r-=2.0*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r+=2.0*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	r1-=(2.0/3.0)*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r1+=2.0*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r1-=2.0*staFunc->ScalProd(Pi,dd);
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	r+=staFunc->ScalProd(bb,dd)/6;
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r-=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r+=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	r1+=staFunc->ScalProd(bb,dd)/6;
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r1-=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r1+=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	r+=staFunc->ScalProd(bb,dd)/6;
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r-=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r+=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        ns2_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);
        ns2_DefRHS(x1,dd);
	r1-=staFunc->ScalProd(bb,dd)/6;
        staFunc->AddScaledVector(x,v0,dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r1+=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r1-=staFunc->ScalProd(bb,dd)/2;
        staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
        ns2_DefRHS(x1,dd);
	r1+=staFunc->ScalProd(bb,dd)/6;
        r1/=(dddx*dddx*dddx);                                         
    }
    *res=(lambda*(s+r)+omega*(s1+r1))/2/(1.0-lambda);
  }
end:	       
  test[0]=*res;
  return 0;
}


Entry(Int2) ns2_ProcDegenNS(Int2 Ind, Vector x, Vector v, CharPtr *msg) {
  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Degenerate Neimark-Sacker (?)");  
  else {
/*  provide bifurcation data for DN -> NS starter */
/*  BifVector[0]:			 _itmap */
/*  BifVector[1] to [nphase]:  	         NS continuation vector */
/*  BifVector[nphase+1] to [2*nphase]: 	 global vector L */
/*  BifVector[2*nphase+1]:               kappa=cos(theta)  */

    BifVector[0]=(VAR)_itmap;
    MemCopy(BifVector+1,VB1,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+nphase+1,VB2,sizeof(FloatHi)*nphase);
    BifVector[2*nphase+1]=x1[nphase+2+nappend];
    stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+2);
    sprintf(buf,"Degenerate Neimark-Sacker: cos(theta)=%g",BifVector[2*nphase+1]);    
  }  
  *msg=buf;
  return 0;
}

/* kappa-1 */
Entry(Int2) ns2_TestR1(Vector x, Vector v, FloatHiPtr res) {
  *res=1.0;
  if (nphase < 2) goto end;
  *res=x[nphase+2+nappend]-1.0;
end:
  test[1]=*res;  
  return 0;
}

Entry(Int2) ns2_ProcR1(Int2 Ind, Vector x, Vector v, CharPtr *msg) {
  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Resonance 1:1 (?)");
  else {
Vector v1,W0,W1,F1;
Matrix D;
Int2 i,j,k;
double a20,b20,b11,f0,f1,f2,f3,f4,f5,d;
    v1=staFunc->CreateVector(2*nphase);
    W0=staFunc->CreateVector(nphase);
    W1=staFunc->CreateVector(nphase);
    D=staFunc->CreateMatrix(2*nphase,2*nphase);
    F1=staFunc->CreateVector(nphase);

/*   compute (generalized)eigenvectors   */
    if (JacRHS (x,C1)) return 1;
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) A[i][i]-=1.0;
    
    staFunc->CopyMatrixBlock(A,0,0,D,0,0,nphase,nphase);
    staFunc->CopyMatrixBlock(A,0,0,D,nphase,nphase,nphase,nphase);
    for (i=0; i<nphase; i++) D[i+nphase][i]=-1;
    i=staFunc->SngSolve(D,_eps,v1);
    if (i != 2*nphase-2) {
      sprintf(buf,"1:1 Resonance: eigenvectors (?)");
      goto final;
    }
    staFunc->CopyVectorElements(v1,0,V,0,nphase);       /* eigenvector V0 */
    d=sqrt(staFunc->ScalProd(V,V));
    staFunc->NormalizeVector(V);
    staFunc->CopyVectorElements(v1,nphase,W,0,nphase);  /* generalized eigenvector V1 */
    staFunc->MultiplyVectorScalar(W,1.0/d,W);
    staFunc->AddScaledVector(W,V,-staFunc->ScalProd(V,W),W);

/*  provide bifurcation data for R1->* starters */
/*  BifVector[0]:			 _itmap */
/*  BifVector[1] to [nphase]:            eigenvector V (A*V=V, <V,V>=1)*/
/*  BifVector[nphase+1]:                 kappa=1  */
    BifVector[0]=(VAR)_itmap;
    MemCopy(BifVector+1,V,sizeof(FloatHi)*nphase);
    BifVector[nphase+1]=1.0;
    stagenData->BifDataOutLen=sizeof(FloatHi)*(nphase+2);
    
/* compute the adjoint eigenvectors */
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) A[i][i]-=1.0;
    staFunc->TransposeMatrix(A,B);
    staFunc->ClearMatrix(D);
    staFunc->CopyMatrixBlock(B,0,0,D,0,0,nphase,nphase);
    staFunc->CopyMatrixBlock(B,0,0,D,nphase,nphase,nphase,nphase);
    for (i=0; i<nphase; i++) D[i+nphase][i]=-1;
    i=staFunc->SngSolve(D,_eps,v1);
    if (i != 2*nphase-2) {
      sprintf(buf,"1:1 Resonance: adjoint eigenvectors (?)");
      goto final;
    }
    staFunc->CopyVectorElements(v1,0,W1,0,nphase);       /* adjoint eigenvector W1 */
    staFunc->CopyVectorElements(v1,nphase,W0,0,nphase);  /* adjoint generalized eigenvector W0 */
    d=staFunc->ScalProd(V,W0);
    staFunc->MultiplyVectorScalar(W1,1.0/d,W1);
    staFunc->AddScaledVector(W0,W1,-staFunc->ScalProd(W0,W),W0);
    staFunc->MultiplyVectorScalar(W0,1.0/d,W0);

/* compute the coefficients of the normal form map: x'=x+y; y'=Y+ax^2+bxy */
/*  a20=<W0,B(V,V)>, b20=<W1,B(V,V)>,  b11=<W1,B(V,W)> */
     if (stagenData->Der2)  {
        H=SymHess(x);
        for (a20=0,b20=0,b11=0,i=0; i<nphase; i++) {
           for (j=0; j<nphase; j++) {
              for (k=0; k<nphase; k++) { 
	         d=HessianElem(i,j,k);
                 a20+=W0[i]*d*V[j]*V[k];
                 b20+=W1[i]*d*V[j]*V[k];
                 b11+=W1[i]*d*V[j]*W[k];
	      }
           }
        }

     } else {
	staFunc->ClearVector(v0);
/*  a20=<W0,B(V,V)>, b20=<W1,B(V,V)> */
	staFunc->CopyVectorElements(V,0,v0,0,nphase);
        ns2_DefRHS(x,F1);
        f0=staFunc->ScalProd(W0,F1);
        f3=staFunc->ScalProd(W1,F1);
        staFunc->AddScaledVector(x,v0,ddx,x1);
        ns2_DefRHS(x1,F1);
        f1=staFunc->ScalProd(W0,F1);
        f4=staFunc->ScalProd(W1,F1);
        staFunc->AddScaledVector(x,v0,-ddx,x1);
        ns2_DefRHS(x1,F1);
        f2=staFunc->ScalProd(W0,F1);
        f5=staFunc->ScalProd(W1,F1);
	a20=(f1-2*f0+f2)/ddx/ddx;
	b20=(f4-2*f3+f5)/ddx/ddx;

/*  b11=<W1,B(V,W)>  */
        staFunc->AddScaledVector(V,W,1.0,F1);
        staFunc->CopyVectorElements(F1,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        ns2_DefRHS(x1,F1);
	b11=staFunc->ScalProd(W1,F1);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_DefRHS(x1,F1);
	b11+=staFunc->ScalProd(W1,F1);
        staFunc->AddScaledVector(V,W,-1.0,F1);
        staFunc->CopyVectorElements(F1,0,v0,0,nphase);
        staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
        ns2_DefRHS(x1,F1);
	b11-=staFunc->ScalProd(W1,F1);
        staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
        ns2_DefRHS(x1,F1);
	b11-=staFunc->ScalProd(W1,F1);
	b11/=(ddx*ddx); 
     }
    if (b20>=0)
      sprintf(buf,"1:1 Resonance: a=%g, b=%g",b20/2,a20+b11-b20);
    else
      sprintf(buf,"1:1 Resonance: a=%g, b=%g",-b20/2,-(a20+b11-b20));
final:
    D=staFunc->DestroyMatrix(D);
    v1=staFunc->DestroyVector(v1);
    W0=staFunc->DestroyVector(W0);
    W1=staFunc->DestroyVector(W1);
    F1=staFunc->DestroyVector(F1);
  }
  *msg=buf;
  return 0;
}

/* det(A-I) */
Entry(Int2) ns2_TestFoldNS(Vector x, Vector v, FloatHiPtr res) {
Int2 i;
  *res=1.0;
  if (nphase < 2) goto end;
  if (JacRHS (x,C1)) return 1;
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  for (i=0; i<nphase; i++) A[i][i]-=1.0;
  if (staFunc->Decomp(A)) {
    *res=0.0;
  } else {
    staFunc->GetMatrixData(A,res,NULL);
  } 
end:
  test[2]=*res;  
  return 0;
}

Entry(Int2) ns2_ProcFoldNS(Int2  Ind, Vector x, Vector v, CharPtr *msg) {
Int2 i;
  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Fold - Neimark-Sacker (?)"); 
  else {
/* reevaluate Jacobian */
   if (JacRHS (x,C1)) return 1;
   staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);

   staFunc->CopyVectorElements(VB1,0,V,0,nphase);
   staFunc->CopyVectorElements(VB2,0,W,0,nphase);
  
/* orthogonalize and normalize eigenvectors */
    staFunc->NormalizeVector(V);
    staFunc->AddScaledVector(W,V,-staFunc->ScalProd(V,W),L);
    staFunc->NormalizeVector(L);
  
/* provide bifurcation data for FN -> * starters */
    if (JacRHS (x,C1)) return 1;
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) A[i][i]-=1.0;
    if (staFunc->SngSolve(A,_eps,W) != nphase-1) {
      sprintf(buf,"Fold + Neimark-Sacker: fold vector (?)"); 
      *msg=buf;
      return 0;
    }   
    staFunc->NormalizeVector(W);
/*  BifVector[0]:                        _itmap */
/*  BifVector[1] to [nphase]:            fold-vector (AW=W, <W,W>=1) */
/*  BifVector[nphase+1] to [2*nphase]:   NS continuation vector */
/*  BifVector[2*nphase+1] to [3*nphase]: global vector L */
/*  BifVector[3*nphase+1]:               kappa */
    BifVector[0]=(VAR)_itmap;
    MemCopy(BifVector+1,W,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+nphase+1,V,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+2*nphase+1,L,sizeof(FloatHi)*nphase);
    BifVector[3*nphase+1]=x[nphase+2+nappend];
    stagenData->BifDataOutLen=sizeof(FloatHi)*(3*nphase+2);
    
    if (BifVector[3*nphase+1]>0) {
      sprintf(buf,"Fold + Neimark-Sacker: cos(theta)=%g",BifVector[3*nphase+1]);
    } else {
      sprintf(buf,"Fold + neutral saddle: kappa=%g",BifVector[3*nphase+1]);
    }      
  }  
  *msg=buf;
  return 0;
}

/* det(A+I) */
Entry(Int2) ns2_TestFlipNS(Vector x, Vector v, FloatHiPtr res) {
Int2 i;
  *res=1.0;
  if (nphase < 2) goto end;
  if (JacRHS (x,C1)) return 1;
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  for (i=0; i<nphase; i++) A[i][i]+=1.0;
  if (staFunc->Decomp(A)) {
    *res=0.0;
  } else {
    staFunc->GetMatrixData(A,res,NULL);
  } 
end:
  test[3]=*res;  
  return 0;
}

Entry(Int2) ns2_ProcFlipNS(Int2  Ind, Vector x, Vector v, CharPtr *msg) {
Int2 i;
  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Flip - Neimark-Sacker (?)"); 
  else {
/* reevaluate Jacobian */
   if (JacRHS (x,C1)) return 1;
   staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);

   staFunc->CopyVectorElements(VB1,0,V,0,nphase);
   staFunc->CopyVectorElements(VB2,0,W,0,nphase);
  
/* orthogonalize and normalize eigenvectors */
    staFunc->NormalizeVector(V);
    staFunc->AddScaledVector(W,V,-staFunc->ScalProd(V,W),L);
    staFunc->NormalizeVector(L);
  
/* provide bifurcation data for FN -> * starters */
    if (JacRHS (x,C1)) return 1;
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) A[i][i]+=1.0;
    if (staFunc->SngSolve(A,_eps,W) != nphase-1) {
      sprintf(buf,"Flip + Neimark-Sacker: fold vector (?)"); 
      *msg=buf;
      return 0;
    }   
    staFunc->NormalizeVector(W);
/*  BifVector[0]:                        _itmap */
/*  BifVector[1] to [nphase]:            flip-vector (AW=-W, <W,W>=1) */
/*  BifVector[nphase+1] to [2*nphase]:   NS continuation vector */
/*  BifVector[2*nphase+1] to [3*nphase]: global vector L */
/*  BifVector[3*nphase+1]:               kappa */
    BifVector[0]=(VAR)_itmap;
    MemCopy(BifVector+1,W,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+nphase+1,V,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+2*nphase+1,L,sizeof(FloatHi)*nphase);
    BifVector[3*nphase+1]=x[nphase+2+nappend];
    stagenData->BifDataOutLen=sizeof(FloatHi)*(3*nphase+2);
    
    if (BifVector[3*nphase+1]>0) {
      sprintf(buf,"Flip + Neimark-Sacker: cos(theta)=%g",BifVector[3*nphase+1]);
    } else {
      sprintf(buf,"Flip + neutral saddle: kappa=%g",BifVector[3*nphase+1]);
    }      
  }  
  *msg=buf;
  return 0;
}

/*  Double NS test function */
Entry(Int2) ns2_TestDoubleNS(Vector x, Vector v, FloatHiPtr res) {
VAR beta,kappa,lambda,lambda1,phi,omega,d,s,r,maxV,maxW;
Int2 i,ii,j,jj,Ret,imaxV,imaxW;  
  *res=1.0;
  if (nphase < 4) goto end; 
  if (JacRHS(x,C1)) { 	    /* C1,A,AA,B,BB,V,W - global */   
     test[4]=*res;
     return 1;
  }
  if (fabs(x[nphase+2+nappend]) >= 1.0) {
/* neutral saddle */
     kappa=x[nphase+2+nappend];
     lambda=kappa+sqrt(kappa*kappa-1.0);
     lambda1=kappa-sqrt(kappa*kappa-1.0);

/* compute real adjoint eigenvectors associated with lambda and lambda1=1/lambda */
     staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
     staFunc->TransposeMatrix(AA,AT);
     for (i=0; i<nphase; i++) AT[i][i]-=lambda;
     if (staFunc->SngSolve(AT,_eps,V) != nphase-1) {
        *res=1.0;
        test[4]=*res;
        return 1;
     } 
     staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
     staFunc->TransposeMatrix(AA,AT);
     for (i=0; i<nphase; i++) AT[i][i]-=lambda1;
     if (staFunc->SngSolve(AT,_eps,W) != nphase-1) {
        *res=1.0;
        test[4]=*res;
        return 1;
     } 
     staFunc->NormalizeVector(V);   
     staFunc->AddScaledVector(W,V,-staFunc->ScalProd(V,W),W);    /* ortogonalize */
     staFunc->NormalizeVector(W);
  } else {	
/* Neimark-Sacker */
    lambda=x[nphase+2+nappend];         /* cos(theta) = kappa */
    omega=sqrt(1.0-lambda*lambda);      /* sin(theta) */

/* compute real and imaginary part of the adjoint eigenvector */
    staFunc->CopyMatrixBlock(C1,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,VV) != 2*nphase-2) {
       *res=1.0;
       test[4]=*res;
       return 1;
    }
    staFunc->CopyVectorElements(VV,0,V,0,nphase);
    staFunc->CopyVectorElements(VV,nphase,W,0,nphase);

    /* normalize real and imaginary parts <V,V>=1, <W,W>=1, <V,W>=0 */
    d=staFunc->ScalProd(V,V);
    s=staFunc->ScalProd(W,W);
    r=staFunc->ScalProd(V,W);
    if (r!=0) {
      beta=sqrt(fabs(4.0*r*r+(s-d)*(s-d)));
      if (2.0*fabs(2.0*r) > beta) {
        phi=atan(fabs((s-d)/(2.0*r)));
      } else {
        phi=PID2-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++) {
        VV[i]=V[i]*cos(phi)-W[i]*sin(phi);
        VV[i+nphase]=V[i]*sin(phi)+W[i]*cos(phi);
      }
      staFunc->CopyVectorElements(VV,0,V,0,nphase);
      staFunc->CopyVectorElements(VV,nphase,W,0,nphase);
    }
    staFunc->NormalizeVector(V);
    staFunc->NormalizeVector(W);
  }

/* find two basis vectors with max projection to Span(V,W) */
  imaxV=0;
  maxV=fabs(V[0]);
  imaxW=0;
  maxW=fabs(W[0]);
  for (i=0; i<nphase; i++) {
    if (fabs(V[i])>maxV) {
      imaxV=i;
      maxV=fabs(V[i]);
    }
    if (fabs(W[i])>maxW) {
      imaxW=i;
      maxW=fabs(W[i]);
    }
  }
  if (imaxV==imaxW) {
    imaxW++;
    if (imaxW==nphase) imaxW=0;
  }
 
/* fill in the matrix B with vectors to be orthogonalized as rows */
  staFunc->ClearMatrix(B);
  staFunc->CopyVectorMatrixRow(V,B,0);
  staFunc->CopyVectorMatrixRow(W,B,1);
  for (i=0,ii=2,jj=0; i<nphase; i++,jj++) {
    if (i==imaxV || i==imaxW) continue;
    B[ii++][jj]=1;
  }
  
/* Gram-Schmidt orthogonalization of rows of B */
  for (i=2; i<nphase; i++) {
    staFunc->CopyMatrixRowVector(B,i,V);
    staFunc->CopyMatrixRowVector(B,i,F1);
    for (j=0; j<i; j++) {
      staFunc->CopyMatrixRowVector(B,j,W);
      s=staFunc->ScalProd(V,W);
      staFunc->AddScaledVector(F1,W,-s,F1);
    }
    staFunc->NormalizeVector(F1);
    staFunc->CopyVectorMatrixRow(F1,B,i);
  }
 
/* compute elements of the restriction matrix AH */ 
  staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
  staFunc->ClearMatrix(AH); 
  for (i=2; i<nphase; i++) {
    staFunc->CopyMatrixRowVector(B,i,V);
    staFunc->MultiplyMatrixVector(AA,V,W);
    for (j=2; j<nphase; j++) {
      staFunc->CopyMatrixRowVector(B,j,V);
      AH[j-2][i-2]=staFunc->ScalProd(V,W);
    }
  }
/* compute det of bi-product matrix BP=AoA-I */
  staFunc->ClearMatrix(DE);
{
  Int2 p,q,r,s;
  i=-1;
  for (p=1; p<nphase-2; p++) {
    for (q=0; q<p; q++) {
      i++;
      j=-1;
      for (r=1; r<nphase-2; r++) {
        for (s=0; s<r; s++) {
          j++;
          DE[i][j]=AH[p][r]*AH[q][s]-AH[q][r]*AH[p][s];
	  if (i==j) DE[i][j]-=1.0;
        }
      }
    }
  }
}
  if (staFunc->Decomp(DE)) {
    *res=0.0;
  } else {
    staFunc->GetMatrixData(DE,res,NULL);
  }
end:
  test[4]=*res;
  return 0;
}

Entry(Int2) ns2_ProcDoubleNS(Int2 Ind, Vector x, Vector v, CharPtr *msg) {
  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Double NS (?)"); 
  else 
{int i,j,imin,jmin;
 double r,s,d,phi,beta,gamma,lambda,lambda1,omega,theta;
 Vector Re,Im;
    Re=staFunc->CreateVector(nphase-2); 
    Im=staFunc->CreateVector(nphase-2); 
    staFunc->EigenValues(AH,Re,Im);

    /* find the pair of extra 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;
        }
      }
    }
    JacRHS(x,C1);
    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(C1,0,0,A,0,0,nphase,nphase);
      for (i=0; i<nphase; i++) A[i][i]-=lambda; 
      if (staFunc->SngSolve(A,_eps,V) != nphase-1) {
        sprintf(buf,"Double NS: lambda=%g, eigenvectors(?)",lambda);
        goto polufinal;
     } 
      staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
      for (i=0; i<nphase; i++) A[i][i]-=lambda1;
      if (staFunc->SngSolve(A,_eps,W) != nphase-1) {
        sprintf(buf,"Double NS: lambda1=%g, eigenvectors(?)",lambda1);
        goto polufinal;
      } 
      staFunc->NormalizeVector(V);   
      staFunc->NormalizeVector(W);
      staFunc->AddScaledVector(V,W,-1.0,F1);             
      staFunc->NormalizeVector(F1);
      staFunc->AddScaledVector(V,W,1.0,V);             
      staFunc->NormalizeVector(V);
/* provide bifurcation data for NN -> NS starter */
/*  Secondary branch: */
/*  BifVector[2*nphase+2] to [3*nphase+1]: 2-NS continuation vector*/
/*  BifVector[3*nphase+2] to [4*nphase+1]: 2-global vector L */
/*  BifVector[4*nphase+2]:                 2-kappa  */
      MemCopy(BifVector+2*nphase+2,V,sizeof(FloatHi)*nphase);
      MemCopy(BifVector+3*nphase+2,F1,sizeof(FloatHi)*nphase);
      BifVector[4*nphase+2]=(lambda+lambda1)/2.0;   /* kappa */
      stagenData->BifDataOutLen=sizeof(FloatHi)*(4*nphase+3);
      if (x[nphase+2+nappend] <= 0.0) {    
        sprintf(buf,"Double neutral saddle");
      } else {
        sprintf(buf,"Neimark-Sacker + neutral saddle");
      }
    } 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);
      
      /* compute real and imaginary part of the eigenvector */
      staFunc->ClearMatrix(BB);
      staFunc->CopyMatrixBlock(C1,0,0,BB,0,0,nphase,nphase);
      staFunc->CopyMatrixBlock(C1,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,VV) != 2*nphase-2) {
        sprintf(buf,"Double NS: sin_2=%g, cos_2=%g, eigenvectors(?)", omega,lambda);
        goto polufinal;
      }
      staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
      staFunc->CopyVectorElements(VV,nphase,Qi,0,nphase);

      /* normalize real and imaginary parts */
      d=staFunc->ScalProd(Qr,Qr);
      s=staFunc->ScalProd(Qi,Qi);
      r=staFunc->ScalProd(Qr,Qi);
      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++) {
          VV[i]=Qr[i]*cos(phi)-Qi[i]*sin(phi);
          VV[i+nphase]=Qr[i]*sin(phi)+Qi[i]*cos(phi);
        }
      }

      staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
      d=staFunc->ScalProd(Qr,Qr);
      gamma=1/sqrt(d);
      staFunc->MultiplyVectorScalar(VV,gamma,VV);
      staFunc->CopyVectorElements(VV,nphase,F1,0,nphase);
      staFunc->NormalizeVector(F1); 

/* provide bifurcation data for NN -> NS starter */
/*  Secondary branch: */
/*  BifVector[2*nphase+2] to [3*nphase+1]: 2-NS continuation vector*/
/*  BifVector[3*nphase+2] to [4*nphase+1]: 2-global vector L */
/*  BifVector[4*nphase+2]:                 2-kappa  */
        MemCopy(BifVector+2*nphase+2,VV,sizeof(FloatHi)*nphase);
        MemCopy(BifVector+3*nphase+2,F1,sizeof(FloatHi)*nphase);
        BifVector[4*nphase+2]=cos(theta);
        stagenData->BifDataOutLen=sizeof(FloatHi)*(4*nphase+3);

        if (x[nphase+2+nappend] <= 0.0) {    
          sprintf(buf,"Neutral saddle + Neimark-Sacker");
        } else {
          sprintf(buf,"Double Neimark-Sacker");
        }
    }
/* provide bifurcation data for NN -> NS starter */
/* Primary branch: */
/*  BifVector[0]:			 _itmap */
/*  BifVector[1] to [nphase]:            1-NS continuation vector */
/*  BifVector[nphase+1] to [2*nphase]:   1-global vector L */
/*  BifVector[2*nphase+1]:               1-kappa  */

    BifVector[0]=_itmap;
    MemCopy(BifVector+1,VB1,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+1+nphase,VB2,sizeof(FloatHi)*nphase); 
    BifVector[2*nphase+1]=x[nphase+2+nappend];

polufinal:
    Re=staFunc->DestroyVector(Re); 
    Im=staFunc->DestroyVector(Im); 
}  
  *msg=buf;
  return 0;
}

/* kappa+1 */
Entry(Int2) ns2_TestR2(Vector x, Vector v, FloatHiPtr res) {
  *res=1.0;
  if (nphase < 2) goto end;
  *res=x[nphase+2+nappend]+1.0;
end:
  test[5]=*res; 
  return 0;
}

Entry(Int2) ns2_ProcR2(Int2  Ind, Vector x, Vector v, CharPtr *msg) {
  stagenData->BifDataOutLen=0;
  if (Ind) 
    sprintf(buf,"Resonance 1:2 (?)");
  else {
Vector v1,W0,W1,F1;
Matrix D;
Int2 i,j,k;
double a20,b20,b11,f0,f1,f2,f3,f4,f5,d;
    v1=staFunc->CreateVector(2*nphase);
    W0=staFunc->CreateVector(nphase);
    W1=staFunc->CreateVector(nphase);
    D=staFunc->CreateMatrix(2*nphase,2*nphase);
    F1=staFunc->CreateVector(nphase);

/*   compute (generalized)eigenvectors   */
    if (JacRHS (x,C1)) return 1;
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) A[i][i]+=1.0;
    
    staFunc->CopyMatrixBlock(A,0,0,D,0,0,nphase,nphase);
    staFunc->CopyMatrixBlock(A,0,0,D,nphase,nphase,nphase,nphase);
    for (i=0; i<nphase; i++) D[i+nphase][i]=-1;
    i=staFunc->SngSolve(D,_eps,v1);
    if (i != 2*nphase-2) {
      sprintf(buf,"1:2 Resonance: eigenvectors (?)");
      goto final;
    }
    staFunc->CopyVectorElements(v1,0,V,0,nphase);       /* eigenvector V0 */
    d=sqrt(staFunc->ScalProd(V,V));
    staFunc->NormalizeVector(V);
    staFunc->CopyVectorElements(v1,nphase,W,0,nphase);  /* generalized eigenvector V1 */
    staFunc->MultiplyVectorScalar(W,1.0/d,W);
    staFunc->AddScaledVector(W,V,-staFunc->ScalProd(V,W),W);

/*  provide bifurcation data for R2->* starters */
/*  BifVector[0]:			 _itmap */
/*  BifVector[1] to [nphase]:            null-vector V */
/*  BifVector[nphase+1]:               kappa=-1  */
    BifVector[0]=(VAR)_itmap;
    MemCopy(BifVector+1,V,sizeof(FloatHi)*nphase);
    BifVector[nphase+1]=-1.0;
    stagenData->BifDataOutLen=sizeof(FloatHi)*(nphase+2);
    
/* compute the adjoint eigenvectors */
    staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
    for (i=0; i<nphase; i++) A[i][i]+=1.0;
    staFunc->TransposeMatrix(A,B);
    staFunc->ClearMatrix(D);
    staFunc->CopyMatrixBlock(B,0,0,D,0,0,nphase,nphase);
    staFunc->CopyMatrixBlock(B,0,0,D,nphase,nphase,nphase,nphase);
    for (i=0; i<nphase; i++) D[i+nphase][i]=-1;
    i=staFunc->SngSolve(D,_eps,v1);
    if (i != 2*nphase-2) {
      sprintf(buf,"1:2 Resonance: adjoint eigenvectors (?)");
      goto final;
    }
    staFunc->CopyVectorElements(v1,0,W1,0,nphase);       /* adjoint eigenvector W1 */
    staFunc->CopyVectorElements(v1,nphase,W0,0,nphase);  /* adjoint generalized eigenvector W0 */
    d=staFunc->ScalProd(V,W0);
    staFunc->MultiplyVectorScalar(W1,1.0/d,W1);
    staFunc->AddScaledVector(W0,W1,-staFunc->ScalProd(W0,W),W0);
    staFunc->MultiplyVectorScalar(W0,1.0/d,W0);

    sprintf(buf,"1:2 Resonance");

final:
    D=staFunc->DestroyMatrix(D);
    v1=staFunc->DestroyVector(v1);
    W0=staFunc->DestroyVector(W0);
    W1=staFunc->DestroyVector(W1);
    F1=staFunc->DestroyVector(F1);
  }
  *msg=buf;
  return 0;
}

/* kappa+0.5 */
Entry(Int2) ns2_TestR3(Vector x, Vector v, FloatHiPtr res) {
  *res=1.0;
  if (nphase < 2) goto end;
  *res=x[nphase+2+nappend]+0.5;
end:
  test[6]=*res; 
  return 0;
}

Entry(Int2) ns2_ProcR3(Int2  Ind, Vector x, Vector v, CharPtr *msg) {
  if (Ind) {
    stagenData->BifDataOutLen=0;
    sprintf(buf,"1:3 Resonance (?)"); 
  } else {
    BifVector[0]=(VAR)_itmap; 
    MemCopy(BifVector+1,VB1,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+nphase+1,VB2,sizeof(FloatHi)*nphase);
    BifVector[2*nphase+1]=x[nphase+2+nappend];
    stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+2);
    sprintf(buf,"1:3 Resonance"); 
  }
  *msg=buf;
  return 0; 
}

/* kappa */
Entry(Int2) ns2_TestR4(Vector x, Vector v, FloatHiPtr res) {
  *res=1.0;
  if (nphase < 2) goto end;
  *res=x[nphase+2+nappend];
end:
  test[7]=*res; 
  return 0;
}

Entry(Int2) ns2_ProcR4(Int2  Ind, Vector x, Vector v, CharPtr *msg) {
  if (Ind) {
    stagenData->BifDataOutLen=0;
    sprintf(buf,"1:4 Resonance (?)"); 
  } else {
    BifVector[0]=(VAR)_itmap; 
    MemCopy(BifVector+1,VB1,sizeof(FloatHi)*nphase);
    MemCopy(BifVector+nphase+1,VB2,sizeof(FloatHi)*nphase);
    BifVector[2*nphase+1]=x[nphase+2+nappend];
    stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+2);
    sprintf(buf,"1:4 Resonance"); 
  }
  *msg=buf;
  return 0; 
}

/****************************************************/
/*  Adaptation of the bordering vectors V1,V2,W1,W2 */
Entry(Int2) ns2_Adapter(Vector x, Vector v) {
int i,j,l,Ret;
VAR scal;
Matrix S;

  if (JacRHS(x,C1)) return 1;
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
/* update borders */
  staFunc->TransposeMatrix(BP,BT);
  Ret=staFunc->Decomp(BT);
  if (Ret) return 1;
  staFunc->ClearVector(PP1);
  PP1[nphase]=1.0;
  staFunc->Solve(BT,PP1);
  staFunc->CopyVector(PP1,W1);
  staFunc->ClearVector(PP2);
  PP2[nphase+1]=1.0;
  staFunc->Solve(BT,PP2);
  staFunc->CopyVector(PP2,W2);
  staFunc->CopyVector(Q1,V1);
  staFunc->CopyVector(Q2,V2); 

  staFunc->NormalizeVector(V1);
  staFunc->NormalizeVector(W1);
  staFunc->NormalizeVector(V2);
  staFunc->NormalizeVector(W2);

/* update indices */
  if (stagenData->Der2) {   
    for (i=0;i<nphase;i++) {
      BP[nphase][i]=V1[i];
      BP[i][nphase]=W1[i];
    }
    for (i=0;i<nphase+2;i++) {
      BP[nphase+1][i]=V2[i];
      BP[i][nphase+1]=W2[i];
    }
    H=SymHess(x);
    staFunc->ClearVector(HV1);
    staFunc->ClearVector(HV2);
    staFunc->ClearVector(HW1);
    staFunc->ClearVector(HW2);     
    staFunc->CopyVectorElements(VB1,0,VB,0,nphase);
    staFunc->MultiplyMatrixVector(A,VB,HV1); 
    staFunc->CopyVectorElements(VB2,0,VB,0,nphase);
    staFunc->MultiplyMatrixVector(A,VB,HV2);
    staFunc->TransposeMatrix(A,B);
    staFunc->CopyVectorElements(PP1,0,VB,0,nphase);
    staFunc->MultiplyMatrixVector(B,VB,HW1);
    staFunc->CopyVectorElements(PP2,0,VB,0,nphase);
    staFunc->MultiplyMatrixVector(B,VB,HW2);  
    S=staFunc->CreateMatrix(2,2); 
    for (l=0; l<nphase+2+nappend; l++) {
      for (i=0,S[0][0]=0.0,S[0][1]=0.0,S[1][0]=0.0,S[1][1]=0.0;i<nphase;i++) {
        for (j=0; j<nphase; j++) {
    	  S[0][0]-=HW1[i]*VB1[j]*HessianElem(i,j,l)+PP1[i]*HV1[j]*HessianElem(i,j,l);
          S[0][0]-=2.0*x[nphase+2+nappend]*PP1[i]*VB1[j]*HessianElem(i,j,l);
    	  S[0][1]-=HW1[i]*VB2[j]*HessianElem(i,j,l)+PP1[i]*HV2[j]*HessianElem(i,j,l);
          S[0][1]-=2.0*x[nphase+2+nappend]*PP1[i]*VB2[j]*HessianElem(i,j,l);
    	  S[1][0]-=HW2[i]*VB1[j]*HessianElem(i,j,l)+PP2[i]*HV1[j]*HessianElem(i,j,l);
          S[1][0]-=2.0*x[nphase+2+nappend]*PP2[i]*VB1[j]*HessianElem(i,j,l);
    	  S[1][1]-=HW2[i]*VB2[j]*HessianElem(i,j,l)+PP2[i]*HV2[j]*HessianElem(i,j,l);            
          S[1][1]-=2.0*x[nphase+2+nappend]*PP2[i]*VB2[j]*HessianElem(i,j,l);
        }
      }
      G0[l]=S[0][0];
      G1[l]=S[0][1];
      G2[l]=S[1][0];
      G3[l]=S[1][1];
    }
    S=staFunc->DestroyMatrix(S);
    G0[nphase+2+nappend]=2.0*staFunc->ScalProd(PP1,HV1);
    G1[nphase+2+nappend]=2.0*staFunc->ScalProd(PP1,HV2);
    G2[nphase+2+nappend]=2.0*staFunc->ScalProd(PP2,HV1);
    G3[nphase+2+nappend]=2.0*staFunc->ScalProd(PP2,HV2);        
  } else {
      staFunc->ClearVector(HV1);
      staFunc->ClearVector(HV2);
      staFunc->ClearVector(HW1);
      staFunc->ClearVector(HW2);
      staFunc->CopyVectorElements(VB1,0,HV1,0,nphase);
      staFunc->CopyVectorElements(VB2,0,HV2,0,nphase);
      staFunc->MultiplyMatrixVector(A,HV1,HV1);
      staFunc->MultiplyMatrixVector(A,HV2,HV2);
      staFunc->TransposeMatrix(A,B);
      staFunc->CopyVectorElements(PP1,0,HW1,0,nphase);
      staFunc->CopyVectorElements(PP2,0,HW2,0,nphase);
      staFunc->MultiplyMatrixVector(B,HW1,HW1);
      staFunc->MultiplyMatrixVector(B,HW2,HW2);  
      Jac1=staFunc->CreateMatrix(nphase,nphase);
      Jac2=staFunc->CreateMatrix(nphase,nphase);
      for (l=0; l<nphase+2+nappend; l++) {
        staFunc->ClearVector(xp);
        staFunc->CopyVectorElements(x0,0,xp,0,nphase+2);
        xp[l]+=_dx;
        staFunc->ClearVector(xp1);
        staFunc->CopyVectorElements(x0,0,xp1,0,nphase+2);
        xp1[l]-=_dx;
        JacRHS(xp,C1);
        staFunc->CopyMatrixBlock(C1,0,0,Jac1,0,0,nphase,nphase);
        JacRHS(xp1,C1);
        staFunc->CopyMatrixBlock(C1,0,0,Jac2,0,0,nphase,nphase);
     /* g11 */
        staFunc->CopyVectorElements(VB1,0,K3,0,nphase);
        staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
        staFunc->MultiplyMatrixVector(Jac2,K3,K2); 
        staFunc->AddScaledVector(K1,K2,-1.0,K3);
        staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
        G0[l]=-staFunc->ScalProd(HW1,K3);
        G0[l]-=2.0*x[nphase+2+nappend]*staFunc->ScalProd(PP1,K3);
        staFunc->TransposeMatrix(Jac1,B);
        staFunc->CopyVectorElements(PP1,0,K3,0,nphase);
        staFunc->MultiplyMatrixVector(B,K3,K1);   
        staFunc->TransposeMatrix(Jac2,B);
        staFunc->MultiplyMatrixVector(B,K3,K2);  
        staFunc->AddScaledVector(K1,K2,-1.0,K3);
        staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
        G0[l]-=staFunc->ScalProd(K3,HV1);                                   
     /* g12 */
        staFunc->CopyVectorElements(VB2,0,K3,0,nphase);
        staFunc->MultiplyMatrixVector(Jac1,K3,K1);
        staFunc->MultiplyMatrixVector(Jac2,K3,K2); 
        staFunc->AddScaledVector(K1,K2,-1.0,K3);
        staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
        G1[l]=-staFunc->ScalProd(HW1,K3);
        G1[l]-=2.0*x[nphase+2+nappend]*staFunc->ScalProd(PP1,K3);
        staFunc->TransposeMatrix(Jac1,B);
        staFunc->CopyVectorElements(PP1,0,K3,0,nphase);
        staFunc->MultiplyMatrixVector(B,K3,K1);   
        staFunc->TransposeMatrix(Jac2,B);
        staFunc->MultiplyMatrixVector(B,K3,K2);
        staFunc->AddScaledVector(K1,K2,-1.0,K3);
        staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
        G1[l]-=staFunc->ScalProd(K3,HV2);                                   
     /* g21 */
        staFunc->CopyVectorElements(VB1,0,K3,0,nphase);
        staFunc->MultiplyMatrixVector(Jac1,K3,K1);  
        staFunc->MultiplyMatrixVector(Jac2,K3,K2); 
        staFunc->AddScaledVector(K1,K2,-1.0,K3);
        staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
        G2[l]=-staFunc->ScalProd(HW2,K3);
        G2[l]-=2.0*x[nphase+2+nappend]*staFunc->ScalProd(PP2,K3);
        staFunc->TransposeMatrix(Jac1,B);
        staFunc->CopyVectorElements(PP2,0,K3,0,nphase);
        staFunc->MultiplyMatrixVector(B,K3,K1);     
        staFunc->TransposeMatrix(Jac2,B);
        staFunc->MultiplyMatrixVector(B,K3,K2);  
        staFunc->AddScaledVector(K1,K2,-1.0,K3);
        staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
        G2[l]-=staFunc->ScalProd(K3,HV1);                                   
     /* g22 */
        staFunc->CopyVectorElements(VB2,0,K3,0,nphase);
        staFunc->MultiplyMatrixVector(Jac1,K3,K1); 
        staFunc->MultiplyMatrixVector(Jac2,K3,K2);  
        staFunc->AddScaledVector(K1,K2,-1.0,K3);
        staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
        G3[l]=-staFunc->ScalProd(HW2,K3);
        G3[l]-=2.0*x[nphase+2+nappend]*staFunc->ScalProd(PP2,K3);
        staFunc->TransposeMatrix(Jac1,B);
        staFunc->CopyVectorElements(PP2,0,K3,0,nphase);
        staFunc->MultiplyMatrixVector(B,K3,K1);   
        staFunc->TransposeMatrix(Jac2,B);
        staFunc->MultiplyMatrixVector(B,K3,K2); 
        staFunc->AddScaledVector(K1,K2,-1.0,K3);
        staFunc->MultiplyVectorScalar(K3,1.0/(_dx+_dx),K3);
        G3[l]-=staFunc->ScalProd(K3,HV2);                                             
        }
        Jac1=staFunc->DestroyMatrix(Jac1);
        Jac2=staFunc->DestroyMatrix(Jac2);
        G0[nphase+2+nappend]=2.0*staFunc->ScalProd(PP1,HV1);
        G1[nphase+2+nappend]=2.0*staFunc->ScalProd(PP1,HV2);
        G2[nphase+2+nappend]=2.0*staFunc->ScalProd(PP2,HV1);
        G3[nphase+2+nappend]=2.0*staFunc->ScalProd(PP2,HV2);           
    } 
  for (i=0; i<nphase+3+nappend; i++) {
    GM1[i][0]=G0[i];
    GM1[i][1]=G1[i];
    GM1[i][2]=G2[i];
    GM1[i][3]=G3[i];
  }
/* reselection of index pairs (m1,n1), (m2,n2) */
  staFunc->ClearMatrix(JM2);
  staFunc->CopyMatrixBlock(C1,0,0,JM2,0,0,nphase,nphase+2);    
  for (i=0; i<nphase; i++) JM2[i][i]-=1.0;
  for (i=0; i<nphase+3+nappend; i++) { 
     JM2[nphase][i]=GM1[i][2*m1+n1];
     JM2[nphase+1][i]=GM1[i][2*m2+n2];
     JM2[nphase+2][i]=v[i];
  }  
  staFunc->ClearVector(TV2);
  TV2[nphase]=1.0;
  staFunc->Decomp(JM2);
  staFunc->Solve(JM2,TV2);
  staFunc->ClearMatrix(TVM);
  for (i=0;i<nphase+2;i++) TVM[i][0] = TV2[i];
  staFunc->ClearVector(TV2);
  TV2[nphase+1]=1.0;  
  staFunc->Solve(JM2,TV2);
  for (i=0;i<nphase+2;i++) TVM[i][1] = TV2[i];    
  TVM[nphase+2][2]=1.0;

  staFunc->TransposeMatrix(TVM,TVMT);
  staFunc->MultiplyMatrixMatrix(TVMT,TVM,HM); 
  Ret=staFunc->Decomp(HM);
  if (Ret) return 1;
  staFunc->MultiplyMatrixVector(TVMT,G0,HV);  
  staFunc->Solve(HM,HV);
  staFunc->MultiplyMatrixVector(TVM,HV,G0);  
  staFunc->MultiplyMatrixVector(TVMT,G1,HV);  
  staFunc->Solve(HM,HV);
  staFunc->MultiplyMatrixVector(TVM,HV,G1);      
  staFunc->MultiplyMatrixVector(TVMT,G2,HV);  
  staFunc->Solve(HM,HV);
  staFunc->MultiplyMatrixVector(TVM,HV,G2);
  staFunc->MultiplyMatrixVector(TVMT,G3,HV);  
  staFunc->Solve(HM,HV);
  staFunc->MultiplyMatrixVector(TVM,HV,G3);         

/* Determine component with maximum norm => gives m1,n1 */
  G[0][0] = staFunc->NormOfVector(G0);
  G[0][1] = staFunc->NormOfVector(G1);
  G[1][0] = staFunc->NormOfVector(G2);
  G[1][1] = staFunc->NormOfVector(G3);    
  m1=0;
  n1=0;      
  if (G[0][1] > G[m1][n1]) n1=1;  
  if (G[1][0] > G[m1][n1]) { m1=1; n1=0; }
  if (G[1][1] > G[m1][n1]) { m1=1; n1=1; }
  for (i=0;i<nphase+3;i++) {
    GM1[i][0] = G0[i]; 
    GM1[i][1] = G1[i];
    GM1[i][2] = G2[i];
    GM1[i][3] = G3[i];
  }
/* G0 has biggest norm ; necessary for projection */ 
  for (i=0;i<nphase+3;i++) {
    G0[i] = GM1[i][2*m1+n1];
  }      
  staFunc->NormalizeVector(G0);
/* G1,G2,G3 are projected on G0 */ 
  for (j=0;j<4;j++) {
    if (j != (2*m1+n1)) {
     for (i=0;i<nphase+3;i++) { /* neem een vector, niet de grootste */
        G1[i] = GM1[i][j];
     }
     scal = staFunc->ScalProd(G1,G0);
     staFunc->AddScaledVector(G1,G0,-scal,G1);
     for (i=0;i<nphase+3;i++) { /* zet nieuwe vector terug */
       GM1[i][j] = G1[i];
     }
    }
  }      
/* compute biggest norm between G1,G2 and G3  => gives m2,n2 */
  for (i=0;i<nphase+3;i++) {
    G0[i] = GM1[i][0];
    G1[i] = GM1[i][1];
    G2[i] = GM1[i][2];
    G3[i] = GM1[i][3];
  }
  G[0][0] = staFunc->NormOfVector(G0);
  G[0][1] = staFunc->NormOfVector(G1);
  G[1][0] = staFunc->NormOfVector(G2);
  G[1][1] = staFunc->NormOfVector(G3);    
  G[m1][n1]=-1.0;
  m2=0;
  n2=0;
  if (G[0][1] > G[m2][n2]) n2=1;  
  if (G[1][0] > G[m2][n2]) { m2=1; n2=0; }
  if (G[1][1] > G[m2][n2]) { m2=1; n2=1; }  

  return 0;
}