# Copyright (C) 2001 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Header$
## $Log$
## Revision 1.358  2002/08/07 11:04:40  geraint
## Now updates timestamp on lbl.txt after processing in make. This eliminates the delay where MTT spends a long period of time uselessly recursing through the model on subsequent passes if the label file is older than the bondgraph.
##
## Revision 1.357  2002/06/28 10:13:40  geraint
## Includes fix_c.r in ese2rdae and def2write_r to eliminate occurrances of x**y.
##
## Revision 1.356  2002/06/17 13:14:50  geraint
## Fixed some rules to reduce number of implicit rule searches.
##
## Revision 1.355  2002/06/10 23:18:44  geraint

$1_ode2odes_implicit_%.stamp: $1_cseo.% $1_csex.% $1_smxa.% $1_smxax.%
	touch \$@
$1_ode2odes_implicit.o: $1_cseo.o $1_csex.o $1_smxa.o $1_smxax.o mtt_implicit.o
	@echo "Creating $1_ode2odes_implicit.o"
	ar -cr \$@ \$^

$1_ode2odes_dassl_oct.stamp: $1_ode.oct $1_odeo.oct mtt_dassl.oct
	touch \$@
$1_ode2odes_dassl_o.stamp: $1_ode.o $1_odeo.o mtt_dassl.o
	touch \$@
$1_ode2odes_dassl_%.stamp: $1_ode.% $1_odeo.%
	touch \$@
$1_ode2odes_dassl.o: $1_ode.o $1_odeo.o mtt_dassl.o
	@echo "Creating \$@"
	ar -cr \$@ \$^

$1_ode2odes_${algebraic_solver}.cc:	mtt_Solver.cc mtt_AlgebraicSolver.cc mtt_${algebraic_solver}.hh mtt_${algebraic_solver}.cc

#! /bin/sh

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

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.82  2002/07/24 14:00:12  geraint
## Corrected arguments passed to mtt_write when dumping data (sigint).
##
## Revision 1.81  2002/07/11 13:00:23  geraint
## Declared more function arguments to be "const" - improves compiler optimisation.
##
## Revision 1.80  2002/05/22 09:35:49  geraint
## Added insertor variable to stop sh-mode font-lock from getting hopelessly confused by embedded C++.
##
## Revision 1.79  2002/05/20 13:42:31  gawthrop

    f = feval ("mtt_dassl", args, 1);
    return f(0).${vector_value} ();
#endif // ((CODEGENTARGET == STANDALONE) || (CODEGENTARGET == MATLABMEX))
}

#if ((CODEGENTARGET == STANDALONE) || (CODEGENTARGET == MATLABMEX))
ColumnVector
Fmtt_residual (const ColumnVector &X,
	      const ColumnVector &DX,
	      double t,
	      int &ires)
{
#elif (CODEGENTARGET == OCTAVEDLD)
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 // ((CODEGENTARGET == STANDALONE) || (CODEGENTARGET == MATLABMEX))

    row = 0;
    fcputime.close();
    }
}

void
${sys}_ode2odes (const ColumnVector &state0,
		 const ColumnVector &par,
		 const ColumnVector &simpar)
{
  static double first, dt, last, stepfactor;
  first		= simpar (0);
  last		= simpar (1);
  dt		= simpar (2);
  stepfactor	= simpar (3);

  for (register int i = 0; i < MTTNX; i++)
    {
      x (i) = state0 (i);
    }

  for (register int j = 0, i = 1; i <= ilast; i++)
    {
      u	= mtt_input (x, y, t, par);
      y	= mtt_output (x, u, t, par);
      if (0 == j)
	{
           mtt_write (t, x, y, first, nrows);
	}
EOF
case "$method" in
    "rk4")

          t2 = t + ddt;

        ColumnVector
          x1 (x),
          x2 (x),
          x3 (x);

        k1 = ddt * mtt_rate (x , u, t , par); x1 += k1 * 0.5;
        k2 = ddt * mtt_rate (x1, u, t1, par); x2 += k2 * 0.5;
        k3 = ddt * mtt_rate (x2, u, t1, par); x3 += k3;
        k4 = ddt * mtt_rate (x3, u, t2, par);
        dx = (k1 + 2.0 * (k2 + k3) + k4) / (6.0 * ddt);
      }
EOF
    ;;
    "dassl")
    ;;
    "implicit")
	cat << EOF >> $filename
      dx = mtt_rate (x, u, t, par);
      AA = mtt_smxa (x, u, ddt, par);
      AAx = mtt_smxax (x, u, ddt, par);
EOF
    ;;
    "euler" | *)
	cat << EOF >> $filename
      dx = mtt_rate (x, u, t, par);
EOF
    ;;
esac

## Common stuff
cat <<EOF >> $filename
      open_switches = mtt_logic (x, u, t, par);
      x =  $algorithm; 
      t += ddt;
      j++;
      j = (j == static_cast<int> (stepfactor)) ? 0 : j;
    }

}

#if (CODEGENTARGET == STANDALONE)
  return Fmtt_residual (X, DX, t, ires);
#elif (CODEGENTARGET == OCTAVEDLD)
  static octave_value_list args, f;
  args(0) = octave_value (X);
  args(1) = octave_value (DX);
  args(2) = octave_value (t);
  args(3) = octave_value (static_cast<double>(ires));
  f = feval ("mtt_residual", args, 1);
  return f(0).VECTOR_VALUE ();
#endif // (CODEGENTARGET == STANDALONE)
}


#if (CODEGENTARGET == STANDALONE)