/*
   1D Reaction-Diffusion
   rd__ss
   Any distribution-->StadyStateCurve
*/

#define SKO_EPS		1.e-10

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

#include "stagen.h"

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

#include "rd__ss.h"
#include "_conti.h"	/* continuator's (i.e. generator's) data */

#define sqr(u) ((u)*(u))	/* u**2 */
#define Print(v) staFunc->PrintVector(v,#v"=") 

/* Numbers of visible name classes;
   the order must be the same as in /names section */
#define CL_TIME		1
#define CL_SPACE	2
#define CL_COMPONENTS	3
#define CL_PAR		4
#define CL_EIGV		5
#define CL_TEST		6
#define CL_UTEST	7
#define CL_STF		8


enum {		/* useed as subscripts to stagenData->AuxFunc[] */
      RHS_I=-1,
      RHSD1_I,
      RHSD2_I,
      RHSD3_I,
      RHSD4_I,
      RHSD5_I,
      BC0_I,
      BC0D1_I,
      BC0D2_I,
      BC0D3_I,
      BC0D4_I,
      BC0D5_I,
      BC1_I,
      BC1D1_I,
      BC1D2_I,
      BC1D3_I,
      BC1D4_I,
      BC1D5_I
};

typedef EntryPtr(int,BCPtr,(VARPTR u, VARPTR ux, VARPTR p, VARPTR t, VARPTR f));
typedef EntryPtr(VARPTR ,BCDPtr,(VARPTR u, VARPTR ux, VARPTR p, VARPTR t));


Local(StagenDataPtr) stagenData=NULL;	/* global copy of Starter's parameter passed by content */
Local(rd__ss_DataPtr) staData=NULL;	/* global copy of our own data area */
Local(StdLinAlgDataPtr) staFunc=NULL;	/* global copy of pointer to standard linear algebra functions */
Local(RdLinAlgDataPtr) rdFunc=NULL;	/* global copy of pointer to specific for RD linear algebra functions */
Local(UtilitiesDataPtr) staUtil=NULL;	/* global copy of pointer to utility functions */
Local(Int2Ptr) active;	        	/* index of active parameter(s) */
Local(Int2) nappend;			/* number of appended user functions */ 
Local(Int2) npar,nphase,nspace;		/* dims of par and phase spaces */
Local(Int2) nmesh;
Local(Vector) x1=NULL;			/* ptr to vector - the first point produced by starter */
Local(Vector) x0=NULL;			/* approximate first point */
Local(Vector) v0=NULL;			/* approximate tangent vector at first point */
Local(FloatHiPtr) P=NULL;		/* current set of params used by Decoder */
Local(FloatHiPtr) E=NULL;		/* current set of eigenvalues */
Local(FloatHiPtr) J;			/* current set of Jacobian vals */
Local(FloatHiPtr) H;                   	/* current set of Hessian vals */
Local(Int2) Der2num;			/* number of 2nd order derivatives */
Local(VAR) zero=0.0;			/* dummy for time */
Local(Boolean) NumDer1=TRUE;		/* TRUE-compute derivatives numericaly */

/********** DEBUG **********/
#define ROWS	30
Local(void) PrintFunc(Vector v) {
  CharPtr m=MemNew(ROWS*(nmesh+1)+1);
  VAR vmax,vmin;
  Int2 i;
  MemFill(m,' ',ROWS*(nmesh+1));
  for (i=0; i<ROWS; i++) m[i*(nmesh+1)+nmesh]='\n';
  vmax=vmin=v[0];
  for (i=1; i<nmesh; i++) {
    if (vmax<v[i]) vmax=v[i];
    if (vmin>v[i]) vmin=v[i];
  }
  for (i=0; i<nmesh; i++) {
    m[((int)((vmax-v[i])*(ROWS-1)/(vmax-vmin)))*(nmesh+1)+i]='*';
  }
  printf(m);
  m=MemFree(m);
}
#undef ROWS
/********** DEBUG **********/

/* Working space for test functions  */

Local(rdMatrix) A0=NULL,A=NULL,B=NULL,C=NULL;
Local(Matrix) JacW,JacV,AppMat;		/* NOT rdMatrix! */
Local(Vector) Re=NULL,Im=NULL;		/* Eigenvalues */
Local(VAR) L2Norm;
Local(Vector) STF=NULL;
Local(VAR) test[3];	/* last computed value for each test function */
Local(Int2) i;  

/* workspace ARS description */
Local(Vector) X,XN,U_0,U_X,U_XX;
Local(Vector) FF,BB,CC,DD; 
Local(VAR) d3fmin=1.0e-5;
/* end workspace */            

/* work vectors for Newtonstep */
Vector alf, bet, gam, FP, xn, gg, dgg, aa, bb, z;

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

/* Compute interpolation coefficients. */
/* staData->Jinterpol: 0 - linear interpolation, 1 - cubic spline */
Local(Int2) Spline (Vector x, Vector f, Vector b, Vector c, Vector d) {
Int2 i,n,n1,n2;
VAR t;
  n=staFunc->GetVectorDim(x);
  n1=n-1;
  n2=n-2;
  switch (staData->Jinterpol) {
    case 1:		/* cubic spline */
      d[0]=x[1]-x[0];
      c[1]=(f[1]-f[0])/d[0];
      for (i=1; i<n1; i++) {	
	d[i]=x[i+1]-x[i];
	b[i]=2.0*(x[i+1]-x[i-1]);
	c[i+1]=(f[i+1]-f[i])/d[i];
	c[i]=c[i+1]-c[i];
      }
      b[0]=-d[0];
      b[n1]=-d[n2];
      c[0]=c[2]/(x[3]-x[1])-c[1]/(x[2]-x[0]);
      c[n1]=c[n2]/(x[n1]-x[n1-2])-c[n1-2]/(x[n2]-x[n1-3]);
      c[0]=c[0]*d[0]*(d[0]/(x[3]-x[0]));  
      c[n1]=-c[n1]*d[n2]*(d[n2]/(x[n1]-x[n1-3]));
      for (i=1; i<n; i++) {
	t=d[i-1]/b[i-1];
	b[i]=b[i]-t*d[i-1];
	c[i]=c[i]-t*c[i-1];
      }
      c[n1]=c[n1]/b[n1];
      for (i=n2; i>-1; i--)
	c[i]=(c[i]-d[i]*c[i+1])/b[i];
      b[n1]=(f[n1]-f[n2])/d[n2]+d[n2]*(c[n2]+2.0*c[n1]);
      for (i=0; i<n1; i++) {
	b[i]=(f[i+1]-f[i])/d[i]-d[i]*(c[i+1]+2.0*c[i]);
	d[i]=(c[i+1]-c[i])/d[i];
	c[i]=3.0*c[i];
      }
      c[n1]=3.0*c[n1];
      d[n1]=d[n2];
      break;
    case 0:		/* linear interpolation */
      for (i=0; i<n1; i++) b[i]=(f[i+1]-f[i])/(x[i+1]-x[i]);
      break;
    }
  return 0;
}

/* Raznostnaya f''' . START */
Local(Int2) F3num (Vector x, Vector f, Vector dddf) {
Int2 i,n,n1,n2,n3,n4,n5;
VAR t,z1,z2,z3,z4,z5;

  n=staFunc->GetVectorDim(x);
  n1=n-1;
  n2=n-2;
  n3=n-3;
  n4=n-4;
  n5=n-5;  
  z1=-(x[1]+x[2]+x[3]+x[4]-4.0*x[0])/(x[1]-x[0])/(x[2]-x[0])/(x[3]-x[0])/(x[4]-x[0]);
  z2= (x[2]+x[3]+x[4]-3.0*x[0])/(x[1]-x[0])/(x[2]-x[1])/(x[3]-x[1])/(x[4]-x[1]);
  z3=-(x[1]+x[3]+x[4]-3.0*x[0])/(x[2]-x[1])/(x[3]-x[2])/(x[2]-x[0])/(x[4]-x[2]);
  z4= (x[1]+x[2]+x[4]-3.0*x[0])/(x[3]-x[2])/(x[4]-x[3])/(x[3]-x[1])/(x[3]-x[0]);
  z5=-(x[1]+x[2]+x[3]-3.0*x[0])/(x[4]-x[3])/(x[4]-x[2])/(x[4]-x[1])/(x[4]-x[0]);
  dddf[0]=6.0*(z1*f[0]+z2*f[1]+z3*f[2]+z4*f[3]+z5*f[4]);
  
    z1=-(x[2]+x[3]+x[4]-3.0*x[1])/(x[1]-x[0])/(x[2]-x[0])/(x[3]-x[0])/(x[4]-x[0]);
    z2= (x[0]+x[2]+x[3]+x[4]-4.0*x[1])/(x[1]-x[0])/(x[2]-x[1])/(x[3]-x[1])/(x[4]-x[1]);
    z3=-(x[0]+x[3]+x[4]-3.0*x[1])/(x[2]-x[1])/(x[3]-x[2])/(x[2]-x[0])/(x[4]-x[2]);
    z4= (x[0]+x[2]+x[4]-3.0*x[1])/(x[3]-x[2])/(x[4]-x[3])/(x[3]-x[1])/(x[3]-x[0]);
    z5=-(x[0]+x[2]+x[3]-3.0*x[1])/(x[4]-x[3])/(x[4]-x[2])/(x[4]-x[1])/(x[4]-x[0]);
    dddf[1]=6.0*(z1*f[0]+z2*f[1]+z3*f[2]+z4*f[3]+z5*f[4]);
  
      for (i=2; i<n2; i++) {
      z1=-(x[i-1]+x[i+1]+x[i+2]-3.0*x[i])/(x[i-1]-x[i-2])/(x[i]-x[i-2])/(x[i+1]-x[i-2])/(x[i+2]-x[i-2]);
      z2= (x[i-2]+x[i+1]+x[i+2]-3.0*x[i])/(x[i-1]-x[i-2])/(x[i]-x[i-1])/(x[i+1]-x[i-1])/(x[i+2]-x[i-1]);
      z3=-(x[i-2]+x[i-1]+x[i+1]+x[i+2]-4.0*x[i])/(x[i]-x[i-1])/(x[i+1]-x[i])/(x[i]-x[i-2])/(x[i+2]-x[i]);
      z4= (x[i-2]+x[i-1]+x[i+2]-3.0*x[i])/(x[i+1]-x[i])/(x[i+2]-x[i+1])/(x[i+1]-x[i-2])/(x[i+1]-x[i-1]);
      z5=-(x[i-2]+x[i-1]+x[i+1]-3.0*x[i])/(x[i+2]-x[i+1])/(x[i+2]-x[i])/(x[i+2]-x[i-1])/(x[i+2]-x[i-2]);
      dddf[i]=6.0*(z1*f[i-2]+z2*f[i-1]+z3*f[i]+z4*f[i+1]+z5*f[i+2]);
      }
  
    z1=-(x[n3]+x[n4]+x[n5]-3.0*x[n2])/(x[n2]-x[n1])/(x[n3]-x[n1])/(x[n4]-x[n1])/(x[n5]-x[n1]);
    z2= (x[n1]+x[n3]+x[n4]+x[n5]-4.0*x[n2])/(x[n2]-x[n1])/(x[n3]-x[n2])/(x[n4]-x[n2])/(x[n5]-x[n2]);
    z3=-(x[n1]+x[n4]+x[n5]-3.0*x[n2])/(x[n3]-x[n2])/(x[n4]-x[n3])/(x[n3]-x[n1])/(x[n5]-x[n3]);
    z4= (x[n1]+x[n3]+x[n5]-3.0*x[n2])/(x[n4]-x[n3])/(x[n5]-x[n4])/(x[n4]-x[n2])/(x[n4]-x[n1]);
    z5=-(x[n1]+x[n3]+x[n4]-3.0*x[n2])/(x[n5]-x[n4])/(x[n5]-x[n3])/(x[n5]-x[n2])/(x[n5]-x[n1]);
    dddf[n2]=6.0*(z1*f[n1]+z2*f[n2]+z3*f[n3]+z4*f[n4]+z5*f[n5]);
  
  z1=-(x[n2]+x[n3]+x[n4]+x[n5]-4.0*x[n1])/(x[n2]-x[n1])/(x[n3]-x[n1])/(x[n4]-x[n1])/(x[n5]-x[n1]);
  z2= (x[n3]+x[n4]+x[n5]-3.0*x[n1])/(x[n2]-x[n1])/(x[n3]-x[n2])/(x[n4]-x[n2])/(x[n5]-x[n2]);
  z3=-(x[n2]+x[n4]+x[n5]-3.0*x[n1])/(x[n3]-x[n2])/(x[n4]-x[n3])/(x[n3]-x[n1])/(x[n5]-x[n3]);
  z4= (x[n2]+x[n3]+x[n5]-3.0*x[n1])/(x[n4]-x[n3])/(x[n5]-x[n4])/(x[n4]-x[n2])/(x[n4]-x[n1]);
  z5=-(x[n2]+x[n3]+x[n4]-3.0*x[n1])/(x[n5]-x[n4])/(x[n5]-x[n3])/(x[n5]-x[n2])/(x[n5]-x[n1]);
  dddf[n1]=6.0*(z1*f[n1]+z2*f[n2]+z3*f[n3]+z4*f[n4]+z5*f[n5]);
  
  return 0;
}

/* Vosstanovlenie funktsii. START */
Local(Int2) Fspline (Vector x, Vector f, Vector b, Vector c, Vector d, VAR z, VAR* ff) {
Int2 i,n;
VAR t;

  n=staFunc->GetVectorDim(x);

  for (i=0; i<n && z>=x[i++]; );	/* the body is empty here! */
  /* ASSERT z<x[i] */
  i-=2;

  t=z-x[i];

  switch (staData->Jinterpol) {
    case 1:	/* cubic spline */
      *ff=f[i]+t*(b[i]+t*(c[i]+t*d[i]));
      break;
    case 0:	/* linear interpolation */
      *ff=f[i]+t*b[i];
      break;
  }
  
  return 0;
}

/* Lineinoe vosstanovlenie funktsii i proizvodnoi. START */
Local(Int2) FdFspline (Vector x, Vector f, Vector a, VAR z, VAR* ff, VAR* dff) {
Int2 i,n;
VAR t;

  n=staFunc->GetVectorDim(x);

  for (i=0; i<n && z>=x[i++]; );	/* the body is empty here! */
  /* ASSERT z<x[i] */
  i-=2;

  t=z-x[i];

  *ff=f[i]+t*a[i];
  *dff=a[i];
  
  return 0;
}

/* Rasstanovka uzlov po integralu. START */
Local(Int2) Integral (VAR C, Vector x, Vector g, Vector a, Vector xint, Vector Xint) {
          
Int2 i,j,n;
VAR CC, Cm, xx;

  n=staFunc->GetVectorDim(x)-1;
  xint[0]=0.0;
  Cm=C;
  xx=0.0;
    for (i=j=0; j<n-1;) {
    CC=(g[i]+g[i+1]+a[i]*(xx-x[i]))*(x[i+1]-xx)/2.;
    if(CC<=Cm) {
      Cm=Cm-CC;
      i=i+1;
      xx=x[i];
    } else {
      if(fabs(a[i])==0.0) xint[j+1]=xx+Cm/g[i];
      else xint[j+1]=x[i]+(sqrt((g[i]+a[i]*(xx-x[i]))*(g[i]+a[i]*(xx-x[i]))+2.0*Cm*a[i])-g[i])/a[i];
      xx=xint[j+1];
      j=j+1;
      Cm=C;
      }
    }
  xint[n]=1.0;
 
 
  Xint[n]=1.0;
  Cm=C;
  xx=1.0;
  for (i=j=n; j>1;) {
    CC=(g[i-1]+a[i-1]*(xx-x[i-1])/2.0)*(xx-x[i-1]);
    if(CC<=Cm) {
      Cm=Cm-CC;
      i=i-1;
      xx=x[i];
    } else {
      if(fabs(a[i-1])==0.0) Xint[j-1]=xx-Cm/g[i-1];
      else Xint[j-1]=x[i-1]+(sqrt((g[i-1]+a[i-1]*(xx-x[i-1]))*(g[i-1]+a[i-1]*(xx-x[i-1]))-2.0*Cm*a[i-1])-g[i-1])/a[i-1];
      xx=Xint[j-1];
      j=j-1;
      Cm=C;
      }
    }
  Xint[0]=0.0;

  for (i=1; i<n; i++) xint[i]=(xint[i]+Xint[i])/2;

  return 0;
}


/* Newton. START */
Local(Int2) Newtonstep (Vector x, Vector x0, Vector g, Vector a, VARPTR xshift, Int2 *INDEX) {

Int2 i,n,n1,n2;
VAR ss,ssold,varf,varfm;

  n=staFunc->GetVectorDim(x)-1;
  n1=n-1;
  n2=n-2;

  for (i=0; i<n; i++) FdFspline (x0, g, a, x[i], &gg[i], &dgg[i]);

  ssold=0.0;
  FP[1]=gg[0]*x[1]-gg[1]*(x[2]-x[1]);
  bet[1]=gg[0]+gg[1]-dgg[1]*(x[2]-x[1]);
  gam[1]=-gg[1];
  ssold=ssold+fabs(FP[1]);
  
  for (i=2; i<n1; i++) {FP[i]=gg[i-1]*(x[i]-x[i-1])-gg[i]*(x[i+1]-x[i]);
                        alf[i]=-gg[i-1]+dgg[i-1]*(x[i]-x[i-1]);
                        bet[i]=gg[i-1]+gg[i]-dgg[i]*(x[i+1]-x[i]);
                        gam[i]=-gg[i];
			if (fabs(bet[i])<fabs(alf[i])+fabs(gam[i])) {
			  *INDEX=1;
			  return 0;
			}
                        ssold=ssold+fabs(FP[i]);
                       }
		  
  FP[n1]=gg[n2]*(x[n1]-x[n2])-gg[n1]*(1.0-x[n1]);
  alf[n1]=-gg[n2]+dgg[n2]*(x[n1]-x[n2]);
  bet[n1]=gg[n2]+gg[n1]-dgg[n1]*(1.0-x[n1]);
  ssold=ssold+fabs(FP[n1]);
  
  ssold=ssold/n1;
  varfm=fabs(FP[1]);
  for (i=2; i<n; i++) {varf=fabs(FP[i]);
                       if(varf>=varfm) varfm=varf;
                      }
  
  aa[1]=-gam[1]/bet[1];
  bb[1]=-FP[1]/bet[1];
  for(i=2; i<n1; i++) {aa[i]=-gam[i]/(alf[i]*aa[i-1]+bet[i]);
                       bb[i]=-(FP[i]+alf[i]*bb[i-1])/(alf[i]*aa[i-1]+bet[i]);
                      }
  z[n1]=-(FP[n1]*aa[n2]+alf[n1]*bb[n2])/(alf[n1]*aa[n2]+bet[n1]);
  
  for(i=n2; i>0; i--) z[i]=aa[i]*z[i+1]+bb[i];		      
  
  xn[0]=0.0;
  for (i=1,*xshift=0.0; i<n; i++) {xn[i]=x[i]+z[i];
                             (*xshift)+=fabs(z[i]);}
  xn[n]=1.0;
  *xshift=*xshift/n1;
  
  for (i=0; i<n; i++) FdFspline (x0, g, a, xn[i], &gg[i], &dgg[i]);

  ss=0.0;
  FP[1]=gg[0]*xn[1]-gg[1]*(xn[2]-xn[1]);
  ss=ss+fabs(FP[1]);
  
  for (i=2; i<n1; i++) {FP[i]=gg[i-1]*(xn[i]-xn[i-1])-gg[i]*(xn[i+1]-xn[i]);
                        ss=ss+fabs(FP[i]);
                       }
		  
  FP[n1]=gg[n2]*(xn[n1]-xn[n2])-gg[n1]*(1.0-xn[n1]);
  ss=ss+fabs(FP[n1]);
  
  ss=ss/n1;
  varfm=fabs(FP[1]);
  for (i=2; i<n; i++) {varf=fabs(FP[i]);
                       if(varf>=varfm) varfm=varf;
                      }
  *INDEX=0;
  if (ss>=ssold)
    *INDEX=1; 
  else
    for(i=0; i<n; i++)
      if (xn[i+1]<=xn[i]) {
	*INDEX=1;
	break;
      }          
  if(*INDEX == 0) staFunc->CopyVector(xn,x);
  
  return 0;
}

/* Rasstanovka uzlov. START */
Local(Int2) Points (Vector x, Vector U,  Vector xn) {
Int2 i,j,it,INDEX,n;
Int2 maxit=5;   /*  bylo 5 */
VAR dxmin=1.0e-10;
VAR C,xshift,gs;
Vector xnewt=NULL, f, dddf, g, a, xint;

  f=staFunc->CreateVector(nmesh);
  dddf=staFunc->CreateVector(nmesh);
  g=staFunc->CreateVector(nmesh);
  a=staFunc->CreateVector(nmesh);
  xint=staFunc->CreateVector(nmesh);
  xnewt=staFunc->CreateVector(nmesh);
  
  n=staFunc->GetVectorDim(x)-1;
  
  for (i=0; i<nphase; i++) {
                           for (j=0; j<nmesh; j++) f[j]=U[i+j*nphase];
                           F3num (x, f, dddf);
                           for (j=0; j<nmesh; j++) g[j]+=dddf[j]*dddf[j];
                           }
        
  for (j=0,gs=0.0; j<nmesh; j++) {
                                 g[j]=sqrt(sqrt(g[j]/nphase));
                                 gs+=g[j];
                                 }
				 gs/=nmesh;
				 			 
  if (gs<=sqrt(d3fmin)) {for (j=0; j<nmesh; j++) xnewt[j]=j*(1.0/n);} 
  else {
       for (j=0; j<nmesh; j++) if (g[j]<=d3fmin) g[j]=d3fmin;
       for (j=0, C=0.0; j<n; j++) C+=(g[j+1]+g[j])*(x[j+1]-x[j])/2.0;
                                  C=C/n;
       for (j=0; j<n; j++) a[j]=(g[j+1]-g[j])/(x[j+1]-x[j]);
       Integral (C, x, g, a, xint,xnewt);

       staFunc->CopyVector(xint,xnewt);
       for(it=1; it<=maxit; it++) {
                                   Newtonstep ( xnewt, x, g, a, &xshift, &INDEX);
				   if(INDEX != 0) break;
				   if(xshift <= dxmin) break;}
       }  
/*				   staFunc->CopyVector(xnewt,xn); */
  for (i=1; i<nmesh; i++) xn[i]=(x[i]+xnewt[i])/2;

  f=staFunc->DestroyVector(f);
  dddf=staFunc->DestroyVector(dddf);
  g=staFunc->DestroyVector(g);
  a=staFunc->DestroyVector(a);
  xint=staFunc->DestroyVector(xint);
  xnewt=staFunc->DestroyVector(xnewt);
		  
  return 0;
}

#ifdef WIN32
static int _USERENTRY ss_SortCmp(const void *f, const void *s) {
#else
Entry(int) ss_SortCmp(const void *f, const void *s) {
#endif
  if (*(VARPTR )f==*(VARPTR )s) return 0;
  return (*(VARPTR )f<*(VARPTR )s) ? -1 : 1;
}

Entry(Int2) ss_JacDefStadyState (Vector x, rdMatrix Jac);

/**************/
/* Adaptation */
Entry(Int2) ss_Adapter(Vector U, Vector V) {
Int2 i,j,Ret;
  Points (X, U, XN);
  /* componentwise spline */
  for (i=0; i<nphase; i++) {
    for (j=0; j<nmesh; j++) FF[j]=U[i+j*nphase];
    Spline(X,FF,BB,CC,DD);
    for (j=0; j<nmesh; j++) Fspline (X, FF, BB, CC, DD, XN[j], &U[i+j*nphase]);  
  }
  staFunc->CopyVector(XN,X);

/* recomputation of the tangent vector  

  ss_JacDefStadyState (U,A);
  rdFunc->AppendMatrixVector(A,V,B);
  Ret=rdFunc->Decomp(B);
  if (Ret == 0) {
    staFunc->ClearVector(v0);
    v0[nmesh*nphase+nappend]=1.0;
    rdFunc->Solve(B,v0);
    staFunc->NormalizeVector(v0);
    staFunc->AddScaledVector(v0,V,-1,x1);
    printf("s=%g\n",staFunc->NormOfVector(x1));
    staFunc->CopyVector(v0,V);
  }  
*/  
  return (0);
}


#define _Dx       staData->dx		/* user's increment */
Local(VAR) _dx;				/* effective increment */

/* Initialize function common for Starter and Decoder */
Local(Int2) Init(StagenDataPtr stagenptr) {
  Int2 i,j;
  /* Set local pointers to data area and functions array */
  stagenData=stagenptr;
  staData=(rd__ss_DataPtr)stagenptr->StaPtr;
  /* For access to standard linea algebra */
  staFunc=(StdLinAlgDataPtr)stagenptr->StdLinAlg;
  staUtil=(UtilitiesDataPtr)stagenptr->Utilities;
  /* For access to RD specific linea algebra (from .con file) */
  rdFunc=(RdLinAlgDataPtr)stagenptr->StaFunc;
  /* user's increment -> effective increment (nearest exact binary value) */
  _dx=pow(2.0,-(int)(-log10(_Dx)/log10(2.0))-1);
  /* Get dimensions of Phase and Par spaces */
  npar=*stagenptr->ClassDim[CL_PAR];
  nphase=*stagenptr->ClassDim[CL_COMPONENTS];
  nspace=*stagenptr->ClassDim[CL_SPACE];	/* ==1 */
  /* Call udf SetInitPoint function, if any */
  if (!stagenData->ContDataGenLen && FuncParActive(staData->SetInitPoint)) {	/* there is UF to set InitPoint */
    Vector tmp;
    VAR PNTR work=NULL;
    Int2 Index_P;
    Boolean error=FALSE;
    Index_P=stagenptr->ParIndex(staData,P);
    staData->nmeshold=0;
    staData->space0=MemGet(sizeof(VAR)*nspace,FALSE);
    staData->coord0=MemGet(sizeof(VAR)*nphase,FALSE);
    if (FuncParCall(staData->SetInitPoint)(-1,staData)) {	/* init (e.g. fopen) */
	error=TRUE;
	staUtil->ErrMessage("Nonzero return code from function SetInitPoint (when ip=-1)"); 
    } else {       
      while (FuncParCall(staData->SetInitPoint)(staData->nmeshold,staData)==0) { /* while there is a point */
        (staData->nmeshold)++;
        if (work) work=MemMore(work,sizeof(VAR)*(nphase+nspace)*staData->nmeshold);
        else work=MemGet(sizeof(VAR)*(nphase+1),FALSE);
        MemCopy(work+(nphase+nspace)*(staData->nmeshold-1),staData->space0,sizeof(VAR)*nspace);
        MemCopy(work+(nphase+nspace)*(staData->nmeshold-1)+nspace,staData->coord0,sizeof(VAR)*nphase);
      }
      FuncParCall(staData->SetInitPoint)(-2,staData);	/* term (e.g. fclose) */
    }
    staData->coord0=MemFree(staData->coord0);
    staData->space0=MemFree(staData->space0);
    qsort(work,staData->nmeshold,sizeof(VAR)*(nphase+1),ss_SortCmp);
    staData->X=MemFree(staData->X);
    staData->U=MemFree(staData->U);
    staData->lenX=sizeof(VAR)*staData->nmeshold;
    staData->X=MemGet(staData->lenX,FALSE);		/* old mesh points */
    staData->lenU=sizeof(VAR)*nphase*staData->nmeshold;
    staData->U=MemGet(staData->lenU,FALSE);	/* old components */
    for (i=0; i<staData->nmeshold; i++) {
      staData->X[i]=work[i*(nphase+1)];
      MemCopy(staData->U+i*nphase,work+i*(nphase+1)+1,sizeof(VAR)*nphase);
    }
    work=MemFree(work);
    if (!error && staData->nmeshold<7) {
      error=TRUE;
      staUtil->ErrMessage("The number of mesh points must be more than 6");
    }
    if (!error && staData->X[0]!=0.0) {
      error=TRUE;
      staUtil->ErrMessage("The first mesh point must be 0.0");
    }
    if (!error && staData->X[staData->nmeshold-1]!=1.0) {
      error=TRUE;
      staUtil->ErrMessage("The last mesh point must be 1.0");
    }
    if (!error && staData->nmesh<7) {
      error=TRUE;
      staUtil->ErrMessage("The value of MeshPoints parameter must be more than 6");
    }
    for (i=1; i<staData->nmeshold; i++)
      if (!error && staData->X[i-1]==staData->X[i]) {
	error=TRUE;
	staUtil->ErrMessage("Two identical mesh points");
	break;
      }
    if (error) {
      staData->X=MemFree(staData->X);
      staData->U=MemFree(staData->U);
      staData->lenX=0;
      staData->lenU=0;
      return 10;
    }
    stagenptr->UpdatePar(Index_P);
  } else {
    if (staData->X==NULL) {
      staUtil->ErrMessage("The initial point was not selected.\n"
                          "Specify InitPointProc or pick up a computed point");
      return 1;
    }
  }
  /* interpolate to reqested number of mesh points */
  if (staData->nmeshold!=staData->nmesh) {
    VAR PNTR Xnew;
    VAR PNTR Unew;
    staData->lenX=sizeof(VAR)*staData->nmesh;
    staData->lenU=sizeof(VAR)*staData->nmesh*nphase;
    Xnew=MemGet(staData->lenX,FALSE);
    Unew=MemGet(staData->lenU,FALSE);
    FF=staFunc->CreateVector(staData->nmeshold);
    X=staFunc->CreateVector(staData->nmeshold);
    BB=staFunc->CreateVector(staData->nmeshold);
    CC=staFunc->CreateVector(staData->nmeshold);
    DD=staFunc->CreateVector(staData->nmeshold);
    MemCopy(X,staData->X,sizeof(VAR)*staData->nmeshold);
    for (i=0; i<staData->nmesh; i++)
      Xnew[i]=i*1.0/(staData->nmesh-1);
    /* componentwise spline */
    for (i=0; i<nphase; i++) {
      for (j=0; j<staData->nmeshold; j++)
	FF[j]=staData->U[i+j*nphase];
      Spline(X,FF,BB,CC,DD);
      for (j=0; j<staData->nmesh; j++)
	Fspline (X, FF, BB, CC, DD, Xnew[j], &Unew[i+j*nphase]);  
    }
    FF=staFunc->DestroyVector(FF);
    X=staFunc->DestroyVector(X);
    BB=staFunc->DestroyVector(BB);
    CC=staFunc->DestroyVector(CC);
    DD=staFunc->DestroyVector(DD);
    staData->X=MemFree(staData->X);
    staData->U=MemFree(staData->U);
    staData->X=Xnew;
    staData->U=Unew;
    staData->nmeshold=staData->nmesh;
  }
  /* Do all the checks related to active parameters and find them */
  for (i=0,nappend=0; i<nUTest; i++) if (ActiveUTestFuns[i]==UDF_APPEND) nappend++;
  active=MemNew(sizeof(Int2)*(nappend+1));    
  if (staUtil->CheckParameters(staData->Active,npar,active,1+nappend)) {
    active=MemFree(active);
    return 1;
  }
  /* Checks related to analytical/numerical differentiation of RHS */
  if (!stagenData->Der1 && staData->dx <= 0) {
    staUtil->ErrMessage("Analytical Jacoby matrix is undefined.\nRecompile system or set positive 'Increment'");
    return 2;
  }
  NumDer1=stagenData->Der1==NULL;
  if ((staData->ActiveSing[0] || staData->ActiveSing[1] || staData->ActiveSing[2]) && 
      !stagenData->Der2 && staData->dx <= 0) {
    staUtil->ErrMessage("Analytical Hessian matrix is undefined.\nRecompile system or set positive 'Increment'");
    active=MemFree(active);
    return 3;
  }
  if (staData->ActiveSing[2] && !stagenData->Der3 && staData->dx <= 0) {
    staUtil->ErrMessage("Analytical matrix of third derivatives is undefined.\nRecompile system or set positive 'Increment'");
    active=MemFree(active);
    return 4;
  }
  if (staData->dx <= 0 && nappend > 0) {
    staUtil->ErrMessage("Set positive 'Increment'");
    active=MemFree(active);
    return 5;
  }
  nmesh=staData->nmesh;
  /* Create a copy parameters */
  if (!P) {
    P=MemNew(sizeof(FloatHi)*npar);
    MemCopy(P,staData->P,sizeof(FloatHi)*npar);
    /* Allocate memory */
    x0=staFunc->CreateVector(nmesh*nphase+1+nappend);
    x1=staFunc->CreateVector(nmesh*nphase+1+nappend);
    v0=staFunc->CreateVector(nmesh*nphase+1+nappend);
    /* MinMax */
    STF=staFunc->CreateVector(1+2*nphase);    /* L2 norm, Min & Max */
    /* Initialize RD algebra package */
    rdFunc->Init(staFunc,nphase,nmesh,nappend);
    /* Allocate memory for matrices used in test functions and in the default processor */
    A0=rdFunc->CreateMatrix(nmesh*nphase,nmesh*nphase+1+nappend);
    A=rdFunc->CreateMatrix(nmesh*nphase+nappend,nmesh*nphase+1+nappend);
    B=rdFunc->CreateMatrix(nmesh*nphase+1+nappend,nmesh*nphase+1+nappend);
    BifVector=MemNew(sizeof(FloatHi)*(4*nmesh*nphase+3+3*nappend));   
    /* global */
    X=staFunc->CreateVector(nmesh);		/* mesh points (grid) */
    MemCopy(X,staData->X,sizeof(VAR)*nmesh);
    U_0=staFunc->CreateVector(nphase);		/* u */
    U_X=staFunc->CreateVector(nphase);		/* u' */
    U_XX=staFunc->CreateVector(nphase);		/* u'' */
    /* Adapter's workspace */
    XN=staFunc->CreateVector(nmesh);
    FF=staFunc->CreateVector(nmesh);
    BB=staFunc->CreateVector(nmesh);
    CC=staFunc->CreateVector(nmesh);
    DD=staFunc->CreateVector(nmesh);
    /* Newtonstep's work space */
    xn=staFunc->CreateVector(nmesh);
    gg=staFunc->CreateVector(nmesh);
    dgg=staFunc->CreateVector(nmesh);
    alf=staFunc->CreateVector(nmesh-1);
    bet=staFunc->CreateVector(nmesh-1);
    gam=staFunc->CreateVector(nmesh-1);
    FP=staFunc->CreateVector(nmesh-1);
    aa=staFunc->CreateVector(nmesh);
    bb=staFunc->CreateVector(nmesh);
    z=staFunc->CreateVector(nmesh);
    /* DefFunc work space */
    AppMat=staFunc->CreateMatrix(nappend,nphase*nmesh+1+nappend);
    JacW=staFunc->CreateMatrix(nphase,nphase);
    JacV=staFunc->CreateMatrix(nphase,nphase);
  }

  Der2num=ind2_(nphase+npar-1,nphase+npar-1)+1;

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

  return 0;
}

Local(void) Term(Boolean OK) {
    stagenData->BifDataOutLen=0;
    if (OK&&P&&(stagenData->Reason==RF_CLEAN)) {
      stagenData->ContDataStaLen=nmesh*sizeof(VAR);
      stagenData->ContDataStaPtr=MemNew(stagenData->ContDataStaLen);
      MemCopy(stagenData->ContDataStaPtr,X,nmesh*sizeof(VAR));
      stagenData->ProcessPoint(stagenData->ProcFuncNum, x1, "Last point");
    }
    active=MemFree(active);
    if (P) {
      P=MemFree(P);
      E=MemFree(E);
      x0=staFunc->DestroyVector(x0);
      x1=staFunc->DestroyVector(x1);
      v0=staFunc->DestroyVector(v0);
      STF=staFunc->DestroyVector(STF);
      Re=staFunc->DestroyVector(Re);
      Im=staFunc->DestroyVector(Im);
      A=rdFunc->DestroyMatrix(A);
      A0=rdFunc->DestroyMatrix(A0);
      B=rdFunc->DestroyMatrix(B);
      BifVector=MemFree(BifVector);
      /* free working space */
      XN=staFunc->DestroyVector(XN);
      FF=staFunc->DestroyVector(FF);
      BB=staFunc->DestroyVector(BB);
      CC=staFunc->DestroyVector(CC);
      DD=staFunc->DestroyVector(DD);
      X=staFunc->DestroyVector(X);
      U_0=staFunc->DestroyVector(U_0);
      U_X=staFunc->DestroyVector(U_X);
      U_XX=staFunc->DestroyVector(U_XX);
      alf=staFunc->DestroyVector(alf);
      bet=staFunc->DestroyVector(bet);
      gam=staFunc->DestroyVector(gam);
      FP=staFunc->DestroyVector(FP);
      xn=staFunc->DestroyVector(xn);
      aa=staFunc->DestroyVector(aa);
      bb=staFunc->DestroyVector(bb);
      z=staFunc->DestroyVector(z);
      gg=staFunc->DestroyVector(gg);
      dgg=staFunc->DestroyVector(dgg);
      AppMat=staFunc->DestroyMatrix(AppMat);
      JacW=staFunc->DestroyMatrix(JacW);
      JacV=staFunc->DestroyMatrix(JacV);
    }
  /* Clean up RD algebra package */
  rdFunc->Term();
}

/**********************************************/
/* Starters to continue the stady state curve */

Local(Int2) SimpleInit(Vector X1, Vector X0, Vector V0);
Local(Int2) DefRHS (Vector u, Vector f);

/* dd->ss Starter */
Entry(Int2) d_ss_Starter(StagenDataPtr stagenptr) {
Int2 i;
  if (!stagenptr) {	/* this is request to terminate */
    Term(TRUE);
    return 0;
  }
  /* initialize local variables */
  if (Init(stagenptr)) return 1;
  /* Create the first point for generator */

  if (!stagenptr->ContDataGenLen) {   /* not "Continuation" */
    MemCopy(x1,staData->U,sizeof(FloatHi)*nmesh*nphase);
    for (i=0; i<=nappend; i++) x1[nmesh*nphase+i]=staData->P[active[i]];
    /* Initialize x0,v0 */
    if (SimpleInit(x1,x0,v0)) {
      staUtil->ErrMessage("No convergence to a stady state from the initial distribution"); 
      Term(FALSE);
      return 1;  
    }
    if (((ContDataPtr)stagenptr->GenPtr)->nadapt) ss_Adapter(x0,v0);
  } else {
    MemCopy(X,stagenptr->ContDataStaPtr,sizeof(VAR)*nmesh);
    stagenptr->ContDataStaPtr=MemFree(stagenptr->ContDataStaPtr);
    stagenptr->ContDataStaLen=0;
  }

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

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

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

Entry(Int2) ss_Starter(StagenDataPtr stagenptr) {
  Int2 i,Ret;
  Vector v1;
  rdMatrix D,DE;
  if (!stagenptr) {	/* this is request to terminate */
    Term(TRUE);
    return 0;
  }
  if (Init(stagenptr)) return 1;	/* initialize local variables */

  /* Create the first point for generator */

  if (!stagenptr->ContDataGenLen) {   /* not "Continuation" */
    v1=staFunc->CreateVector(nmesh*nphase+1+nappend);
    D=rdFunc->CreateMatrix(nmesh*nphase+nappend,nmesh*nphase+1+nappend);
    DE=rdFunc->CreateMatrix(nmesh*nphase+1+nappend,nmesh*nphase+1+nappend);
    MemCopy(x0,staData->U,sizeof(FloatHi)*nmesh*nphase);
    for (i=0; i<=nappend; i++) x0[nmesh*nphase+i]=staData->P[active[i]];
    Ret=ss_JacDefStadyState (x0,D);
    for (i=0; i<nmesh*nphase+1+nappend; i++) {
      v1[i]=1.0;
      rdFunc->AppendMatrixVector(D,v1,DE);
      Ret=rdFunc->Decomp(DE);
      if (Ret == 0) {
        v0[nmesh*nphase+nappend]=stagenData->Forward ? 1 : -1;
        rdFunc->Solve(DE,v0);
        break;
      };
      v1[i]=0.0;
    };
    if (((ContDataPtr)stagenptr->GenPtr)->nadapt) ss_Adapter(x0,v0);
    v1=staFunc->DestroyVector(v1);
    D=rdFunc->DestroyMatrix(D);
    DE=rdFunc->DestroyMatrix(DE);
    if (Ret) {
       staUtil->ErrMessage(
       "Unable to find a tangent vector to the stady state curve at the initial point"); 
       Term(FALSE);
       return 1;
    };
    staFunc->NormalizeVector(v0);
  } else {
    MemCopy(X,stagenptr->ContDataStaPtr,sizeof(VAR)*nmesh);
    stagenptr->ContDataStaPtr=MemFree(stagenptr->ContDataStaPtr);
    stagenptr->ContDataStaLen=0;
  }
  /* Initialize continuation */
  ((ContDataPtr)stagenptr->GenPtr)->X0=x0;
  ((ContDataPtr)stagenptr->GenPtr)->V0=v0;
  /* That's all */
  return 0;
}

/* bp-ss Starter  */

#define _branch   staData->branch     /* switch data  */
/* _dx: see above */

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

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

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

  /* Initialize the first point */

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

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

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

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


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

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

    Vector F,W=NULL,Z=NULL;
    rdMatrix A, B0, B=NULL;
    VAR errx, errf;
    
    ret=2;
    F = staFunc->CreateVector(n-1);
    A = rdFunc->CreateMatrix (n-1,n);
    B0 = rdFunc->CreateMatrix (n,n);

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

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

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

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

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

    return (ret);
}


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

/************************/
/* Some common routines */
Local(void) CopyToX(Vector vp) {
}

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

Local(void) CopyToPa(VARPTR  p) {
  Int2 i;
  for (i=0; i<1+nappend; i++) P[active[i]]=p[i];
}

Local(void) CopyToE(Vector Re, Vector Im) {
}

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

Entry(Int2) ss_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 nmesh;
    default:		/* extract next M-point from G-point pointed by vpoint */
      /* vpoint is a Vector */
      indirect[CL_TIME]=&zero;
      indirect[CL_SPACE]=X+oper-1;
      indirect[CL_COMPONENTS]=(VARPTR )vpoint+(oper-1)*nphase;
      indirect[CL_PAR]=P;
      CopyToP(vpoint);
      indirect[CL_EIGV]=&zero;
      indirect[CL_TEST]=test;
      indirect[CL_STF]=STF;
  }
  return 0;
}

Local(Vector) U;	/* vector came from generator */
Local(Int2) ipoint;

/***************************/
/* User-defined functions  */
Entry(Int2) ss_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;
}

/************************/
/* Caclulate u, u' and u'' */

Local(VAR) phil0,phil1,phil2;
Local(VAR) phim,phi,phip,psim,psi,psip,etam,eta,etap;
Local(VAR) phirN,phirN1,phirN2;
Local(VAR) Xi;

Local(void) Der1ULeft(Vector u) {
  Int2 j;
  phil0=-(X[1]+X[2])/(X[1]*X[2]);
  phil1= X[2]/(X[1]*(X[2]-X[1]));
  phil2=-X[1]/(X[2]*(X[2]-X[1]));
  for (j=0; j<nphase; j++)
    U_X[j]=phil0*u[j]+phil1*u[j+nphase]+phil2*u[j+2*nphase];
}

Local(void) Der12UIntern(Vector u, Int2 i) {
  Int2 j,b=nphase*i;
  VAR h,H,s;
  h=X[i]-X[i-1];
  H=X[i+1]-X[i];
  s=h+H;

/* Doedel's formulas */

  phim=-(h+s)/(3*h*s);
  phi=(H-h)/(3*h*H);
  phip=(H+s)/(3*H*s);
  psim=2/(h*s);
  psi=-2/(h*H);
  psip=2/(H*s);
  etam=(h-H)*(H+s)/(9*h*s);
  eta=(h+s)*(H+s)/(9*h*H);
  etap=(H-h)*(h+s)/(9*H*s);
  Xi=X[i]+(H-h)/3;
  for (j=0; j<nphase; j++) {
    U_0[j]= etam*u[b-nphase+j]+eta*u[b+j]+etap*u[b+nphase+j];
    U_X[j]= phim*u[b-nphase+j]+phi*u[b+j]+phip*u[b+nphase+j];
    U_XX[j]=psim*u[b-nphase+j]+psi*u[b+j]+psip*u[b+nphase+j];
  }
}

Local(void) Der1URight(Vector u) {
  Int2 N=nmesh-1;
  Int2 j,b=nphase*N;
  phirN=(2-X[N-1]-X[N-2])/((1-X[N-1])*(1-X[N-2]));
  phirN1=(X[N-2]-1)/((1-X[N-1])*(X[N-1]-X[N-2]));
  phirN2=(1-X[N-1])/((X[N-1]-X[N-2])*(1-X[N-2]));
  for (j=0; j<nphase; j++)
    U_X[j]=phirN2*u[b-2*nphase+j]+phirN1*u[b-nphase+j]+phirN*u[b+j];
}

/*******/
/* RHS */
Local(Int2) DefRHS (Vector u, Vector f) {
  Int2 N=nmesh-1;
  Int2 i,j,b;
  CopyToP(u);
  b=0;
  /* u'(0) */
  Der1ULeft(u);
  if (((BCPtr)(stagenData->AuxFunc[BC0_I]))(u+b,U_X,P,&zero,f+b))	/* left boundary condition */
    return 1;
  b+=nphase;
  /* u(i),   i=1,...,N-1 */
  /* u'(i),  i=1,...,N-1 */
  /* u''(i), i=1,...,N-1 */
  for (i=1; i<N; i++, b+=nphase) {
    Der12UIntern(u,i);
    if (stagenData->Rhs(U_0,U_X,U_XX,&Xi,P,&zero,f+b))	/* equation */
      return 2;
  }
  /* u'(N) */
  Der1URight(u);
  if (((BCPtr)(stagenData->AuxFunc[BC1_I]))(u+b,U_X,P,&zero,f+b))	/* right boundary condition */
    return 3;
  return 0;
}

/**************************************/
/* Defining function for stady states */
Entry(Int2) ss_DefStadyState (Vector x, Vector y) {
  DefRHS(x,y);
  if (nappend) ss_DefUser(x,y+nmesh*nphase);
  return 0;
}

#include "rd_symd.c"

/* f0(u,.,.,.) */
Entry(Int2) ss_BC0u(VARPTR  u0, Vector f) {
  if (((BCPtr)(stagenData->AuxFunc[BC0_I]))(u0,U_X,P,&zero,f))	/* left boundary condition */
    return 1;
  return 0;
}

/* f0(.,u_x,.,.) */
Entry(Int2) ss_BC0u_x(VARPTR  u_x0, Vector f) {
  if (((BCPtr)(stagenData->AuxFunc[BC0_I]))(U+0,u_x0,P,&zero,f))	/* left boundary condition */
    return 1;
  return 0;
}

/* f0(.,.,p,.) */
Entry(Int2) ss_BC0p(VARPTR  p, Vector f) {
  CopyToPa(p);
  if (((BCPtr)(stagenData->AuxFunc[BC0_I]))(U+0,U_X,P,&zero,f))	/* left boundary condition */
    return 1;
  return 0;
}

/* F(u,.,.,.,.,.) */
Entry(Int2)  ss_RHSu(VARPTR  u, Vector f) {
  if (stagenData->Rhs(u,U_X,U_XX,&Xi,P,&zero,f))		/* equation */
    return 1;
  return 0;
}

/* F(.,u_x,.,.,.,.) */
Entry(Int2)  ss_RHSu_x(VARPTR  u_x, Vector f) {
  if (stagenData->Rhs(U_0,u_x,U_XX,&Xi,P,&zero,f))		/* d(equation)/d(u') */
    return 1;
  return 0;
}

/* F(.,.,u_xx,.,.,.) */
Entry(Int2)  ss_RHSu_xx(VARPTR  u_xx, Vector f) {
  if (stagenData->Rhs(U_0,U_X,u_xx,&Xi,P,&zero,f))		/* d(equation)/d(u'') */
    return 1;
  return 0;
}

/* F(.,.,.,.,p,.) */
Entry(Int2)  ss_RHSp(VARPTR  p, Vector f) {
  CopyToPa(p);
  if (stagenData->Rhs(U_0,U_X,U_XX,&Xi,P,&zero,f))		/* d(equation)/d(p) */
    return 1;
  return 0;
}

/* f1(u,.,.,.) */
Entry(Int2) ss_BC1u(VARPTR  u1, Vector f) {
  if (((BCPtr)(stagenData->AuxFunc[BC1_I]))(u1,U_X,P,&zero,f))	/* right boundary condition */
    return 1;
  return 0;
}

/* f1(.,u_x,.,.) */
Entry(Int2) ss_BC1u_x(VARPTR  u_x1, Vector f) {
  if (((BCPtr)(stagenData->AuxFunc[BC1_I]))(U+(nmesh-1)*nphase,u_x1,P,&zero,f))	/* left boundary condition */
    return 1;
  return 0;
}

/* f1(.,.,p,.) */
Entry(Int2) ss_BC1p(VARPTR  p, Vector f) {
  CopyToPa(p);
  if (((BCPtr)(stagenData->AuxFunc[BC1_I]))(U+(nmesh-1)*nphase,U_X,P,&zero,f))	/* left boundary condition */
    return 1;
  return 0;
}

/* Copy values of analitical derivatives to a[i],b[i],c[i] */
Local(void) CopyToMatrixU(Matrix tm, VARPTR sa, Int2 firstcol, Int2 rowlen) {
  Int2 i;
  for (i=0; i<nphase; i++)
    MemCopy(tm[i],sa+i*rowlen+firstcol,sizeof(VAR)*nphase);
}

/* Copy values of analitical derivatives to delta[i] */
Local(void) CopyToMatrixP(Matrix tm, VARPTR sa, Int2 firstcol, Int2 rowlen) {
  Int2 i,j;
  for (i=0; i<nphase; i++)
    for (j=0; j<=nappend; j++)
      tm[i][j]=sa[i*rowlen+firstcol+active[j]];
}

/********************/
/* Jacobian of RHS  */
Local(Int2) JacRHS (Vector u, rdMatrix Jac) {
#define N (nmesh-1)
  Matrix ma,mb,mc,md;
  VARPTR dm;
  Int2 i;
  U=u;
  CopyToP(u);
  /* a[0] b[0] c[0] delta[0] */
  ma=rdFunc->GetMatrixA(Jac,0);
  mb=rdFunc->GetMatrixB(Jac,0);
  mc=rdFunc->GetMatrixC(Jac,0);
  md=rdFunc->GetMatrixDelta(Jac,0);
  Der1ULeft(u);
  if (NumDer1) {
    if (staFunc->DerA((DiffFun)ss_BC0u, nphase, u+0, _dx, JacW)) return 1;
    if (staFunc->DerA((DiffFun)ss_BC0u_x, nphase, U_X, _dx, ma)) return 1;
  } else {
    dm=((BCDPtr)(stagenData->AuxFunc[BC0D1_I]))(u+0,U_X,P,&zero);
    if (!dm) return 11;
    CopyToMatrixU(JacW,dm,0,2*nphase+npar);
    CopyToMatrixU(ma,dm,nphase,2*nphase+npar);
  }
  staFunc->CopyMatrix(ma,mb);
  staFunc->CopyMatrix(ma,mc);
  staFunc->AddScaledMatrix(JacW,ma,phil0,ma);
  staFunc->AddScaledMatrix(NULL,mb,phil1,mb);
  staFunc->AddScaledMatrix(NULL,mc,phil2,mc);
  if (NumDer1) {
    if (staFunc->DerA((DiffFun)ss_BC0p, nphase, u+nmesh*nphase, _dx, md)) return 1;
  } else
    CopyToMatrixP(md,dm,2*nphase,2*nphase+npar);
  CopyToP(u);
  /* a[i] b[i] c[i] delta[i], i=1,...,N-1 */
  for (ipoint=1; ipoint<N; ipoint++) {
    ma=rdFunc->GetMatrixA(Jac,ipoint);
    mb=rdFunc->GetMatrixB(Jac,ipoint);
    mc=rdFunc->GetMatrixC(Jac,ipoint);
    md=rdFunc->GetMatrixDelta(Jac,ipoint);
    Der12UIntern(u,ipoint);
    if (NumDer1) {
      if (staFunc->DerA((DiffFun)ss_RHSu, nphase, U_0, _dx, JacV)) return 1;
      if (staFunc->DerA((DiffFun)ss_RHSu_x, nphase, U_X, _dx, JacW)) return 1;
      if (staFunc->DerA((DiffFun)ss_RHSu_xx, nphase, U_XX, _dx, ma)) return 1;
    } else {
      dm=stagenData->Der1(U_0,U_X,U_XX,&Xi,P,&zero);
      if (!dm) return 12;
      CopyToMatrixU(JacV,dm,0,3*nphase+nspace+npar);
      CopyToMatrixU(JacW,dm,nphase,3*nphase+nspace+npar);
      CopyToMatrixU(ma,dm,2*nphase,3*nphase+nspace+npar);
    }
    staFunc->AddScaledMatrix(NULL,ma,psi,mb);
    staFunc->AddScaledMatrix(NULL,ma,psip,mc);
    staFunc->AddScaledMatrix(NULL,ma,psim,ma);
    staFunc->AddScaledMatrix(ma,JacW,phim,ma);
    staFunc->AddScaledMatrix(ma,JacV,etam,ma);
    staFunc->AddScaledMatrix(mb,JacW,phi,mb);
    staFunc->AddScaledMatrix(mb,JacV,eta,mb);
    staFunc->AddScaledMatrix(mc,JacW,phip,mc);
    staFunc->AddScaledMatrix(mc,JacV,etap,mc);
    if (NumDer1) {
      if (staFunc->DerA((DiffFun)ss_RHSp, nphase, u+nmesh*nphase, _dx, md)) return 1;
    } else
      CopyToMatrixP(md,dm,3*nphase+nspace,3*nphase+nspace+npar);
    CopyToP(u);
  }
  /* a[N] b[N] c[N] delta[N] */
  ma=rdFunc->GetMatrixA(Jac,N);
  mb=rdFunc->GetMatrixB(Jac,N);
  mc=rdFunc->GetMatrixC(Jac,N);
  md=rdFunc->GetMatrixDelta(Jac,N);
  Der1URight(u);
  if (NumDer1) {
    if (staFunc->DerA((DiffFun)ss_BC1u, nphase, u+N*nphase, _dx, JacW)) return 1;
    if (staFunc->DerA((DiffFun)ss_BC1u_x, nphase, U_X, _dx, mc)) return 1;
  } else {
    dm=((BCDPtr)(stagenData->AuxFunc[BC1D1_I]))(u+N*nphase,U_X,P,&zero);
    if (!dm) return 13;
    CopyToMatrixU(JacW,dm,0,2*nphase+npar);
    CopyToMatrixU(mc,dm,nphase,2*nphase+npar);
  }
  staFunc->CopyMatrix(mc,mb);
  staFunc->CopyMatrix(mc,ma);
  staFunc->AddScaledMatrix(JacW,mc,phirN,mc);
  staFunc->AddScaledMatrix(NULL,mb,phirN1,mb);
  staFunc->AddScaledMatrix(NULL,ma,phirN2,ma);
  if (NumDer1) {
    if (staFunc->DerA((DiffFun)ss_BC1p, nphase, u+nmesh*nphase, _dx, md)) return 1;
  } else
    CopyToMatrixP(md,dm,2*nphase,2*nphase+npar);
  CopyToP(u);
  return 0;
#undef N
}
/*****************************************************/
/* Jacobian of defining conditions for stady states  */
Entry(Int2) ss_JacDefStadyState (Vector u, rdMatrix Jac) {
  if (JacRHS(u, Jac)) return 1;
  if (nappend) {
    staFunc->DerA(ss_DefUser, nappend, u, _dx, AppMat);
    for (ipoint=0; ipoint<nmesh; ipoint++) {
      staFunc->CopyMatrixBlock(AppMat, 0,nphase*ipoint, 
                             rdFunc->GetMatrixP(Jac,ipoint),0,0,nappend,nphase);
    }
    staFunc->CopyMatrixBlock(AppMat,0,nphase*nmesh, 
                             rdFunc->GetMatrixFi(Jac),0,0,nappend,nappend+1);
  }
  return 0;
}
#define GessElem(i,j,k) (*(H+Der2num*(i)+ind2_(j,k)))
#define MixedGessElem(i,j,k) GessElem(i,j,nphase+active[k])
#define GessParElem(i,j,k) GessElem(i,nphase+active[j],nphase+active[k])

/***********/
/* L2 Norm */

/* see appropriate header file for values of staData->Jinterpol */
Local(VAR) L2_Norm(Vector U) {
Int2 i,j;
VAR h,fi,s;
  /* componentwise spline */
  for (j=0,s=0; j<nphase; j++) {
    for (i=0; i<nmesh; i++) FF[i]=U[j+i*nphase];
    Spline(X,FF,BB,CC,DD);
    switch (staData->Jinterpol) {
      case 1:		/* cubic spline */
	for (i=0,fi=0; i<nmesh-1; i++) {
	  h=X[i+1]-X[i];
	  fi+=h*(FF[i]*FF[i]+
	      h*(FF[i]*BB[i]+
	      h*((BB[i]*BB[i]+2*FF[i]*CC[i])/3+
	      h*((FF[i]*DD[i]+BB[i]*CC[i])/2+
	      h*((2*BB[i]*DD[i]+CC[i]*CC[i])/5+
	      h*DD[i]*(CC[i]/3+h*DD[i]/7))))));  
	}
	s+=fi;
	break;
      case 0:		/* linear interpolation */
	for (i=0,fi=0; i<nmesh-1; i++) {
	  h=X[i+1]-X[i];
	  fi+=h*(FF[i]*FF[i]+
	      h*(FF[i]*BB[i]+
	      h*(BB[i]*BB[i]/3)));  
	}
	s+=fi;
	break;
    }
  }
  s=sqrt(s/nphase);
  return s;

}

/**********/
/* MinMax */
Local(void) MinMax(Vector U) {
  VAR z,Uz,z1,z2;
  Int2 i,j;
#define fmin	(STF+1)
#define fmax	(STF+1+nphase)
  for (i=0; i<nphase; i++) {
    for (j=0; j<nmesh; j++) FF[j]=U[i+j*nphase];
    fmin[i]=FF[0];
    fmax[i]=FF[0];
    switch (staData->Jinterpol) {
      case 0:		/* linear interpolation */
	for (j=1; j<nmesh; j++){
	  if (FF[j]>fmax[i]) fmax[i]=FF[j];
	  if (FF[j]<fmin[i]) fmin[i]=FF[j];
	}
	break;
      case 1:		/* cubic spline */
	Spline (X,FF,BB,CC,DD);
	for (j=0; j<nmesh-1; j++) {
	  if (fabs(DD[j])<SKO_EPS) {
	    if (fabs(CC[j])<SKO_EPS) {
	      if (FF[j]>fmax[i]) fmax[i]=FF[j];
	      if (FF[j]<fmin[i]) fmin[i]=FF[j];
	      if (FF[j+1]>fmax[i]) fmax[i]=FF[j+1];
	      if (FF[j+1]<fmin[i]) fmin[i]=FF[j+1];
	    } else {
	      z=-BB[j]/(2*CC[j]);
	      if ((z>X[j]) && (z<X[j+1])) {
		if (FF[j]>fmax[i]) fmax[i]=FF[j];
		if (FF[j]<fmin[i]) fmin[i]=FF[j];
		if (FF[j+1]>fmax[i]) fmax[i]=FF[j+1];
		if (FF[j+1]<fmin[i]) fmin[i]=FF[j+1];
		Uz=FF[j]-BB[j]*BB[j]/(4*CC[j]);
		if (Uz>fmax[i]) fmax[i]=Uz;
		if (Uz<fmin[i]) fmin[i]=Uz;
	      } else {
		if (FF[j]>fmax[i]) fmax[i]=FF[j];
		if (FF[j]<fmin[i]) fmin[i]=FF[j];
		if (FF[j+1]>fmax[i]) fmax[i]=FF[j+1];
		if (FF[j+1]<fmin[i]) fmin[i]=FF[j+1];
	      }
	    }
	  } else {
	    z=CC[j]*CC[j]-3*DD[j]*BB[j];
	    if (z<0) {
	      if (FF[j]>fmax[i]) fmax[i]=FF[j];
	      if (FF[j]<fmin[i]) fmin[i]=FF[j];
	      if (FF[j+1]>fmax[i]) fmax[i]=FF[j+1];
	      if (FF[j+1]<fmin[i]) fmin[i]=FF[j+1];
	      continue;
	    }
	    z=sqrt(z);
	    z1=(-CC[j]-z)/(3*DD[j]);
	    z2=(-CC[j]+z)/(3*DD[j]);
	    if (z2<X[j]) {
	      if (FF[j]>fmax[i]) fmax[i]=FF[j];
	      if (FF[j]<fmin[i]) fmin[i]=FF[j];
	      if (FF[j+1]>fmax[i]) fmax[i]=FF[j+1];
	      if (FF[j+1]<fmin[i]) fmin[i]=FF[j+1];
	      continue;
	    }
	    if (z1<X[j] && X[j]<=z2 && z2<=X[j+1]) {
	      if (FF[j]>fmax[i]) fmax[i]=FF[j];
	      if (FF[j]<fmin[i]) fmin[i]=FF[j];
	      if (FF[j+1]>fmax[i]) fmax[i]=FF[j+1];
	      if (FF[j+1]<fmin[i]) fmin[i]=FF[j+1];
	      z=z2-X[j];
	      Uz=FF[j]+z*(BB[j]+z*(CC[j]+z*DD[j]));
	      if (Uz>fmax[i]) fmax[i]=Uz;
	      if (Uz<fmin[i]) fmin[i]=Uz;
	      continue;
	    }
	    if (z1<X[j] && z2>X[j+1]) {
	      if (FF[j]>fmax[i]) fmax[i]=FF[j];
	      if (FF[j]<fmin[i]) fmin[i]=FF[j];
	      if (FF[j+1]>fmax[i]) fmax[i]=FF[j+1];
	      if (FF[j+1]<fmin[i]) fmin[i]=FF[j+1];
	      continue;
	    }
	    if ((X[j]<=z1 && z2<=X[j+1])) {
	      if (FF[j]>fmax[i]) fmax[i]=FF[j];
	      if (FF[j]<fmin[i]) fmin[i]=FF[j];
	      if (FF[j+1]>fmax[i]) fmax[i]=FF[j+1];
	      if (FF[j+1]<fmin[i]) fmin[i]=FF[j+1];
	      z=z1-X[j];
	      Uz=FF[j]+z*(BB[j]+z*(CC[j]+z*DD[j]));
	      if (Uz>fmax[i]) fmax[i]=Uz;
	      if (Uz<fmin[i]) fmin[i]=Uz;
	      z=z2-X[j];
	      Uz=FF[j]+z*(BB[j]+z*(CC[j]+z*DD[j]));
	      if (Uz>fmax[i]) fmax[i]=Uz;
	      if (Uz<fmin[i]) fmin[i]=Uz;
	      continue;
	    }
	    if (X[j]<=z1 && z1<=X[j+1] && X[j+1]<=z2) {
	      if (FF[j]>fmax[i]) fmax[i]=FF[j];
	      if (FF[j]<fmin[i]) fmin[i]=FF[j];
	      if (FF[j+1]>fmax[i]) fmax[i]=FF[j+1];
	      if (FF[j+1]<fmin[i]) fmin[i]=FF[j+1];
	      z=z1-X[j];
	      Uz=FF[j]+z*(BB[j]+z*(CC[j]+z*DD[j]));
	      if (Uz>fmax[i]) fmax[i]=Uz;
	      if (Uz<fmin[i]) fmin[i]=Uz;
	      continue;
	    }
	    if (FF[j]>fmax[i]) fmax[i]=FF[j];
	    if (FF[j]<fmin[i]) fmin[i]=FF[j];
	    if (FF[j+1]>fmax[i]) fmax[i]=FF[j+1];
	    if (FF[j+1]<fmin[i]) fmin[i]=FF[j+1];
	  }
	}
	break;
    }		/* switch */
  }		/* for nphase */
#undef fmin
#undef fmax
}


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

Entry(Int2) ss_ProcDefault(Int2 Ind, Vector x, Vector v1, CharPtr *msg) {
  if (Ind == -1) {  /* zero of user-defined function */
    stagenData->BifDataOutLen=0;     
    return 0;
  }
  /* Assert: v1==NULL, msg==NULL, Ind==0 */
  MemCopy(x1,x,sizeof(VAR)*(nmesh*nphase+1+nappend));     /* for the last point */
  if (staData->eigcomp) {	/* eigenvalues */
    if (JacRHS(x, A0)) return (1);      
    staFunc->EigenValues(A0,Re,Im);
  }
  if (staData->MinMaxComp) {	/* Min&Max */
    MinMax(x);
  }
  if (staData->L2NormComp) {	/* L2 norm */
    STF[0]=L2_Norm(x);
  }
  return 0;
}

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

Local(Char) buf[80];

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

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

  FloatHi reseq;
  if (ss_JacDefStadyState(x, A)) return (1);
  B = rdFunc->AppendMatrixVector(A, v, B);
  rdFunc->GetMatrixData(B,&reseq,res);
  test[0]=*res;
  return 0;
}

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

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

  stagenData->BifDataOutLen=0;

  if (Ind != 0) {
    sprintf(buf,"Branching (?)");
    *msg=buf;
    return 0;
  } else {
    sprintf(buf,"Branching");
    *msg=buf;
    return 0;

    q=staFunc->CreateVector(nphase+1+nappend);
    p=staFunc->CreateVector(nphase+1+nappend);
    phi=staFunc->CreateVector(nphase+nappend);
    if (nappend) psi=staFunc->CreateVector(nappend);
    D=rdFunc->CreateMatrix(nphase+1+nappend,nphase+1+nappend);
    DT=rdFunc->CreateMatrix(nphase+1+nappend,nphase+1+nappend);
    v2=staFunc->CreateVector(nphase+1+nappend);
    v=staFunc->CreateVector(nphase+1+nappend);
    w=staFunc->CreateVector(nphase+1+nappend);
    x1=staFunc->CreateVector(nphase+1+nappend);
    F=staFunc->CreateVector(nphase+nappend);
    F1=staFunc->CreateVector(nappend);

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

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

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

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

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

/* Limit point (LP): p[active[0]]-component of the tangent vector to the curve or det */
Entry(Int2) ss_TestFold(Vector x, Vector v, FloatHiPtr res) {
int i;
FloatHi reseq;
  *res=1.0;
  if (nappend==0) {
    *res=v[nmesh*nphase];
  } else {
    if (ss_JacDefStadyState(x, A)) return (1);
    B = rdFunc->AppendMatrixVector(A, v, B);
    rdFunc->GetMatrixData(B,res,NULL);
  }
  test[1]=*res;
  return 0;
}

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

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

  stagenData->BifDataOutLen=0;

  if (Ind == 0) {

sprintf(buf,"Limit point");
*msg=buf;
return 0;

    V=staFunc->CreateVector(nphase);
    W=staFunc->CreateVector(nphase);
    A=rdFunc->CreateMatrix(nphase,nphase);
    B=rdFunc->CreateMatrix(nphase,nphase);
    J=rdFunc->CreateMatrix(nphase,nphase+1+nappend);
    x1=staFunc->CreateVector(nphase+1+nappend);
    F=staFunc->CreateVector(nphase);

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

/*  provide bifurcation data for lp->lp starter */
    staFunc->CopyVectorElements(v,0,V,0,nphase);
    staFunc->NormalizeVector(V);
    MemCopy(BifVector,V,sizeof(FloatHi)*nphase);
    stagenData->BifDataOutLen=sizeof(FloatHi)*nphase;

/*  compute adjoint null-vector */
    rdFunc->TransposeMatrix(A,B);
    if (staFunc->SngSolve(B,eps,W) != nphase-1) {
      sprintf(buf,"Limit point: a=?");
    } else {
      staFunc->NormalizeVectorRelative(W,V);

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

      } else {
        for (i=0; i<=nappend; i++) v[nphase+i]=0.0;
        DefRHS(x,F);
        f0=staFunc->ScalProd(W,F);
        staFunc->AddScaledVector(x,v,_dx,x1);
        DefRHS(x1,F);
        f1=staFunc->ScalProd(W,F);
        staFunc->AddScaledVector(x,v,-_dx,x1);
        DefRHS(x1,F);
        f2=staFunc->ScalProd(W,F);
        sprintf(buf,"Limit point: a=%g",(f1-2*f0+f2)/_dx/_dx);
      }
    }
    rdFunc->DestroyMatrix(J);
    rdFunc->DestroyMatrix(B);
    rdFunc->DestroyMatrix(A);
    staFunc->DestroyVector(W);
    staFunc->DestroyVector(W);
    staFunc->DestroyVector(x1);
    staFunc->DestroyVector(F);
  } else {
    sprintf(buf,"Limit point (?)");
  }
  *msg=buf;
  return 0;
}
/* Hopf (H): Determinant of the bialternate product 2A*I */
Entry(Int2) ss_TestHopf(Vector x, Vector v, FloatHiPtr res) {
  int i,j,p,q,r,s;

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

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

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

int i,j,imin,jmin;
double d,s,r,h,beta,gamma,phi,lambda,omega,eps=0.01;
rdMatrix AA,AT,BB;
Vector V,Qr,Qi,Pr,Pi;
   
  AA=rdFunc->CreateMatrix(nphase,nphase);
  AT=rdFunc->CreateMatrix(nphase,nphase);
  BB=rdFunc->CreateMatrix(2*nphase,2*nphase);
  V=staFunc->CreateVector(2*nphase);
  Qr=staFunc->CreateVector(nphase);
  Qi=staFunc->CreateVector(nphase);
  Pr=staFunc->CreateVector(nphase);
  Pi=staFunc->CreateVector(nphase);

  stagenData->BifDataOutLen=0;
  if (Ind == 0) {
    ss_JacDefStadyState(x,A);
    rdFunc->CopyMatrixBlock(A,0,0,AA,0,0,nphase,nphase);
    rdFunc->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 */
      rdFunc->CopyMatrixBlock(A,0,0,AA,0,0,nphase,nphase);
      for (i=0; i<nphase; i++) AA[i][i]-=lambda;
      if (staFunc->SngSolve(AA,eps,Qr) != nphase-1) {
        sprintf(buf,"Neutral saddle (?)");
        goto final;
      } 
      rdFunc->CopyMatrixBlock(A,0,0,AA,0,0,nphase,nphase);
      for (i=0; i<nphase; i++) AA[i][i]+=lambda;
      if (staFunc->SngSolve(AA,eps,Qi) != nphase-1) {
        sprintf(buf,"Neutral saddle (?)");
        goto final;
      } 
      staFunc->NormalizeVector(Qr);   
      staFunc->NormalizeVectorRelative(Qi,Qr);    /* ensures smallest possible angle */
      staFunc->NormalizeVector(Qi);

      /* provide bifurcation data for h->h starter */
      staFunc->AddScaledVector(Qr,Qi,1.0,Pr);
      staFunc->NormalizeVector(Pr);
      MemCopy(BifVector,Pr,sizeof(FloatHi)*nphase);
      staFunc->AddScaledVector(Qr,Qi,-1.0,Pr);
      staFunc->NormalizeVector(Pr);
      MemCopy(BifVector+nphase,Pr,sizeof(FloatHi)*nphase);
      BifVector[2*nphase]=-lambda*lambda;
      stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+1);      
      sprintf(buf,"Neutral saddle: lambda=%g",lambda);
    } else {
/* Hopf */
      lambda=Re[imin];
      omega=fabs(Im[imin]);
      /* compute real and imaginary part of the eigenvector */
      rdFunc->CopyMatrixBlock(A,0,0,BB,0,0,nphase,nphase);
      rdFunc->CopyMatrixBlock(A,0,0,BB,nphase,nphase,nphase,nphase);
      for (i=0; i<nphase; i++) {
        BB[i][i+nphase]=omega;
        BB[i][i]-=lambda;
        BB[i+nphase][i]=-omega;
        BB[i+nphase][i+nphase]-=lambda;
      }
      if (staFunc->SngSolve(BB,eps,V) != 2*nphase-2) {
        sprintf(buf,"Hopf (?)");
        goto final;
      }
      staFunc->CopyVectorElements(V,0,Qr,0,nphase);
      staFunc->CopyVectorElements(V,nphase,Qi,0,nphase);
      /* normalize real and imaginary parts <Qr,Qr>=1, <Qr,Qi>=0 */
        d=staFunc->ScalProd(Qr,Qr);
        s=staFunc->ScalProd(Qi,Qi);
        r=staFunc->ScalProd(Qr,Qi);
	gamma=1.0/sqrt(d+s);
        if (r==0) {
          staFunc->MultiplyVectorScalar(V,gamma,V);
        } else {    
          beta=0.5*(d-s)/r;
          if (beta == 0.0) {
            phi=0.5*sqrt(2.0);
            for (i=0; i<nphase; i++) {
              V[i]=gamma*phi*(Qr[i]-Qi[i]);
              V[i+nphase]=gamma*phi*(Qr[i]+Qi[i]);
            }
          } else {
            phi=0.5*atan(-1.0/beta);
            for (i=0; i<nphase; i++) {
              V[i]=gamma*(Qr[i]*cos(phi)-Qi[i]*sin(phi));
              V[i+nphase]=gamma*(Qr[i]*sin(phi)+Qi[i]*cos(phi));
            }
          }
        }
        staFunc->CopyVectorElements(V,0,Qr,0,nphase);
        staFunc->CopyVectorElements(V,nphase,Qi,0,nphase);
        d=staFunc->ScalProd(Qr,Qr);
        gamma=1/sqrt(d);
        staFunc->MultiplyVectorScalar(V,gamma,V);
        staFunc->CopyVectorElements(V,0,Qr,0,nphase);
        staFunc->CopyVectorElements(V,nphase,Qi,0,nphase);

      /* provide bifurcation data for h->h starter */
       MemCopy(BifVector,V,sizeof(FloatHi)*2*nphase);
       BifVector[2*nphase]=omega*omega;
       stagenData->BifDataOutLen=sizeof(FloatHi)*(2*nphase+1);

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

      /* normalize the adjoint vector <Pr,Qr>+<Pi,Qi>=1, <Pr,Qi>-<Pi,Qr>=0 */
       d=staFunc->ScalProd(Pr,Qr);
       s=staFunc->ScalProd(Pi,Qi);
       r=staFunc->ScalProd(Pr,Qi);
       h=staFunc->ScalProd(Pi,Qr);
       phi=atan((r-h)/(d+s));

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

       sprintf(buf,"Hopf: omega=%g",omega);
    }
  } else {
    sprintf(buf,"Hopf or neutral saddle (?)");
  }
  *msg=buf;
final:
  AA=rdFunc->DestroyMatrix(AA);
  AT=rdFunc->DestroyMatrix(AT);
  BB=rdFunc->DestroyMatrix(BB);
  V=staFunc->DestroyVector(V);
  Qr=staFunc->DestroyVector(Qr);
  Qi=staFunc->DestroyVector(Qi);
  Pr=staFunc->DestroyVector(Pr);
  Pi=staFunc->DestroyVector(Pi);
  return 0;
}