#! /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.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 " int itime;"$
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 " itime = 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,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 " $1_ode();"$
write " }"$
write " /* Set up system inputs */"$
write " $1_input(time,itime);"$
write " time = time + DT;"$
write " itime = itime + 1;"$
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