%A1DCEC.m
%Edwin Kite

%A 1-D Catchment Evolution Code

%modified from moundgrowthform.m

%Purpose: Simulate the evolution of a catchment.
%TO DO: Add flexure.

if exist('wwf','var') == 0;            wwf = 5;         end%wall width ratio 
if exist('nx','var') == 0;             nx = 150;        end%number of finite differences 
if exist('nt','var') == 0;             nt = 5e2;        end%2e3; %Max. timestep
if exist('Uo_prime','var') == 0;       Uo_prime = 0.0;  end%Background (synoptic) wind speed as a fraction of maximum wind speed
if exist('alpha','var') == 0;          alpha = 3;       end%1;     end%Exponent in expression dz/dt \propto U^alpha, where U is wind speed.
if exist('RoverL','var') == 0;         RoverL = 4*(1-2/wwf);    end%Gaussian falloff distance for wind
if exist('d_prime','var') == 0;        d_prime = 0.4;   end%;1.1;   end%Nondimensional sediment accumulation rate
if exist('dt_scale','var') == 0;       dt_scale = 0.05; end%dt adjustment (relative to anticipated (d-e) for a run)
if exist('nplot','var') == 0;          nplot = 0.0;     end%Background (synoptic) wind speed as a fraction of maximum wind speed
if exist('k_e','var') == 0;            k_e = 1e-2;      end %Erodibility coefficient; *does not matter* because d' divides through by k_e
if exist('num_diff','var') == 0;       num_diff = 2.5e-4; end%1e-3; end %Numerical diffusivity (for model stability; not modeling a physical process)
if exist('wallslopeangle','var') == 0; wallslopeangle = 1; end %Wall slope angle in degrees
if exist('oscillated','var') == 0;     oscillated = 0;  end %toggle for sinusoidal oscillations in deposition rate (in time; still uniform in space)
if exist('milankamp','var') == 0;      milankamp = 0.5; end %amplitude of sinusoidal oscillations if present
if exist('dryspell','var') == 0;       dryspell = 0;    end
if exist('picturesque','var') == 0;    picturesque = 0; end
if exist('showvelocities','var') ==0;  showvelocities = 0; end 
if exist('sumerosion','var') ==0;      sumerosion = 0;  end %(Mike Lamb suggestion) - erosion = (sum of upslope and downslope), rather than maximum of the two.

eod = round(nx/4); %end of deposition (by aeolian processes)

current_plot_title=strcat( ...
    'alpha',num2str(alpha),'_R/L_',num2str(RoverL),'_d`_',num2str(d_prime),'_Uo`_',num2str(Uo_prime));%'k_e',num2str(k_e),

L = (1-2/wwf)*(nx/2)/RoverL; %The (1-2/wwf) is so L is defined relative to the floor, not the total size
nplot = nplot + 1;
mdz = 0; %maximum change in z
cfl = 0; %worst value of stability criterion during run
eti = 0; %evolutionary track index

%*********SET BASEMENT TOPOGRAPHY**********
if exist('keepgoing','var') == 1; 
    z(1,:) = z(end,:);
    figure;plot(z(1,:))
    z(2:end,:) = [];
else
    z = zeros(nx,1)';
    
    ws = tan(wallslopeangle*pi/180);%wallslope
    
    z(1,1)  = (nx/wwf)*(ws);
    z(1,nx) = (nx/wwf)*(ws);    
    
    basement = [(nx/wwf)*(ws):-1*(ws):1*(ws),zeros((1-2/wwf)*nx,1)',1*ws:ws:ws*(nx/wwf)]; %will be vector
    basementi = interp1([1:nx],basement,[1:nx-1]+0.5); %basement at intermediate points in Lax-Wendroff scheme

    disp('WARNING: The 2 points at either end of the model domain must stay fixed to define slopes in the Lax-Wendroff scheme')

    z(z<basement) = basement(z<basement);
end

%*********SET PARAMETERS WHICH DEPEND ON BASEMENT TOPOGRAPHY**********

lsw = zeros(nx,nx-1);rsw=zeros(nx,nx-1);

for l = 1:nx
         dist(l,:) = exp(-abs(l - [1.5:1:(nx-0.5)])/L); %Distance from a point to a slope. For integer points
         weight(l,:) = dist(l,:)./sum(dist(l,:)); %Normalize to unity weight
         if l~=1;lsw(l,1:(l-1)) = weight(l,1:(l-1))./sum(weight(l,1:[l-1]));end %Left-slope weights
         if l~=nx;rsw(l,l:end) = weight(l,l:end)./sum(weight(l,l:end));end %Right-slope weights
end

%Initial condition slope:
    icslope = diff(z(1,:));
    for l = 1:nx
        if l~=1
            ls(l) = sum(icslope(1:(l-1)) .* weight(l,1:(l-1)))./sum(weight(l,1:[l-1])); %Left slope, i.e. wind blowing from high x to low x
        end
        if l~=nx
            rs(l) = sum(icslope(l:end) .* weight(l,l:end))./sum(weight(l,l:end)); %Left slope, i.e. wind blowing from high x to low x
        end
            %ls(l) = sum(icslope .*weight(l,:)); %Left slope, i.e. wind blowing from high x to low x
    end
        ls(1) = ls(2); rs(nx) = rs(nx-1); %The first point has no slopes to the left of it, and conversely for the last point.
    
    if sumerosion == 0;    
        icesi = max(abs(ls),abs(rs));%,abs(rs)); %Effective slope *at integer points* (erosion propto effective_slope^alpha)
        Uo = Uo_prime*max(icesi);%Background wind as function of max slope wind
        icei = k_e*(Uo + 0.5*(icesi(1:nx-1)+icesi(2:nx))).^alpha; %erosion rate *at halfstep points*    
        %d = d_prime*max(icei);%max erosion rate, e.g. on crater/canyon walls
    elseif sumerosion == 1
        icesi = sum([abs(ls).^alpha;abs(rs).^alpha]);%,abs(rs)); %Effective slope *at integer points* (erosion propto effective_slope^alpha)
        Uo = Uo_prime*max(max(abs(ls),abs(rs)));%Background wind as function of max slope wind - note this is not a function of the max *summed* wind *erosion*
        icei = k_e*(Uo + 0.5*(icesi(1:nx-1)+icesi(2:nx))); %erosion rate *at halfstep points*    
    end
    

    d = d_prime*mean(icei(z(1,1:nx-1)==0 & z(1,2:nx)==0));%erosion rate on crater floor
    disp(strcat('d=',num2str(d)))
    dt = dt_scale/(1e-4+abs(d-mean(icei(z(1,1:nx-1)==0 & z(1,2:nx)==0))));%Adjust dt so growth and decay are initially slow
    dt = min(dt,1000); %Prevent excessively large dt
    disp(strcat('dt=',num2str(dt)))

    lswi = zeros(nx-1,nx-2);rswi=zeros(nx-1,nx-2);

for l = 1:nx-1
        disti(l,:) = exp(-abs(l + 0.5 - [2:(nx-1)])/L);
        weighti(l,:) = disti(l,:)./sum(disti(l,:)); %Normalize to unity weight
        if l~=1; lswi(l,1:(l-1)) = weighti(l,1:(l-1))./sum(weighti(l,1:[l-1]));end %left slope weights, intermediate (precalc for speed)
        if l~=nx-1;rswi(l,l:end) = weighti(l,l:end)./sum(weighti(l,l:end));end %right slope weights, intermediate
end

d = d .* [zeros(eod,1); ones(nx-eod,1)]'; %Deposition suppressed at x values less than 'end of deposition'
dbase = d;
n = 0;
exitcond = -1;
clear et


dt = dt
%*********MAIN LOOP***********************
while exitcond == -1% -1 --> do not exit
    n = n + 1;
    
    %FOR TESTING ONLY! Deposition-rate oscillations
    if oscillated==1
     d = dbase + milankamp*dbase*cos(2*pi*n/200);
     dtest(n) = d;
    end
    if dryspell==1
       if n>220
           d = dbase.*0;
       end
    end
    
    if rem(n-2,100)==1; %The '-2' is to allow {rsi, lsi} display at small n 
        disp(n);
        if picturesque == 1 %track and plot progress
            hold on
                hold off
                if exist('oldmaxshearspeed','var')==1;plot(oldmaxshearspeed(1:nx/2),'m','LineWidth',2);end
                hold on
                if exist('firstmaxshearspeed','var')==1;plot(firstmaxshearspeed(1:nx/2),'g','LineWidth',2);end
                h1 = plot((rsi(1:nx/2)),'r');              
                h2 = plot((lsi(1:nx/2)),'k');
                h3 = plot(abs(rsi(1:nx/2)),'r','LineWidth',2);
                h4 = plot(abs(lsi(1:nx/2)),'k','LineWidth',2);
                grid on
                h5 = plot(z(end,1:nx/2)./max(z(end,1:nx/2)),'b--');
                set(gca,'FontSize',18)
                legend('Previous max. wind','Initial max. wind','Wind blowing from right', 'Wind blowing from left', ...
                       'Magnitude of right wind','Magnitude of left wind','Normalized topography')
                xlabel('x (no units)','FontSize',18)
                ylabel('shear velocity OR dimensionless elevation','FontSize',18)
                title(current_plot_title,'Interpreter','none')
                oldmaxshearspeed = max(abs(rsi),abs(lsi));
                if exist('firstmaxshearspeed','var')==0; 
                  firstmaxshearspeed = max(abs(rsi),abs(lsi));
                end
                xlim([1,nx/2])
                drawnow expose
                pause
        end
    end
      
    slope = diff(z(n,:));
    
    %Two-step Lax-Wendroff
    %Step 1: Calculate effective slopes at each integer point.
    for l = 1:nx  %slope 1 is between elevations 1 and 2, slope 2 is between 2 and 3, e.t.c.
        
        if l~=1
            ls(l) = sum(slope(1:(l-1)) .* lsw(l,1:(l-1))); %Left slope, i.e. wind blowing from high x to low x
        end
        if l~=nx
            rs(l) = sum(slope(l:end) .* rsw(l,l:end)); %Left slope, i.e. wind blowing from high x to low x
        end
    end
            ls(1) = ls(2); rs(nx) = rs(nx-1); %The first point has no slopes to the left of it, and conversely for the last point.
    
    if sumerosion == 0
        esi = max(abs(ls),abs(rs)); %Effective slope *at integer points* (erosion propto effective_slope^alpha)
        %Calculate flux at *half-step points*
        ei = k_e*(Uo + 0.5*(esi(1:nx-1)+esi(2:nx))).^alpha; %erosion rate *at halfstep points*    
    elseif sumerosion == 1        
        esi = sum([abs(ls).^alpha;abs(rs).^alpha]);%,abs(rs)); %Effective slope *at integer points* (erosion propto effective_slope^alpha)
        ei = k_e*(Uo + 0.5*(esi(1:nx-1)+esi(2:nx))); %erosion rate *at halfstep points*    
    end
        
        ui = 0.5*(z(n,1:nx-1)+z(n,2:nx)) + dt*(0.5*(d(1:nx-1)+d(2:nx)) - ei);
        ui(ui<basementi) = basementi(ui<basementi); %impose 'no basalt erosion by wind' condition 
    
    %Step 1.5: Calculate fluxes at the new, halfstep points. Notice that 'edge locations' (basalt) are held fixed
    %At the halfstep, we can only define slopes at points 2, 3, 4, ... nx-1.
    slopei = diff(ui); %([[z(n,1)-ui(1),diff(ui),z(n,end)-ui(end)]);
    
    %Using the slopes defined with these integers, define 'effective slopes' IN BETWEEN the *halfstep* points (1.5, 2.5,
    %e.t.c.)
    for l = 2:nx-2
        lsi(l) = sum(slopei(1:(l-1)) .* lswi(l,1:(l-1))); %Left slope, i.e. wind blowing from high x to low x
        rsi(l) = sum(slopei(l:end) .* rswi(l,l:end)); %Left slope, i.e. wind blowing from high x to low x        
     end
        l=nx-1; lsi(l) = sum(slopei(1:(l-1)) .* lswi(l,1:(l-1))); %Left slope, i.e. wind blowing from high x to low x
        l=1; rsi(l) = sum(slopei(l:end) .* rswi(l,l:end)); %Left slope, i.e. wind blowing from high x to low x
        lsi(1) = lsi(2); rsi(nx-1) = rsi(nx-2); %The first point has no slopes to the left of it, and conversely for the last point.
    
    if sumerosion == 0
        es = max(abs(lsi),abs(rsi)); 
        %Calculate flux at *integer points* (only defined for 2:nx-1)
        %Note that es(1) maps to x=1.5; es(end) maps to x=99.5
        e = k_e*(Uo + 0.5*(es(1:nx-2)+es(2:nx-1))).^alpha; %erosion rate *at integer points*    
    elseif sumerosion == 1
        es = sum([abs(lsi).^alpha;abs(rsi).^alpha]);%,abs(rs)); %Effective slope *at integer points* (erosion propto effective_slope^alpha)
        e = k_e*(Uo + 0.5*(es(1:nx-2)+es(2:nx-1))); %erosion rate *at halfstep points*    
    end
    
    %Step 2: Using erosive flux calculated at integer points, calculate updated values        
    
    nd = num_diff*((z(n,1:nx-2)-basement(1:nx-2)) ...
               - 2*(z(n,2:nx-1)-basement(2:nx-1)) ...
                  +(z(n,3:nx)-basement(3:nx))); %numerical diffusivity calculated *for integer points*

    z(n+1,2:nx-1) = z(n,2:nx-1) + dt*(d(2:nx-1) - e  + nd); %Add numerical diffusivity only in the final steps
    z(n+1,1) = z(n); z(n+1,end) = z(n,end);
    
    z(n+1,z(n+1,:)<basement) = basement(z(n+1,:)<basement);
    oldz = z(n+1,:);
    
    a(n+1,2:nx-1) = d(2:nx-1)-e; %net accumulation;
    
    if max(abs(z(n+1,:) - z(n,:))) > 1
        error('Timestep is too long! Large changes in a single timestep.')
    else
        mdz = max(mdz,max(z(n+1,:) - z(n,:)));
        cfl(n) = 0;
    end
    
    %Find edge of moat 
    if rem(n-1,10) == 0;
        eti = eti + 1; %Evolutionary track index
     
        %Exit conditions
        if n>nt;exitcond=1; %max timestep exceeded
        elseif max(abs(z(n+1,:) - z(n,:))) < 1e-8;exitcond=2; disp('Converged'); %converged
        else %no peaks on floor if no accumulation
           [pks,locs] = findpeaks(z(n+1,:));
            peaksonfloor = find(locs>nx/5 & locs<=(nx*4/5));
            if max(size(peaksonfloor))>=2; %Do multiple peaks exist on the crater floor?
                for peaknumber=peaksonfloor(1:end-1)%Is there "exposed basement" between any of them - i.e., are
                %they seperate mounds?
                 if min(z(n+1,locs(peaknumber):locs(peaknumber+1))) > 1e-2%ignoring tiny mounds
                     if oscillated == 0; %Suppress this exit condition when deposition-rate oscillations are turned on
                        exitcond = 3; disp('Multiple seperated mounds in crater interior')
                     end
                 end
                end            
            elseif min(z(n+1,10:nx-9) - basement(10:nx-9)) > 1/2*(max(basement)-min(basement)); exitcond=5; disp('Will very likely fill crater'); end 
            %multiple mounds in moat; further calcs. are excessively sensitive
            %to value of num_diff (probably related to choice of dx; I am using dx=100)
        end
               
    end
    
end

disp({'exit condition: ',exitcond})
disp({'mdz',mdz,'max(cfl)',max(cfl)})

hold on