% Marissa E. Cameron, Bridget Smith-Konter
% Code to calculate Coulomb and Sense of Shear

% change increment at the bottom (current z=z+100 for 100 meter interval)
z=0
% observation depth (meters)
while z<=2000

% change increment at the bottom (current m=m+10 for 10 degree interval)
m=0
% orbital position 
while m<=360


% This is not particularly concise, but it is how I use the filename as variables and to output the resulting file.
namea='m'
name1=num2str(m)
name2='_midptlonlatbetaTTPTPP'
name3='TsTn.dat'

filename=[namea, name1, name2, name3]

% Change coef and H/d here.
coef=0.2;		% coef. friction
Hd=1.0;			% depth averaged parameter

% comptue depth averaged stresses

% load a file
f=load(filename)

% extract long, lat, shear and normal vectors

long=f(:,1);
lat=f(:,2);
Ts=f(:,7);
Tn=-f(:,8);  % make tensile negative, compressional positive

da_Tn=Hd*Tn;
da_Ts=Hd*Ts;

% separate tension and compression regions

for i=1:length(da_Tn)
	if da_Tn(i) > 0	% set up compressional vector segments
	Tn_comp(i)=da_Tn(i);
	Tn_tens(i)=nan;
	elseif da_Tn(i) <= 0	% set up tensional vector of segments
	Tn_tens(i)=da_Tn(i);
	Tn_comp(i)=nan;
	end
end

Tn_tens_max=min(Tn_tens)  % compute max negative of tensile stress, compare to lithostatic

 
% set up density and gravity
% CURRENT VALUES FOR GANYMEDE

rho = 1000;		% kg/m^3
g = 1.428;		% m/s^2

% calculate lithostatic stress

ls=rho*g*z./(10^3)	% kPa


% calculate shear failure stress for compressional segments and for tensile segments < lithostatic
% note that previous program had (-1)*Tn - but sign change has been adopted when stresses loaded in

tsf_comp=coef*(ls + Tn_comp);

for i=1:length(Tn_tens)
if abs(Tn_tens(i)) < ls
tsf_tens(i)=coef*(ls + Tn_tens(i));
else 
tsf_tens(i)=nan;
end

end


% calculate coulomb stress

tc_comp=abs(Ts)-tsf_comp';
tc_tens=abs(Ts)-tsf_tens';

% extract signs of shear stress and coulomb stress?

%tss=sign(da_Ts);
%tcs=sign(tc);

% output coulomb stress and shear
% note that tensile stress negative accounted for when loading file, therefore da_Tn below does not have a sign change applied. In other versions, there may be a negative in front of da_Tn but this is only if the file did not have a negative applied when read in and the above loop has a negative (-1)Tn


% Roundabout way of using filename variables to create new filenames.
output=[long lat da_Ts da_Tn tc_comp tc_tens]';
namex=[namea name1 name2]
namez_a='z'
namez_b=num2str(z)
namez=[namez_a namez_b]
fullname=['failure_', namex, '_', namez, '.out']

fid=fopen(fullname,'w');
% 6.2f describes rounding to different decimals. Look up if you need to change significant figures and rounding scheme.
fprintf(fid,'%6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n',output);
st=fclose(fid);

m=m+10
end

z=z+100
end

