function u_end = forward_1d(u_start, stps, drw)
% advance the last column of  U_START forward through STPS
%
% INPUTS 
%   U_START: The last column is the set of starting values
%   STPS:    If scaler- the number of steps to move forward, 
%            with time increments of PARAMS.dt
%            If vector:  the specific time steps to run.
% OUTPUTS
%    U_END:  output has all the information in U_START as well 
%            as the new increments computed here.
%
%--------------------------------------------------------------------------
% the state (T,S) of the member is advanced
% by solving numerically the stochastic reaction convection diffusion equation
%
% dT/dt = c1 * div (grad T) + v' * grad T - c2*(T-Ta) - c4*dS/dt + SS
% dS/dt = -f(S,T,Ti)
%
% where
% T is the fire temperature
% S is the fuel supply density
% c1 is the diffusion coefficient
% v incorporates ground slope and wind effect
% dS/dt is the rate of burning
% Ta is the ambient temperature
% c2*(T-Ta) is the rate of the heat escaping from the fire
% f(S,T,Ti)=S.*max((T/Ti).^3-1,0)
%           is the burn intensity (I do not know any better)
% c4 is the specific heat generated by burning
% SS are additional stochastic terms described in Mandel et al., ICCS04
% 
% to initialize ensemble member correctly, run fire_init.m
%
% note: this is a very simplified model, but it ehibits a nice moving flame
% front behavior
%
% method: the equation is written as
% u'=F(u,t)+R(u,t)
% where R are random terms.
% Time is advanced by the Crank-Nicholson method
% u{n+1}=u{n}+0.5*dt*(F(u{n+1},t{n+1})+F(u{n},t{n}))+R
% where u=[T;S]
% The new value n{n+1} is found using the quasi-Newton method
% u{n+1} <- u{n+1} - J'(u{n},t{n})\J(u{n},u{n+1}.t{n},t{n+1})
%
% Modified:  Craig Johns, August 2004
% Based on:  Jan Mandel,  August 2004 
%
% import variables from global PARAMS
global PARAMS 

if(nargin==2)
  drw=1;
end

% state
dm = size(u_start);
if(length(stps)==1); 
    t0   = dm(2);
    stps = t0 + (1:stps); 
elseif(length(stps)>1)
    t0   = min(stps)-1;
    stps = t0 + 1 + stps-min(stps);
end

% params
dt   = PARAMS.dt;
Smax = PARAMS.Smax;
Ti   = PARAMS.Ti;
m    = PARAMS.mesh.m;

% setup
u = u_start(:,end);
T = u_start(1:m,end);              % temperature
S = u_start((m+1):(2*m),end);      % fuel supply

% output
u_end = [u_start, zeros(dm(1),length(stps))];





kk =0;    
for k = stps;
    k;
    kk   = kk+1;
    Sold = S;
    u    = odetrap(@F1d,dt*(t0+kk-1),dt,u); 
    u_end(:,k) = u;
    T    = u(1:m);
    S    = u(m+1:2*m);
    %warnings = 0;
    %if warnings
    if any(S<0),
        warning('S negative')
    end
    if any (S>Sold)
       warning('S increased')
    end
    if any(T<0 | Ti >Ti *4),
        warning('T out of range')
    end
    %end
    ns=0;
    if(drw)
    plot(1:m,T./Ti,'r',1:m,S,'g');
       title(sprintf('Time = %1.3g',dt*k))
       axis([1,m,-.05,7]);
       pause(.001);
    end
end





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v=F1d(t,u,mn)
% the function F for fire_advance_1d
% u     
global PARAMS
% import dimensions and parameters
m  = PARAMS.mesh.m;    % meshsize in x direction
c1 = PARAMS.c1;        % diffusion coefficient 
c2 = PARAMS.c2;        % heat escape coefficient
c3 = PARAMS.c3;        % burning speed coefficient
Ti = PARAMS.Ti;        % ignition temperature, for burn rate and  display scale
Ta = PARAMS.Ta;        % ambient temperature
c4 = PARAMS.c4;        % specific heat coefficient
dx = PARAMS.mesh.dx;         % meshstep in the x direction
vx = PARAMS.mesh.vx;         % the wind/slope vector


% split state 
if length(u) ~= 2*m,
    error('bad state length')
end
T = u(1:m);
S = u(m+1:2*m);
v = zeros(2*m,1);

dm = 2:m-1; % interior mode where the difference formulas are evaluated

% dS/dt = -c3*S.*max((T/Ti).^3-1,0)
burn = c3.*S(dm).*max((T(dm)/Ti).^3-1,0); 

% dT/dt = c1 * div (grad T) + v' * grad T - c2*(T-Ta)+ c4*dS/dt
diffusion  = c1*(1/(dx*dx)).*(T(dm-1)-2*T(dm)+T(dm+1));
convection = vx(dm).*(0.5/dx).*(T(dm-1)-2*T(dm+1));
outflow    = -c2*(T(dm)-Ta);
heat       = diffusion + convection + outflow + c4*burn;

v(dm)   = heat;
v(m+dm) = -burn; % insert in output