#! /bin/sh
######################################
##### Model Transformation Tools #####
######################################
# Bourne shell script: ode2obs_r
# Odrinary differential equations to observer function equations
# P.J.Gawthrop 14 June 1991, 8 Aug 1991, 2 April 1992, 14 April 1994, 28 Dec 94,
# 12th July 1995, April 1996
# Copyright (c) P.J.Gawthrop 1991, 1992, 1994, 1995, 1996.
###############################################################
## Version control history
###############################################################
## $Id$
## $Log$
## Revision 1.5 1998/04/07 08:12:12 peterg
## Added affine form.
##
## Revision 1.4 1998/04/07 05:45:12 peterg
## Reverted to an older version due to untraceable bug in new version
##
## Revision 1.3 1996/08/25 10:07:05 peter
## Remove a du state ment causaing touble
## - but needs more work.
##
## Revision 1.2 1996/08/25 08:38:14 peter
## Error handling added.
##
## Revision 1.1 1996/08/25 08:37:44 peter
## Initial revision
##
###############################################################
#Inform user
echo Creating $1_obs.r
# Remove the old log file
rm -f ode2obs_r.log
# Use reduce to accomplish the transformation
reduce >ode2obs_r.log << EOF
%Read the formatting function
in "$MTTPATH/trans/reduce_matrix.r";
OFF Echo;
OFF Nat;
ON NERO;
in "$1_def.r";
MTTdxs := MTTdX; %Save the symbolic form of dX
%Set default values - reset by obspar file.
MTTGPCNy := 2;
MTTGPCNu := 0;
%%in "$1_sympar.r";
in "$1_ode.r";
%%in "$1_simp.r";
% Read the obs form parameters.
in "$1_obspar.r";
%Create the U vector of input derivatives.
MTTNuu := (MTTGPCNu+1)*MTTNu;
MTTNuu1 := MTTGPCNu*MTTNu;
%FOR i := 1:MTTNu DO
% MTTUU(i,1) := MTTU(i,1);
%END;
%
%IF MTTGPCNu>0 THEN
%BEGIN
% FOR i := 1:MTTNu DO
% MTTUU(i+MTTNu,1) := MTTdU(i,1);
% END;
%END;
MTTU := MTTU;
MTTdU := MTTdU; % ---- removed temporarily, needs def file change?
%Create the Y vector of output derivatives.
MTTNyy := (MTTGPCNy+1)*MTTNy;
Matrix MTTYY(MTTNyy,1);
FOR i := 1:MTTNy DO
MTTYY(i,1) := MTTY(i,1);
END;
l := MTTNy;
FOR i := 1:MTTGPCNy DO
FOR j := 1:MTTNy DO
BEGIN
l := l+1;
MTTYY(l,1) := 0;
FOR k := 1:MTTNx DO %Derivatives wrt x
BEGIN
xk := MTTX(k,1);
MTTYY(l,1) := MTTYY(l,1) + df(MTTYY(l-MTTNy,1), xk, 1)*MTTdX(k,1);
END;
IF MTTGPCNu>0 THEN
FOR k := 1:MTTGPCNu DO %Non-zero derivatives of u
BEGIN
FOR kk := 1:MTTNu DO
BEGIN
uk := MTTUU(kk,k);
MTTYY(l,1) := MTTYY(l,1) + df(MTTYY(l,1), uk, 1)*MTTUU(1,k+1);
END;
END;
END;
END;
%%Create O_x - derivative of YY wrt x
%Matrix MTTO_x(MTTNyy,MTTNx);
%FOR j := 1:MTTNx DO
% BEGIN
% xj := MTTX(j,1);
% FOR i := 1:MTTNyy DO
% BEGIN
% MTTO_x(i,j) := df(MTTYY(i,1), xj);
% END;
% END;
%%Create O_u - derivative of YY wrt u (Assumes GPC Nu = 0)
%MTTNNu := (MTTGPCNu+1)*MTTNu;
%Matrix MTTO_u(MTTNyy,MTTNNu);
%FOR j := 1:MTTNNu DO
% BEGIN
% uj := MTTu(j,1);
% FOR i := 1:MTTNyy DO
% BEGIN
% MTTO_u(i,j) := df(MTTYY(i,1), uj);
% END;
% END;
%%Create O_uu - derivative of O_u wrt u (Assumes GPC Nu = 0)
%%This is a multi-dimensional matrix kth elements stacked sideways.
%Matrix MTTO_uu(MTTNyy,MTTNNu*MTTNNu);
%FOR k := 1:MTTNNu DO
% BEGIN
% uk := MTTu(k,1);
% FOR j := 1:MTTNNu DO
% BEGIN
% FOR i := 1:MTTNyy DO
% BEGIN
% jk := j+(k-1)*MTTNu;
% MTTO_uu(i,jk) := df(MTTO_u(i,j), uk);
% END;
% END;
% END;
%%Create O_ux - derivative of O_u wrt x
%%This is a multi-dimensional matrix kth elements stacked sideways.
%Matrix MTTO_ux(MTTNyy,MTTNu*MTTNx);
%FOR k := 1:MTTNx DO
% BEGIN
% xk := MTTx(k,1);
% FOR j := 1:MTTNu DO
% BEGIN
% FOR i := 1:MTTNyy DO
% BEGIN
% jk := j+(k-1)*MTTNu;
% MTTO_ux(i,jk) := df(MTTO_u(i,j), xk);
% END;
% END;
% END;
% Try and split the obs function into affine form (O(x,u) = o(x) + h(x)u)
affine := 1;
MTTNuu := (MTTGPCNu+1)*MTTNu;
write "1 affine := ", affine, ";";
%IF (affine=1) THEN
%BEGIN
MATRIX MTTObs_o(MTTNyy,1);
MATRIX MTTObs_h(MTTNyy,MTTNuu);
write "2 affine := ", affine, ";";
FOR i := 1:MTTNyy DO
BEGIN
write "3 affine := ", affine, ";";
MTTObs_o(i,1) := MTTYY(i,1);
k := 0;
FOR j := 1:MTTNu DO
FOR jj := 0:MTTGPCNu DO
BEGIN
k := k+1;
%Expand as polynomial in u_j^[jj]
coeffs := coeff(MTTObs_o(i,1), MTTUU(j,jj+1));
MTTObs_o(i,1) := first(coeffs);
IF length(coeffs)>1 THEN
MTTObs_h(i,k) := second(coeffs);
IF length(coeffs)>2 THEN
affine := -1;
END;
END;
%END;
%%Create the _obs.r file
OUT "$1_obs.r";
IF affine=1 THEN
write "affine := 1; % The O function is affine in u"
ELSE
write "affine := 0; % The O function is not affine in u";
write "% The matrix sizes";
write "MTTNyy := ", MTTNyy, ";";
write "MTTNuu := ", MTTNuu, ";";
mtt_matrix := MTTYY$
mtt_matrix_n := MTTNyy$
mtt_matrix_m := 1$
mtt_matrix_name := "MTTYY"$
reduce_matrix();
IF affine=1 THEN
BEGIN
mtt_matrix := MTTObs_o$
mtt_matrix_n := MTTNyy$
mtt_matrix_m := 1$
mtt_matrix_name := "MTTObs_o"$
reduce_matrix();
mtt_matrix := MTTObs_h$
mtt_matrix_n := MTTNyy$
mtt_matrix_m := MTTNuu$
mtt_matrix_name := "MTTObs_h"$
reduce_matrix();
END;
%MTTO_x := MTTO_x;
%MTTO_u := MTTO_u;
%MTTO_uu := MTTO_uu;
%MTTO_ux := MTTO_ux;
write ";END;";
SHUT "$1_obs.r";
quit;
EOF
# Now invoke the standard error handling.
mtt_error_r ode2obs_r.log