

clear
format long

diary deform.txt

% Matlab solution of Deformation exercise

% deformation

vp=6000;
vs=3500;
ro=2700;

% a)

disp(' a) ')

lam=ro*vp^2-2*ro*vs^2;
mu=ro*vs^2;

disp(sprintf(' lam = %g N/m2 ',lam));
disp(sprintf(' mu = %g N/m2 ',mu));


% b)

disp(' b) ');
e11=-.26*1e-6;
e12=-.69*1e-6;
e22= .92*1e-6;
e=[e11 e12; e12 e22];

s11=lam*(e11+e22)+ 2*mu*e11;
s22=lam*(e11+e22)+ 2*mu*e22;
s12=2*mu*e12;
s=[s11 s12; s12 s22];

disp( ' deformation ')
disp(e)
disp(' stress changes (Nm) ')
disp(s)


disp(' c) ')

[v,d]=eig(s);
disp(' Eigenvectors ')
disp(v)
disp(' Eigenvalues ')
diag(d)

ang=atan(v(1,1)/v(1,2))/pi*180;

disp(sprintf(' Angle from North: %g deg ',ang));


disp(' d)  ')

e11=.101*1e-6;
e12=.005*1e-6;
e22=-.02*1e-6;
e=1000*[e11 e12; e12 e22];

s11=lam*(e11+e22)+ 2*mu*e11;
s22=lam*(e11+e22)+ 2*mu*e22;
s12=2*mu*e12;
s=[s11 s12; s12 s22];

disp( ' deformation ')
disp(e)
disp(' stress changes (Nm) ')
disp(s)


disp(' e) ' )

a=[1000 0 0];na=a;
b=[0 1000 0];nb=b;

% before
da=e*a(1:2)'; na(1:2)=na(1:2)+da';
db=e*b(1:2)'; nb(1:2)=nb(1:2)+db';

area=cross(a,b);
narea=cross(na,nb);

disp(sprintf(' Area before : %g ',area(3)))
disp(sprintf(' Area after : %g ',narea(3)))
disp(sprintf(' Difference %g m^2 ',-area(3)+narea(3)))

% alternative

[v,d]=eig(e);
disp(' Eigenvectors ')
disp(v)
disp(' Eigenvalues ')
diag(d)

disp(' Change after 1000 years ')
disp(sprintf(' Area before : %g ',area(3)))
disp(sprintf(' Area after : %g ',(1000*(1+d(1,1)))*(1000*(1+d(2,2)))))
disp(sprintf(' Alt Difference %g m^2 ',-area(3)+(1000*(1+d(1,1)))*(1000*(1+d(2,2)))))

% Landers 

e11=-.26*1e-6;
e12=-.69*1e-6;
e22= .92*1e-6;
e=[e11 e12; e12 e22];

a=[1000 0 0];na=a;
b=[0 1000 0];nb=b;

% before
da=e*a(1:2)'; na(1:2)=na(1:2)+da';
db=e*b(1:2)'; nb(1:2)=nb(1:2)+db';

area=cross(a,b);
narea=cross(na,nb);

disp(sprintf(' Area before Landers: %g ',area(3)))
disp(sprintf(' Area after  Landers: %g ',narea(3)))
disp(sprintf(' Difference %g m^2 ',-area(3)+narea(3)))

% alternative

[v,d]=eig(e);
disp(' Eigenvectors ')
disp(v)
disp(' Eigenvalues ')
diag(d)

disp(sprintf(' Area before Landers: %g ',area(3)))
disp(sprintf(' Area after  Landers: %g ',(1000*(1+d(1,1)))*(1000*(1+d(2,2)))))
disp(sprintf(' Alt Difference %g m^2 ',-area(3)+(1000*(1+d(1,1)))*(1000*(1+d(2,2)))))


diary off
