# Copyright (C) 2001 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Header$
## $Log$
## Revision 1.309.2.1  2001/05/04 04:07:24  geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## Revision 1.309  2001/04/28 03:15:03  geraint
## Fixed comment (interfered with "mtt help representations").
##
## Revision 1.308  2001/04/25 22:17:45  geraint
## Fixed icd.txt2 dependency.
##
## Revision 1.307  2001/04/15 21:15:41  geraint
## Added interface definition rep: _ICD.(txt|c|cc|m).
##
## Revision 1.305  2001/04/11 09:44:26  gawthrop

# Default not verbose
verbose=' -s'

# Default integration method
integration_method=implicit;

# Default algebraic equation solver
algebraic_solver=Reduce_Solver;

# Default no info
info_switch=''

# Default use m, not oct files
m='m';

# Default use ps files

			integration_method=rk4;
			mtt_switches="$mtt_switches rk4";
			;;
		    *)
			echo $1 is an unknown integration method - use euler, rk4 or implicit;
                        exit;;
		esac;;
	-ae )
                mtt_switches="$mtt_switches $1";
		case $2 in
		    fsolve | hybrd)
			mtt_switches="$mtt_switches $2";
			algebraic_solver=Hybrd_Solver;
			shift;
			;;
		    hooke)
			mtt_switches="$mtt_switches $2";
			algebraic_solver=HJ_Solver;
			shift;
			;;
		    reduce)
			mtt_switches="$mtt_switches $2";
			algebraic_solver=Reduce_Solver;			
			Solve='-A';
			shift;
			;;
		    *)
			echo $1 is an unknown solver - use hybrd, hooke or reduce;
			exit;;
	        esac;;
	-s )
                mtt_switches="$mtt_switches $1";
		sensitivity=sensitivity ;;
	-ss )
                mtt_switches="$mtt_switches $1";
		steadystate_computation=yes ;;
	-d )

    echo '         -I  prints more information'
    echo '         -abg start at abg.m representation'
    echo '         -c  c-code generation'
    echo '         -d  <dir>  use directory <dir>'
    echo '         -dc Maximise derivative (not integral) causality'
    echo '         -dc Maximise derivative (not integral) causality'
    echo '         -i <implicit|euler|rk4>   Use implicit, euler or rk4 integration'
    echo '         -ae <reduce|hybrd|hooke>   Solve algebraic equations with Reduce, hybrd (fsolve) or "Hooke and Jeeves"'
    echo '         -o ode is same as dae'
    echo '         -oct use oct files in place of m files where appropriate'
    echo '         -opt optimise code generation'
    echo '         -p  print environment variables'
    echo '         -partition partition hierachical system'
    echo '         -pdf generate pdf in place of ps'
    echo '         -r  reset time stamp on representation'

  rm -f *_csex.cc *_numpar.cc *_smxa.cc *_state.cc 
  rm -f *_ode.cc *_odeo.cc
  rm -f *_logic.m *_logic.cc *_logic.oct
  rm -f *_state.m *_state.cc *_state.oct
  rm -f *_ode2odes.* *.dat2 MTT.core
  rm -f *_modpar.txt *_modpar.r
  rm -f *_ICD.txt *_ICD.c *_ICD.cc *_ICD.m
  rm -f *_ae.r *_ae.m *_ae.cc *_ae.oct
  rm -fR *_rep MTT_work
  exit
fi

# Clean up named system
if [ "$2" = "Clean" ] && [ "$3" = "" ]; then
  echo 'Removing all generated files for system ' $1

  rm -f $1_csex.cc $1_numpar.cc $1_smxa.cc $1_state.cc 
  rm -f $1_ode.cc $1_odeo.cc
  rm -f $1_logic.m $1_logic.cc $1_logic.oct 
  rm -f $1_state.m $1_state.cc $1_state.oct
  rm -f $1_ode2odes.* $1.dat2
  rm -f $1_modpar.txt $1_modpar.r
  rm -f $1_ICD.txt $1_ICD.c $1_ICD.cc $1_ICD.m
  rm -f $1_ae.r $1_ae.m $1_ae.cc $1_ae.oct
  rm -fR $1_rep MTT_work
  exit
fi

if [ "$2" = "rep" ]; then
     documenttype=article
      # See if we are making a book -- ie representation rep on a directory

.PRECIOUS: %.cc # Don't let mtt delete them
$1_%.cc:  $1_%.m
	mtt_m2cc.sh  $1 \$* cat 

mtt_%.cc: ${MTT_LIB}/cc/mtt_%.cc
	cp ${MTT_LIB}/cc/mtt_\$*.cc mtt_\$*.cc

mtt_%.hh: ${MTT_LIB}/cc/mtt_%.hh
	cp ${MTT_LIB}/cc/mtt_\$*.hh mtt_\$*.hh

## .o files
.PRECIOUS: $1_%.o
$1_%.o: $1_%.cc $1_def.h $1_sympar.h $1_cr.h
	echo Compiling $1_\$*.cc
	${MTT_CXX} ${MTT_CXXFLAGS} ${MTT_CXXINCS} -c $< -DSTANDALONE

.PRECIOUS: mtt_%.o

else
	mtt_update $1 logic
endif

# Dummy target
FORCE:

$1_ode2odes_common.m    : $1_ae.m   $1_input.m   $1_logic.m   $1_numpar.m   $1_simpar.m   $1_state.m
	@echo > /dev/null
$1_ode2odes_common.cc   : $1_ae.cc  $1_input.cc  $1_logic.cc  $1_numpar.cc  $1_simpar.cc  $1_state.cc
	@echo > /dev/null
$1_ode2odes_common.o    : $1_ae.o   $1_input.o   $1_logic.o   $1_numpar.o   $1_simpar.o   $1_state.o
	@echo "Creating $1_ode2odes_common.o"
	ar -cr \$@ \$^
$1_ode2odes_common.oct  : $1_ae.oct $1_input.oct $1_logic.oct $1_numpar.oct $1_simpar.oct $1_state.oct
	@echo > /dev/null

$1_ode2odes_euler.m  $1_ode2odes_rk4.m    : $1_ode.m     $1_odeo.m
	@echo > /dev/null
$1_ode2odes_euler.cc $1_ode2odes_rk4.cc   : $1_ode.cc    $1_odeo.cc
	@echo > /dev/null
$1_ode2odes_euler.o  $1_ode2odes_rk4.o    : $1_ode.o     $1_odeo.o   mtt_euler.o

	@echo > /dev/null
$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_implicit.oct: $1_cseo.oct  $1_csex.oct $1_smxa.oct $1_smxax.oct mtt_implicit.oct
	@echo > /dev/null

mtt_Solver.cc:	mtt_Solver.hh
$1_ode2odes_${algebraic_solver}.cc:	mtt_Solver.cc mtt_${algebraic_solver}.hh mtt_${algebraic_solver}.cc
$1_ode2odes_${algebraic_solver}.o:	mtt_Solver.o mtt_${algebraic_solver}.o
	@echo "Creating $1_ode2odes_${algebraic_solver}.o"
	ar -cr \$@ \$^

#SUMMARY numpar	numerical parameter declaration (m) 
$1_numpar.m:  $1_numpar.txt $1_sympars.txt
	mtt_txt2m $1 numpar

#SUMMARY numpar	numerical parameter declaration (c) 
#SUMMARY numpar	numerical parameter declaration (view) 

$1_dae.m: $1_def.r $1_dae.r  $1_sympars.txt
	dae_r2m $1; matlab_tidy $1_dae.m; matlab_tidy $1_daeo.m
$1_dae.c: $1_def.r $1_dae.r  $1_sympar.r
	dae_r2c $1; c_tidy $1_dae.c
$1_dae.tex: $1_dae.r $1_simp.r
	dae_r2tex $partition $1; latex_tidy $1_dae.tex

#SUMMARY ae	algebraic equations - unknown inputs (r)
#SUMMARY ae	algebraic equations - unknown inputs (m)
#SUMMARY ae	algebraic equations - unknown inputs (cc)

$1_ae.r: $1_cse.r
	touch $1_ae.r
$1_ae.m: $1_ae.r
	mtt_r2m $1 ae

#SUMMARY cse	constrained-state equations (r)
#SUMMARY cse*	constrained-state equations (m)
#SUMMARY cse*	constrained-state equations (oct)
#SUMMARY cse	constrained-state equations (tex)
#SUMMARY cse*	constrained-state equations (view)
#SUMMARY cse	constrained-state equations (ps)

                $1_smxa.m $1_smxax.m\
                $1_simpar.m $1_numpar.m $1_state.m $1_input.m \
                $1_csex.m $1_cseo.m  $1_logic.m
ifeq ($using_oct,yes)
	touch $1_ode2odes.m # Create a dummy which wont' be used
	mtt $mtt_switches -q -u $1 ode2odes oct
else
	make_ode2odes $1 m $integration_method $algebraic_solver
endif
endif
ifneq ($integration_method,implicit)
$1_ode2odes.m : $1_def.r $1_sympars.txt\
		$1_simpar.m $1_numpar.m $1_state.m $1_input.m \
		$1_ode.m $1_odeo.m  $1_logic.m
ifeq ($using_oct,yes)
	echo "*** Warning: Shouldn't be here! Creating dummy $1_ode2odes.m"
	touch $1_ode2odes.m # Create a dummy which wont' be used
	mtt $mtt_switches -q -u $1 ode2odes oct
else
	make_ode2odes $1 m $integration_method $algebraic_solver
endif

endif

#SUMMARY ode2odes Simulation function (m)
#SUMMARY ode2odes Simulation function (cc)
#SUMMARY ode2odes Simulation function (oct)
#SUMMARY ode2odes Simulation function (exe)
$1_ode2odes.exe: $1_def.h $1_sympar.h\
		 $1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o $1_ode2odes_${algebraic_solver}.o
	echo Creating $1_ode2odes.exe
	${MTT_CXX} ${MTT_CXXFLAGS} -o $1_ode2odes.exe\
	$1_ode2odes.o $1_ode2odes_common.o $1_ode2odes_${integration_method}.o $1_ode2odes_${algebraic_solver}.o\
	${MTT_LDFLAGS} ${MTT_CXXLIBS}

$1_ode2odes.o: $1_ode2odes.cc $1_ode2odes_common.o $1_ode2odes_${integration_method}.o $1_ode2odes_${algebraic_solver}.o
	echo Creating $1_ode2odes.o
	${MTT_CXX} ${MTT_CXXFLAGS} ${MTT_CXXINCS} -c $1_ode2odes.cc -DSTANDALONE

$1_ode2odes.oct: $1_ode2odes.cc $1_ode2odes_common.oct $1_ode2odes_${integration_method}.oct
	touch $1_ode2odes.m
	echo Creating $1_ode2odes.oct
	$MKOCTFILE $1_ode2odes.cc

$1_ode2odes.cc: $1_def.r $1_sympars.txt\
		$1_ode2odes_common.m $1_ode2odes_common.cc\
		$1_ode2odes_${integration_method}.m $1_ode2odes_${integration_method}.cc mtt_Solver.cc mtt_${algebraic_solver}.cc mtt_${algebraic_solver}.hh
	touch $1_ode2odes.m
	make_ode2odes $1 cc $integration_method $algebraic_solver

#Conversion of m to p to c
#SUMMARY ode	ordinary differential equations (c)
#SUMMARY ode	ordinary differential equations (p)
#SUMMARY state	state declaration (c) 
#SUMMARY state	state declaration (p) 
$1_simpar.p : $1_def.r $1_simpar.m

# Copyright (c) P.J.Gawthrop, 1991, 1994.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.1  2000/12/28 12:21:31  peterg
## Put under RCS
##
## Revision 1.2  1997/01/06 21:17:10  peterg
## Removed:  OFF Exp; OFF GCD;
##
## Revision 1.1  1996/08/25 10:05:45  peter
## Initial revision
##
###############################################################

in "$1_def.r";

%Read the substitution file
in "$1_subs.r";

%Read the constrained-state equations file
in "$1_cse.r";

%Read the algebraic equations file
in "$1_ae.r";

IF MTTNx>0 THEN
 IF MTTNz>0 THEN
   MTTdXX := MTTE^(-1)*MTTEdX
 ELSE
   MTTdXX := MTTEdX;

# Copyright (c) P.J.Gawthrop 1991, 1992, 1994.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.15.2.2  2001/05/05 20:50:16  geraint
## Fixed errors when MTTNx=0.
##
## Revision 1.15.2.1  2001/05/04 04:07:24  geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## Revision 1.15  2001/03/19 02:28:52  geraint
## Branch merge: merging-ode2odes-exe back to MAIN.
##
## Revision 1.14.2.1  2001/03/19 00:29:08  geraint
## Parse switches (-A) before calling def2write_r.
## Update $1_def.* instead of removing.
##
## Revision 1.14  2000/12/28 12:24:35  peterg
## *** empty log message ***
##

  BEGIN
    MTT_sol_i := first(MTT_sol); MTT_sol := rest(MTT_sol);
    set(lhs(MTT_sol_i),rhs(MTT_sol_i));
  END;

  % No algebraic variables left!
  MTTNYz := 0;
END; % IF MTTNyz>0 and $solve

OUT "$1_ae.r";
IF (MTTNyz>0) THEN % not $solve (or perhaps solution failed?)
BEGIN
    WRITE "MATRIX MTTyz(",MTTNyz,",1)";
    WRITE "%File: $1_ae.r";
    FOR i := 1:MTTNyz DO
	WRITE "MTTyz(",i,",1) := ",MTTyz(i,1);
END; % if MTTNyz>0 (and !$solve)
WRITE ";END;";
SHUT "$1_ae.r";

% Create the matrix declarations
OUT "$1_cse.r1";
write "%";
IF (MTTNx > 0) THEN
BEGIN
    write "MATRIX MTTEdx(", MTTNx, ",", 1, ")$"; 
    write "MATRIX MTTE(", MTTNx, ",", MTTNx, ")$"; 
END;
SHUT "$1_cse.r1";

OUT "$1_csex.r1";
write "%File:$1_csex.r1";
IF (MTTNx > 0) THEN
    write "MATRIX MTTEdx(", MTTNx, ",", 1, ")$"; 
SHUT "$1_csex.r1";

IF MTTNy>0 THEN
BEGIN
  OUT "$1_cseo.r1";
  write "MATRIX MTTY(", MTTNy, ",", 1, ")$"; 
  SHUT "$1_cseo.r1";
END;

%%Create the _cse.r file
OUT "$1_cse.r2";
write "%File: $1_cse.r";
in ("$1_cse_write.r");

# Copyright (C) 2000 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.7  2001/04/11 09:44:26  gawthrop
## Fixed cc and c problems to do with pow(x,y) and integers
## mtt/lib/reduce/fix_c.r is included in rdae2dae and cse2smx_lang for
## -c, -cc and -oct options
##
## Revision 1.6  2000/11/29 20:48:53  peterg
## Zapped unnecessary Npar creation
##
## Revision 1.5  2000/11/09 10:12:24  peterg
## Removed debugging line
##
## Revision 1.4  2000/10/14 16:19:54  peterg

	;;
    ode)
        matrices='dX'
        ns="$Nx"
        ms="1"
	;;
    odeo)
        matrices='Y'
        ns="$Ny"
        ms="1 1"
	;;
    sm)
        matrices='A B C D'
        ns="$Nx $Nx $Ny $Ny"
        ms="$Nx $Nu $Nx $Nu"
	;;

% %%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Field separator
if nargin<3
  FS = ';';
end;

arg = '';
if strcmp(args, '')==0
  L = length(args);
  args_count = 0;
  for i=1:n
    arg_count = 0;
    arg = '';
    if args_count == L
      break;
    end;  
    while args_count < L
      args_count = args_count+1;
      arg_count = arg_count+1;
      ch = str2ch(args,args_count);
      if ch==FS
	break;
      end;
      arg = [arg ch];
    end;
  end;
end;

#! /bin/sh

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

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.57.2.1  2001/05/04 04:07:24  geraint
## Numerical solution of algebraic equations.
## sys_ae.cc written for unsolved inputs.
## Solution of equations using hybrd from MINPACK (as used by Octave fsolve).
##
## Revision 1.57  2001/04/01 03:38:54  geraint
## Reset row to zero after write to file, ready for subsequent runs.
## Eliminates SIGSEGV in Octave when _ode2odes called multiple times.
##
## Revision 1.56  2001/03/30 15:13:58  gawthrop
## Rationalised simulation modes to each return mtt_data
##
## Revision 1.55  2001/03/27 13:21:59  geraint
## Octave version compatibility for save_ascii_data(_for_plotting).
##
## Revision 1.54  2001/03/27 01:14:27  geraint

filename=${sys}_ode2odes.${lang}

if [ -n "$3" ]; then
  method=$3    
else
  method=implicit  
fi

if [ -n "$4" ]; then
    algebraic_solver=$4
else
  algebraic_solver="Reduce_Solver"
fi

echo Creating $filename with $method integration method

# Find system constants
Nx=`mtt_getsize $sys x` # States
Nu=`mtt_getsize $sys u` # Inputs 
Ny=`mtt_getsize $sys y` # Inputs  

	feval_header=parse.h
	save_ascii_data_function=save_ascii_data_for_plotting
	;;
esac

cat <<EOF  > $filename
#include <octave/oct.h>
#include <octave/load-save.h>
#include <octave/lo-mappers.h>
#include <octave/ov-struct.h>
#include <octave/variables.h>

#ifndef STANDALONE
#include <octave/${feval_header}>
#endif

#include "${sys}_def.h"
#include "${sys}_sympar.h"

#ifdef STANDALONE
#include <csignal>
#include <fstream>

#include "mtt_${algebraic_solver}.hh"

extern ColumnVector F${sys}_ae (
	ColumnVector &x,
	ColumnVector &u,
	const double &t,
	ColumnVector &par);

extern ColumnVector F${sys}_input (
	ColumnVector &x,
	ColumnVector &y,
	const double &t,
	ColumnVector &par);

	ColumnVector &u,
	const double &t,
	ColumnVector &par);
 
EOF
fi
cat <<EOF >> $filename

void set_signal_handlers (void);

#endif // STANDALONE

inline ColumnVector
mtt_input (ColumnVector &x,
	   ColumnVector &y,
	   const double &t,
	   ColumnVector &par)
{
#ifdef STANDALONE
  static ColumnVector u  (MTTNU);
  static ColumnVector ui (MTTNYZ);
  static ColumnVector U  (MTTNU+MTTNYZ);
 
  static ${algebraic_solver} ae(F${sys}_ae,MTTNPAR,MTTNU,MTTNX,MTTNY,MTTNYZ);

  u = F${sys}_input (x, y, t, par);
  if (MTTNYZ == 0)
    {
      return u;
    }
  else
    {
      return ae.solve(x,u,t,par);
    }
#else
  static octave_value_list args, f;
  args (0) = octave_value (x);
  args (1) = octave_value (y);
  args (2) = octave_value (t);
  args (3) = octave_value (par);
  f = feval ("${sys}_input", args, 1);

	   ostream &file = cout)
{
  static Matrix data;
  static int row;

  if (dump_data)
    {
       if (row > 0)
         {
           Matrix written_data = data.extract (0, 0, row-1, data.cols ()-1);
           $save_ascii_data_function (file, written_data, "MTT_data");
	 }
       return;
    }

  const int nx = x.length (), ny = y.length ();
  register int col = 0;

  if (0 == row)

#ifdef STANDALONE
void dump_data (ostream &file)
{
  ColumnVector null (0.0);
  mtt_write (0.0, null, null, 0, true, file);
}



void handle_signal (int signum)
{
  // handle some signals to ensure data is written.
  cerr << "# Writing data to MTT.core (signal " << signum << ")" << endl;
  ofstream corefile ("MTT.core");
  dump_data (corefile);
  switch (signum)
    {
      case SIGFPE:
        // Intel chips do not raise SIGFPE for DIVZERO :-(
//        raise (SIGABRT);
        break;
      case SIGINT:
	break;
      case SIGQUIT:
	signal (SIGQUIT, SIG_DFL);
	raise (SIGQUIT);
	break;

# Copyright (C) 2000 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.31  2001/04/03 14:49:42  gawthrop
## Revised to incorporate new ssim (sensitivity simulation)
## representation (m only just now).
##
## Revision 1.30  2001/03/30 15:13:58  gawthrop
## Rationalised simulation modes to each return mtt_data
##
## Revision 1.29  2001/03/27 13:10:23  geraint
## Improved determination of Octave version.
##
## Revision 1.28  2001/02/17 03:48:17  geraint

# Representation-specific stuff

eqnargs='mttx,mttu,mttt,mttpar'
inputeqnargs='mttx,mtty,mttt,mttpar'
case $rep in
    ae)
	states=yes;
	inputs=yes;
	parameters=yes;
	output=mttyz;
	args=$eqnargs;
	;;
    cse)
	states=yes;
        inputs=yes;
	parameters=yes;
        output='mttdx,mtte'
        args=$eqnargs
	;;

	  states=no;
	  parameters=no;
          declareinputs=no;
        else
	  states=yes;
	  parameters=yes;
          declareinputs=yes
	  declarestates=yes
        fi
	;;
    logic)
	states=no;
	inputs=no;
	parameters=yes;
        output=mttopen

$Lc Inputs $Rc
EOF

N=`n2m 1 $Nu`
for i in $N; do
  echo $constant_declaration 'mttu'$i' = mttu('$i$minusone');'
done
cat <<EOF

$Lc Unknown Inputs $Rc
EOF

Ni=`n2m 1 $Nyz`
for i in $Ni; do
  echo $constant_declaration 'mttui'$i' = mttu('$Nu+$i$minusone');'
done
fi

if [ "$declareinputs" = "yes" ]; then
  declare_vars input
fi

if [ "$declaredummies" = "yes" ]; then

	    arg_type="Octave_map"
	    arg_size="";
	    ;;
	mttu)
	    arg_type="ColumnVector"
	    arg_size="(MTTNU)"
	    ;;
	mttyz)
	    arg_type="ColumnVector"
	    arg_size="(MTTNYZ)"
	    ;;
	mtty)
	    arg_type="ColumnVector"
	    arg_size="(MTTNY)"
	    ;;
	mttt)
	    arg_type="const double"
	    arg_size=""

# Copyright (C) 2000 by Peter J. Gawthrop

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.27  2001/02/08 00:39:56  geraint
## Removed clear from code generated by: mtt -c -i euler -o
##
## Revision 1.26  2000/11/27 11:52:10  peterg
## Changed some greps to be case insensitive
##
## Revision 1.25  2000/10/14 08:50:07  peterg
## Use new mtt_header
##
## Revision 1.24  2000/10/11 09:07:53  peterg

echo Creating $out

# Remove the old log file
rm -f mtt_r2m.log

#Header
case $rep in
    ae|csex|cseo|ode|odeo)
    vectorise=yes
       ;;
    *)
esac

    mtt_header  $1 $rep $ext  > mtt_junk

#define pi 3.14159 // Predefine pi

#ifndef HAVE_USEFUL_FUNCTIONS_HH
#define HAVE_USEFUL_FUNCTIONS_HH


#ifndef __cplusplus
#define inline			/* strip */



#endif // ! __cplusplus


static inline double
max (const double x1, const double x2)
{
  return static_cast<double>((x1 >= x2) ? x1 : (x1 < x2) ? x2 : 0);
}

static inline double
min (const double x1, const double x2)
{
  return static_cast<double>((x1 <= x2) ? x1 : (x1 > x2) ? x2 : 0);
}

static inline double
nonsingular (const double x)
{
  return static_cast<double>((x == 0) ? 1.0e-30 : x);
}

static inline double
sign (const double x)
{
  return static_cast<double>((x > 0) ? +1 : (x < 0) ? -1 : 0);
}


// Octave functions
#ifdef __cplusplus