#! /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.17 1998/05/15 07:47:23 peterg
## Includes sign.c by default
##
## Revision 1.16 1998/05/14 08:25:52 peterg
## Corrected time skew.
##
## Revision 1.15 1998/05/14 08:20:29 peterg
## Start time at DT - the result at time zero is computed outside the
## loop
##
## Revision 1.14 1998/05/13 08:57:27 peterg
## Now uses simpar.h in place of odes.h
##
## Revision 1.13 1998/02/24 13:34:45 peterg
## Back under RCS
##
# Revision 1.12 1997/05/15 08:39:56 peterg
# Don't initialise states - now done in numpar file.
#
# Revision 1.11 1997/05/12 16:00:54 peterg
# Removed itime again.,
#
# 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 <math.h>"$
write "#include ""sign.c"" "$
write "#include ""$1_simpar.h"" "$
write "#include ""$1_ode.c"" "$
write "#include ""$1_input.c"" "$
write "#include ""$1_numpar.c"" "$
write "#include ""$1_state.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();"$
%Set up initial state
write "/* Set up initial state */"$
write " $1_state();"$
%Initialise main (Euler) integration loop
write "/* Initialise main (Euler) integration loop */"$
write " time = 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-DT))"$
write " {"$
%Write to output to file
write "/* Write to output file */"$
write " fprintf(fpso, ""%5.4g "",time);"$
write " for (i=1; i<=MTTNY; i++)"$
write " fprintf(fpso, ""%5.4g "", y[i]);"$
write " fprintf(fpso, ""\n"");"$
%Write to state to file
write "/* Write to state file */"$
write " fprintf(fps, ""%5.4g "",time);"$
write " for (i=1; i<=MTTNX; i++)"$
write " fprintf(fps, ""%5.4g "", x[i]);"$
write " fprintf(fps, ""\n"");"$
write " time = time + DT;"$
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 " }"$
%Write to files
%Write to output to file
write "/* Write to output file */"$
write " fprintf(fpso, ""%5.4g "",time);"$
write " for (i=1; i<=MTTNY; i++)"$
write " fprintf(fpso, ""%5.4g "", y[i]);"$
write " fprintf(fpso, ""\n"");"$
%Write to state to file
write "/* Write to state file */"$
write " fprintf(fps, ""%5.4g "",time);"$
write " for (i=1; i<=MTTNX; i++)"$
write " fprintf(fps, ""%5.4g "", x[i]);"$
write " fprintf(fps, ""\n"");"$
write " fprintf(fps, ""];\n"");"$
write " fprintf(fpso, ""];\n"");"$
write " return 0;"$
write "}"$
SHUT "$1_odes.c";
EOF