/*
   __gh1.c
   Generalized Hopf curve by means of maximally extended system
   December 8, 1998
*/

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

#include "stagen.h"

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

#include "__gh1.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_EIGV		4
#define CL_TEST	        5
#define CL_UTEST	6

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


Local(StagenDataPtr) stagenData=NULL;	/* global copy of Starter's parameter passed by content */
Local(__gh_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) E=NULL;		/* current set of eigenvalues */
Local(VAR) zero=0;

#include "_symdiff.c"

/* Working space for defining and test functions  */

Local(Matrix) A,B,C,C1,C2,AA,AT,BB;
Local(Vector) vn,xp,xp1,V,W,F,Vp,L,VV,VB,Q,Qr,Qi,PP1,PP2,Pr,Pi,Q0r,Q0i,aa,bb,cc,dd,rr;
Local(Vector) V1,V2,V3,OLD,RL,K1,K2,K3,OK1,OK2;
Local(Vector) Re=NULL,Im=NULL;
Local(VAR) test[1],dx,ddx,dddx;

#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=(__gh_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];
  /* 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+3));    
  if (staUtil->CheckParameters(staData->Active,npar,active,3+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;
  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 */
    E=MemNew(sizeof(FloatHi)*2*nphase);
    Re=staFunc->CreateVector(nphase);
    Im=staFunc->CreateVector(nphase);
    /* Allocate memory for working space used in defining and test functions */
    A=staFunc->CreateMatrix(nphase,nphase);
    B=staFunc->CreateMatrix(nphase,nphase);
    
    V1=staFunc->CreateVector(nphase);
    V2=staFunc->CreateVector(nphase);
    V3=staFunc->CreateVector(nphase);
    VB=staFunc->CreateVector(nphase);
    VV=staFunc->CreateVector(2*nphase);
    Q=staFunc->CreateVector(m);
    PP1=staFunc->CreateVector(m+2);
    PP2=staFunc->CreateVector(m+2);
    OLD=staFunc->CreateVector(m+1);
    RL=staFunc->CreateVector(nphase);
    K1=staFunc->CreateVector(m+2);
    K2=staFunc->CreateVector(m+2);
    K3=staFunc->CreateVector(m+2);
    OK1=staFunc->CreateVector(m+2);
    OK2=staFunc->CreateVector(m+2);
  
    AA=staFunc->CreateMatrix(nphase,nphase);
    BB=staFunc->CreateMatrix(2*nphase,2*nphase);
    C=staFunc->CreateMatrix(nphase+nappend,8*nphase+6+nappend);
    C1=staFunc->CreateMatrix(nphase,8*nphase+6+nappend);
    C2=staFunc->CreateMatrix(7*nphase+5,8*nphase+6+nappend);
    AT=staFunc->CreateMatrix(nphase,nphase);
    Qr=staFunc->CreateVector(nphase);
    Qi=staFunc->CreateVector(nphase);
    Q0r=staFunc->CreateVector(nphase);
    Q0i=staFunc->CreateVector(nphase);    
    Pr=staFunc->CreateVector(nphase);
    Pi=staFunc->CreateVector(nphase);
    aa=staFunc->CreateVector(nphase);
    bb=staFunc->CreateVector(nphase);
    cc=staFunc->CreateVector(nphase);
    dd=staFunc->CreateVector(nphase);
    rr=staFunc->CreateVector(nphase);
    x0=staFunc->CreateVector(8*nphase+6+nappend);
    x1=staFunc->CreateVector(8*nphase+6+nappend);
    v0=staFunc->CreateVector(8*nphase+6+nappend);
    vn=staFunc->CreateVector(nphase);
    xp=staFunc->CreateVector(nphase+3+nappend);
    xp1=staFunc->CreateVector(8*nphase+6+nappend);
    V=staFunc->CreateVector(nphase);
    W=staFunc->CreateVector(nphase);
    F=staFunc->CreateVector(nphase);
    L=staFunc->CreateVector(nphase);
    Vp=staFunc->CreateVector(nphase+3+nappend);
    BifVector=MemNew(sizeof(FloatHi)*(4*nphase+2));
  }
  Der2num=ind2_(nphase+npar-1,nphase+npar-1)+1;
  Der3num=ind3_(nphase+npar-1,nphase+npar-1,nphase+npar-1)+1;
  Der4num=ind4_(nphase+npar-1,nphase+npar-1,nphase+npar-1,nphase+npar-1)+1;
  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)) {
/* no continuation data */
    stagenData->ContDataStaLen=2*nphase*sizeof(FloatHi);
    stagenData->ContDataStaPtr=MemNew(stagenData->ContDataStaLen);
    MemCopy(stagenData->ContDataStaPtr,Q0r,(nphase)*sizeof(FloatHi));
    MemCopy(stagenData->ContDataStaPtr+(nphase),Q0i,(nphase)*sizeof(FloatHi));
    stagenData->BifDataOutLen=0;
    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);
    E=MemFree(E);
    Re=staFunc->DestroyVector(Re);
    Im=staFunc->DestroyVector(Im);
    A=staFunc->DestroyMatrix(A);
    B=staFunc->DestroyMatrix(B);
    V1=staFunc->DestroyVector(V1);
    V2=staFunc->DestroyVector(V2);
    V3=staFunc->DestroyVector(V3);
    Q=staFunc->DestroyVector(Q);
    PP1=staFunc->DestroyVector(PP1);
    PP2=staFunc->DestroyVector(PP2);
    AA=staFunc->DestroyMatrix(AA);
    BB=staFunc->DestroyMatrix(BB);
    VV=staFunc->DestroyVector(VV);
    VB=staFunc->DestroyVector(VB);
    OLD=staFunc->DestroyVector(OLD);
    RL=staFunc->DestroyVector(RL);    
    K1=staFunc->DestroyVector(K1);
    K2=staFunc->DestroyVector(K2);
    K3=staFunc->DestroyVector(K3);
    OK1=staFunc->DestroyVector(OK1);
    OK2=staFunc->DestroyVector(OK2); 
    C=staFunc->DestroyMatrix(C);
    C1=staFunc->DestroyMatrix(C1);
    C2=staFunc->DestroyMatrix(C2);
    AT=staFunc->DestroyMatrix(AT);
    Qr=staFunc->DestroyVector(Qr);
    Qi=staFunc->DestroyVector(Qi);
    Q0r=staFunc->DestroyVector(Q0r);
    Q0i=staFunc->DestroyVector(Q0i);    
    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);
    vn=staFunc->DestroyVector(vn);
    xp=staFunc->DestroyVector(xp);
    xp1=staFunc->DestroyVector(xp1);
    V=staFunc->DestroyVector(V);
    W=staFunc->DestroyVector(W);
    F=staFunc->DestroyVector(F);
    Vp=staFunc->DestroyVector(Vp);
    L=staFunc->DestroyVector(L);
    BifVector=MemFree(BifVector);  
  }
} 

Local(Int2) JacDefEquilib (Vector x, Matrix J);
Entry(Int2) gh1_JacDefGeneralizedHopf (Vector x, Matrix J);
Entry(Int2) gh1_DefGeneralizedHopf (Vector x, Vector J);
Entry(Int2) gh1_DefBord (Vector x, Vector y);
Entry(Int2) gh1_DefRHS (Vector x, Vector y);

/************/
/* Starters */
/* gh -> gh Starter */
Entry(Int2) gh1_Starter(StagenDataPtr stagenptr) {
  Int2 i,j,k,Ret;
  Vector v1,v2,v3;
  Matrix D,DE;
  VAR beta, phi, omega, d, h, r, s;
 
  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 (stagenptr->ContDataGenLen) {       /* "Extension" */
   /* no continuation data */ 
    MemCopy(Q0r,stagenptr->ContDataStaPtr,sizeof(FloatHi)*(nphase));
    MemCopy(Q0i,stagenptr->ContDataStaPtr+(nphase),sizeof(FloatHi)*(nphase));
    MemFree(stagenptr->ContDataStaPtr);
    stagenptr->ContDataStaLen=0;
    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<3+nappend; i++) x0[nphase+i]=staData->P[active[i]];            
    JacRHS(x0,C1);    
  } else {
    v1=staFunc->CreateVector(8*nphase+6+nappend);
    v2=staFunc->CreateVector(8*nphase+6+nappend);
    D=staFunc->CreateMatrix(8*nphase+5+nappend,8*nphase+6+nappend);
    DE=staFunc->CreateMatrix(8*nphase+6+nappend,8*nphase+6+nappend);

    MemCopy(x0,staData->X,sizeof(FloatHi)*nphase);
    for (i=0; i<3+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] to [nphase-1]:        Hopf continuation vector	*/
/*  BifVector[nphase] to [2*nphase-1]: global vector L 		*/
/*  BifVector[2*nphase]:	       -lambda*lambda   	*/

/*      MemCopy(V,stagenData->BifDataInPtr,sizeof(FloatHi)*nphase); */
/*      MemCopy(L,stagenData->BifDataInPtr+nphase,sizeof(FloatHi)*nphase); */
 
/*    } else { */
      Ret=0;
{int i,j,imin,jmin;
 double r,s,d,phi,beta,gamma,lambda;
    staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
    staFunc->EigenValues(AA,Re,Im);
    /* find the pair of eigenvalues with zero sum */
    s=fabs(Re[0]+Re[1])+fabs(Im[0]+Im[1]);
    imin=0;
    jmin=1;
    for (i=0; i<nphase; i++) {
      for (j=i+1; j<nphase; j++) {
        d=fabs(Re[i]+Re[j])+fabs(Im[i]+Im[j]);
        if (d<s) {
          imin=i;
          jmin=j;
          s=d;
        }
      }
    }
    if (fabs(Im[imin])+fabs(Im[jmin])==0) {
/* Neutral saddle */
      lambda=fabs(Re[imin]);
      /* compute real eigenvectors associated with lambda and -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,V) != 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]+=lambda;
      if (staFunc->SngSolve(AA,_eps,W) != nphase-1) {
        staUtil->ErrMessage("Unable to find (-)eigenvector at the initial point");
        Ret=1;
        goto polufinal;
      } 
      staFunc->NormalizeVector(V);   
      staFunc->NormalizeVector(W);
      staFunc->AddScaledVector(V,W,-1.0,L);             
      staFunc->NormalizeVector(L);
      staFunc->AddScaledVector(V,W,1.0,V);             
      staFunc->NormalizeVector(V);
      x0[8*nphase+3]=-lambda*lambda;
    } else {
/* Hopf */
      lambda=Re[imin];
      omega=fabs(Im[imin]);
      /* 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,V,0,nphase);
      staFunc->CopyVectorElements(VV,nphase,W,0,nphase);
      /* normalize real and imaginary parts */
        d=staFunc->ScalProd(V,V);
        s=staFunc->ScalProd(W,W);
        r=staFunc->ScalProd(V,W);
        if (r!=0) {
          beta=0.5*(d-s)/r;
          if (beta == 0.0) {
            for (i=0; i<nphase; i++) {
              VV[i]=V[i]-W[i];
              VV[i+nphase]=V[i]+W[i];
            }
          } else {
            phi=0.5*atan(-1.0/beta);
            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);
        d=staFunc->ScalProd(V,V);
        gamma=1/sqrt(d);
        staFunc->MultiplyVectorScalar(VV,gamma,VV);

        staFunc->CopyVectorElements(VV,0,V,0,nphase);
        staFunc->CopyVectorElements(VV,nphase,L,0,nphase);
        staFunc->NormalizeVector(L);
        x0[8*nphase+3]=omega;
    }
polufinal:
      if (Ret) {
        v1=staFunc->DestroyVector(v1);
        v2=staFunc->DestroyVector(v2);
        D=staFunc->DestroyMatrix(D);
        DE=staFunc->DestroyMatrix(DE);
        Term(FALSE);
        return 1;
      }
  /*  } else bij BifData */
}       
/* Compute Qr, Qi */
    staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
    staFunc->ClearMatrix(BB);
    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][i+nphase]=omega;
       BB[i+nphase][i]=-omega;
    }
    staFunc->ClearVector(VV);
    if (staFunc->SngSolve(BB,_eps*100,VV) < 2*nphase-2) {
       return 1;
    }
    staFunc->CopyVectorElements(VV,0,Qr,0,nphase);
    staFunc->CopyVectorElements(VV,nphase,Qi,0,nphase);

/* Compute Pr, Pi */
    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+nphase][i]=omega;
    }
    staFunc->ClearVector(VV);
    if (staFunc->SngSolve(BB,_eps*100,VV) < 2*nphase-2) {
       return 1;
    }
    staFunc->CopyVectorElements(VV,0,Pr,0,nphase);
    staFunc->CopyVectorElements(VV,nphase,Pi,0,nphase);

/* Normalization conditions */
    
    d=staFunc->ScalProd(Qr,Qr);
    s=staFunc->ScalProd(Qi,Qi);
    staFunc->MultiplyVectorScalar(Qr,1.0/(sqrt(d+s)),Qr);
    staFunc->MultiplyVectorScalar(Qi,1.0/(sqrt(d+s)),Qi);
    staFunc->CopyVector(Qr,Q0r);
    staFunc->CopyVector(Qi,Q0i);    
    staFunc->CopyVectorElements(Qr,0,x0,nphase+3,nphase);
    staFunc->CopyVectorElements(Qi,0,x0,2*nphase+3,nphase);

    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));

    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(Pr,beta,Pr);
    staFunc->MultiplyVectorScalar(Pi,beta,Pi);
    staFunc->CopyVectorElements(Pr,0,x0,3*nphase+3,nphase);
    staFunc->CopyVectorElements(Pi,0,x0,4*nphase+3,nphase);
    
/* Compute V1,V2,V3 */
    if (stagenData->Der2) {
      CopyToX(x0);
      CopyToP(x0);
      H=stagenData->Der2(X,P,&zero);
      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]+Qr[k]*Qi[j]);
          }
        }
        aa[i]=s;
        bb[i]=r;
        cc[i]=-h;
      }        
    } else { /* numerical part */
        v3=staFunc->CreateVector(8*nphase+6+nappend);
	staFunc->ClearVector(v3);
	staFunc->CopyVectorElements(Qr,0,v3,0,nphase);
	gh1_DefRHS(x0,aa);
        staFunc->AddScaledVector(x0,v3,ddx,x1);
	gh1_DefRHS(x1,bb);
	staFunc->AddScaledVector(x0,v3,-ddx,x1);
	gh1_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,v3,0,nphase);
	gh1_DefRHS(x0,rr);
        staFunc->AddScaledVector(x0,v3,ddx,x1);
	gh1_DefRHS(x1,bb);
	staFunc->AddScaledVector(x0,v3,-ddx,x1);
	gh1_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,v3,0,nphase);
	staFunc->AddScaledVector(x0,v3,0.5*ddx,x1);
	gh1_DefRHS(x1,rr);
	staFunc->AddScaledVector(x0,v3,-0.5*ddx,x1);
	gh1_DefRHS(x1,cc);
	staFunc->AddScaledVector(rr,cc,1.0,dd);
	
        staFunc->AddScaledVector(Qr,Qi,-1.0,rr);
	staFunc->CopyVectorElements(rr,0,v3,0,nphase);
	staFunc->AddScaledVector(x0,v3,0.5*ddx,x1);
	gh1_DefRHS(x1,rr);
	staFunc->AddScaledVector(x0,v3,-0.5*ddx,x1);
	gh1_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);
        v3=staFunc->DestroyVector(v3);        
    }
    staFunc->ClearVector(V1);
    if (JacRHS(x0,C1)) return 1;
    staFunc->CopyMatrixBlock(C1,0,0,AA,0,0,nphase,nphase);
    staFunc->AddScaledVector(aa,bb,1.0,V1);
    Ret=staFunc->Decomp(AA);
    if (Ret) { /* Zero-Hopf case */
      return 1;      
    }
    staFunc->Solve(AA,V1);
    staFunc->CopyVectorElements(V1,0,x0,5*nphase+3,nphase);
    
    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]=2.0*omega;
     BB[i+nphase][i]=-2.0*omega;
    }
    Ret=staFunc->Decomp(BB);
    if (Ret) return 3;
    staFunc->Solve(BB,VV);
    staFunc->CopyVectorElements(VV,0,x0,6*nphase+3,nphase);
    staFunc->CopyVectorElements(VV,0,V2,0,nphase);
    staFunc->CopyVectorElements(VV,nphase,x0,7*nphase+3,nphase);
    staFunc->CopyVectorElements(VV,nphase,V3,0,nphase);  

/* set initial values for omega1 and omega2 */
    x0[8*nphase+4]=0.0;
    x0[8*nphase+5]=omega;

/* Compute a tangent vector to the Generalized Hopf curve */
    gh1_DefGeneralizedHopf(x0,v0);
    Ret=gh1_JacDefGeneralizedHopf(x0,D);  

    for (i=0; i<8*nphase+6+nappend; i++) {
      v1[i]=1.0;
      staFunc->AppendMatrixVector(D,v1,DE);
      Ret=staFunc->Decomp(DE);
      if (Ret == 0) {
        v2[8*nphase+5+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 Generalized Hopf 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<3+nappend; i++) P[active[i]]=vp[nphase+i];
}

Local(void) CopyToE(Vector Re, Vector Im) {
  MemCopy(E,Re,sizeof(FloatHi)*nphase);
  MemCopy(E+nphase,Im,sizeof(FloatHi)*nphase);
}


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

Entry(Int2) gh1_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_EIGV]=E;
      CopyToE(Re,Im);
      indirect[CL_TEST]=test;
  }
  return 0;
}
/***************************/
/* User-defined functions  */
Entry(Int2) gh1_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) gh1_DefRHS (Vector x, Vector y) {
  CopyToX(x);
  CopyToP(x);
  stagenData->Rhs(X,P,&zero,y);	
  return 0;
}

/*************************************/
/* Defining function for equilibrium */
Local(Int2) DefEquilib (Vector x, Vector y) {
  gh1_DefRHS(x,y);
  if (nappend) gh1_DefUser(x,y+nphase);
  return 0;
}

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

/***************************************************/
/* Jacobian of defining conditions for equilibria  */
 Local(Int2) JacDefEquilib (Vector x, Matrix Jac) {
  staFunc->ClearMatrix(Jac);
  JacRHS(x, Jac);
  if (nappend) staFunc->Der(gh1_DefUser, nappend, x, dx, Jac+nphase);
  return 0;
}

/***************************************/
/* Generalized Hopf defining condition */
Entry(Int2) gh1_DefBord (Vector x, Vector y) {
Int2 i,j,k,l;
VAR d,h,r,s,omega,s2;

  if (JacRHS(x,C1)) return 1;	  /* A,C1 - global */
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  staFunc->CopyVectorElements(x,nphase+3,Qr,0,nphase);
  staFunc->CopyVectorElements(x,2*nphase+3,Qi,0,nphase);  
  staFunc->CopyVectorElements(x,3*nphase+3,Pr,0,nphase);  
  staFunc->CopyVectorElements(x,4*nphase+3,Pi,0,nphase);  
  staFunc->CopyVectorElements(x,5*nphase+3,V1,0,nphase);  
  staFunc->CopyVectorElements(x,6*nphase+3,V2,0,nphase);  
  staFunc->CopyVectorElements(x,7*nphase+3,V3,0,nphase);  
  omega=x[8*nphase+3];

/* A Qr + omega0 Qi = 0 ; Qr and Qi are assumed global */ 
  staFunc->MultiplyMatrixVector(A,Qr,VB);
  staFunc->AddScaledVector(VB,Qi,omega,y);

/* A Qi - omega0 Qr = 0 */  
  staFunc->MultiplyMatrixVector(A,Qi,VB);
  staFunc->AddScaledVector(VB,Qr,-x[8*nphase+3],y+nphase);

/* AT Pr - omega1 Pr - omega2 Pi = 0 */
  staFunc->TransposeMatrix(A,B);
  staFunc->ClearVector(RL);
  staFunc->MultiplyMatrixVector(B,Pr,VB);
  staFunc->AddScaledVector(VB,Pr,-x[8*nphase+4],RL);
  staFunc->AddScaledVector(RL,Pi,-x[8*nphase+5],y+2*nphase);
/*  staFunc->CopyVectorElements(VB,0,y,2*nphase,nphase); */

/* -AT Pi + omega1 Pi - omega2 Pr = 0 */
  staFunc->ClearVector(RL);
  staFunc->MultiplyMatrixVector(B,Pi,VB);
  staFunc->MultiplyVectorScalar(VB,-1.0,VB);
  staFunc->AddScaledVector(VB,Pi,x[8*nphase+4],RL);
  staFunc->AddScaledVector(RL,Pr,-x[8*nphase+5],y+3*nphase);
/*  staFunc->CopyVectorElements(VB,0,y,3*nphase,nphase); */

/* <Q0r,Qr> + <Q0i,Qi> - 1 = 0 */  
  r=staFunc->ScalProd(Q0r,Qr);
  s=staFunc->ScalProd(Q0i,Qi);
  y[4*nphase]=r+s-1;

/* -<Q0i,Qr> + <Q0r,Qi> = 0 */
  r=staFunc->ScalProd(Q0i,Qr);
  s=staFunc->ScalProd(Q0r,Qi);
  y[4*nphase+1]=-r+s;

/* <Pr,Qr> + <Pi,Qi> - 1 = 0 ; <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);
  y[4*nphase+2]=d+s-1;
  y[4*nphase+3]=r-h;

/* computation of Fu V1 - B(q,Conjugate(q)) = 0 */   
  if (stagenData->Der2) {
    CopyToX(x);
    CopyToP(x);
    H=stagenData->Der2(X,P,&zero);
    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]+Qi[j]*Qi[k]);
         r+=HessianElem(i,j,k)*(Qr[j]*Qr[k]-Qi[j]*Qi[k]);
         h+=HessianElem(i,j,k)*(Qr[j]*Qi[k]+Qr[k]*Qi[j]);
        }
      }
      aa[i]=s; /* aa contains B(q,conjuguate(q)) */
      bb[i]=r; /* bb contains real part of B(q,q) */
      cc[i]=h; /* cc contains imaginary part of B(q,q) */
    }
  } else { /* numerical part */
      staFunc->ClearVector(v0);
      staFunc->CopyVectorElements(Qr,0,v0,0,nphase);
      gh1_DefRHS(x,aa);
      staFunc->AddScaledVector(x,v0,ddx,x1);
      gh1_DefRHS(x1,bb);
      staFunc->AddScaledVector(x,v0,-ddx,x1);
      gh1_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);
      gh1_DefRHS(x,rr);
      staFunc->AddScaledVector(x,v0,ddx,x1);
      gh1_DefRHS(x1,bb);
      staFunc->AddScaledVector(x,v0,-ddx,x1);
      gh1_DefRHS(x1,cc);
      staFunc->AddScaledVector(bb,rr,-2.0,rr);
      staFunc->AddScaledVector(rr,cc,1.0,dd);
      staFunc->MultiplyVectorScalar(dd,1.0/ddx/ddx,dd);

      staFunc->AddScaledVector(aa,dd,1.0,aa);
      staFunc->AddScaledVector(aa,dd,-2.0,bb);

      staFunc->AddScaledVector(Qr,Qi,1.0,rr);
      staFunc->CopyVectorElements(rr,0,v0,0,nphase);
      staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
      gh1_DefRHS(x1,rr);
      staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
      gh1_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);
      gh1_DefRHS(x1,rr);
      staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
      gh1_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);      
  }  

  staFunc->MultiplyMatrixVector(A,V1,RL);
  staFunc->AddScaledVector(RL,aa,-1.0,y+4*nphase+4);

/* computation of (2i omega0 I - Fu)(V2 + i V3) - B(q,q) = 0 */
  staFunc->MultiplyMatrixVector(A,V2,RL);
  staFunc->MultiplyVectorScalar(RL,-1.0,RL);
  staFunc->AddScaledVector(RL,V3,-2.0*x[8*nphase+3],VB);
  staFunc->AddScaledVector(VB,bb,-1.0,y+5*nphase+4);
  staFunc->MultiplyMatrixVector(A,V3,RL);
  staFunc->MultiplyVectorScalar(RL,-1.0,RL);
  staFunc->AddScaledVector(RL,V2,2.0*x[8*nphase+3],VB);
  staFunc->AddScaledVector(VB,cc,-1.0,y+6*nphase+4);

/* computation of first Lyapunov-coefficient ; use V1,V2,V3 */
  if (stagenData->Der3) {
    CopyToX(x);
    CopyToP(x);
    D3=stagenData->Der3(X,P,&zero);	/* The 1st elem must be time! */
    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]));
	  }
	}
      }
    }	
  } else { /* numerical part */
       staFunc->ClearVector(v0);
       staFunc->CopyVectorElements(Qr,0,v0,0,nphase);
       staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r=(2.0/3.0)*staFunc->ScalProd(Pr,dd);
       staFunc->AddScaledVector(x,v0,dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r-=2.0*staFunc->ScalProd(Pr,dd);
       staFunc->AddScaledVector(x,v0,-dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r+=2.0*staFunc->ScalProd(Pr,dd);
       staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r-=(2.0/3.0)*staFunc->ScalProd(Pr,dd);
	
       staFunc->CopyVectorElements(Qi,0,v0,0,nphase);
       staFunc->AddScaledVector(x,v0,3*dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r+=(2.0/3.0)*staFunc->ScalProd(Pi,dd);
       staFunc->AddScaledVector(x,v0,dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r-=2.0*staFunc->ScalProd(Pi,dd);
       staFunc->AddScaledVector(x,v0,-dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r+=2.0*staFunc->ScalProd(Pi,dd);
       staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r-=(2.0/3.0)*staFunc->ScalProd(Pi,dd);
	
       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);
       gh1_DefRHS(x1,dd);
       r+=staFunc->ScalProd(bb,dd)/6;
       staFunc->AddScaledVector(x,v0,dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r-=staFunc->ScalProd(bb,dd)/2;
       staFunc->AddScaledVector(x,v0,-dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r+=staFunc->ScalProd(bb,dd)/2;
       staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r-=staFunc->ScalProd(bb,dd)/6;

       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);
       gh1_DefRHS(x1,dd);
       r+=staFunc->ScalProd(bb,dd)/6;
       staFunc->AddScaledVector(x,v0,dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r-=staFunc->ScalProd(bb,dd)/2;
       staFunc->AddScaledVector(x,v0,-dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r+=staFunc->ScalProd(bb,dd)/2;
       staFunc->AddScaledVector(x,v0,-3*dddx/2,x1);
       gh1_DefRHS(x1,dd);
       r-=staFunc->ScalProd(bb,dd)/6;
       r/=(dddx*dddx*dddx);                     
  } 

  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]*(-2.0*Qr[j]*V1[k]+Qr[j]*V2[k]+Qi[j]*V3[k]) +
            Pi[i]*(-2.0*Qi[j]*V1[k]+Qr[j]*V3[k]-Qi[j]*V2[k]));
     }
   }
  } else {
      staFunc->ClearVector(v0);
/*  -2*<Pr,B(Qr,V1)>  */
      staFunc->AddScaledVector(Qr,V1,1.0,dd);
      staFunc->CopyVectorElements(dd,0,v0,0,nphase);
      staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s=-2.0*staFunc->ScalProd(Pr,dd);
      staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s-=2.0*staFunc->ScalProd(Pr,dd);

      staFunc->AddScaledVector(Qr,V1,-1.0,dd);
      staFunc->CopyVectorElements(dd,0,v0,0,nphase);
      staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s+=2.0*staFunc->ScalProd(Pr,dd);
      staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s+=2.0*staFunc->ScalProd(Pr,dd);
/*  -2*<Pi,B(Qi,V1)>  */
      staFunc->AddScaledVector(Qi,V1,1.0,dd);
      staFunc->CopyVectorElements(dd,0,v0,0,nphase);
      staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s-=2.0*staFunc->ScalProd(Pi,dd);
      staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s-=2.0*staFunc->ScalProd(Pi,dd);

      staFunc->AddScaledVector(Qi,V1,-1.0,dd);
      staFunc->CopyVectorElements(dd,0,v0,0,nphase);
      staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s+=2.0*staFunc->ScalProd(Pi,dd);
      staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s+=2.0*staFunc->ScalProd(Pi,dd);
/*  <Pr,B(Qr,V2)>  */
      staFunc->AddScaledVector(Qr,V2,1.0,dd);
      staFunc->CopyVectorElements(dd,0,v0,0,nphase);
      staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s+=staFunc->ScalProd(Pr,dd);
      staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s+=staFunc->ScalProd(Pr,dd);

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

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

      staFunc->AddScaledVector(Qr,V3,-1.0,dd);
      staFunc->CopyVectorElements(dd,0,v0,0,nphase);
      staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s-=staFunc->ScalProd(Pi,dd);
      staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s-=staFunc->ScalProd(Pi,dd);
/*  <Pi,B(Qi,V2)>  */
      staFunc->AddScaledVector(Qi,V2,1.0,dd);
      staFunc->CopyVectorElements(dd,0,v0,0,nphase);
      staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s-=staFunc->ScalProd(Pi,dd);
      staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s-=staFunc->ScalProd(Pi,dd);

      staFunc->AddScaledVector(Qi,V2,-1.0,dd);
      staFunc->CopyVectorElements(dd,0,v0,0,nphase);
      staFunc->AddScaledVector(x,v0,0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s+=staFunc->ScalProd(Pi,dd);
      staFunc->AddScaledVector(x,v0,-0.5*ddx,x1);
      gh1_DefRHS(x1,dd);
      s+=staFunc->ScalProd(Pi,dd);
      s/=(ddx*ddx);      
  }
  y[7*nphase+4]=r+s;
  return 0;
}

/****************************************************/
/* Jacobian of Generalized Hopf defining conditions */
Local(Int2) JacDefBord (Vector x, Matrix Jac) {
int i,j,k,l,m;
VAR s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,s15,s16,s17,s18,s19,s20,
    s21,s22,s23,s24,s25,s26,s27,s28,s29,s30,s31,s32,s33,s34,s35,s36,s77;

  JacRHS(x,C1);
  staFunc->CopyMatrixBlock(C1,0,0,A,0,0,nphase,nphase);
  staFunc->TransposeMatrix(A,B);
  staFunc->ClearMatrix(Jac);
  if (stagenData->Der2 && stagenData->Der3 && stagenData->Der4) {
    CopyToX(x);
    CopyToP(x);
    H=stagenData->Der2(X,P,&zero);
    D3=stagenData->Der3(X,P,&zero);
    D4=stagenData->Der4(X,P,&zero);
    for (i=0; i<nphase; i++) {
      for (k=0; k<nphase+3; k++) {
        for (j=0;j<nphase; j++) {
          Jac[i][k]+=HessianElem(i,j,k)*Qr[j];
          Jac[i+nphase][k]+=HessianElem(i,j,k)*Qi[j];
          Jac[i+2*nphase][k]+=HessianElem(j,i,k)*Pr[j];
          Jac[i+3*nphase][k]-=HessianElem(j,i,k)*Pi[j];
          Jac[i+4*nphase+4][k]+=HessianElem(i,j,k)*x[j+5*nphase+3];
          Jac[i+5*nphase+4][k]-=HessianElem(i,j,k)*x[j+6*nphase+3];          
          Jac[i+6*nphase+4][k]-=HessianElem(i,j,k)*x[j+7*nphase+3];
        }  
      }
    }
    for (i=0; i<nphase; i++) {
      for (k=0; k<nphase; k++) {
        for (j=0;j<nphase; j++) {
          Jac[i+4*nphase+4][k+nphase+3]-=2.0*HessianElem(i,j,k)*Qr[j];
          Jac[i+4*nphase+4][k+2*nphase+3]-=2.0*HessianElem(i,j,k)*Qi[j];
          Jac[i+5*nphase+4][k+nphase+3]-=2.0*HessianElem(i,j,k)*Qr[j];
          Jac[i+5*nphase+4][k+2*nphase+3]+=2.0*HessianElem(i,j,k)*Qi[j];
          Jac[i+6*nphase+4][k+nphase+3]-=2.0*HessianElem(i,j,k)*Qi[j];
          Jac[i+6*nphase+4][k+2*nphase+3]-=2.0*HessianElem(i,j,k)*Qr[j];          
        }  
      }
    }  
    staFunc->CopyMatrixBlock(C1,0,0,Jac,0,nphase+3,nphase,nphase);
    staFunc->CopyMatrixBlock(C1,0,0,Jac,nphase,2*nphase+3,nphase,nphase);
    staFunc->CopyMatrixBlock(C1,0,0,Jac,4*nphase+4,5*nphase+3,nphase,nphase);
    for (i=0; i<nphase; i++) {
      Jac[i][i+2*nphase+3]=x[8*nphase+3];
      Jac[i][8*nphase+3]=Qi[i];
      Jac[i+nphase][i+nphase+3]=-x[8*nphase+3];
      Jac[i+nphase][8*nphase+3]=-Qr[i];
      Jac[i+2*nphase][i+4*nphase+3]=-x[8*nphase+5];
      Jac[i+2*nphase][8*nphase+4]=-Pr[i];
      Jac[i+2*nphase][8*nphase+5]=-Pi[i];
      Jac[i+3*nphase][i+3*nphase+3]=-x[8*nphase+5];
      Jac[i+3*nphase][8*nphase+4]=Pi[i];
      Jac[i+3*nphase][8*nphase+5]=-Pr[i];
      Jac[4*nphase][i+nphase+3]=Q0r[i];
      Jac[4*nphase][i+2*nphase+3]=Q0i[i];
      Jac[4*nphase+1][i+nphase+3]=-Q0i[i];
      Jac[4*nphase+1][i+2*nphase+3]=Q0r[i];      
      Jac[4*nphase+2][i+nphase+3]=Pr[i];
      Jac[4*nphase+2][i+2*nphase+3]=Pi[i];
      Jac[4*nphase+2][i+3*nphase+3]=Qr[i];
      Jac[4*nphase+2][i+4*nphase+3]=Qi[i];
      Jac[4*nphase+3][i+nphase+3]=-Pi[i];
      Jac[4*nphase+3][i+2*nphase+3]=Pr[i];
      Jac[4*nphase+3][i+3*nphase+3]=Qi[i];
      Jac[4*nphase+3][i+4*nphase+3]=-Qr[i];      
      Jac[i+5*nphase+4][i+7*nphase+3]=-2.0*x[8*nphase+3];
      Jac[i+5*nphase+4][8*nphase+3]=-2.0*x[i+7*nphase+3];
      Jac[i+6*nphase+4][i+6*nphase+3]=2.0*x[8*nphase+3];
      Jac[i+6*nphase+4][8*nphase+3]=2.0*x[i+6*nphase+3];
      s1=s3=s5=s7=s9=s11=s13=s14=s15=s16=s17=s19=0.0;
      s21=s22=s23=s77=0.0;
      for (j=0; j<nphase; j++) {
        Jac[i+2*nphase][j+3*nphase+3]=B[i][j];
        Jac[i+3*nphase][j+4*nphase+3]=-B[i][j];
        Jac[i+5*nphase+4][j+6*nphase+3]=-A[i][j];
        Jac[i+6*nphase+4][j+7*nphase+3]=-A[i][j];
        s2=s4=s6=s8=s10=s12=s18=s20=0.0;
        for (k=0; k<nphase; k++) {
          for (l=0; l<nphase; l++) {
            s2+=D3Elem(j,k,l,i)*Qr[k]*Qr[l];
            s4+=D3Elem(j,k,l,i)*Qi[k]*Qi[l];
            s6+=D3Elem(j,k,l,i)*Qr[k]*Qi[l];
            Jac[7*nphase+4][i+3*nphase+3]+=D3Elem(i,j,k,l)*Qr[j]*Qr[k]*Qr[l]+D3Elem(i,j,k,l)*Qr[j]*Qi[k]*Qi[l];
            Jac[7*nphase+4][i+4*nphase+3]+=D3Elem(i,j,k,l)*Qi[j]*Qi[k]*Qi[l]+D3Elem(i,j,k,l)*Qr[j]*Qr[k]*Qi[l];
          }
          s8+=HessianElem(j,k,i)*V1[k];
          s10+=HessianElem(j,k,i)*V2[k];
          s12+=HessianElem(j,k,i)*V3[k];          
          Jac[7*nphase+4][i+3*nphase+3]-=2.0*HessianElem(i,j,k)*Qr[j]*V1[k]
          -HessianElem(i,j,k)*Qr[j]*V2[k]-HessianElem(i,j,k)*Qi[j]*V3[k];
          Jac[7*nphase+4][i+4*nphase+3]-=2.0*HessianElem(i,j,k)*Qi[j]*V1[k]
          +HessianElem(i,j,k)*Qi[j]*V2[k]-HessianElem(i,j,k)*Qr[j]*V3[k];
          s18+=HessianElem(j,k,i)*Qr[k];
          s20+=HessianElem(j,k,i)*Qi[k];
        }
        s1+=Pr[j]*s2;
        s3+=Pr[j]*s4;
        s5+=Pi[j]*s6; /* changed from s5+=Pr[j]*s6 */
        s7+=Pr[j]*s8;
        s9+=Pr[j]*s10;
        s11+=Pi[j]*s12;
        s13+=Pr[j]*s6;
        s14+=Pi[j]*s4;
        s15+=Pi[j]*s2;
        s16+=Pr[j]*s12;
        s17+=Pi[j]*s10;
        s19+=Pr[j]*s18;
        s21+=Pi[j]*s20;
        s22+=Pr[j]*s20;
        s23+=Pi[j]*s18;
        s77+=Pi[j]*s8;
      }
      Jac[i+2*nphase][i+3*nphase+3]-=x[8*nphase+4];
      Jac[i+3*nphase][i+4*nphase+3]+=x[8*nphase+4];      
      Jac[7*nphase+4][i+nphase+3]=3.0*s1+s3+2.0*s5-2.0*s7+s9+s11;
      Jac[7*nphase+4][i+2*nphase+3]=2.0*s13+3.0*s14+s15-2.0*s77+s16-s17;
      Jac[7*nphase+4][i+5*nphase+3]=-2.0*s19-2.0*s21;
      Jac[7*nphase+4][i+6*nphase+3]=s19-s21;
      Jac[7*nphase+4][i+7*nphase+3]=s22+s23;
    }
    for (i=0; i<nphase; i++) {
      s25=s30=s33=0.0;
      for (j=0; j<nphase; j++) {
        s24=s26=s27=s28=s29=s31=s32=s34=s35=s36=0.0;
        for (k=0; k<nphase; k++) {
          for (l=0; l<nphase; l++) {
            for (m=0; m<nphase; m++) {
             s24+=D4Elem(j,k,l,m,i)*Qr[k]*Qr[l]*Qr[m];
             s26+=D4Elem(j,k,l,m,i)*Qr[k]*Qi[l]*Qi[m];
             s27+=D4Elem(j,k,l,m,i)*Qi[k]*Qi[l]*Qi[m];
             s28+=D4Elem(j,k,l,m,i)*Qr[k]*Qr[l]*Qi[m];             
            }          
            s29+=D3Elem(j,k,l,i)*Qr[k]*x[l+5*nphase+3];
            s31+=D3Elem(j,k,l,i)*Qi[k]*x[l+5*nphase+3];
            s32+=D3Elem(j,k,l,i)*Qr[k]*x[l+6*nphase+3];
            s34+=D3Elem(j,k,l,i)*Qi[k]*x[l+7*nphase+3];
            s35+=D3Elem(j,k,l,i)*Qr[k]*x[l+7*nphase+3];
            s36+=D3Elem(j,k,l,i)*Qi[k]*x[l+6*nphase+3];
          }
        }
        s25+=Pr[j]*s24+Pr[j]*s26+Pi[j]*s27+Pi[j]*s28;
        s30+=Pr[j]*s29+Pi[j]*s31;
        s33+=Pr[j]*(s32+s34)+Pi[j]*(s35-s36);        
      }
      Jac[7*nphase+4][i]=s25-2.0*s30+s33;
    }
    for (i=0; i<3; i++) {
      s25=s30=s33=0.0;
      for (j=0; j<nphase; j++) {
        s24=s26=s27=s28=s29=s31=s32=s34=s35=s36=0.0;
        for (k=0; k<nphase; k++) {
          for (l=0; l<nphase; l++) {
            for (m=0; m<nphase; m++) {
             s24+=D4Elem(j,k,l,m,nphase+active[i])*Qr[k]*Qr[l]*Qr[m];
             s26+=D4Elem(j,k,l,m,nphase+active[i])*Qr[k]*Qi[l]*Qi[m];
             s27+=D4Elem(j,k,l,m,nphase+active[i])*Qi[k]*Qi[l]*Qi[m];
             s28+=D4Elem(j,k,l,m,nphase+active[i])*Qr[k]*Qr[l]*Qi[m];             
            }          
            s29+=D3Elem(j,k,l,nphase+active[i])*Qr[k]*x[l+5*nphase+3];
            s31+=D3Elem(j,k,l,nphase+active[i])*Qi[k]*x[l+5*nphase+3];
            s32+=D3Elem(j,k,l,nphase+active[i])*Qr[k]*x[l+6*nphase+3];
            s34+=D3Elem(j,k,l,nphase+active[i])*Qi[k]*x[l+7*nphase+3];
            s35+=D3Elem(j,k,l,nphase+active[i])*Qr[k]*x[l+7*nphase+3];
            s36+=D3Elem(j,k,l,nphase+active[i])*Qi[k]*x[l+6*nphase+3];
          }
        }
        s25+=Pr[j]*s24+Pr[j]*s26+Pi[j]*s27+Pi[j]*s28;
        s30+=Pr[j]*s29+Pi[j]*s31;
        s33+=Pr[j]*(s32+s34)+Pi[j]*(s35-s36);        
      }
      Jac[7*nphase+4][i+nphase]=s25-2.0*s30+s33;
    }    
    for (i=0; i<nphase; i++) {
      for (l=0; l<nphase; l++) {
        for (j=0; j<nphase; j++) {
          for (k=0; k<nphase; k++) {
            Jac[i+4*nphase+4][l]-=D3Elem(i,j,k,l)*Qr[j]*Qr[k]
            +D3Elem(i,j,k,l)*Qi[j]*Qi[k];
            Jac[i+5*nphase+4][l]-=D3Elem(i,j,k,l)*Qr[j]*Qr[k]
            -D3Elem(i,j,k,l)*Qi[j]*Qi[k];
            Jac[i+6*nphase+4][l]-=2.0*D3Elem(i,j,k,l)*Qr[j]*Qi[k];
          }
        }
      }
      for (l=0; l<3; l++) {
        for (j=0; j<nphase; j++) {
          for (k=0; k<nphase; k++) {
            Jac[i+4*nphase+4][l+nphase]-=D3Elem(i,j,k,nphase+active[l])*Qr[j]*Qr[k]
            +D3Elem(i,j,k,nphase+active[l])*Qi[j]*Qi[k];
            Jac[i+5*nphase+4][l+nphase]-=D3Elem(i,j,k,nphase+active[l])*Qr[j]*Qr[k]
            -D3Elem(i,j,k,nphase+active[l])*Qi[j]*Qi[k];
            Jac[i+6*nphase+4][l+nphase]-=2.0*D3Elem(i,j,k,nphase+active[l])*Qr[j]*Qi[k];
          }
        }
      }

    }
  } else {
     if (staFunc->Der(gh1_DefBord, 7*nphase+5, x, dx, Jac)) return 1;
    }  
  return 0;
}

/*************************************************/
/* Defining functions for Generalized Hopf curve */

Entry(Int2) gh1_DefGeneralizedHopf (Vector x, Vector y) {
  if (DefEquilib(x,y)) return 1;
  if (gh1_DefBord(x,y+nphase+nappend)) return 2;
  return 0;
}

/*************************************************/
/* Jacobian of defining functions for Generalized Hopf curve */  
Entry(Int2) gh1_JacDefGeneralizedHopf (Vector x, Matrix Jac) {
  staFunc->ClearMatrix(Jac);
  if (JacDefEquilib(x,C)) return 1;
  staFunc->CopyMatrixBlock(C,0,0,Jac,0,0,nphase+nappend,nphase+3+nappend);
  if (JacDefBord(x,C2)) return 2;
  staFunc->CopyMatrixBlock(C2,0,0,Jac,nphase+nappend,0,7*nphase+5,8*nphase+6+nappend);
  return 0;
}

/*********************/
/* Default processor */

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

  /* Assert: v1==NULL, msg==NULL, Ind==0 */
  gh1_DefGeneralizedHopf (x, xp1);  /* compute global Q for test functions */
  MemCopy(x1,x,sizeof(FloatHi)*(8*nphase+6+nappend));     /* for the last point */
  if (!staData->eigcomp) return 0;
  if (JacRHS(x,C1)) return (1); 

  staFunc->EigenValues(C1,Re,Im);
  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));

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

  *res=1.0;
  test[0]=*res;  
  return 0;
}

Entry(Int2) gh1_ProcDummy(Int2 Ind, Vector x, Vector v, CharPtr *msg) {
 
  return 0;
}

Entry(Int2) gh1_Adapter(Vector x, Vector v) {
VAR r,s;

 staFunc->CopyVector(Qr,Q0r);
 staFunc->CopyVector(Qi,Q0i);
 r=staFunc->ScalProd(Q0r,Q0r);
 s=staFunc->ScalProd(Q0i,Q0i);
 staFunc->MultiplyVectorScalar(Q0r,1.0/(sqrt(r+s)),Q0r);
 staFunc->MultiplyVectorScalar(Q0i,1.0/(sqrt(r+s)),Q0i); 

 return 0;
}