#! /bin/sh

     ###################################### 
     ##### Model Transformation Tools #####
    ######################################

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.64  2001/08/08 02:15:00  geraint
## Rationalisation of solver code, beginning with algebraic solvers.
##
## Revision 1.63  2001/08/07 04:39:24  geraint
## Consolidated dassl and residual functions.
##
## Revision 1.62  2001/08/01 22:14:32  geraint
## Bug fix for dassl.
##
## Revision 1.61  2001/08/01 04:06:07  geraint

#endif
}

#ifdef STANDALONE
ColumnVector
Fmtt_residual (const ColumnVector &X,
	      const ColumnVector &DX,
	      double t)
	      double t,
	      int &ires)
{
#else // !STANDALONE
DEFUN_DLD (mtt_residual, args, , "mtt_residual")
{
    static ColumnVector	X  (MTTNX+MTTNYZ);
    static ColumnVector	DX (MTTNX+MTTNYZ);
    static double      	t;
    static int		&ires;

    X  = args(0).${vector_value} ();
    DX = args(1).${vector_value} ();
    t  = args(2).double_value ();
    ires = static_cast<int>(args(3).double_value ());
#endif // STANDALONE

    static ColumnVector residual (MTTNX+MTTNYZ);
    static ColumnVector U (MTTNU+MTTNYZ);
    static ColumnVector u (MTTNU);
    static ColumnVector y (MTTNY,0.0);
    static ColumnVector par (MTTNPAR);

    yz = feval ("${sys}_ae", new_args, 1)(0).${vector_value} ();
#endif

    for (register int i = 0; i < MTTNX; i++)
      residual (i) = dx(i) - DX(i);
    
    if (MTTNYZ > 0)
      {
      residual.insert (yz,MTTNX);
	residual.insert (yz,MTTNX);

	// XXX:
	// DASSL needs residual to be dependent on Ui and Uidot
	// mtt_dassl always sets the initial Ui to zero, so
	// Ui - h*Uidot should equal zero BUT, we don't know h
	static double t_old;
	double step;
	if (t != t_old)
	  {
	    step = t - t_old;
	    t = t_old;
	  }
	else
	  step = t;	
	for (register int i = MTTNX; i < MTTNX+MTTNYZ; i++)
	  residual(i) += X(i) - DX(i)*step;
      }

#ifdef STANDALONE
    return residual;
#else // !STANDALONE
    return octave_value (residual);
#endif // STANDALONE
}