#! /bin/sh

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

# Bourne shell script: ode2odes_r2c
# Reduce ordinary differential equations to  differential-algebraic 
# equations solution in the form of a c program.

# Euler integration of the state is included.

# NB Arrays should are  defined to be one larger than expected 
# - the 0 element is not used.

# Copyright (c) P.J.Gawthrop 1997.

###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
# Revision 1.10  1997/05/10  10:05:15  peterg
# Put _input in inner loop in front of call to _ode
#
## Revision 1.9  1997/05/10 09:54:34  peterg
## Moved _input call to after the inner integration loop.
##
## Revision 1.8  1997/05/10 08:12:23  peterg
## Put second argument into _input.
##
## Revision 1.7  1997/05/10 07:01:15  peterg
## _input called in outer loop only.
## time updated in outer loop only.
## Integer time (itime) introduced and updated in outer loop - maybe
## useful for discrete events.
## Integer time (itime) passed to _input.
##
## Revision 1.6  1997/05/06 13:53:32  peterg
## Now uses the preprocessor to declare sizes -- MTTNX etc
##
# Revision 1.5  1997/05/01  13:50:11  peterg
# Replaced float by double.
#
# Revision 1.4  1997/05/01  13:43:44  peterg
# Changed double to float.
#
# Revision 1.3  1997/05/01  11:15:33  peterg
# Back under RCS
#
# Revision 1.2  1997/03/20  14:36:56  peterg
# Includes the sympar.h file
#
## Revision 1.1  1997/01/21 22:54:54  peterg
## Initial revision
##
###############################################################


# Inform user
echo Creating $1_odes.c

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

# Use reduce to accomplish the transformation
reduce >ode2odes_r2c.log << EOF

%Read the reduce definitions file
in "$1_def.r";

%Set up the number of argument variables to zero in case the user has forgotten
MTTNVar := 0;

%Read the symbolic parameters file
in "$1_sympar.r";

ON BigFloat, NumVal;
PRECISION 16; %Compatible with Matlab
%OFF Nat;

ON NERO;        % Suppress zero elements

%Generate the Header part
OUT "$1_odes.c";

write "/*"$
write "Program to solve ode for system $1"$
write "NB Arrays are defined to be one larger than expected"$
write " - the 0 element is not used."$

write "File $1_odes.c"$
write "Generated by MTT"$
write "*/"$
write " "$


%Program heading
write "#define MTTNX ", MTTNx $
write "#define MTTNY ", MTTNy $
write "#define MTTNU ", MTTNu $
write "#define MTTNX1 ", MTTNx+1 $
write "#define MTTNY1 ", MTTNy+1 $
write "#define MTTNU1 ", MTTNu+1 $
write "#define MTTNX2 ", MTTNx+2 $
write "#define MTTNY2 ", MTTNy+2 $
write "#define MTTNU2 ", MTTNu+2 $

write "#include <stdio.h>"$
write "#include ""$1_odes.h"" "$
write "#include ""$1_ode.c"" "$
write "#include ""$1_input.c"" "$
write "#include ""$1_numpar.c"" "$

%External (global) variable list
write "#include ""$1_sympar.h"" "$


write "/* Declare standard arrays */"$
write "double y[MTTNY1]; /* $1_ode output */"$
write "double dx[MTTNX1]; /* $1_ode state derivative */"$ 
write "double x[MTTNX1]; /* $1_ode state */"$ 
write "double u[MTTNU1]; /* $1_ode input */"$


write "/* Files */ "$
write "  FILE *fopen(), *fps, *fpso;"$

write "main()"$
write "  "$
write "{"$



write "/* Counters etc*/ "$
write "  double time;"$
write "  double dt;"$
write "  int i;"$
write "  int k;"$

write "/*functions */ "$
write "  extern  $1_numpar();"$

%Open the output files
write "/* %Open the output file */"$
write "fps = fopen(""$1_odes.m"", ""w"");  "$
write "fpso = fopen(""$1_odeso.m"", ""w"");  "$

%Set up user-defined constants
write "/* Set up user-defined constants */"$
write "  $1_numpar();"$

%Initialise main (Euler) integration loop
write "/* Initialise main (Euler) integration loop */"$
write "  time = 0.0;"$
write "  dt = DT/STEPFACTOR;"$
write "  for (i=1; i<=MTTNX; i++)"$
write "      x[i] = 0.0;"$

%Set up system inputs
write "/* Set up system inputs */"$
write "  for (i=1; i<=", MTTNu, "; i++)"$
write "      u[i] = 1.0;"$

write "  $1_input(0.0);"$

write "    fprintf(fps, ""function data = ", "$1_odes \n"");"$
write "    fprintf(fps, ""data = [\n"");"$
write "    fprintf(fpso, ""function data = ", "$1_odeso \n"");"$
write "    fprintf(fpso, ""data = [\n"");"$

% Compute the first output
write "/* Compute the first output */"$
write "  $1_ode(y,dx,x,u);"$

%Main (Euler) integration loop
write "/* Main (Euler) integration loop */"$

write "  while (time<LAST)"$
write "  {"$

%Write to output to file
write "/* Write to output file */"$
write "    fprintf(fpso, ""%5.4f "",time);"$
write "    for (i=1; i<=MTTNY; i++)"$
write "      fprintf(fpso, ""%5.4f "", y[i]);"$
write "    fprintf(fpso, ""\n"");"$

%Write to state to file
write "/* Write to state file */"$
write "    fprintf(fps, ""%5.4f "",time);"$
write "    for (i=1; i<=MTTNX; i++)"$
write "      fprintf(fps, ""%5.4f "", x[i]);"$
write "    fprintf(fps, ""\n"");"$


write "/* Inner integration loop */"$
write "    for (k=1; k<=STEPFACTOR; k++)"$
write "    {"$
write "      for (i=1; i<=MTTNX; i++)"$
write "        x[i] = x[i] + dx[i]*dt;"$

write "      /* Set up system inputs */"$
write "      $1_input(time);"$
write "      $1_ode();"$
write "    }"$


write "      time = time + DT;"$
write "  }"$

%Write to files
%Write to output to file
write "/* Write to output file */"$
write "    fprintf(fpso, ""%5.4f "",time);"$
write "    for (i=1; i<=MTTNY; i++)"$
write "      fprintf(fpso, ""%5.4f "", y[i]);"$
write "    fprintf(fpso, ""\n"");"$

%Write to state to file
write "/* Write to state file */"$
write "    fprintf(fps, ""%5.4f "",time);"$
write "    for (i=1; i<=MTTNX; i++)"$
write "      fprintf(fps, ""%5.4f "", x[i]);"$
write "    fprintf(fps, ""\n"");"$


write "  fprintf(fps, ""];\n"");"$
write "  fprintf(fpso, ""];\n"");"$

write "  return 0;"$
write "}"$


SHUT "$1_odes.c";


EOF



