%script to bring the stress analysis together for the example of a fault in
%the South Mountains
%Note this assumes tension positive
%JRA 4/27/09

clear all
figure(1)
clf

%--------------------------------------------------------------------------
%Set up Receiver Fault
%We have strike and dip for the fault and then trend and plunge of the
%slickenline

%Haddad station 142
strike = 340; measured_dip = 32;
slickenlinetrend = 075; slickenlineplunge = 24;

%compute pole to plane and get Normal vector
poletrend = strike-90;
if poletrend<0
    poletrend=poletrend+360;
end
poleplunge=90-measured_dip;

[l,m,n] =plunge_trend_to_dir_cosines(poleplunge,poletrend);
ld1 = -l; md1 = -m; nd1 = cosd(poleplunge);
[dip, dipdir] = dir_cosines_to_plunge_trend(ld1, md1, nd1);
fprintf(1,'dip is %3.1f and dip dir is %3.1f\n',dip,dipdir)

N=[l;m;n];

%now for the slicks and to correct for the fact that they are not exactly
%on the plane as reported this way
%see problem on page 13 and 14 from Ragan Chapter 7

slickenlinepoletrend = slickenlinetrend - 90;
if slickenlinepoletrend <0;
    slickenlinepoletrend=slickenlinepoletrend+360;
end
slickenlinepoleplunge = 0;
[ls,ms,ns] =plunge_trend_to_dir_cosines(slickenlinepoleplunge,slickenlinepoletrend);
slickenlinepole = [ls,ms,ns];

A = cross(slickenlinepole,N');
mag = sqrt(A(1)^2 + A(2)^2 + A(3)^2);
A = A./mag;

[intersectplunge, intersecttrend] = dir_cosines_to_plunge_trend(A(1), A(2), A(3));
[lslick,mslick,nslick] =plunge_trend_to_dir_cosines(slickenlineplunge, slickenlinetrend);
slick=[lslick,mslick,nslick];
angular_difference = acosd(dot(A,slick));
if angular_difference>90
    angular_difference=180-angular_difference;
end

fprintf(1,'Measured slickenline plunge = %3.1f and trend = %3.1f\n', slickenlineplunge, slickenlinetrend)
fprintf(1,'Intersection plunge = %3.1f and trend = %3.1f\n', intersectplunge, intersecttrend)
fprintf(1,'Angular difference = %3.1f\n', angular_difference)

%--------------------------------------------------------------------------
%Set up stress tensor
%assume that the principal stresses are appropriate for normal faulting
%conditions so maximum stress is the vertical stress
sv = -26.7.*12; %assume 26.5 MPa per km and 12 km depth
shmin = sv.*0.1; %assume the 1 direction is the minimum horizontal stress and is 10% sv
shmax = sv.*0.5; %assume the 2 direction is intermediate
S = [shmin 0 0; 
         0 shmax 0; 
         0   0 sv] 
%buildrotationmatrix2(xprimetrend, xprimeplunge, yprimetrend,yprimeplunge,zprimetrend,zprimeplunge, talkandplot)
%Note that the plotting starts in here
R = buildrotationmatrix2(      30,            0,          120,           0,          0,          90, 1);

rotatedS = R'*S*R

%-------------------------------------------------------------------------
%Now resolve the stresses
[T, B, sigma, Ts_mag, Ts] = resolvestresses(rotatedS, N);
Ts_normalized = Ts./Ts_mag;
%how good is the fit?
angular_misfit = acosd(dot(A,Ts_normalized));
if angular_misfit>90
    angular_misfit=180-angular_misfit;
end

%--------------------------------------------------------------------------
%Now for a quick bit on the coulomb
%I am going to switch signs now!
friction = 0.65;
figure(2)
clf
maxstress = -1.*min(min(S));
sigma_n = 0:maxstress;
tau_strength = friction.*sigma_n;
plot(sigma_n,tau_strength,'k-')
axis equal
hold on
plot(-sigma,Ts_mag, 'r.')
xlabel('sigma, MPa')
ylabel('tau, MPa')
cff = Ts_mag+sigma.*friction; %compute Coulomb Failure Function
s = sprintf('CFF = %3.1f MPa\n',cff);
text(-sigma,Ts_mag,s)


%-------------------------------------------------------------------------
%Now report the results
T_mag = sqrt(sum(T.^2));
fprintf(1,'traction magnitude %3.4f\n', T_mag)

%normalize components of T to get its direction cosines
lt=T(1)./T_mag; mt = T(2)./T_mag; nt = T(3)./T_mag;
[traction_plunge, traction_trend] = dir_cosines_to_plunge_trend(lt, mt, nt);
fprintf(1,'traction plunge = %3.1f  trend = %3.1f\n',traction_plunge, traction_trend)

fprintf(1,'normal traction mag %3.2f shear traction mag %3.2f\n', sigma, Ts_mag)
%let's check that the normal and shear are components of the traction
testmag = sqrt(sum(sigma.^2 + Ts_mag.^2));
fprintf(1,'check that components make same length as traction: %3.4f =?= %3.4f\n',testmag, T_mag)
fprintf(1,'Angular misfit between resolved traction and observed slickenline = %3.1f\n', angular_misfit)

%--------------------------------------------------------------------------
%Now let's plot everything else on a stereonet
figure(1)
hold on
plotpointcolor(poleplunge,poletrend,'b.');
greatcircle(dip,dipdir);
plotpointcolor(traction_plunge, traction_trend, 'm.');
B_mag = sqrt(B(1)^2 + B(2)^2 + B(3)^2);
B = B./B_mag;
[B_plunge, B_trend] = dir_cosines_to_plunge_trend(B(1), B(2), B(3));
plotpointcolor(B_plunge, B_trend, 'ms');
[plunge, trend] = dir_cosines_to_plunge_trend(Ts_normalized(1), Ts_normalized(2), Ts_normalized(3));
plotpointcolor(plunge,trend,'md');
plotpointcolor(intersectplunge, intersecttrend, 'g.');
plotpointcolor(slickenlineplunge, slickenlinetrend, 'gx');

legend('primitive', 'sigma1', 'sigma2','sigma3', 'xyz stress1', 'xyz stress2', 'xyz stress3', 'pole', 'fault plane', 'T', 'B', 'Ts', 'slick', 'correct slick', 'Location','BestOutside')



