function output = starmap_solver(par)
%STARMAP_SOLVER
%   A second order staggered grid finite difference solver for
%   linear hyperbolic moment approximations (such as P_N and
%   SP_N) to the equations of radiative transfer in a 2D
%   geometry. Needs to be called from an example file by
%   STARMAP_SOLVER(PAR), where PAR defines the problem
%   parameters and system matrices.
%   OUTPUT = STARMAP_SOLVER(PAR) writes the final state of the
%   zeroth moment and its (x,y) coordinates into the struct
%   OUTPUT.
%
%   Version 1.9
%   Copyright (c) 03/13/2015 Benjamin Seibold and Martin Frank
%   http://www.math.temple.edu/~seibold
%   http://www.mathcces.rwth-aachen.de/5people/frank/start
%   Contributers: Edgar Olbrant (v1.0), Kerstin Kuepper (v1.5)
%
%   StaRMAP project website:
%   http://math.temple.edu/~seibold/research/starmap/

%========================================================================
% License
%========================================================================
%   Permission is hereby granted, free of charge, to any person obtaining
%   a copy of this software and associated documentation files (the
%   "Software"), to deal in the Software without restriction, including
%   without limitation the rights to use, copy, modify, merge, publish,
%   distribute, sublicense, and/or sell copies of the Software, and to
%   permit persons to whom the Software is furnished to do so, subject to
%   the following conditions:
%   * The above copyright notice and this permission notice shall be
%     included in all copies or substantial portions of the Software.
%   * Appropriate references shall be given to the name of the Software
%     (StaRMAP), the website of the Software
%     (http://www.math.temple.edu/~seibold/research/starmap)
%     and the appropriate research article(s) corresponding to the
%     Software, in particular:
%     B. Seibold, M. Frank, ACM Transactions on Mathematical Software,
%     41(1) [4]. DOI: 10.1145/2590808 (http://arxiv.org/abs/1211.2205)
%
%   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
%   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
%   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
%   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
%   BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
%   ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
%   CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
%   SOFTWARE.

%========================================================================
% Parameter Defaults
%========================================================================
if ~isfield(par,'Mx')||~isfield(par,'My')         % If no system matrices
    par.Mx = [0 1 0;1 0 0;0 0 0]/sqrt(3); par.n_mom = 1;     % specified,
    par.My = par.Mx([1 3 2],[1 3 2]); par.closure = 'P';        % use P1.
end
n_sys = length(par.Mx);                    % Number of system components.
if isfield(par,'mom_order')
    f_mom_order = par.mom_order;          % Create copy of par.mom_order.
    highest_mom = find(f_mom_order==max(f_mom_order)); % Highest moments.
else highest_mom = [];
end
if ~isfield(par,'diffusion'),par.diffusion = 0; end   % Default: no diff.
if par.diffusion&&~isfield(par,'mom_order') % If diff., but no mom_order.
    disp('Diffusion correction requires par.mom_order.')
    par.diffusion = 0;
end
flag_filter = isfield(par,'filterfunction')&&isfield(par,'f_position');
if flag_filter&&~isfield(par,'mom_order')%If filtering, but no mom_order.
    disp('Filtering requires par.mom_order.'), flag_filter = false;
end
if ~flag_filter, par.filterfunction = @ones; par.f_position = ''; end
if isfield(par,'sigma_sm')&&~isfield(par,'mom_order')
    disp('Anisotropic scattering requires par.mom_order.')
end
if par.diffusion                        % If diffusion correction active.
    par.closure = sprintf('D%s%d',par.closure,par.n_mom-1);%Closure name.
    if par.diffusion<0, par.closure = ['c',par.closure]; end % Crescendo.
    dfac = abs(par.diffusion);                     % Acceleration factor.
    if dfac==1, dfac = 100; end          % Default value of acceleration.
    par.Mx(highest_mom,:) = par.Mx(highest_mom,:)*dfac; % Accelerated ad-
    par.My(highest_mom,:) = par.My(highest_mom,:)*dfac;%vection matrices.
else par.closure = sprintf('%s%d',par.closure,par.n_mom); dfac = 1;
end
flag_ani_scatter = isfield(par,'sigma_sm')&&isfield(par,'mom_order');
if ~flag_ani_scatter % If no anisotropic part of scattering kernel given,
    par.mom_order = [0,ones(1,n_sys-1-length(highest_mom)),...   % set up
        ones(1,length(highest_mom))*((par.diffusion~=0)+1)];       % data 
    par.sigma_sm = @zero; % structures for isotropic scattering, possibly
end                  % indexing highest moments for diffusion correction.
if ~isfield(par,'name'),     par.name = 'Output'; end     % Default name.
if ~isfield(par,'ic'),       par.ic = @zero; end% Default: no init. cond.
if ~isfield(par,'sigma_a'),  par.sigma_a = @zero; end%Default: no absorp.
if ~isfield(par,'sigma_s0'), par.sigma_s0 = @zero; end%Defa.: no scatter.
if ~isfield(par,'source'),  par.source = @zero; end % Default: no source.
if ~isfield(par,'source_ind') % If no source moment vector specified, use
    par.source_ind = 1:((nargin(par.source)>=0)+...   % zeroth moment for
        (nargin(par.source)>3)*(n_sys-1));   % isotropic, and all moments
end                                             % for anisotropic source.
if ~isfield(par,'ax'),       par.ax = [0 1 0 1]; end    % Problem domain.
if ~isfield(par,'n'),        par.n = [100 100]; end%Number of grid cells.
if ~isfield(par,'timestep'), par.timestep = 0; end %0=hyp; 1=par; 2=par0.
if ~isfield(par,'cn'),       par.cn = 0.99; end     % Default CFL number.
if ~isfield(par,'bc'),       par.bc = [0 0]; end  % 0=periodic; 1=extrap.
if ~isfield(par,'t_plot'),   par.t_plot = []; end       % Plotting times.
if ~isfield(par,'tfinal'),   par.tfinal = par.t_plot(end); end
if ~isfield(par,'output'),   par.output = @default_output; end
if ~isfield(par,'mom_output'), par.mom_output = 1; end % Moments plotted.
if strcmp(par.mom_output,'all'), par.mom_output = 1:n_sys; end

%========================================================================
% Initialization
%========================================================================
% Create index matrices and vectors for grid type
Ix = cellfun(@find,num2cell(par.Mx',1),'Un',0);         % Grid dependency
Iy = cellfun(@find,num2cell(par.My',1),'Un',0);                % indices.
c11 = 1; c22 = []; c21 = []; c12 = [];% First component is on (1,1)-grid.
while length([c11,c22,c21,c12])<n_sys                         % Determine
    c21 = unique(vertcat(c21',Ix{c11},Iy{c22})');            % indices of
    c12 = unique(vertcat(c12',Iy{c11},Ix{c22})');            % grid types
    c11 = unique(vertcat(c11',Ix{c21},Iy{c12})');           % from system
    c22 = unique(vertcat(c22',Iy{c21},Ix{c12})');             % matrices.
end
if isempty(c22), c22 = []; end     % Create empty matrix of correct size.
gtx = ones(1,n_sys); gty = gtx;                              % Grid types
gtx([c21,c22]) = 2; gty([c12,c22]) = 2;                     % in x and y.

% Grid size and time step control
h = (par.ax([2 4])-par.ax([1 3]))./par.n;                    % Grid size.
dt_hyp = min(h)/abs(eigs(par.Mx,1,'lm'))/2;  % Max. hyperbolic time step.
M = -1i*(par.Mx/h(1)+par.My/h(2));% Scaled matrices (for par. time step).
time_dependent = [nargin(par.sigma_a)>2,nargin(par.sigma_s0)>2,...
            nargin(par.sigma_sm)>3]&[1 1 flag_ani_scatter]; 

% Create staggered grids and coefficients
x{1} = par.ax(1)+h(1)/2:h(1):par.ax(2)-h(1)/2;              % Cell center
y{1} = par.ax(3)+h(2)/2:h(2):par.ax(4)-h(2)/2;             % coordinates.
x{2} = par.ax(1)+h(1)*(1-par.bc(1)):h(1):par.ax(2);           % Staggered
y{2} = par.ax(3)+h(2)*(1-par.bc(2)):h(2):par.ax(4);        % coordinates.
[X,Y] = cellfun(@ndgrid,[x;x]',[y;y],'Un',0);          % Staggered grids.
n1 = par.n; n2 = size(X{2,2});                % Sizes of staggered grids.
extendx = {[1:par.bc(1),1:n1(1),par.bc(1)*(n1(1)-1)+1],...    % Extension
    [n2(1)*(1:1-par.bc(1)),1:n2(1)]};                       % indices for
extendy = {[1:par.bc(2),1:n1(2),par.bc(2)*(n1(2)-1)+1],...     % boundary
    [n2(2)*(1:1-par.bc(2)),1:n2(2)]};                       % conditions.
sg = [gtx;gty;par.mom_order]';                 % Construct staggered grid 
sg = unique(sg(2:end,:),'rows')';            % indices with moment order.

% Initialize field variables
sS = cell(2,2,max(par.mom_order)); sT = sS; ET = sS; % Scattering fields.
U = cell(1,n_sys); dxU = U; dyU = U; Q = U;   % Unknowns and derivatives.
for j = 1:n_sys                           % If initial condition only two 
    U{j} = X{gtx(j),gty(j)}*0+...                        % arguments, set
        capargs(par.ic,X{gtx(j),gty(j)},Y{gtx(j),gty(j)},j)*...     % all
        (nargin(par.ic)>2||j==1);               % higher moments to zero.
    Q{j} = X{gtx(j),gty(j)}*0;      % Initialize source terms with zeros.
end

% Initialize filtering variables
if flag_filter
    if strcmp(par.f_position,'full')
        filterTerm = cell(2,2,max(f_mom_order));      % Full filter term.
    elseif strcmp(par.f_position,'substep')
        filterTerm = cell(2,2,max(f_mom_order),2); % Substep filter term.
    end
    f_sg = [gtx;gty;f_mom_order]';         % Staggered grid indices with 
    f_sg = unique(f_sg(2:end,:),'rows')';                % moment order.
end

%========================================================================
% Time Loop
%========================================================================
par.t_plot = [par.t_plot par.tfinal inf]; plot_count = 1;
cputime = zeros(1,3);
t = 0; dt = 0; step_count = 0;
while t<par.tfinal      % Loop until final time is reached (or exceeded).
    % Evaluate material parameters and determine maximal time step
    tic
    tleft = (par.tfinal-t)/par.cn;  % Time remaining, scaled with 1/CFL#.
    for k = double(t>0):1      % Loop twice for t=0, otherwise only once.
        if ~any(time_dependent)&&t>0&&~par.diffusion, break, end
        flag_check_stab = k>0&&par.timestep>0&&...% If not first time and
            ~(par.diffusion<0&&t==0); % parab. step, and not cresc. @t=0.
        if flag_check_stab
            % Determine maximum admissible time step
            dt = max(dt,dt_hyp);  % Initial guess at least hyperbolic dt.
            decay = ones(1,n_sys)*min(min(sA{1,1}));  % Smallest sigma_a.
            for j = 2:n_sys        % Smallest total cross-section values.
                decay(j) = min(min(sT{gtx(j),gty(j),par.mom_order(j)}));
            end
            la = growth_factor(par,M,[c11,c22],[c21,c12],decay,dt);
            dt_bracket = [1 1]*dt; la_bracket = [1 1]*la;
            la0 = 1+1e-8;  % Allow eigenv. bound of 1+a bit (->no decay).
            while ~diff(la_bracket<la0)...%While eigs on same side of la0
                &&(dt_bracket(1)<=tleft...    % and (t_final not exceeded
                ||la_bracket(1)>la0)                % or dt(1) unstable).
                ind = (la_bracket(1)<la0)+1;% dt(1) unstable:1; stable:2.
                dt_bracket = dt_bracket(ind)*[1 2]*2^(ind-2); % Scale dt.
                la_bracket = la_bracket([2 1]); % Overwrite new time step
                la_bracket(ind) = growth_factor(par,M,...  % with new ei-
                    [c11,c22],[c21,c12],decay,dt_bracket(ind));%genvalue.
                if min(la_bracket)<1, la0 = 1; end% Lower eig bound to 1.
            end
            while dt_bracket(1)/dt_bracket(2)<0.98...   % Bisection iter-
                &&dt_bracket(1)<=tleft   % ation, if not reaching tfinal.
                dt_mid = sqrt(prod(dt_bracket)); % Geom. mid point of dt.
                la_mid = growth_factor(par,M,...          % Eigenvalue at
                    [c11,c22],[c21,c12],decay,dt_mid);     % midpoint dt.
                ind = 1+(la_mid>=la0); % Index in bracket to be replaced.
                dt_bracket(ind) = dt_mid; la_bracket(ind) = la_mid;
            end
        else dt_bracket = dt_hyp;             % Use hyperbolic time step.
        end
        if k>0                                  % If not very first time.
            dt = min(dt_bracket(1),tleft)*par.cn;      %Actual time step.
            sA0 = sA; sT0 = sT;   % Copies of parameters at current time.
            if ~any(time_dependent)&&~par.diffusion% If autonom. or diff.
                sA1 = sA{1,1}; sT1 = sT; break   % use current parameters
            end                                          % and exit loop.
        end
        % Evaluate material parameters at t+dt
        evaluate = time_dependent|(dt==0&[1 1 flag_ani_scatter]);
        if evaluate(1)                            % Absorption: evaluate, 
            sA = cellfun(@(x,y)capargs(...            % if time-dependent
                par.sigma_a,x,y,t+dt),X,Y,'Un',0);       % or first time.
        end
        if evaluate(2)       % Scattering, isotropic component: evaluate,
            for j = 1:size(sg,2)       % if time-dependent or first time.
                sS{sg(1,j),sg(2,j),sg(3,j)} = capargs(par.sigma_s0,...
                    X{sg(1,j),sg(2,j)},Y{sg(1,j),sg(2,j)},t+dt);
            end
        end
        if evaluate(3)    % Scattering, anisotropic components: evaluate,
            for j = 1:size(sg,2)       % if time-dependent or first time,
                sS{sg(1,j),sg(2,j),sg(3,j)} = ...% if anisotropic at all.
                    sS{sg(1,j),sg(2,j),sg(3,j)}-capargs(par.sigma_sm,...
                    X{sg(1,j),sg(2,j)},Y{sg(1,j),sg(2,j)},sg(3,j),t+dt);
            end
        end
        if any(evaluate)||par.diffusion% If any above evaluated, or diff.
            for j = 1:size(sg,2)                    % sigma_t^m = sigma_a
                sT{sg(1,j),sg(2,j),sg(3,j)} = ... % +sigma_s^0-sigma_s^m.
                    sA{sg(1,j),sg(2,j)}+sS{sg(1,j),sg(2,j),sg(3,j)};
            end
        end
        if par.diffusion               % If diffusion correction closure.
            for j = find(sg(3,:)==max(sg(3,:)))%For highest moment order:
                tm = t+dt+eps+(par.diffusion>0)*1e12;  % accelerate decay
                sT{sg(1,j),sg(2,j),sg(3,j)} = ...  % and apply crescendo.
                    dfac./(expm1div(-sT{sg(1,j),sg(2,j),sg(3,j)}*tm)*tm);
            end
        end
        % Determine material parameters used to conduct time step
        if k>0                                  % If not very first time.
            sA1 = (sA0{1,1}+sA{1,1})/2; sT1 = sT;               % Average
            for j = 1:size(sg,2)                               % material
                sT1{sg(1,j),sg(2,j),sg(3,j)} = ...           % parameters
                    (sT0{sg(1,j),sg(2,j),sg(3,j)}+....             % at t
                    sT{sg(1,j),sg(2,j),sg(3,j)})/2;           % and t+dt.
            end
        end
        if flag_check_stab % If stability investigations were done above.
            % Check whether avg. material parameters are stable with dt
            decay = ones(1,n_sys)*min(min(sA1));      % Smallest sigma_a.
            for j = 2:n_sys        % Smallest total cross-section values.
                decay(j) = min(min(sT1{gtx(j),gty(j),par.mom_order(j)}));
            end
            if ~growth_factor(par,M,...% If step with averaged parameters
                [c11,c22],[c21,c12],decay,dt)<la0;    % is unstable, fall
                sA1 = sA0{1,1}; sT1 = sT0;% back to parameters at time t.
            end
        end
    end
    % Compute decay terms
    if t==0||any(time_dependent)||par.diffusion<0% Recompute decay terms.
        EA = expm1div(-sA1*dt/2);         % Decay term for zeroth moment.
        for j = 1:size(sg,2)
            ET{sg(1,j),sg(2,j),sg(3,j)} = ...         % Decay of sigma_t.
                expm1div(-sT1{sg(1,j),sg(2,j),sg(3,j)}*dt/2);
        end
    end
    % Time increment
    t = t+dt; step_count = step_count+1;
    % Evaluate source
    for j = par.source_ind     % Source: evaluate only active components.
        Q{j} = capargs(par.source,X{gtx(j),gty(j)},...  % Evaluate source
            Y{gtx(j),gty(j)},t-dt/2,j)*...    % at half-time of step, and
            (1+any(j==highest_mom)*(dfac-1)); % scale with accel. factor.
    end
    % Evaluate filter function
    if (t==dt||nargin(par.filterfunction)>5)&&flag_filter  % Evaluate, if
        for j = 1:size(f_sg,2)       % first time step or time-dependent.
            if strcmp(par.f_position,'full')
                filterTerm{f_sg(1,j),f_sg(2,j),f_sg(3,j)} = ...
                    capargs(par.filterfunction,par,dt,f_sg(3,j),...
                    X{f_sg(1,j),f_sg(2,j)},Y{f_sg(1,j),f_sg(2,j)},t);
            elseif strcmp(par.f_position,'substep')
                filterTerm{f_sg(1,j),f_sg(2,j),f_sg(3,j),1} = ...
                    capargs(par.filterfunction,par,dt,f_sg(3,j),...
                    X{f_sg(1,j),f_sg(2,j)},Y{f_sg(1,j),f_sg(2,j)},...
                    t-dt/2);
                filterTerm{f_sg(1,j),f_sg(2,j),f_sg(3,j),2} = ...
                    capargs(par.filterfunction,par,dt,f_sg(3,j),...
                    X{f_sg(1,j),f_sg(2,j)},Y{f_sg(1,j),f_sg(2,j)},t);
            end
        end
    end
    cputime(2) = cputime(2)+toc;
    % Update solution
    U0 = U(par.mom_output);                     % Store current solution.
    k = 1; tic
    for step = [1 2 1]                   % Strang splitting of sub-steps.
        switch step
        case 1                    % Half step with (2,1) and (1,2) grids.
        for j = c11       % Compute update at (2,1) and (1,2) from (1,1).
            dxU{j} = diff(U{j}(extendx{1},:),[],1)/h(1);
            dyU{j} = diff(U{j}(:,extendy{1}),[],2)/h(2);
        end
        for j = c22       % Compute update at (1,2) and (2,1) from (2,2).
            dxU{j} = diff(U{j}(extendx{2},:),[],1)/h(1);
            dyU{j} = diff(U{j}(:,extendy{2}),[],2)/h(2);
        end
        for j = [c21 c12]   % Update components on (2,1) and (1,2) grids.
            W = -sumcell([dxU(Ix{j}),dyU(Iy{j})],...
                [par.Mx(j,Ix{j}),par.My(j,Iy{j})]);
            U{j} = U{j}+dt/2*(W+Q{j}-...
                sT1{gtx(j),gty(j),par.mom_order(j)}.*U{j}).*...
                ET{gtx(j),gty(j),par.mom_order(j)};
            if strcmp(par.f_position,'substep')       % Multiplication by
                U{j} = U{j}.*filterTerm{gtx(j),gty(j),f_mom_order(j),k};
            end                                    % substep filter term.
        end
        case 2               % Two half steps with (1,1) and (2,2) grids.
        for j = c21       % Compute update at (1,1) and (2,2) from (2,1).
            dxU{j} = diff(U{j}(extendx{2},:),[],1)/h(1);
            dyU{j} = diff(U{j}(:,extendy{1}),[],2)/h(2);
        end
        for j = c12       % Compute update at (2,2) and (1,1) from (2,1).
            dxU{j} = diff(U{j}(extendx{1},:),[],1)/h(1);
            dyU{j} = diff(U{j}(:,extendy{2}),[],2)/h(2);
        end
        for j = [c11 c22]   % Update components on (1,1) and (2,2) grids.
            W = -sumcell([dxU(Ix{j}),dyU(Iy{j})],...
                [par.Mx(j,Ix{j}),par.My(j,Iy{j})]);
            for k = 1:2                         % Perform two half-steps.
                if j==1        % Zeroth moment decays by absorption only.
                    U{j} = U{j}+dt/2*(W+Q{j}-sA1.*U{j}).*EA;
                else                  % All other moments decay normally.
                    U{j} = U{j}+dt/2*(W+Q{j}-...
                        sT1{gtx(j),gty(j),par.mom_order(j)}.*U{j}).*...
                        ET{gtx(j),gty(j),par.mom_order(j)};
                    if strcmp(par.f_position,'substep')  % Multiplication
                        U{j} = U{j}.*...        % by substep filter term.
                            filterTerm{gtx(j),gty(j),f_mom_order(j),k};
                    end
                end
            end
        end
        end
    end
    if strcmp(par.f_position,'full') % Multiplication by full filter term 
        for j= 2:n_sys                         % at the end of time step.
            U{j} = U{j}.*filterTerm{gtx(j),gty(j),f_mom_order(j)};
        end
    end
    cputime(1) = cputime(1)+toc;
    % Plotting
    tic
    while t>=par.t_plot(plot_count)-1e-14  % If current time has exceeded
        lambda = (par.t_plot(plot_count)-t+dt)/dt; %plotting time, define
        Uplot = cellfun(@(x,y)(1-lambda)*x+lambda*y,...    % solution via
            U0,U(par.mom_output),'Un',0); % linear interpolation in time.
        if plot_count==length(par.t_plot)-1      % Is final time reached?
            if nargout                   % If yes, save solution at final
                output = struct('x',x(gtx),'y',y(gty),'U',U);     % time.
            end
        else                           % If not, invoke plotting routine.
            xplot = x(gtx(par.mom_output));             % Assign grids at
            yplot = y(gty(par.mom_output));          % outputted moments.
            if length(par.mom_output)==1        % If only a single moment
                xplot = xplot{:}; yplot = yplot{:}; Uplot = Uplot{:};% is
            end                         % plotted, remove cell structure.
            if nargout(par.output)                     % Call output rou,
                par = par.output(par,xplot,yplot,Uplot,plot_count);%tine-
            else                          % allowing for it to modify the
                par.output(par,xplot,yplot,Uplot,plot_count)% struct par.
            end
        end
        plot_count = plot_count+1;
    end
    cputime(3) = cputime(3)+toc;
end
% Output problem overview and statistics
fprintf('%s with %s\n %d moments, %dx%d grid, %0.0f time steps\n',... 
    par.name,par.closure,n_sys,par.n,step_count)           % Display test
cputime = reshape([cputime;cputime/sum(cputime)*1e2],1,[]);   % case info
fprintf(['CPU-times\n advection:%15.2fs%5.0f%%\n ab',... % and CPU times.
    'sorb/scatter/dt:%7.2fs%5.0f%%\n plotting:%16.2fs%5.0f%%\n'],cputime)


%========================================================================
% Technical Functions
%========================================================================
function la = growth_factor(par,M,c_even,c_odd,decay,dt)
% Calculate largest eigenvalue of von-Neumann growth factor matrix.
G = diag(dt*expm1div(-decay*dt/2))*(M+diag(-decay/2));
I = eye(length(decay));
G_e = G; G_e(c_odd,:) = 0; G_e = G_e+I;    %Growth factor of even update.
G_o = G; G_o(c_even,:) = 0; G_o = G_o+I;   % Growth factor of odd update.
G = G_o*G_e*G_e*G_o;            % Growth factor matrix of full time step.
la = max(abs(eig(G)));      % Largest eigenvalue of growth factor matrix.
if par.timestep==2, la = la*exp(decay(1)*dt); end   % Relative stability.

function S = sumcell(A,w)
% Add vector of cells A, weighted with vector w.
S = A{1}*w(1); for j = 2:length(w), S = S+A{j}*w(j); end

function z = capargs(fct,varargin)
% Call function fct with as many arguments as it requires (at least 1),
% and ignore further arguments.
narg = max(nargin(fct),1);
z = fct(varargin{1:narg});

function f = zero(varargin)
% Zero function.
f = zeros(size(varargin{1}));

function y = expm1div(x)
% Function (exp(x)-1)/x that is accurate for x close to zero.
y = 1+x*.5+x.^2/6;
ind = abs(x)>2e-4;
y(ind) = (exp(x(ind))-1)./x(ind);

function default_output(par,x,y,U,step)
% Default plotting routine.
clf, imagesc(x,y,U'), axis xy equal tight
title(sprintf('%s at t = %0.2f',par.name,par.t_plot(step))), drawnow
