Monday, August 10, 2015

Matlab Program - Geodetic response due to a tensile fault (i.e., Iceland Mi dAtlantic Ridge)

function [Ux,Uy,Uz,Ur]=Dyke_Model_red(extension,depth);

% Program Written by: Stephanie Scheiber
% Date: 2007-2008
% Reference: Okada, Y. (1985). Surface Deformation due to Shear and
% Tensile Faults in a Half-Space. Bull. Seis. Soc. Amer. 75, 4, 1135-
% 1154pp.
% Description: Models vertical and horizontal ground deformation on the
% surface of the Earth surface (2-D grid) due to the implacement of a
% dyke. Locking depth, length, width, extension and dip of dyke may be
% varied.


% User inputted values
L=1000; % Length of dyke (km). Large value is used to reduce edge effects
W=1000; % Width of dyke (km). Large value is used to reduce edge effects
d=1000+depth; % Depth to the bottom of dyke (km). "depth" is the locking depth. Large value is used to reduce edge effects
extens=extension; % Dyke-normal rate  (mm/yr)
dip=89.9; % Dip of dyke (degrees). Cannot be vertical as causes division by zero and therefore an undefined solution
rdist2=500; % Length of grid area (km)
sspace=1; % Grid spacing(km)
mu=1; % Medium constant assumed to be equal to one
lamda=1; % Medium constant assumed to be equal to one

% Converting from kilometres to metres
L = L*10^3;            
W = W*10^3;
d = d*10^3;
extens = extens/10^3;
rdist2 = rdist2*10^3;
sspace = sspace*10^3;
% Converting from degrees to radians
dip2 = dip*pi/180;        

% MESHGRID generates X and Y matrices for the grid area. Here the
% x and y values for one profile through the grid area (over the
% centre of the dyke) are generated
% in order to save computation time
[x,y]=meshgrid(1,(-1*rdist2/2):sspace:(rdist2/2)+1); 

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

% NEAR VERTICAL DYKE - Equations from Okada (1985) for a tensile fault:
    
% For E=x-L and n=p (Chinnery's notation)
Ea=x-L;
pa=(y*cos(dip2))+(d*sin(dip2));
na=pa;                                 
qa=(y*sin(dip2))-(d*cos(dip2));
Ya=(na*cos(dip2))+(qa*sin(dip2));
Da=(na*sin(dip2))-(qa*cos(dip2));
Ra=sqrt((Ea.*Ea)+(na.*na)+(qa.*qa));
Xa=sqrt((Ea.*Ea)+(qa.*qa));

const1a=(mu/(lamda+mu));
denoma=Ra+Da;
num4a=1/cos(dip2);
num5a=sin(dip2)/cos(dip2);
num6a=log(Ra+na);
num7a=log(Ra+Da);
num8a=na.*(Xa+(qa*cos(dip2)));
num9a=Xa.*(Ra+Xa)*sin(dip2);
num10a=(Ea.*(Ra+Xa)*cos(dip2))+0.000000000000000001;
I4a=const1a*num4a*(num7a-(sin(dip2)*num6a));
I5a=const1a*2*num4a*atan((num8a+num9a)./(num10a));
I1a=((-1*const1a*num4a*Ea)./denoma)-(num5a*I5a);
I3a=(const1a*(((num4a*Ya)./denoma)-num6a))+(num5a*I4a);

const2a=extens/(2*pi);
num1a=qa./(Ra.*(Ra+na));
num2a=qa./(Ra.*(Ra+Ea));
num3a=(Ea.*na)./(qa.*Ra);
Uxa=const2a*((qa.*num1a)-(I3a*((sin(dip2))^2)));
Uya=const2a*((-Da.*num2a)-(sin(dip2)*((Ea.*num1a)-atan(num3a)))-(I1a*((sin(dip2))^2)));
Uza=const2a*((Ya.*num2a)+(cos(dip2)*((Ea.*num1a)-atan(num3a)))-(I5a*((sin(dip2))^2)));

% For E=x-L and n=p-W (Chinnery's notation)
Eb=x-L;
pb=(y*cos(dip2))+(d*sin(dip2));
nb=pb-W;                                 
qb=(y*sin(dip2))-(d*cos(dip2));
Yb=(nb*cos(dip2))+(qb*sin(dip2));
Db=(nb*sin(dip2))-(qb*cos(dip2));
Rb=sqrt((Eb.*Eb)+(nb.*nb)+(qb.*qb));
Xb=sqrt((Eb.*Eb)+(qb.*qb));

const1b=(mu/(lamda+mu));
denomb=Rb+Db;
num4b=1/cos(dip2);
num5b=sin(dip2)/cos(dip2);
num6b=log(Rb+nb);
num7b=log(Rb+Db);
num8b=nb.*(Xb+(qb*cos(dip2)));
num9b=Xb.*(Rb+Xb)*sin(dip2);
num10b=(Eb.*(Rb+Xb)*cos(dip2))+0.000000000000000001;
I4b=const1b*num4b*(num7b-(sin(dip2)*num6b));
I5b=const1b*2*num4b*atan((num8b+num9b)./(num10b));
I1b=((-1*const1b*num4b*Eb)./denomb)-(num5b*I5b);
I3b=(const1b*(((num4b*Yb)./denomb)-num6b))+(num5b*I4b);

const2b=extens/(2*pi);
num1b=qb./(Rb.*(Rb+nb));
num2b=qb./(Rb.*(Rb+Eb));
num3b=(Eb.*nb)./(qb.*Rb);
Uxb=const2b*((qb.*num1b)-(I3b*((sin(dip2))^2)));
Uyb=const2b*((-Db.*num2b)-(sin(dip2)*((Eb.*num1b)-atan(num3b)))-(I1b*((sin(dip2))^2)));
Uzb=const2b*((Yb.*num2b)+(cos(dip2)*((Eb.*num1b)-atan(num3b)))-(I5b*((sin(dip2))^2)));

% For E=x+L and n=p (Chinnery's notation)
Ec=x+L;
pc=(y*cos(dip2))+(d*sin(dip2));
nc=pc;                                 
qc=(y*sin(dip2))-(d*cos(dip2));
Yc=(nc*cos(dip2))+(qc*sin(dip2));
Dc=(nc*sin(dip2))-(qc*cos(dip2));
Rc=sqrt((Ec.*Ec)+(nc.*nc)+(qc.*qc));
Xc=sqrt((Ec.*Ec)+(qc.*qc));
const1c=(mu/(lamda+mu));
denomc=Rc+Dc;
num4c=1/cos(dip2);
num5c=sin(dip2)/cos(dip2);
num6c=log(Rc+nc);
num7c=log(Rc+Dc);
num8c=nc.*(Xc+(qc*cos(dip2)));
num9c=Xc.*(Rc+Xc)*sin(dip2);
num10c=(Ec.*(Rc+Xc)*cos(dip2))+0.000000000000000001;
I4c=const1c*num4c*(num7c-(sin(dip2)*num6c));
I5c=const1c*2*num4c*atan((num8c+num9c)./(num10c));
I1c=((-1*const1c*num4c*Ec)./denomc)-(num5c*I5c);
I3c=(const1c*(((num4c*Yc)./denomc)-num6c))+(num5c*I4c);

const2c=extens/(2*pi);
num1c=qc./(Rc.*(Rc+nc));
num2c=qc./(Rc.*(Rc+Ec));
num3c=(Ec.*nc)./(qc.*Rc);
Uxc=const2c*((qc.*num1c)-(I3c*((sin(dip2))^2)));
Uyc=const2c*((-Dc.*num2c)-(sin(dip2)*((Ec.*num1c)-atan(num3c)))-(I1c*((sin(dip2))^2)));
Uzc=const2c*((Yc.*num2c)+(cos(dip2)*((Ec.*num1c)-atan(num3c)))-(I5c*((sin(dip2))^2)));

% For E=x+L and n=p-W (Chinnery's notation)
Ed=x+L;
pd=(y*cos(dip2))+(d*sin(dip2));
nd=pd-W;                                 
qd=(y*sin(dip2))-(d*cos(dip2));
Yd=(nd*cos(dip2))+(qd*sin(dip2));
Dd=(nd*sin(dip2))-(qd*cos(dip2));
Rd=sqrt((Ed.*Ed)+(nd.*nd)+(qd.*qd));
Xd=sqrt((Ed.*Ed)+(qd.*qd));

const1d=(mu/(lamda+mu));
denomd=Rd+Dd;
num4d=1/cos(dip2);
num5d=sin(dip2)/cos(dip2);
num6d=log(Rd+nd);
num7d=log(Rd+Dd);
num8d=nd.*(Xd+(qd*cos(dip2)));
num9d=Xd.*(Rd+Xd)*sin(dip2);
num10d=(Ed.*(Rd+Xd)*cos(dip2))+0.000000000000000001;
I4d=const1d*num4d*(num7d-(sin(dip2)*num6d));
I5d=const1d*2*num4d*atan((num8d+num9d)./(num10d));
I1d=((-1*const1d*num4d*Ed)./denomd)-(num5d*I5d);
I3d=(const1d*(((num4d*Yd)./denomd)-num6d))+(num5d*I4d);

const2d=extens/(2*pi);
num1d=qd./(Rd.*(Rd+nd));
num2d=qd./(Rd.*(Rd+Ed));
num3d=(Ed.*nd)./(qd.*Rd);
Uxd=const2d*((qd.*num1d)-(I3d*((sin(dip2))^2)));
Uyd=const2d*((-Dd.*num2d)-(sin(dip2)*((Ed.*num1d)-atan(num3d)))-(I1d*((sin(dip2))^2)));
Uzd=const2d*((Yd.*num2d)+(cos(dip2)*((Ed.*num1d)-atan(num3d)))-(I5d*((sin(dip2))^2)));

% Calculating the Total Vertical (Uz) and Horizontal (Ur) Displacement
% along the profile
Ux=zeros(501,1); % Ux component is zero along the profile
Uy=Uyc-Uyd-Uya+Uyb;
Uz=Uzc-Uzd-Uza+Uzb;
Ur=sqrt((Ux.*Ux)+(Uy.*Uy));



APPENDIX 3:   Response due to Hekla Mogi source
between 2001 and 2006

function [TOTAL0106]=Mogi_Hekla_01_06(VCh01_06,dpth01_06,q);

% Program written by: Stephanie Scheiber
% Date: 2007-2008
% Input and output values: Volume change (VCh01_06) and depth (dpth01_06)
% of Hekla source; Time period over which inflation/deflation of Mogi
% source has occurred (q); Rate of deformation at each site (TOTAL0106,
% horizontal deformation is the component parallel to spreading (N78W))
% Reference: Mogi, K. (1958). "Relations between the eruptions of various
% volcanoes and the deformations of the ground surfaces around them."Bull % of Earthquake Research Institute 36: 99-134.
% Description: Horizontal and vertical displacements and rates (for 2001-
% 2006) over a Mogi source

% User inputted values
VChange=VCh01_06; % Volume Change (km3)
depth=dpth01_06; % Chamber depth (km)
rdist=108; % Length of grid area (km)
sspace=200; % Grid spacing (km)

% Converting from kilometres to metres
depth2 = depth*10^3;
VChange2 = VChange*10^9;
rdist2 = rdist*10^3;

% MESHGRID generates X and Y matrices for the grid area
[x,y]=meshgrid((-1*rdist2/2):sspace:(rdist2/2));
                                 
% Radial Distance
r = sqrt(x.^2+y.^2);

% Vertical displacement due to a Mogi source
const = 3/(4*pi);
Uz = (const*VChange2*depth2)./((depth2^2+r.^2).^(3/2));

% Horizontal displacement due to a Mogi source
Ur = (const*VChange2*r)./((depth2^2+r.^2).^(3/2));

% Site and Mogi source locations in WGS84, UTM Zone 27N
mogix=565000.4372-000;
mogiy=7096999.581-000;

ISAKx=560874.8819;
KROKx=578057.3671;
LIHOx=585129.5376;
BLAUx=580352.5266;
LAUFx=580005.0291;
HRAFx=587179.4827;
HRASx=588221.7342;
DOMAx=591324.9819;
FROSx=595822.0785;
LAHRx=594912.3594;
SATUx=585955.0103;
THRAx=588598.9273;
KKLOx=595844.3353;
KGILx=599775.4237;

ISAKy=7110983.683;
KROKy=7105794.749;
LIHOy=7097061.677;
BLAUy=7092475.799;
LAUFy=7086893.04;
HRAFy=7093773.73;
HRASy=7091573.266; 
DOMAy=7102811.594;
FROSy=7099601.906;
LAHRy=7097344.121;
SATUy=7084821.087;
THRAy=7078205.599;
KKLOy=7082878.089;
KGILy=7083000.709;

% Determining the calculated displacement at each site due to the Mogi
% source deformation.

%UTM co-ordinates are transformed into grid co-ordinates (column and row values).
ISAKcoly=int16(271-((ISAKy-mogiy)/200));
ISAKrowx=int16(271-((ISAKx-mogix)/200));
KROKcoly=int16(271-((KROKy-mogiy)/200));
KROKrowx=int16(271-((KROKx-mogix)/200));
FROScoly=int16(271-((FROSy-mogiy)/200));
FROSrowx=int16(271-((FROSx-mogix)/200));
LAHRcoly=int16(271-((LAHRy-mogiy)/200));
LAHRrowx=int16(271-((LAHRx-mogix)/200));
SAcoly=int16(271-((SATUy-mogiy)/200));
SArowx=int16(271-((SATUx-mogix)/200));
TAcoly=int16(271-((THRAy-mogiy)/200));
TArowx=int16(271-((THRAx-mogix)/200));
KKcoly=int16(271-((KKLOy-mogiy)/200));
KKrowx=int16(271-((KKLOx-mogix)/200));
KGcoly=int16(271-((KGILy-mogiy)/200));
KGrowx=int16(271-((KGILx-mogix)/200));
Lcoly=int16(271-((LIHOy-mogiy)/200));
Lrowx=int16(271-((LIHOx-mogix)/200));
Bcoly=int16(271-((BLAUy-mogiy)/200));
Browx=int16(271-((BLAUx-mogix)/200));
LFcoly=int16(271-((LAUFy-mogiy)/200));
LFrowx=int16(271-((LAUFx-mogix)/200));
HFcoly=int16(271-((HRAFy-mogiy)/200));
HFrowx=int16(271-((HRAFx-mogix)/200));
HScoly=int16(271-((HRASy-mogiy)/200));
HSrowx=int16(271-((HRASx-mogix)/200));
Dcoly=int16(271-((DOMAy-mogiy)/200));
Drowx=int16(271-((DOMAx-mogix)/200));

% Angle between the site and Mogi source relative to an east-west axis
ISAKangle=atan((ISAKy-mogiy)/(ISAKx-mogix))*180/pi;
KROKangle=atan((KROKy-mogiy)/(KROKx-mogix))*180/pi;
FROSangle=atan((FROSy-mogiy)/(FROSx-mogix))*180/pi;
LAHRangle=atan((LAHRy-mogiy)/(LAHRx-mogix))*180/pi;
SAangle=atan((SATUy-mogiy)/(SATUx-mogix))*180/pi;
TAangle=atan((THRAy-mogiy)/(THRAx-mogix))*180/pi;
KKangle=atan((KKLOy-mogiy)/(KKLOx-mogix))*180/pi;
KGangle=atan((KGILy-mogiy)/(KGILx-mogix))*180/pi;
Langle=atan((LIHOy-mogiy)/(LIHOx-mogix))*180/pi;
Bangle=atan((BLAUy-mogiy)/(BLAUx-mogix))*180/pi;
LAangle=atan((LAUFy-mogiy)/(LAUFx-mogix))*180/pi;
HFangle=atan((HRAFy-mogiy)/(HRAFx-mogix))*180/pi;
HSangle=atan((HRASy-mogiy)/(HRASx-mogix))*180/pi;
Dangle=atan((DOMAy-mogiy)/(DOMAx-mogix))*180/pi;

% Determining the exact location of the site relative to the Mogi source in
% order to determine the signs (positive or negative) of the x and y
% displacement components for each site. Also determining the displacement
% experienced at each site due to Mogi deformation [in m]
if mogix>=560874.8819 && mogiy<=7110983.683;
    ISAK=Ur(ISAKrowx,ISAKcoly);
elseif mogix<=560874.8819 && mogiy>=7110983.683;
    ISAK=Ur(ISAKrowx,ISAKcoly)*-1;
elseif mogix>=560874.8819 && mogiy>=7110983.683;
    ISAK=Ur(ISAKrowx,ISAKcoly);
else
    ISAK=Ur(ISAKrowx,ISAKcoly)*-1;
end;
if mogix>=578057.3671 && mogiy<=7105794.749;
    KROK=Ur(KROKrowx,KROKcoly)*-1;
elseif mogix<=578057.3671 && mogiy>=7105794.749;
    KROK=Ur(KROKrowx,KROKcoly);
elseif mogix>=578057.3671 && mogiy>=7105794.749;
    KROK=Ur(KROKrowx,KROKcoly)*-1;
else
    KROK=Ur(KROKrowx,KROKcoly);
end;
if mogix>=595822.0785 && mogiy<=7099601.906;
    FROS=Ur(FROSrowx,FROScoly)*-1;
elseif mogix<=595822.0785 && mogiy>=7099601.906;
    FROS=Ur(FROSrowx,FROScoly);
elseif mogix>=595822.0785 && mogiy>=7099601.906;
    FROS=Ur(FROSrowx,FROScoly)*-1;
else
    FROS=Ur(FROSrowx,FROScoly);
end;
if mogix>=591324.9819 && mogiy<=7102811.594;
    DOMA=Ur(Drowx,Dcoly)*-1;
elseif mogix<=591324.9819 && mogiy>=7102811.594;
    DOMA=Ur(Drowx,Dcoly);
elseif mogix>=591324.9819 && mogiy>=7102811.594;
    DOMA=Ur(Drowx,Dcoly)*-1;
else
    DOMA=Ur(Drowx,Dcoly);
end;
if mogix>=594912.3594 && mogiy<=7097344.121;
    LAHR=Ur(LAHRrowx,LAHRcoly)*-1;
elseif mogix<=594912.3594 && mogiy>=7097344.121;
    LAHR=Ur(LAHRrowx,LAHRcoly);
elseif mogix>=594912.3594 && mogiy>=7097344.121;
    LAHR=Ur(LAHRrowx,LAHRcoly)*-1;
else
    LAHR=Ur(LAHRrowx,LAHRcoly);
end;
if mogix>=588598.9273 && mogiy<=7078205.599;
    THRA=Ur(TArowx,TAcoly)*-1;
elseif mogix<=588598.9273 && mogiy>=7078205.599;
    THRA=Ur(TArowx,TAcoly);
elseif mogix>=588598.9273 && mogiy>=7078205.599;
    THRA=Ur(TArowx,TAcoly)*-1;
else
    THRA=Ur(TArowx,TAcoly);
end;
if mogix>=595844.3353 && mogiy<=7082878.089;
    KKLO=Ur(KKrowx,KKcoly)*-1;
elseif mogix<=595844.3353 && mogiy>=7082878.089;
    KKLO=Ur(KKrowx,KKcoly);
elseif mogix>=595844.3353 && mogiy>=7082878.089;
    KKLO=Ur(KKrowx,KKcoly)*-1;
else
    KKLO=Ur(KKrowx,KKcoly);
end;
if mogix>=599775.4237 && mogiy<=7083000.709;
    KGIL=Ur(KGrowx,KGcoly)*-1;
elseif mogix<=599775.4237 && mogiy>=7083000.709;
    KGIL=Ur(KGrowx,KGcoly);
elseif mogix>=599775.4237 && mogiy>=7083000.709;
    KGIL=Ur(KGrowx,KGcoly)*-1;
else
    KGIL=Ur(KGrowx,KGcoly);
end;
if mogix>=588221.7342 && mogiy<=7091573.266;
    SATU=Ur(SArowx,SAcoly)*-1;
elseif mogix<=588221.7342 && mogiy>=7091573.266;
    SATU=Ur(SArowx,SAcoly);
elseif mogix>=588221.7342 && mogiy>=7091573.266;
    SATU=Ur(SArowx,SAcoly)*-1;
else
    SATU=Ur(SArowx,SAcoly);
end;
if mogix>=585129.5376 && mogiy<=7097061.677;
    LIHO=Ur(Lrowx,Lcoly)*-1;
elseif mogix<=585129.5376 && mogiy>=7097061.677;
    LIHO=Ur(Lrowx,Lcoly);
elseif mogix>=585129.5376 && mogiy>=7097061.677;
    LIHO=Ur(Lrowx,Lcoly)*-1;
else
    LIHO=Ur(Lrowx,Lcoly);
end;
if mogix>=580352.5266 && mogiy<=7092475.799;
    BLAU=Ur(Browx,Bcoly)*-1;
elseif mogix<=580352.5266 && mogiy>=7092475.799;
    BLAU=Ur(Browx,Bcoly);
elseif mogix>=580352.5266 && mogiy>=7092475.799;
    BLAU=Ur(Browx,Bcoly)*-1;
else
    BLAU=Ur(Browx,Bcoly);
end;
if mogix>=580005.0291 && mogiy<=7086893.04;
    LAUF=Ur(LFrowx,LFcoly)*-1;
elseif mogix<=580005.0291 && mogiy>=7086893.04;
    LAUF=Ur(LFrowx,LFcoly);
elseif mogix>=580005.0291 && mogiy>=7086893.04;
    LAUF=Ur(LFrowx,LFcoly)*-1;
else
    LAUF=Ur(LFrowx,LFcoly);
end;
if mogix>=587179.4827 && mogiy<=7093773.73;
    HRAF=Ur(HFrowx,HFcoly)*-1;
elseif mogix<=587179.4827 && mogiy>=7093773.73;
    HRAF=Ur(HFrowx,HFcoly);
elseif mogix>=587179.4827 && mogiy>=7093773.73;
    HRAF=Ur(HFrowx,HFcoly)*-1;
else
    HRAF=Ur(HFrowx,HFcoly);
end;
if mogix>=588221.7342 && mogiy<=7091573.266;
    HRAS=Ur(HSrowx,HScoly)*-1;
elseif mogix<=588221.7342 && mogiy>=7091573.266;
    HRAS=Ur(HSrowx,HScoly);
elseif mogix>=588221.7342 && mogiy>=7091573.266;
    HRAS=Ur(HSrowx,HScoly)*-1;
else
    HRAS=Ur(HSrowx,HScoly);
end;

% Calculating the x and y components of the horizontal displacements [in m]
ISAKx=ISAK*cos(ISAKangle*pi/180);
ISAKy=ISAK*sin(ISAKangle*pi/180);
KROKx=KROK*cos(KROKangle*pi/180);
KROKy=KROK*sin(KROKangle*pi/180);
FROSx=FROS*cos(FROSangle*pi/180);
FROSy=FROS*sin(FROSangle*pi/180);
LAHRx=LAHR*cos(LAHRangle*pi/180);
LAHRy=LAHR*sin(LAHRangle*pi/180);
SATUx=SATU*cos(SAangle*pi/180);
SATUy=SATU*sin(SAangle*pi/180);
THRAx=THRA*cos(TAangle*pi/180);
THRAy=THRA*sin(TAangle*pi/180);
KKLOx=KKLO*cos(KKangle*pi/180);
KKLOy=KKLO*sin(KKangle*pi/180);
KGILx=KGIL*cos(KGangle*pi/180);
KGILy=KGIL*sin(KGangle*pi/180);
LIHOx=LIHO*cos(Langle*pi/180);
LIHOy=LIHO*sin(Langle*pi/180);
BLAUx=BLAU*cos(Bangle*pi/180);
BLAUy=BLAU*sin(Bangle*pi/180);
LAUFx=LAUF*cos(LAangle*pi/180);
LAUFy=LAUF*sin(LAangle*pi/180);
HRAFx=HRAF*cos(HFangle*pi/180);
HRAFy=HRAF*sin(HFangle*pi/180);
HRASx=HRAS*cos(HSangle*pi/180);
HRASy=HRAS*sin(HSangle*pi/180);
DOMAx=DOMA*cos(Dangle*pi/180);
DOMAy=DOMA*sin(Dangle*pi/180);

% Vertical displacements for each site
ISAKz=Uz(ISAKrowx,ISAKcoly);
KROKz=Uz(KROKrowx,KROKcoly);
FROSz=Uz(FROSrowx,FROScoly);
LAHRz=Uz(LAHRrowx,LAHRcoly);
SATUz=Uz(SArowx,SAcoly);
THRAz=Uz(TArowx,TAcoly);
KKLOz=Uz(KKrowx,KKcoly);
KGILz=Uz(KGrowx,KGcoly);
LIHOz=Uz(Lrowx,Lcoly);
BLAUz=Uz(Browx,Bcoly);
LAUFz=Uz(LFrowx,LFcoly);
HRAFz=Uz(HFrowx,HFcoly);
HRASz=Uz(HSrowx,HScoly);
DOMAz=Uz(Drowx,Dcoly);

% Table of horizontal (x, y, and radial) and vertical displacements
KROK1=[KROK,KROKx,KROKy,KROKz];
ISAK1=[ISAK,ISAKx,ISAKy,ISAKz];
BLAU1=[BLAU,BLAUx,BLAUy,BLAUz];
LIHO1=[LIHO,LIHOx,LIHOy,LIHOz];
LAUF1=[LAUF,LAUFx,LAUFy,LAUFz];
HRAS1=[HRAS,HRASx,HRASy,HRASz];
HRAF1=[HRAF,HRAFx,HRAFy,HRAFz];
DOMA1=[DOMA,DOMAx,DOMAy,DOMAz];
FROS1=[FROS,FROSx,FROSy,FROSz];
LAHR1=[LAHR,LAHRx,LAHRy,LAHRz];
SATU1=[SATU,SATUx,SATUy,SATUz];
THRA1=[THRA,THRAx,THRAy,THRAz];
KKLO1=[KKLO,KKLOx,KKLOy,KKLOz];
KGIL1=[KGIL,KGILx,KGILy,KGILz];
TOTAL=[ISAK1;KROK1;BLAU1;LIHO1;LAUF1;HRAS1;HRAF1;DOMA1;FROS1;LAHR1;SATU1;THRA1;KKLO1;KGIL1];

% Calculating the rate of deformation due to the Mogi source for each site [in m/yr]
% (i.e., Displacement divided by "q" years, inputted at the start)
KROK_v=KROK/q;
KROKx_v=KROKx/q;
KROKy_v=KROKy/q;
ISAK_v=ISAK/q;
ISAKx_v=ISAKx/q;
ISAKy_v=ISAKy/q;
BLAU_v=BLAU/q;
BLAUx_v=BLAUx/q;
BLAUy_v=BLAUy/q;
LIHO_v=LIHO/q;
LIHOx_v=LIHOx/q;
LIHOy_v=LIHOy/q;
LAUF_v=LAUF/q;
LAUFx_v=LAUFx/q;
LAUFy_v=LAUFy/q;
HRAS_v=HRAS/q;
HRASx_v=HRASx/q;
HRASy_v=HRASy/q;
HRAF_v=HRAF/q;
HRAFx_v=HRAFx/q;
HRAFy_v=HRAFy/q;
DOMA_v=DOMA/q;
DOMAx_v=DOMAx/q;
DOMAy_v=DOMAy/q;
FROS_v=FROS/q;
FROSx_v=FROSx/q;
FROSy_v=FROSy/q;
LAHR_v=LAHR/q;
LAHRx_v=LAHRx/q;
LAHRy_v=LAHRy/q;
SATU_v=SATU/q;
SATUx_v=SATUx/q;
SATUy_v=SATUy/q;
THRA_v=THRA/q;
THRAx_v=THRAx/q;
THRAy_v=THRAy/q;
KKLO_v=KKLO/q;
KKLOx_v=KKLOx/q;
KKLOy_v=KKLOy/q;
KGIL_v=KGIL/q;
KGILx_v=KGILx/q;
KGILy_v=KGILy/q;

% Vertical velocities
KROKz_v=KROKz/q;
ISAKz_v=ISAKz/q;
BLAUz_v=BLAUz/q;
LIHOz_v=LIHOz/q;
LAUFz_v=LAUFz/q;
HRASz_v=HRASz/q;
HRAFz_v=HRAFz/q;
DOMAz_v=DOMAz/q;
FROSz_v=FROSz/q;
LAHRz_v=LAHRz/q;
SATUz_v=SATUz/q;
THRAz_v=THRAz/q;
KKLOz_v=KKLOz/q;
KGILz_v=KGILz/q;

% Table of horizontal (x, y, and radial) and vertical velocities
KROK2=[KROK_v,KROKx_v,KROKy_v,KROKz_v];
ISAK2=[ISAK_v,ISAKx_v,ISAKy_v,ISAKz_v];
BLAU2=[BLAU_v,BLAUx_v,BLAUy_v,BLAUz_v];
LIHO2=[LIHO_v,LIHOx_v,LIHOy_v,LIHOz_v];
LAUF2=[LAUF_v,LAUFx_v,LAUFy_v,LAUFz_v];
HRAS2=[HRAS_v,HRASx_v,HRASy_v,HRASz_v];
HRAF2=[HRAF_v,HRAFx_v,HRAFy_v,HRAFz_v];
DOMA2=[DOMA_v,DOMAx_v,DOMAy_v,DOMAz_v];
FROS2=[FROS_v,FROSx_v,FROSy_v,FROSz_v];
LAHR2=[LAHR_v,LAHRx_v,LAHRy_v,LAHRz_v];
SATU2=[SATU_v,SATUx_v,SATUy_v,SATUz_v];
THRA2=[THRA_v,THRAx_v,THRAy_v,THRAz_v];
KKLO2=[KKLO_v,KKLOx_v,KKLOy_v,KKLOz_v];
KGIL2=[KGIL_v,KGILx_v,KGILy_v,KGILz_v];
TOTALv=[ISAK2;KROK2;BLAU2;LIHO2;LAUF2;HRAS2;HRAF2;DOMA2;FROS2;LAHR2;SATU2;THRA2;KKLO2;KGIL2];

% Calculating the component of the horizontal velocity parallel to spreading
% (i.e. N78E) for each site [in m/yr]
angleK=(atan(KROKy_v/KROKx_v))*180/pi;
TangleK=360+(90-angleK);
DiffK=abs(282-TangleK);
comptK=KROK_v*cos(DiffK*pi/180);
comptxK=comptK*cos((282-270)*pi/180)*-1*1000;
comptyK=comptK*sin((282-270)*pi/180)*1000;

angleI=(atan(ISAKy_v/ISAKx_v))*180/pi;
TangleI=270-angleI;
DiffI=abs(282-TangleI);
comptI=ISAK_v*cos(DiffI*pi/180);
comptxI=comptI*cos((282-270)*pi/180)*-1*1000;
comptyI=comptI*sin((282-270)*pi/180)*1000;

angleB=(atan(BLAUy_v/BLAUx_v))*180/pi;
TangleB=90-angleB;
DiffB=abs(282-TangleB);
comptB=BLAU_v*cos(DiffB*pi/180);
comptxB=comptB*cos((282-270)*pi/180)*-1*1000;
comptyB=comptB*sin((282-270)*pi/180)*1000;

angleL=(atan(LIHOy_v/LIHOx_v))*180/pi;
TangleL=360+(90-angleL);
DiffL=abs(282-TangleL);
comptL=LIHO_v*cos(DiffL*pi/180);
comptxL=comptL*cos((282-270)*pi/180)*-1*1000;
comptyL=comptL*sin((282-270)*pi/180)*1000;

angleLA=(atan(LAUFy_v/LAUFx_v))*180/pi;
TangleLA=90-angleLA;
DiffLA=abs(282-TangleLA);
comptLA=LAUF_v*cos(DiffLA*pi/180);
comptxLA=comptLA*cos((282-270)*pi/180)*-1*1000;
comptyLA=comptLA*sin((282-270)*pi/180)*1000;

angleH=(atan(HRASy_v/HRASx_v))*180/pi;
TangleH=90-angleH;
DiffH=abs(282-TangleH);
comptH=HRAS_v*cos(DiffH*pi/180);
comptxH=comptH*cos((282-270)*pi/180)*-1*1000;
comptyH=comptH*sin((282-270)*pi/180)*1000;

angleHF=(atan(HRAFy_v/HRAFx_v))*180/pi;
TangleHF=90-angleHF;
DiffHF=abs(282-TangleHF);
comptHF=HRAF_v*cos(DiffHF*pi/180);
comptxHF=comptHF*cos((282-270)*pi/180)*-1*1000;
comptyHF=comptHF*sin((282-270)*pi/180)*1000;

angleD=(atan(DOMAy_v/DOMAx_v))*180/pi;
TangleD=90-angleD;
DiffD=abs(282-TangleD);
comptD=DOMA_v*cos(DiffD*pi/180);
comptxD=comptD*cos((282-270)*pi/180)*-1*1000;
comptyD=comptD*sin((282-270)*pi/180)*1000;

angleF=(atan(FROSy_v/FROSx_v))*180/pi;
TangleF=90-angleF;
DiffF=abs(282-TangleF);
comptF=FROS_v*cos(DiffF*pi/180);
comptxF=comptF*cos((282-270)*pi/180)*-1*1000;
comptyF=comptF*sin((282-270)*pi/180)*1000;

angleLH=(atan(LAHRy_v/LAHRx_v))*180/pi;
TangleLH=90-angleLH;
DiffLH=abs(282-TangleLH);
comptLH=LAHR_v*cos(DiffLH*pi/180);
comptxLH=comptLH*cos((282-270)*pi/180)*-1*1000;
comptyLH=comptLH*sin((282-270)*pi/180)*1000;

angleS=(atan(SATUy_v/SATUx_v))*180/pi;
TangleS=90-angleS;
DiffS=abs(282-TangleS);
comptS=SATU_v*cos(DiffS*pi/180);
comptxS=comptS*cos((282-270)*pi/180)*-1*1000;
comptyS=comptS*sin((282-270)*pi/180)*1000;

angleT=(atan(THRAy_v/THRAx_v))*180/pi;
TangleT=90-angleT;
DiffT=abs(282-TangleT);
comptT=THRA_v*cos(DiffT*pi/180);
comptxT=comptT*cos((282-270)*pi/180)*-1*1000;
comptyT=comptT*sin((282-270)*pi/180)*1000;

DiffKK=abs(282-(360+(90-KKangle)));
comptKK=KKLO_v*cos(DiffKK*pi/180);
comptxKK=comptKK*cos((282-270)*pi/180)*-1*1000;
comptyKK=comptKK*sin((282-270)*pi/180)*1000;

DiffKG=abs(282-(360+(90-KGangle)));
comptKG=KGIL_v*cos(DiffKG*pi/180);
comptxKG=comptKG*cos((282-270)*pi/180)*-1*1000;
comptyKG=comptKG*sin((282-270)*pi/180)*1000;

% Table of components of the horizontal velocities parallel to plate
% spreading

TOTAL0106=[comptxI,comptyI;comptxK,comptyK;comptxL,comptyL;comptxD,comptyD;comptxF,comptyF;comptxLH,comptyLH;comptxB,comptyB;comptxHF,comptyHF;comptxH,comptyH;comptxLA,comptyLA;comptxS,comptyS;comptxT,comptyT;comptxKK,comptyKK;comptxKG,comptyKG];

No comments:

Post a Comment