%Exercise 3_2, Basin analysis
%% Evaluation of beta
% Sonia Scarselli, ETH-Zurich, sonia.scarselli@erdw.ethz.ch
% Calculation of the Stretch Factor from Thermal Subsidence Data
% This program calculates the stretch factor that best fits the thermal
% subsidence data 
...derived from a decompaction and backstripping analysis of a borehole or surface 
    ...stratigraphic section. 

clear;

% Define some initial parameters
y_c     = 35000;                      % Initial crustal thickness in m              [m]
y_l     = 125000;                     % Initial lithospheric thickness in m         [m]
rho_m0  = 3330;                       % Density of the mantle at 0 degrees celcius  [kg/m^3]
rho_c0  = 2800;                       % Density of the crust at 0 degrees celcius   [kg/m^3]
rho_s   = 1000;                       % Density of the infilling material           [kg/m^3]
alpha_v = 3.28e-5;                    % Thermal expansivity                         [1/K]
Tm      = 1333;                       % Temperature of the mantle                   [C]
kappa   = 1e-6;                       % Thermal diffusivity                         [m^2/s]
time_my = [0, 55, 65, 100];           % Time since end of rifting                   [my]
time_s  = time_my*365*24*3600*1e6;    % Time since end of rifing in seconds
sub     = [-.217 1.031 1.251 1.854];  % Subsidence in km

tau     = y_l^2/(pi^2*kappa);
E0      = 4*y_l*rho_m0*alpha_v*Tm/(pi^2*(rho_m0-rho_s));
 
x       = ((1-exp(-time_s/tau)));  %X-axis 
y       = sub*1000;                %Y-axis

plot(x,y,'o-r')                      %Plot the data points
xlabel('1-exp(-t/\tau)')
ylabel('Thermal subsidence [m]')

% Calculate the best-fit through the data points by using the fitting tool of matlab (see the figure window -> Tools -> basic fitting)
% y = p1*x + p2    -> best fit gives p1 = 2218.6 (slope)

slope_bestfit = 2218.6;

% We know that the slope of the best-fit line through the points is given by the formula 
...slope = E0*(beta/pi)*sin(pi/beta).We know E0 and need to find beta.

% Two methods are provided to calculate the stretch factor. The first is simply by manual trial and error, 
...where you keep modifying your estimate of the stretch factor until the linear best fit slope is correct. 

% The second approach is doing this with a computer. This is called iterations and is shown below.

% Approach 1, do it by hand
beta  = .2;
slope = E0*beta/pi*sin(pi/beta)  %if slope is not equal to slope_bestfit, modify beta

% Approach 2, let the computer do the work
beta            = 1
dbeta           = .1;
iteration_error = 1;
while abs(iteration_error) > 1e-10     %do it until the difference between the best-fitted slope and the calculated slope is smaller than 1e-10
    
   slope           = E0*beta/pi*sin(pi/beta);
   iteration_error = slope-slope_bestfit;
   
   if iteration_error<0
       beta  = beta + dbeta;
   else
       dbeta = dbeta/2;
       beta  = beta-dbeta;
   end
   
end
beta