Local(FloatHiPtr) J;			/* current set of Jacobian vals */
Local(FloatHiPtr) H;                   	/* current set of Hessian vals */
Local(FloatHiPtr) D3;                   /* current set of third-order dervals */
Local(Int2) Der2num;			/* number of 2nd order derivatives */
Local(Int2) Der3num;			/* number of 3nd order derivatives */

Local(void) CopyToX(Vector x);
Local(void) CopyToP(Vector x);

/**************************************************************************/
/* Assign values to the Jacobian matrix of RHS using symbolic derivatives */
Local(void) SymJac(Vector x, Matrix Jac) {
int i,j,nactive;
  CopyToX(x);
  CopyToP(x);
  J=stagenData->Der1(NULL,NULL,NULL,X,P,&zero);
  /*                 U    U_S  U_SS S P T */
  nactive=staFunc->GetVectorDim(x)-nphase;
  for (i=0; i<nphase; i++) {
    MemCopy(Jac[i],J+(nphase+npar)*i,sizeof(FloatHi)*nphase);
    for (j=0; j<nactive; j++) 
         Jac[i][nphase+j]=*(J+(nphase+npar)*i+nphase+active[j]);
  }
}

/* Returns second order derivatives */ 
Local(FloatHi) HessianElem(int i, int j, int k) {

#define HessElem(i,j,k) (*(H+Der2num*(i)+ind2_(j,k)))    
#define MixedHessElem(i,j,k) HessElem(i,j,nphase+active[k])  
  
  if (j<=k) 
    if (k<nphase) 
      return HessElem(i,j,k);
    else 
      return MixedHessElem(i,j,k-nphase);
  else 
    if (j<nphase) 
      return HessElem(i,k,j);
    else
      return MixedHessElem(i,k,j-nphase);

#undef HessElem
#undef MixedHessElem
}

#define Der3Elem(i,j,k,l) (*(D3+Der3num*(i)+ind3_(j,k,l)))
#define Swap(a,b) {t=a; a=b; b=t;}

/* Returns third order derivatives */
Local(FloatHi) D3Elem(int i, int j, int k, int l) {
  int t;
  /* ensure j<=k<=l */
  if (j>k) Swap(j,k);
  if (j>l) Swap(j,l);
  if (k>l) Swap(k,l);
  return Der3Elem(i,j,k,l);
}