// A program to decompact sediments and then backstrip well assuming Airy isostasy
clear
n = 10;

// initial parameters ( actual thicknesses from boreholes)
z = zeros(1,n);
z(1) = 252; // Thickness of the Continental Mesozoic, in m
z(2) = 215; // Thickness of the Post Tasselian in m
z(3) = 378; // Dalma
z(4) = 891; // Binem
z(5) = 44; // Tadrart
z(6) = 78; // Akakus
z(7) = 106; // Tannezuft
z(8) = 201; // Mamuniyat
z(9) = 1534; // PreGlacial Cambro-Ordovician
z(10) = 1; // Bottom meter of C-O so that the column starts at zero

// Surface porosity for each formation from Sclater & Christie
fi = zeros(1,n);
fi(1) = 0.50; // surface porosity of the Continental Mesozoics sandstones facies;
fi(2) = 0.51; // surface porosity of the Post Tasselian sandstones facies;
fi(3) = 0.51; // surface porosity of the Dalma formation;
fi(4) = 0.55; // surface porosity of the Binem Formation;
fi(5) = 0.515; // surface porosity of the Tadrart;
fi(6) = 0.55; // surface porosity of the Akakus;
fi(7) = 0.61; // surface porosity of the Tannezuft;
fi(8) = 0.50; // surface porosity of the Mamuniyat;
fi(9) = 0.49; // surface porosity of the PreGlacial Cambro Ordovician
fi(10)= 0.49; // Bottom metre of C-O

// fi-depth coefficent for each formation = compaction coefficient
c = zeros(1,n);
c(1) = 0.000280; // fi-depth coefficent in m^-1 for sandstones of the Continental Mesozoic;
c(2) = 0.000310; // fi-depth coefficent in m^-1 for sandstones of the Post Tasselian;
c(3) = 0.000310; // fi-depth coefficent in m^-1 for the Dalma Fm;
c(4) = 0.000405; // fi-depth coefficent in m^-1 for the Binem formation;
c(5) = 0.000295; // fi-depth coefficent in m^-1 for the Tadrart formation;

c(6) = 0.000420; // fi-depth coefficent in m^-1 for the Akakus formation;
c(7) = 0.000485; // fi-depth coefficent in m^-1 for the Tannezuft shaley sequence;
c(8) = 0.000280; // fi-depth coefficent in m^-1 for the Mamuniyat formation;
c(9) = 0.000270; // fi-depth coefficent in m^-1 for the PreGlacial Cambro Ordovician formation;
c(10)= 0.00027; // Bottom metre

//depth to the top of the formations
y = zeros(1,n);
y(1) = 0; // Continental Mesozoic Sandstones
y(2) = 252; // Post Tasselian Sandstones
y(3) = 467; // Dalma Fm
y(4) = 845; // Binem Fm
y(5) = 1736; // Tadrart Fm
y(6) = 1780; // Akakus Fm
y(7) = 1858; // Tannezuft shales
y(8) = 1964; // Mamuniyat Fm
y(9) = 2164; // Preglacial Cambro-Ordovician sandstones
y(10)= 3699; // Lowest Metre
yBas = 3700;

// density at the surface conditions
rho_w = 1030; // density of the sea water Kg/m-3;
rho_m = 3300; // density of the mantle
rho = zeros(1,n);
rho(1) = 2655; // grain density of the Continental Mesozoic Sandstones
rho(2) = 2660; // grain density of the Post tasselian Sandstones
rho(3) = 2660; // grain density of the Dalma Fm
rho(4) = 2690; // grain density of the Binem Fm
rho(5) = 2655; // grain density of the Tadrart Fm
rho(6) = 2680; // grain density of the Akakus Fm
rho(7) = 2710; // grain density of the Tannezuft Fm
rho(8) = 2655; // grain density of the Mamuniyat Fm
rho(9) = 2650; // grain density of the Cambro-Ordovician sandstones
rho(10) = 2650; // grain density of the bottom metre

//The paleosea-levels compared to the present day sea-level taken form the
//Exxon curves
sl(1) = 205; // The Cretaceous (Continental Mesozoic)
sl(2) = 20; // The Upper Jurrassic (Post Tasselian)
sl(3) = 20; // Mississippian (Dalma Fm)
sl(4) = 105; // Late Devonian (Binem Fm)
sl(5) = 175; // Early-Mid Devonian (Tadrart Fm)
sl(6) = 105; // Mid-Late Silurian (Akakus Fm)
sl(7) = 145; // Early Silurian (Tannezuft Fm)
sl(8) = 140; // Hirnantian (Mamuniyat Rm)
sl(9) = 165; // Cambro-Ordovician
sl(10) = 100; // Start of the Cambrian
//sl = fliplr(sl);
sl = sl(:,$:-1:1);

//The Paleobathymetry or water depth.
Wd(1) = 0; // The Cretaceous (Continental Mesozoic)
Wd(2) = 0; // The Upper Jurrassic (Post Tasselian)
Wd(3) = 20; // Mississippian (Dalma Fm)
Wd(4) = 40; // Late Devonian (Binem Fm)
Wd(5) = 0; // Early-Mid Devonian (Tadrart Fm)
Wd(6) = 60; // Mid-Late Silurian (Akakus Fm)
Wd(7) = 100; // Early Silurian (Tannezuft Fm)
Wd(8) = 50; // Hirnantian (Mamuniyat Rm)
Wd(9) = 5; // Cambro-Ordovician
Wd(10) = 5; // Start of the Cambrian
//Wd = fliplr(Wd);
Wd = Wd(:,$:-1:1); // Scilab flip L-R

x = zeros(n,n); //create an array for the depths to the top of the units over time

for it = 1:n
    eps = 0.0001 ;
    error = 2*eps ;
        x_top = 0;
        
        for in = n-(it-1):n
        x_iter = 1;
        eps = 0.0001 ;
        error = 2*eps ;
while error > eps
    x(in,it) = z(in) - fi(in)/c(in)*(exp(-c(in)*y(in))-exp(-c(in)..
    *(y(in)+z(in)))) + fi(in)/c(in)*(exp(-c(in)*x_top)-exp(-c(in)..
    *x_iter)) + x_top;
    error = abs(x(in,it)-x_iter)/x_iter;
    x_iter = x(in,it);
        end
    x_top = x(in,it); // reset for the next layer
end
end

//xu = flipud(x); // This creates a matrix with time on the x-axis and each layer on the z-axis
xu = x($:-1:1,:); // Scilab flip U-D
t = [513 445.6 443.7 440.8 416 386 359.2 318 223.4 97]; //prepare for the plot
// carry out the backstripping
//calculating the porosity for each layer through time

fi_depth = zeros(n,n);
rho_i = zeros(1,n);
rho_b = zeros(1,n);
tecsub = zeros(1,n);

for it = 1:n
    x_top = 0;
        for in = n-(it-1):n
        fi_depth(in,it) = (fi(in)/c(in)) * ((exp(-c(in)*x_top)-exp(-c(in)*x(in,it)))/(x(in,it) - x_top));
        rho_i(in) = (fi_depth(in,it)*rho_w + (1-fi_depth(in,it))*rho(in))*(x(in,it)-x_top)/x(n,it);
        x_top = x(in,it);
        end
    rho_i;
    rho_b(it) = sum(rho_i);
    tecsub(it) = x(n,it)*((rho_m - rho_b(it))/(rho_m - rho_w))..
     -sl(it)*(rho_w/(rho_m-rho_w)) + (Wd(it)-sl(it));
end
//Plot of the complete output of the program showing both the decompacted
//subsidence and the tectonic subsidence
figure(1), clf
f=gcf(); // macro version of f=get("current_figure")
f.figure_size=[1024,800];
title ('Subsidence of the Kufrah Basin Well A1-NC198')
subplot (2,1,1)

plot (-t,-xu(1,1:n),'b-o',..
      -t,-xu(2,1:n),'b-o',..
      -t,-xu(3,1:n),'b-o',..
      -t,-xu(4,1:n),'k-x',..
      -t,-xu(5,1:n),'g-x',..
      -t,-xu(6,1:n),'r-d',..
      -t,-xu(7,1:n),'r-d',..
      -t,-xu(8,1:n),'c-^',..
      -t,-xu(9,1:n),'y-^',..
      -t,-xu(10,1:n),'y-^') // set colour and symbology for each line
//    -t,-xu(10,1:n),'y^-','linewidth',2) // Matlab original

xtitle ('The Decompacted Sedimetary thickness for the Kufrah basin')
xlabel ('Time (Ma)')
ylabel ('Depth of each uncompacted layer (m)')
//legend ('Cambro-Ordivician','off','off','off','off','off','Black Shale','Silurian','hide','hide','hide','Devonian','hide','hide','hide','hide','Carboniferous','location','SouthWest')
subplot (2,1,2)
plot (-t,-tecsub,'o-','linewidth',2)
// set (gca,'FontSize',14)
xtitle ('Water loaded tectonic subsidence')
xlabel ('Time (Ma)')
ylabel ('The Tectonic subsidence (m)')
//axis ([-550 0 -2000 0])