% Copyright (C) 2001, 2000 Marc Compere
%
% This program is free software; you can redistribute it and/or modify it
% under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
% General Public License for more details at www.gnu.org/copyleft/gpl.html.
%
% This file is intended for use with Octave.
% --------------------------------------------------------------------

% ode78 (v1.14) Integrates a system of ordinary differential equations using
% 7th order formulas.
%
% This is a 7th-order accurate integrator therefore the local error normally
% expected is O(h^8).  However, because this particular implementation
% uses the 8th-order estimate for xout (i.e. local extrapolation) moving
% forward with the 8th-order estimate will yield errors on the order of O(h^9).
%
% The order of the RK method is the order of the local *truncation* error, d,
% which is the principle error term in the portion of the Taylor series
% expansion that gets dropped, or intentionally truncated.  This is different
% from the local error which is the difference between the estimated solution
% and the actual, or true solution.  The local error is used in stepsize
% selection and may be approximated by the difference between two estimates of
% different order, l(h) = x_(O(h+1)) - x_(O(h)).  With this definition, the
% local error will be as large as the error in the lower order method.
% The local truncation error is within the group of terms that gets multipled
% by h when solving for a solution from the general RK method.  Therefore, the
% order-p solution created by the RK method will be roughly accurate to O(h^(p+1))
% since the local truncation error shows up in the solution as h*d, which is
% h times an O(h^(p)) term, or rather O(h^(p+1)).
% Summary:   For an order-p accurate RK method,
%            - the local truncation error is O(h^p)
%            - the local error used for stepsize adjustment and that
%              is actually realized in a solution is O(h^(p+1))
%
% This requires 13 function evaluations per integration step.
%
% Relevant discussion on step size choice can be found on pp.90,91 in
% U.M. Ascher, L.R. Petzold, Computer Methods for  Ordinary Differential Equations
% and Differential-Agebraic Equations, Society for Industrial and Applied Mathematics
% (SIAM), Philadelphia, 1998
%
% More may be found in the original author's text containing numerous
% applications on ordinary and partial differential equations using Matlab:
%
%     Howard Wilson and Louis Turcotte, 'Advanced Mathematics and 
%     Mechanics Applications Using MATLAB', 2nd Ed, CRC Press, 1997
%
%
% [tout, xout] = ode78(FUN,tspan,x0,ode_fcn_format,tol,trace,count,hmax)
%
% INPUT:
% FUN   - String containing name of user-supplied problem description.
%         Call: xprime = fun(t,x) where FUN = 'fun'.
%         t      - Time (scalar).
%         x      - Solution column-vector.
%         xprime - Returned derivative COLUMN-vector; xprime(i) = dx(i)/dt.
% tspan - [ tstart, tfinal ]
% x0    - Initial value COLUMN-vector.
% ode_fcn_format - this specifies if the user-defined ode function is in
%         the form:     xprime = fun(t,x)   (ode_fcn_format=0, default)
%         or:           xprime = fun(x,t)   (ode_fcn_format=1)
%         Matlab's solvers comply with ode_fcn_format=0 while
%         Octave's lsode() and sdirk4() solvers comply with ode_fcn_format=1.
% tol   - The desired accuracy. (optional, default: tol = 1.e-6).
% trace - If nonzero, each step is printed. (optional, default: trace = 0).
% count - if nonzero, variable 'rhs_counter' is initalized, made global
%         and counts the number of state-dot function evaluations
%         'rhs_counter' is incremented in here, not in the state-dot file
%         simply make 'rhs_counter' global in the file that calls ode78
% hmax  - limit the maximum stepsize to be less than or equal to hmax
%
% OUTPUT:
% tout  - Returned integration time points (row-vector).
% xout  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, xout).

%   Daljeet Singh & Howard Wilson
%   Dept. Of Electrical Engg., The University of Alabama.
%   11-24-1988.
%
% modified by:
% Marc Compere
% CompereM@asme.org
% created : 06 October 1999
% modified: 19 May 2001


function [tout,xout] = ode78(FUN,tspan,x0,ode_fcn_format,tol,trace,count,hmax)

% The Fehlberg coefficients:
alpha_ = [ 2./27., 1/9, 1/6, 5/12, 0.5, 5/6, 1/6, 2/3, 1/3, 1, 0, 1 ]';
beta_  = [ 2/27, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          1/36, 1/12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          1/24, 0, 1/8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          5/12, 0, -25/16, 25/16, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
          0.05, 0, 0, 0.25, 0.2, 0, 0, 0, 0, 0, 0, 0, 0 ;
          -25/108, 0, 0, 125/108, -65/27, 125/54, 0, 0, 0, 0, 0, 0, 0 ;
          31/300, 0, 0, 0, 61/225, -2/9, 13/900, 0, 0, 0, 0, 0, 0 ;
          2, 0, 0, -53/6, 704/45, -107/9, 67/90, 3, 0, 0, 0, 0, 0 ;
          -91/108, 0, 0, 23/108, -976/135, 311/54, -19/60, 17/6, -1/12, 0, 0, 0, 0 ;
          2383/4100, 0, 0, -341/164, 4496/1025, -301/82, 2133/4100, 45/82, 45/164, 18/41, 0, 0, 0 ;
          3/205, 0, 0, 0, 0, -6/41, -3/205, -3/41, 3/41, 6/41, 0, 0, 0 ;
          -1777/4100, 0, 0, -341/164, 4496/1025, -289/82, 2193/4100, 51/82, 33/164, 12/41, 0, 1, 0 ]';
chi_  = [ 0, 0, 0, 0, 0, 34/105, 9/35, 9/35, 9/280, 9/280, 0, 41/840, 41/840]';
psi_  = [ 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, -1 ]';

pow = 1/8; % see p.91 in the Ascher & Petzold reference for more infomation.


if nargin < 8, hmax = (tspan(2) - tspan(1))/2.5; end
if nargin < 7, count = 0; end
if nargin < 6, trace = 0; end
if nargin < 5, tol = 1.e-6; end
if nargin < 4, ode_fcn_format = 0; end

% Initialization
t0 = tspan(1);
tfinal = tspan(2);
t = t0;
% the following step parameters are used in ODE45
% hmax = (tfinal - t)/5;
% hmin = (tfinal - t)/20000;
% h = (tfinal - t)/100;
% The following parameters were taken because the integrator has
% higher order than ODE45. This choice is somewhat subjective.
%hmin = (tfinal - t)/10000;
hmin = (tfinal - t)/1e20;
h = (tfinal - t)/50;
x = x0(:);          % the '(:)' ensures x is initialized as a column vector
f = x*zeros(1,13);  % f needs to be an Nx13 matrix where N=number of rows in x
tout = t;
xout = x.';
tau = tol * max(norm(x,'inf'), 1);

if count==1,
 global rhs_counter
 if ~exist('rhs_counter'),rhs_counter=0;,end
end % if count

if trace
%  clc, t, h, x
%   clc, t, x
   clc, t
end

% The main loop
   while (t < tfinal) & (h >= hmin)
      if t + h > tfinal, h = tfinal - t; end

      % Compute the slopes
      if (ode_fcn_format==0), % (default)
         f(:,1) = feval(FUN,t,x);
         for j = 1: 12,
            f(:,j+1) = feval(FUN, t+alpha_(j)*h, x+h*f*beta_(:,j));
         end
      else, % ode_fcn_format==1
         f(:,1) = feval(FUN,x,t);
         for j = 1: 12,
            f(:,j+1) = feval(FUN, x+h*f*beta_(:,j), t+alpha_(j)*h);
         end
      end % if (ode_fcn_format==1)


      % increment rhs_counter
      if count==1, rhs_counter = rhs_counter + 13; end

      % Truncation error term
      gamma1 = h*41/840*f*psi_;

      % Estimate the error and the acceptable error
      delta = norm(gamma1,'inf');
      tau = tol*max(norm(x,'inf'),1.0);

      % Update the solution only if the error is acceptable
      if delta <= tau
         t = t + h;
         x = x + h*f*chi_;  % this integrator uses local extrapolation
         tout = [tout; t];
         xout = [xout; x.'];
      end
      if trace
%        home, t, h, x
%        home, t, x
         home, t, h
      end

      % Update the step size
      if delta == 0.0
       delta = 1e-16;
      end
      h = min(hmax, 0.8*h*(tau/delta)^pow);

   end;

   if (t < tfinal)
      disp('SINGULARITY LIKELY.')
      t
   end
