Monday, August 10, 2015

Matlab Program - 2D flexure profile fitted to depth data across the Karoo basin, South Africa

function [minALL]=SlabFleX_V2_Eastern_BestFit;

% Input: Width of load and densities

clear

%Import observed Whitehill Fm depth data
load 4_Geosoft5_Flex2D_Models_BV06_2_Noheading_WH.csv;

% Assign variables to data
x1=X4_Geosoft5_Flex2D_Models_BV06_2_Noheading_WH(1:132,1)';
y1=X4_Geosoft5_Flex2D_Models_BV06_2_Noheading_WH(1:132,5)';

% Tectonic regions the profile passes over
Capex=150000; % Edge of Cape sediments (Distance)
NAMx=340000; % Edge of Kaapvaal Craton (Distance)
ENDx=630000;% End of obs profile

% Parameters
    % Beginning position, end position and height of load
stop = 85*10^3; % Edge of load to correlate with northern edge of CFB (m)
width=100*10^3; % Load width (m)
start = stop-width;
    % Density
rho1 = 2550; % Infill (kg/m3);
rho3 = 2550; % Beam (kg/m3);
rho4 = 3300; % Mantle (kg/m3);
    % Young's modulus
E = 100*10^9;
    % Thickness
Tstart=5; % Beginning elastic thickness (km)
Tend=100; % End elastic thickness (km)
Tdiv1=5; % Divisions
Tdiv=Tend/Tdiv1;
for T1=Tstart:Tdiv1:Tend;
T = T1*10^3; 
    % Poisson’s ratio
v = 0.25;
    % Height
for Height1 = 2:1:8; 
Height = Height1*10^3;
    % Acceleration due to gravity
g = 9.81;
rho = [rho1,rho3,rho4];
    % Total load
pb = rho(1)*Height*(stop-start);
    % Linear load density
q = pb/(stop-start);
    %length of the profile (m)
x = 0:5000: ENDx;

% Tectonic regions the profile passes over (cont.)
[rSTOP,cSTOP]=find(x==stop);
[rCape,cCape]=find(x==Capex);
[rNAM,cNAM]=find(x==NAMx);
[rEND,cEND]=find(x==ENDx);

cSTOP2=cSTOP+1;
cCape2=cCape+1;
cEND5=cNAM;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%% - Properties of the beams - %%%%%%%%%%%%%%%%%%%%%%%%%%%%

D = zeros(1,1);
for k = 1 : length(E);
D(k) = (E(k)*T(k)^3)/(12*(1-v(k)));
end
    % Lambda 
if (rho(3)-rho(1))>0
la(1) = (((rho(3)-rho(1))*g)/(4*D(1)))^0.25;
else
    la(1) = (((rho(4)-rho(1))*g)/(4*D(1)))^0.25;
end
    % Winkler foundation
if (rho(3)-rho(1))>0
k = (rho(3)-rho(1));
else
   k = (rho(4)-rho(1));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%% - Flexure of first beam - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

y = zeros(1,1); 
for i = 1 : length(x);
if x(i) < start;                              % To the left of the load
   a = start - x(i);
   b = stop - x(i);   
   y(i) = (q/(2*(k)*g)) * (exp(-la(1)*a)*cos(la(1)*a)-exp(-la(1)*b)...
       *cos(la(1)*b));
end
if x(i) == start;
    x(i) = x(i) + 0.001;
end
if x(i) > start && x(i) < stop;                % Under the load
    a = x(i)-start;
    b = stop - x(i);    
    y(i) = (q/(2*(k)*g)) * (2 -exp(-la(1)*a)*cos(la(1)*a)-exp(-la(1)*b)...
        *cos(la(1)*b));
end
if x(i) == stop;
    x(i) = stop + 0.001;
end
if x(i) > stop;                               % To the right of the load
    a = x(i) - start;
    b = x(i) - stop;    
    y(i) = (-q/(2*(k)*g)) * (exp(-la(1)*a)*cos(la(1)*a)-exp(-la(1)*b)...
        *cos(la(1)*b));
end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - Shift - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

y=(-10*y); % why 10? 
if y(1,cEND5)>y1(1,cEND5); % Match two curves at end of profile (END) 
    ydiff200(T1./5,Height1) = y(1,cEND5)-y1(1,cEND5);
    ydiff2 = y(1,cEND5)-y1(1,cEND5);
    y=y-ydiff2;
elseif  y(1,cEND5)<y1(1,cEND5);
    ydiff200(T1./5,Height1) = y1(1,cEND5)-y(1,cEND5);
    ydiff2 = y1(1,cEND5)-y(1,cEND5);
    y=y+ydiff2;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - Best-fit  - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Difference between observed and calculated data - KVC
ydiff(1,cSTOP2:cEND)=y1(1,cSTOP2:cEND)-y(1,cSTOP2:cEND);
ydiffT(T1./5,Height1)=sum(abs(ydiff(1,cSTOP2:cEND)));
ydiffTN(T1./5,Height1)=sum(abs(ydiff(1,cSTOP2:cEND))); 

Height2(1,Height1)=Height1;
T2(1,T1./5)=T1;

end;
end;

% Minimum difference for profile (best-fit)
min4TN=min(ydiffTN(1:Tdiv,4));% Height of 4 km
min6TN=min(ydiffTN(1:Tdiv,6));% Height of 6 km
min8TN=min(ydiffTN(1:Tdiv,8));% Height of 8 km

% Determine the parameters that give that minimum difference
[r4TN,c4TN]=find(ydiffTN==min4TN);% Height of 4 km
r4TN2=r4TN*5;
[r6TN,c6TN]=find(ydiffTN==min6TN);% Height of 6 km
r6TN2=r6TN*5;
[r8TN,c8TN]=find(ydiffTN==min8TN);% Height of 8 km
r8TN2=r8TN*5;

diff4TN=ydiff200(r4TN,c4TN); % Height of 4 km
diff6TN=ydiff200(r6TN,c6TN); % Height of 6 km
diff8TN=ydiff200(r8TN,c8TN); % Height of 8 km

minALL=[4 r4TN2 min4TN diff4TN; 6 r6TN2 min6TN diff6TN; 8 r8TN2 min8TN diff8TN];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%% - Displaying results - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1);
hold on; colormap jet;
contour(Height2,T2,ydiffTN,25);
colorbar;
axis square;  ylabel('Te (km)','FontSize',11); xlabel('Height (km)','FontSize',11);

axis([2 8 5 80]);hold off;

No comments:

Post a Comment