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