#! /bin/sh
######################################
##### Model Transformation Tools #####
######################################
# Bourne shell script: ode2odesol_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$
###############################################################
# Inform user
echo Creating $1_odesol.c
# Remove the old log file
rm -f ode_r2c.log
# Use reduce to accomplish the transformation
reduce >ode2odesol_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_odesol.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_odesol.c"$
write "Generated by MTT"$
write "*/"$
write " "$
%Program heading
write "#include <stdio.h>"$
%Set up some variables - need to do this better sometime
write "/* Set up some variables - need to do this better sometime */"$
write "#define DT 0.1 /* Time step (for printing) */"$
write "#define LAST 10.0 /* Last time */"$
write "#define STEPFACTOR 1000 /* Integration steps per time step */"$
write "/* Declare standard arrays */"$
write "float y[", MTTNy+1, "]; /* $1_ode output */"$
write "float dx[", MTTNx+1, "]; /* $1_ode state derivative */"$
write "float x[", MTTNx+1, "]; /* $1_ode state */"$
write "float u[", MTTNu+1, "]; /* $1_ode input */"$
%External (global) variable list
write "/* External (global) variable list */ "$
IF MTTNvar>0 THEN
BEGIN
FOR i := 1:MTTNvar DO
IF numberp(MTTVar(i,1))
THEN
BEGIN
% Do nowt
END
ELSE
BEGIN
write "float ", MTTVar(i,1), ";"$
END$
END$
write "/* File */ "$
write " FILE *fopen(), *fp;"$
write "main()"$
write " "$
write "{"$
write "/* Declare standard arrays */"$
write " extern float y[", MTTNy+1, "]; /* $1_ode output */"$
write " extern float dx[", MTTNx+1, "]; /* $1_ode state derivative */"$
write " extern float x[", MTTNx+1, "]; /* $1_ode state */"$
write " extern float u[", MTTNu+1, "]; /* $1_ode input */"$
%External (global) variable list
write "/* External (global) variable list */ "$
IF MTTNvar>0 THEN
BEGIN
FOR i := 1:MTTNvar DO
IF numberp(MTTVar(i,1))
THEN
BEGIN
% Do nowt
END
ELSE
BEGIN
write " extern float ", MTTVar(i,1), ";"$
END$
END$
write "/* Counters etc*/ "$
write " float time;"$
write " float dt;"$
write " int i;"$
write " int k;"$
write "/* functions */ "$
write " extern $1_numpar();"$
%Open the output file
write "/* %Open the output file */"$
write "fp = fopen(""$1_odes.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;"$
write " dt = DT/STEPFACTOR;"$
write " for (i=1; i<=", MTTNx, "; i++)"$
write " x[i] = 0.0;"$
write " for (i=1; i<=", MTTNu, "; i++)"$
write " u[i] = 1.0;"$
write " fprintf(fp, ""data = [\n"");"$
%Main (Euler) integration loop
write "/* Main (Euler) integration loop */"$
write " while (time<LAST)"$
write " {"$
write " for (k=1; k<=STEPFACTOR; k++)"$
write " {"$
write " $1_ode(y,dx,x,u);"$
write " for (i=1; i<=", MTTNx, "; i++)"$
write " x[i] = x[i] + dx[i]*dt;"$
write " }"$
write " fprintf(fp, ""%5.4f "",time);"$
write " for (i=1; i<=", MTTNy, "; i++)"$
write " fprintf(fp, ""%5.4f "", y[i]);"$
write " fprintf(fp, ""\n"");"$
write " time = time + DT;"$
write " }"$
write " fprintf(fp, ""];\n"");"$
write "}"$
SHUT "$1_odesol.c";
EOF