% mountain glacier model
% IJH 18 October 2020
clear;
set(0,'DefaultTextInterpreter','latex');
set(0,'DefaultTextFontSize',12);
set(0,'DefaultAxesFontSize',12);
fn = mfilename;

% dimensionless parameters
pp.n = 3; % flow-law exponent
pp.a_fun = @(x) 1-x; % net accumulation funtion
pp.dt = 1e-2; % timestep

% steady state
x_st = linspace(0,2,101)';
H_st = ((pp.n+2)*x_st.*(1-x_st/2)).^(1/(pp.n+2));

% initial condition
t = 0;
l = 1;
% sigma = l*10.^(linspace(-6,0,101))'; % exponential spacing
n = 400; dx0 = 1e-7; r = fzero(@(r) 1/dx0-(r.^(n+1)-1)/(r-1),1.1); sigma = l*[cumsum(dx0*r.^(0:n)')]; % geometric spacing
% sigma = [0.1*sigma(1:end-1); linspace(0.1,1,201)']; % geometric spacing + constant spacing
x = sigma;
H = max(0,(pp.n+2)*x.*(l/2-x/2)).^(1/(pp.n+2));
x2 = linspace(-4,0,401)'; % ghost characteristics in x<0
H2 = 0*x2.^0; 
x = [x2; x];
H = [H2; H];
s = zeros(size(x)); 

% solve characteristic equations
ts = t;
xs = x;
Hs = H;
ss = s;
i1 = 1:length(x); % indices of characeteristics
i2 = []; % indices of shocks
while t+pp.dt<3
    
    i2 = find(s); % indices of shocks
    i3 = find(isnan(x)); % indices of terminated characteristics
    i4 = find(x<0); % indices of 'ghost' characteristics
    i1 = setdiff(1:length(x),union(i2,[i3;i4])); % indices of characeteristics
    
    a = pp.a_fun(x).*(x>=0); % accumulation rate
    H(i1) = H(i1) + pp.dt*a(i1); % advance H along characteristics
    H(i2) = H(i2) + pp.dt*a(i2); % advance H along shock
    H(H<=0) = 0; % cap H = 0
    x(i1) = x(i1) + pp.dt*H(i1).^(pp.n+1); % advance characteristics
    x(i2) = x(i2) + pp.dt*(H(i2).^(pp.n+1)-H(s(i2)).^(pp.n+1))/(pp.n+2); % advance shocks
    x(i4) = min(x(i4)+pp.dt*1,0); % advance 'ghost' charatacteristics at speed 1 in x<0
    t = t + pp.dt; % advance time
    
    ii = union(i1,i2); % indices of characteristics and shocks
    tmps = find(x(ii(1:end-1))>x(ii(2:end))); % locate intersections
    % label downstream characteristics of shocks
    if ~isempty(tmps)
        for tmp = flip(tmps')
            if s(ii(tmp+1))>0
                s(ii(tmp)) = s(ii(tmp+1)); % if crashed into a shock, take over that shock
            else 
                s(ii(tmp)) = ii(tmp+1); % if crashed into a characteristic, that becomes downstream characteristic
            end
            x(ii(tmp+1)) = NaN; % terminated characteristics
            s(ii(tmp+1)) = 0;
        end
    end
    
    ts = [ts t];
    xs = [xs x];
    Hs = [Hs H];
    ss = [ss s];
end
% blank out various entries
xs(xs<=0) = NaN; 
xs(Hs<=0) = NaN; 
Hs(Hs<=0) = NaN; 
ss(ss<=0) = NaN;

% plot solutions
figure(1); clf; set(gcf,'Paperpositionmode','auto','Units','centimeters','Position',[2 4 20 8]);
width = 1;
ax(1) = subplot(1,2,1);
    xpts = 1:10:size(xs,1);
    plot(xs(xpts,:)',ts(:)','k-','linewidth',width); % characteristics curves
    hold on; plot(xs'.*ss'.^0,ts(:)','r.'); % shocks
    xlabel('$x$')
    ylabel('$t$');
    xlim([0 3]);
    box off;
    
ax(2) = subplot(1,2,2);
    xpts = 1:10:size(xs,1);
    plot3(xs(xpts,:)',ts(:)',0*Hs(xpts,:)','-','color',0.8*[1 1 1],'linewidth',width); % characteristics curves
    hold on; plot3(xs'.*ss'.^0,ts(:)',0*Hs','r.'); % shocks
    tpts = 1:20:size(ts,2);
    hold on; plot3(xs(:,tpts),ones(length(x),1)*ts(tpts),Hs(:,tpts),'b-','linewidth',width); % solution through time
    xlabel('$x$')
    ylabel('$t$');
    zlabel('$H$');  
    xlim([0 3]);
    ylim([0 3]);
    ax(2).Position = ax(2).Position + [0 0.05 0 0];
    
% axes(ax(2)); view(10,40); print(gcf,'-depsc2',[fn]);
% axes(ax(2)); view(10,40); print(gcf,'-depsc2',[fn,'_fig1']); % l = 3;
% axes(ax(2)); view(10,40); print(gcf,'-depsc2',[fn,'_fig2']); % l = 1;