% numerical cole for the paper:
% Probabilistic identification of debris-flow pathways on
% mountain fans within a stochastic framework
% authors: Schiavo, M. et al. (2024)

clc
close all
clear
% GENERAL INFO

% load you DEM data: columns 1 and 2 for X, Y coordinates, columns 3-7 for
% 5 DEMs
% A = 'one of your DEMs';
% no_DEM=5
% Nx : number of pixels along the x-axis
% Ny : number of pixels along the y-axis
% N: number of available DEMs
% NMC : number of Monte Carlo simulations


% PART 1: synthetic DEM terrains
% prepare your N available DEM surveys as column vectors

% temp='[DEM1 DEM2 ... DEM N]';
stdev3 = std(temp,0,2); stdev3(stdev3<0.05)=0;
for k=1:NMC
    random(1,k)=randi([1 NMC-1])./(0.5.*NMC);
end
variance = stdev.*random; clear random;
terrains = nanmean(temp,2);
% perturbate the mean upon local simulated st.devs to have local Gaussian
% pdfs of elevations
terrains = terrains+variance; clear variance;

% now verify hte numerical stability of local elevations
clear temp
for i=1:size(terrains,2)
    temp=terrains(:,1:i);
    first_moment(:,i)=mean(temp,2);
end
% if the first moment at any pixel is stable within your NMC runs, your
% simulations are fine.

% PART 2: Preferential pathways on stochastic terrains
% for j=1:no_DEM

% create a GRID matrix (you can use the X,Y point data) and call it G
A = [G A];
writematrix(A(:,1:3),'DEM%year.txt');
% repeat the A matrix construction for each of columns 3-7
% eventually uncomment row no.5 and 7 and loop the present code (watch out the
% last row END!)

disp('stochastic writing DEM matrix')
bot = A(:,3);

for i = 1:size(bot,2)
    C = bot(:,i);
    x = G(:,1);
    y_delta = G(:,2);
    z = C(:,1);
    xi = unique(x);
    yi = unique(y_delta);
    [X,Y] = meshgrid(xi,yi);
    Z = reshape(z,size(X));
    mat = [x y_delta z];
end

disp('writing text and tiff outputs')
for i = 1:size(bot,2)
    C = bot(:,i) ;
    filename = strcat('mesh',num2str(i),'.txt') ;
    D = [G(:,1:2) C G(:,3:end)];
    filenameV = strcat('meshV',num2str(i),'.txt') ;
    save(filenameV, 'D', '-ascii');
    filenameVV = strcat('meshV',num2str(i),'.mat') ;
    save(filenameVV, 'D', '-mat');
    D1 = D(:,3);
    D2 = reshape(D1,Ny,Nx);
end

% command line instructions
display('command line instructions');

system('cd /');
system('cd Windows\System32\');
system('cd ..');
system('your path');
system('for /L %i in (1,1,NMC) do (pitremove -z meshV%i.txt -fel pitremove%i.tif)');
system('for /L %i in (1,1,NMC) do (d8flowdir -fel pitremove%i.tif -p d8flow%i.tif -sd8 slope%i.tif)');
system('for /L %i in (1,1,NMC) do (aread8 -p d8flow%i.tif -ad8 d8area%i.tif)');

% managing areas
display('managing areas');
for i = 1:size(bot,2)
    filename=strcat('d8area',num2str(i),'.tif');
    R = importdata(filename);    
    n = find(R < 0);
    R(n)= 0;
    n=R;
    L=cast(R,'double');
    L=reshape(L,77500,1);
    filenameR = strcat('d8area',num2str(i),'.mat');
    save(filenameR ,'L','-mat');
end

% managing slopes
display('managing slopes');

for i = 1:size(bot,2)
    filename=strcat('slope',num2str(i),'.tif');
    R = importdata(filename);    
    n = find(R < 0);
    R(n)= 0;
    n=R;
    L2=cast(R,'double');
    L2=reshape(L2,77500,1);
    filenameR = strcat('slope',num2str(i),'.mat');
    save(filenameR ,'L2','-mat');
end

% writing areas and slopes

for i = 1:size(bot,2)
    P = strcat('d8area',num2str(i),'.mat');
    Q = load(P, 'L');
    Q=Q.L;
    D8AREAS(:,i) = Q;
end

for i = 1:size(bot,2)
    P2 = strcat('slope',num2str(i),'.mat');
    Q2 = load(P2, 'L2');
    Q2=Q2.L2;
    SLOPES(:,i) = Q2;
end

cd('..');

SLOPES2 = SLOPES;
SLOPES2(SLOPES2==0)=1e-07;
D8AREAS2 = D8AREAS;
D8AREAS2(D8AREAS2==0)=1;
% r*A_i=Ksat_i*grad(z_i)*w --> Ksat_i=r*A_i/(grad(z_i)*w)
recharge = 3e-08; % m/s
Ksat = recharge.*D8AREAS2./(SLOPES2.*100);

