#! /bin/sh

     ###################################### 
     ##### Model Transformation Tools #####
     ######################################

# Bourne shell script: cse2scse_r
# Reduce constrained-state equations to sensitivity version
# P.J.Gawthrop 10th May 199, 8th August 1991, April 1994, Dec 1994
# Copyright (c) P.J.Gawthrop, 1999


# The sensitivity parameter
theta=$2
thetas=$2"s"

echo $theta $thetas
# Number of states
Nx=`grep "MTTNx " <$1_def.r | awk '{print $3}' | sed 's/;//'`

#Inform user
echo "Creating $1_scse.r (for parameter $theta, $Nx states)"

# Remove the old log file
rm -f cse2scse_r.log

# Use reduce to accomplish the transformation
$SYMBOLIC  << EOF >cse2scse_r.log

%Read the formatting function
in "$MTTPATH/trans/reduce_matrix.r";

%Read the definitions file
in "$1_def.r";

%Read the substitution file
in "$1_subs.r";

%Read the constrained-state equations file
in "$1_cse.r";

% Compute the sensitivity E matrix
matrix MTTssE(MTTNx,MTTNx);
clear MTTx; % Dont use - mkid is better

FOR ii := 1:MTTNx DO
  FOR jj := 1:MTTNx DO
  BEGIN
    % First with respect to theta...
    MTTssE(ii,jj) := df(MTTE(ii,jj), $theta)*$thetas;
    % then with respect to x
    FOR i := 1:MTTNx DO
    BEGIN
      xi := mkid(MTTx,i);
      sxi := mkid(MTTsx,i);
      MTTssE(ii,jj) := MTTssE(ii,jj) + df(MTTE(ii,jj), xi)*sxi;
    END;
  END;

% Compute the sensitivity of the RHS of the cse
matrix MTTssEdx(MTTNx,1);


FOR ii := 1:MTTNx DO
BEGIN 
  % First with respect to theta
  MTTssEdx(ii,1) := $thetas*df(MTTEdx(ii,1), $theta);
  
% Then with respect to x
  FOR i := 1:MTTNx DO
  BEGIN
    xi := mkid(MTTx,i);
    sxi := mkid(MTTsx,i);
    MTTssEdx(ii,1) := MTTssEdx(ii,1) + df(MTTEdx(ii,1), xi)*sxi;
  END; 
END; 

% Now reorganise everything into composite system
%  - odd rows are the system 
%  - even rows are the sensitivity system
% NB at this stage, the states are numbered incorrectly - sorted out below.

%E matrix 
MTTNx2 := 2*MTTNx;
matrix MTTsE(MTTNx2,MTTNx2);
FOR i := 1:MTTNx DO
  FOR j := 1:MTTNx DO
  BEGIN
    MTTsE(2*i-1,2*j-1) := MTTE(i,j);  % System
    MTTsE(2*i,2*j)     := MTTE(i,j);  % Sensitivity system
    MTTsE(2*i,2*j-1)   := MTTssE(i,j);  % Sensitivity system
  END; 

%dX matrix
matrix MTTsEdX(MTTNx2,1);

FOR i := 1:MTTNx DO
BEGIN
  MTTsEdX(2*i-1,1) := MTTEdx(i,1);
  MTTsEdX(2*i,1)   := MTTssEdx(i,1);
END; 

OFF nat;

OUT "$1_scse.r";

write "%File: $1_scse.r";


% E matrix
MTT_Matrix := MTTsE$ 
MTT_Matrix_name := "MTTsE"$
MTT_Matrix_n := MTTNx2$
MTT_Matrix_m := MTTNx2$
Reduce_Matrix()$

% State derivative
MTT_Matrix := MTTsEdX$ 
MTT_Matrix_name := "MTTsEdX"$
MTT_Matrix_n := MTTNx2$
MTT_Matrix_m := 1$
Reduce_Matrix()$

% Output
MTT_Matrix := MTTY$ 
MTT_Matrix_name := "MTTY"$
MTT_Matrix_n := MTTNy$
MTT_Matrix_m := 1$
Reduce_Matrix()$

write "END;";

SHUT "$1_scse.r";
quit;
EOF

# Now invoke the standard error handling.
mtt_error_r cse2scse_r.log


## Now reorganise the  states
mv -f $1_scse.r mtt_junk


##echo "Nx = $Nx"
awk '{
  ## Make sure all MTTn variables are followed by a space
  gsub(/mttx[0-9]*/, "& ");

  for (i=Nx;i>0;i--) {

    state = sprintf("mttx%i ",i);
    newstate = sprintf("mttx%i ",2*i-1);
    sstate = sprintf("mttsx%i",i);
    newsstate = sprintf("mttx%i",2*i);

    gsub(state,newstate);
    gsub(sstate,newsstate);

  }
  print $0
}' Nx=$Nx <mtt_junk > $1_scse.r