% Exercise4_1, Basin Analysis
% Calculation of the Deflection caused by a Distributed Load System
% This program calculates the deflections of an elastic lithosphere from a number of ...
...individual rectangular load blocks, analogous to the loading of the lithosphere by...
    ...a fold-thrust belt. The algorithms follow Jordan (1981). The exercise allows ...
    ...you to see the contributions to the total deflection of the individual load blocks...
    ...It is immediately apparent that when load blocks are widely distributed, their ...
    ...individual deflections mutually interfere.

% The assumptions made in this code are the following:
%    ____ ____
%   |    |    |____
%   |    |    |    |
%   x----x----x----x----x----x-----x-----x-----x
%X= 0    1    2    3    4    5     6     7     8     X-vector describing the coordinates.
%H= 2    2    2    1    0    0     0     0     0     H-vector describing the height of the load.
%
%Sonia Scarselli, ETH-Zentrum, sonia.scarselli@erdw.ethz.ch
clear

% Input  parameters----------------------------------------------------------------------------------------------------
Rho_a      = 3300;            % Density of the asthenosphere               [kg m^3]
Rho_c      = 2400;            % Density of the strata (crust)              [kg m^3]
Rho_air    = 1;               % Density of the air material                [kg m^3]

Te         = 5;               % Effective elastic thickness of lithosphere [km]
E          = 100e9;           % Young's modulus of the crust               [Pa]
v          = 0.25;            % Poisson's ratio of the crust     
D          = E*(Te*1000)^3/(12*(1-v^2));% Flexural rigidity of lithosphere [Nm] --> See the exercise for the equation
g          = 9.81;            % Gravitational acceleration                 [m/s^2]
X          = [0:300];         % Distance in km                             [km]
H          = zeros(size(X));  % Height of the load                         [km]
% ---------------------------------------------------------------------------------------------------------

% Specify the thrust load in km's at this place, It is possible change the Height of the load H
% and evaluate the effect on the deflection.------------------------------------------------------------
%H(1:61)      = 4;        
%H(62:106)    = 2;
H(107:160)    = 5;
%H(161:210)   = 2.5;
%H(120:180)   = 5;
%H(1:61)      = 5;        
%H(62:106)    = 3;
% ---------------------------------------------------------------------------------------------------------


% Performing the calculations--------------------------------
H          = H*1e3;          % Height of load in m
X          = X*1e3;          % Horizontal coordinates in m

w_total    = zeros(size(X)); % Initialise the total deflection
for j=1:length(X)-1
    h            = H(j);                                             % Height of load between x(j) and x(j+1)
    s            = (X(j+1)+X(j))/2;                                  % Distance of center of load to origin (at x=0)
    dX           = (X(j+1)-X(j));                                    % Distance between x(j) and x(j+1)
    a            = dX/2;                                             % Half-width of load
    
    lam          = ((Rho_a-Rho_air)*g/(4*D)).^(1/4);
    w_left       =   h/2.*(Rho_c-Rho_air)./(Rho_a-Rho_air).*(   exp((-lam).*(-X + s -a)).*cos(( lam).*(-X + s -a)) -  exp((-lam).*(-X + s +a)).*cos(( lam).*(-X + s +a)) );
    w_below      =   h/2.*(Rho_c-Rho_air)./(Rho_a-Rho_air).*(2- exp((-lam).*( X - s +a)).*cos(( lam).*( X - s +a)) -  exp((-lam).*(-X + s +a)).*cos(( lam).*(-X + s +a)) );
    w_right      = - h/2.*(Rho_c-Rho_air)./(Rho_a-Rho_air).*(   exp((-lam).*( X - s +a)).*cos(( lam).*( X - s +a)) -  exp((-lam).*( X - s -a)).*cos(( lam).*( X - s -a)) );
    
    % Add all the deflections together, to calculate the deflection for the load at point ind
    w            = zeros(size(X));
    w(1:j-1)     = w_left(1:j-1);
    w(j)         = w_below(j);
    w(j+1:end)   = w_right(j+1:end);
    
    w_total      = w_total + w;
    
end
% ---------------------------------------------------------------------------------------------------------

% Plot the results-----------------------------------------------------------------------------------------
subplot(211)
bar( (X(1:end-1)+X(2:end))./(2*1e3),H(2:end)/1e3,1)
xlabel('Distance from origin [km]')
ylabel('Height of sediment load [km]')

subplot(212)
plot(X/1e3,w_total), axis ij, grid on, hold on

xlabel('Distance from origin [km]')
ylabel('Deflection w [m]')
% ---------------------------------------------------------------------------------------------------------