/*
   1D Reaction-Diffusion
   rd__o     
   Any distribution-->Orbit
*/

#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__o.h"
#include "rd_int.h"  /* integrator'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 {		/* used 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__o_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) x2=NULL;
Local(VAR) t;				/* time */
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;		/* NOT rdMatrix! */
Local(Vector) Re=NULL,Im=NULL;		/* Eigenvalues */
Local(VAR) L2Norm;
Local(Vector) STF=NULL;
Local(VAR) test[1];	/* dummy 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 o_SortCmp(const void *f, const void *s) {
#else
Entry(int) o_SortCmp(const void *f, const void *s) {
#endif
  if (*(VARPTR )f==*(VARPTR )s) return 0;
  return (*(VARPTR )f<*(VARPTR )s) ? -1 : 1;
}

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

/**************/
/* Adaptation */
Entry(Int2) o_Adapter(Vector U, Vector XM) {
Int2 i,j,Ret;
  Points (XM, 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,XM);   
  staFunc->CopyVector(XN,X); 
  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__o_DataPtr)stagenptr->StaPtr;
  /* For access to standard linear algebra */
  staFunc=(StdLinAlgDataPtr)stagenptr->StdLinAlg;
  staUtil=(UtilitiesDataPtr)stagenptr->Utilities;
  /* For access to RD specific linear 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),o_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 (!stagenData->Der2 && staData->dx <= 0) {
    staUtil->ErrMessage("Analytical Hessian matrix is undefined.\nRecompile system or set positive 'Increment'");
    /* active=MemFree(active); */
    return 3;
  }
  if (!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 for eigenvalues */
    x0=staFunc->CreateVector(nmesh*nphase+1+nappend);
    x1=staFunc->CreateVector(nmesh*nphase+1+nappend);
    x2=staFunc->CreateVector(nmesh);
    /* 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 */
    JacW=staFunc->CreateMatrix(nphase,nphase);
    JacV=staFunc->CreateMatrix(nphase,nphase);
  }

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

  ((IntegrDataPtr)stagenptr->GenPtr)->ActiveUserTest=staData->ActiveUserTest;
  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);
      x2=staFunc->DestroyVector(x2);
      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);
      JacW=staFunc->DestroyMatrix(JacW);
      JacV=staFunc->DestroyMatrix(JacV);
    }
  /* Clean up RD algebra package */
  rdFunc->Term();
}

/******************************/
/* Any point -> Orbit Starter */

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

Entry(Int2) o_Starter(StagenDataPtr stagenptr) {
  Int2 i,Ret;
  
  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 */
  x1[nmesh*nphase]=*stagenptr->ClassDim[CL_TIME] ? *staData->T : *staData->Tdef;  
  MemCopy(x1,staData->U,sizeof(FloatHi)*nmesh*nphase);
  MemCopy(x2,staData->X,sizeof(FloatHi)*nmesh);
  /* Initialize continuation */
  ((IntegrDataPtr)stagenptr->GenPtr)->X0=x1;
  ((IntegrDataPtr)stagenptr->GenPtr)->X1=x2;
  /* That's all */
  return 0;
}

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

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
#define _maxnewt  staData->maxnewt
#define _epsx     staData->epsx
#define _epsf     staData->epsf

/************************/
/* 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) o_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]=(VAR *)vpoint+nmesh*nphase;
      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) o_DefUser (Vector x, Vector y) {
Int2 i,j;
  for (i=j=0; i<nUTest; i++) 
    if (ActiveUTestFuns[i]==UDF_APPEND) {
      (stagenData->udFunc)(U,i,&y[j]);
      j++;
    }
  return 0;
}

/***************************/
/* Calculate 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) o_DefStadyState (Vector x, Vector y) {
Int2 i;

  DefRHS(x,y);
  if (nappend) o_DefUser(x,y+nmesh*nphase);
  return 0;
}

#include "rd_symd.c"

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

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

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

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

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

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

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

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

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

/* f1(.,.,p,.) */
Entry(Int2) o_BC1p(VARPTR  p, Vector f) {
  CopyToPa(p);
  if (((BCPtr)(stagenData->AuxFunc[BC1_I]))(U+(nmesh-1)*nphase,U_X,P,&p[0],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)o_BC0u, nphase, u+0, _dx, JacW)) return 1;
    if (staFunc->DerA((DiffFun)o_BC0u_x, nphase, U_X, _dx, ma)) return 1;
  } else {
    dm=((BCDPtr)(stagenData->AuxFunc[BC0D1_I]))(u+0,U_X,P,&u[0]);
    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);
  /* 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)o_RHSu, nphase, U_0, _dx, JacV)) return 1;
      if (staFunc->DerA((DiffFun)o_RHSu_x, nphase, U_X, _dx, JacW)) return 1;
      if (staFunc->DerA((DiffFun)o_RHSu_xx, nphase, U_XX, _dx, ma)) return 1;
    } else {
      dm=stagenData->Der1(U_0,U_X,U_XX,&Xi,P,&u[0]);
      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);
  }
  /* 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)o_BC1u, nphase, u+N*nphase, _dx, JacW)) return 1;
    if (staFunc->DerA((DiffFun)o_BC1u_x, nphase, U_X, _dx, mc)) return 1;
  } else {
    dm=((BCDPtr)(stagenData->AuxFunc[BC1D1_I]))(u+N*nphase,U_X,P,&u[0]);
    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);
  return 0;
#undef N
}
/*****************************************************/
/* Jacobian of defining conditions for stady states  */
Entry(Int2) o_JacDefStadyState (Vector u, rdMatrix Jac) {
  if (JacRHS(u, Jac)) return 1;
  if (nappend) {
    for (ipoint=0; ipoint<nmesh; ipoint++) {
      staFunc->DerA(o_DefUser, nappend, U_0, _dx, rdFunc->GetMatrixP(Jac,ipoint));
    }
    staFunc->DerA(o_DefUser, nappend, U_0, _dx, rdFunc->GetMatrixFi(Jac));
  }
  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)
  printf("MinMax entered\n");
  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
  printf("MinMax left\n");
}


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

Entry(Int2) o_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];

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

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