function starmap_ex_lattice_diffusion
%STARMAP_EX_LATTICE_DIFFUSION
%   Example case for STARMAP_SOLVER, a second order staggered
%   grid finite difference solver for linear hyperbolic moment
%   approximations to radiative transfer in 2D geometry.
%
%   Example: Lattice/Checkerboard geometry, as defined in
%   [Brunner, Holloway, JCP 210 (2005) 386-399].
%   Computation up to time t=2.8. This file is an adaptation
%   of the test case STARMAP_EX_LATTICE to compare the
%   accuracy of various types of diffusion correction closures.
%
%   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/

%   For license, see file starmap_solver.m, as published on
%   http://www.math.temple.edu/~seibold/research/starmap

%========================================================================
% Problem Parameters
%========================================================================
par = struct(...
'name','Lattice Test',... % name of example
'closure','P',... % type of closure (can be 'P' or 'SP')
'timestep',1,... % parabolic time step
'sigma_a',@sigma_a,... % absorption coefficient (defined below)
'sigma_s0',@sigma_s0,... % isotropic scattering coefficient (def. below)
'source',@source,... % source term (defined below)
'ax',[0 7 0 7],... % coordinates of computational domain
'n',[1 1]*120,... % numbers of grid cells in each coordinate direction
'bc',[1 1],... % type of boundary cond. (0 = periodic, 1 = extrapolation)
'tfinal',2.8); % final time

% Define test cases
n_mom = [9 1 1 1 1 3 3 3 3 5 5 5 5];  % List of moment orders to compute.
diffusion = [0 0 1 -1 2 0 1 -1 2 0 1 -1 2];% Types of diffusion to apply.
arrange = [1 2 3 4 5;1 6 7 8 9;1 10 11 12 13];  % How to combine figures.
cs = 3.5;                  % Position of cross-sectional evaluation line.
save_figure = 0; % Flag whether outputs should be saved into image files.

%========================================================================
% Moment System Setup, Solver Execution, and Plotting
%========================================================================
switch par.closure                      % Define closure matrix function.
    case  'P', closurefun = 'starmap_closure_pn';
    case 'SP', closurefun = 'starmap_closure_spn';
end

% Computation
U = cell(1,numel(n_mom)); uy = U; closure = U;
for j = 1:numel(n_mom)
    par.n_mom = n_mom(j); par.diffusion = diffusion(j);
    [par.Mx,par.My,par.mom_order] = feval(closurefun,par.n_mom);
    out = starmap_solver(par);                              % Run solver.
    x = out(1).x; y = out(1).y;                      % Grid axis vectors.
    U{j} = out(1).U; % Numerical solution.
    uy{j} = interp2(x,y,U{j}',y*0+cs,y); % Evaluate sol. along line x=cs.
    uy{j} = log10(max(uy{j},1e-20));% Cap from below and apply log scale.
    closure{j} = sprintf('%s%d',par.closure,n_mom(j)-(diffusion(j)~=0));
    if diffusion(j), closure{j} = ['D',closure{j}]; end    % Modfify name
    if diffusion(j)<0, closure{j} = ['c',closure{j}]; end   % of closure.
    if abs(diffusion(j))>1, closure{j} = [closure{j},'acc']; end
end

% Plot individual results as 2d colormaps
if save_figure
    for j = 1:numel(n_mom)                        % Loop over test cases.
        h = figure;
        plot2d(x,y,U{j}), title(closure{j})
        filename = sprintf('fig_lattice_closure_%s',closure{j});
        set(h,'position',[10 50 500 500])
        set(h,'paperpositionmode','auto')
        print(h,'-dpng',filename,'-r580')
        close(h)
    end
else
    for k = 1:size(arrange,1)         % Loop over number of arrangements.
        figure
        cols = ceil(sqrt(size(arrange,2)));% Number of subfigure columns.
        rows = ceil(size(arrange,2)/cols);    % Number of subfigure rows.
        for j = 1:size(arrange,2)%Loop over all plots in the arrangement.
            subplot(rows,cols,j)
            plot2d(x,y,U{arrange(k,j)}), title(closure{arrange(k,j)})
        end
    end
end

% Plot results together in 1d cross-sections
width = [5,2,3,3,2];                                       % Line widths.
color = [.7 .7 .8;0 .2 0;1 0 0;0 0 1;.6 0 .6];             % Line colors.
style = {'-','-','--','-.',':'};                           % Line styles.
for k = 1:size(arrange,1)             % Loop over number of arrangements.
    h = figure;
    for j = 1:size(arrange,2)   % Loop over all plots in the arrangement.
        plot(y,uy{arrange(k,j)},style{j},'color',color(j,:),...
            'linewidth',width(j))
        if j==1, hold on, end
    end
    hold off, axis([0 7 -7 1])
    title(sprintf('Cut at x=%g',cs))
    xlabel('y'), ylabel('log10(\psi_0)')
    legend(closure(arrange(k,:)))
    if save_figure
        filename = sprintf('fig_lattice_comparison_%d',k);
        set(h,'position',[10 50 600 500])
        set(h,'paperpositionmode','auto')
        print(h,'-dpng',filename,'-r580')
        close(h)
    end
end

%========================================================================
% Problem Specific Functions
%========================================================================
function f = sigma_a(x,y)
% Absorption coefficient.
cx = ceil(x); cy = ceil(y);
g = (ceil((cx+cy)/2)*2==(cx+cy)).*(1<cx&cx<7&1<cy&cy-2*abs(cx-4)<4);
f = (1-g)*0+g*10;

function f = sigma_s0(x,y)
% Isotropic scattering coefficient.
cx = ceil(x); cy = ceil(y);
g = (ceil((cx+cy)/2)*2==(cx+cy)).*(1<cx&cx<7&1<cy&cy-2*abs(cx-4)<4);
f = (1-g)*1+g*0;

function f = source(x,y)
% Radiation source (only for zeroth moment).
f = 3<x&x<4&3<y&y<4;

function plot2d(x,y,U)
cax = [-7 0];                    % Colormap range used for log10 scaling.
vneg = cax(1)-diff(cax)/254;        % Value assigned where U is negative.
V = log10(max(U,1e-20));% Cap U s.t. U>0 and use logarithmic color scale.
Vcm = max(V,cax(1));                      % Cap colormap plot from below.
Vcm(U<0) = vneg;              % Assign special value where U is negative.
imagesc(x,y,Vcm'), axis xy equal tight, caxis([vneg cax(2)])
cm = jet(256); cm(1,:) = [1 1 1]*.5; % Change lowest color entry to gray.
colormap(cm)
xlabel('x'), ylabel('y')
