#! /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