% Generic.m calculates surface displacements for a specified model
%global medium_type source_type idum ux uy uz x y uxtot uytot uztot

if idum == 1
%%%%% Create a 2 dimensional rectangular grid 

xmax    = get(xmaxbut, 'string');	       
xmax    = str2num(xmax);
xmin    = get(xminbut, 'string');	       
xmin    = str2num(xmin);
xspace  = get(xspacebut, 'string');	       
xspace  = str2num(xspace);

ymax    = get(ymaxbut, 'string');	       
ymax    = str2num(ymax);
ymin    = get(yminbut, 'string');	       
ymin    = str2num(ymin);
yspace  = get(yspacebut, 'string');	       
yspace  = str2num(yspace);

xvec = 1e3*[xmin:xspace:xmax];
yvec = 1e3*[ymin:yspace:ymax];

[x0,y0]=meshgrid(xvec,yvec);
XB=[];
YB=[];

%	Medium properties	 
nu    = get(nubut, 'string');	%Poisson's ratio       
nu    = str2num(nu);
mu    = get(mubut, 'string');	%Shear modulus       
mu    = str2num(mu);    
delete(medium_plate);
end

shiftx  = get(shiftxbut, 'string');  %source offset from the origin, x
shiftx  = 1e3*str2num(shiftx);
shifty  = get(shiftybut, 'string');  %source offset from the origin, y
shifty  = 1e3*str2num(shifty);

[x,y]=meshgrid(xvec-shiftx,yvec-shifty);

%initialize data
[rowmax,colmax] = size(x);
ux=zeros(rowmax,colmax);
uy=zeros(rowmax,colmax);
uz=zeros(rowmax,colmax);
if source_type(idum) == 4     %Okada disl.
 ux1=zeros(rowmax,colmax);
 uy1=zeros(rowmax,colmax);
 uz1=zeros(rowmax,colmax);
 ux2=zeros(rowmax,colmax);
 uy2=zeros(rowmax,colmax);
 uz2=zeros(rowmax,colmax);
 ux3=zeros(rowmax,colmax);
 uy3=zeros(rowmax,colmax);
 uz3=zeros(rowmax,colmax);
 tp=zeros(rowmax,colmax);
end
if idum == 1
% load dump;
 uxtot=zeros(rowmax,colmax);
 uytot=zeros(rowmax,colmax);
 uztot=zeros(rowmax,colmax);
end

switch source_type(idum)
case 1	%mogi
p      = get(pressurebut, 'string');  %excess magma pressure
p      = str2num(p);
c      = get(lengthbut, 'string');  %source radius
c      = str2num(c);
d      = get(depthbut, 'string');	%depth, meters 
d      = str2num(d);

save SIC_par p c d

c      = c^3*p*(1-nu)/mu*1e-3;  %leading coefficient
[ux,uy,uz] = calc_mogi(x,y,d,c);

case 2	%penny
d      = get(depthbut, 'string');	%depth, meters 
d      = str2num(d);
h      = get(hbut, 'string');	%ratio d/R (R is crack radius)
h      = str2num(h);
%h      = har(h);
p      = get(pressurebut, 'string');  %excess magma pressure

save SIC_par p d h

p      = str2num(p)*2*(1-nu)/mu*1e-3;
[ux,uy,uz] = calc_penny(x,y,d,h,p);

case 3	%yang&davis spheroid
here1=0
d      = get(depthbut, 'string');	%depth, meters 
d      = str2num(d)
a      = get(abut, 'string');          %major semiaxis
a      = str2num(a)
b      = get(bbut, 'string');          %minor semiaxis
b      = str2num(b)
strike = get(strikebut, 'string');    %strike of a long axis, degrees
strike =  str2num(strike)
theta  = get(thetabut, 'string');     %plunge of a long axis, degrees
theta  =  str2num(theta)
p      = get(pressurebut, 'string');  %excess magma pressure

save SIC_par a b p d theta strike

p      = str2num(p)/mu*1e-3
here=0
[ux,uy,uz] = calc_yang(theta,strike,nu,p,a,b,d,x,y);

case 4	%okada
Un = get(unbut, 'string');		%fault opening, meters
Un = str2num(Un);
Ud = get(udbut, 'string');		%dip slip, meters
Ud = str2num(Ud);
Us = get(usbut, 'string');		%strike slip, meters
Us = str2num(Us);
delta  = get(dipbut, 'string');		%fault dip, degrees
delta  = str2num(delta)*pi/180;
strike = get(strikebut, 'string');	%fault strike, degrees
strike =  str2num(strike)*pi/180;
len    = get(lengthbut, 'string');	%fault length, meters
len    = str2num(len);
wid      = get(widthbut, 'string');     	%fault width, meters
wid      = str2num(wid);
d      = get(depthbut, 'string');%fault depth, meters (top of fault patch)
d      = str2num(d);

save SIC_par Un Ud Us delta strike len wid d 

%fprintf('Fault top is at %d ',(d-wid*sin(delta))/1e3);
%fprintf(' x0= %d ',(shiftx-wid*cos(delta)*cos(strike))/1e3);
%fprintf(' y0= %d\n',(shifty+wid*cos(delta)*sin(strike))/1e3);
fprintf('Fault bottom is at %d ',(d+wid*sin(delta))/1e3);
fprintf(' x0= %d ',shiftx/1e3);
fprintf(' y0= %d\n',shifty/1e3);
sins=sin(strike);
coss=cos(strike);
Lx=0.5*len*sins;
Ly=0.5*len*coss;
W=wid*cos(delta);
xb=[-Lx+W*coss -Lx Lx Lx+W*coss -Lx+W*coss]+shiftx;
yb=[-Ly-W*sins -Ly Ly Ly-W*sins -Ly-W*sins]+shifty; 
XB=[XB xb];
YB=[YB yb];

% fault_type 1=strike,2=dip,3=tensile
% call okada routines

if Us ~= 0
 [ux1,uy1,uz1] = ...
  calc_okada(1,Us,x,y,nu,delta,d,len,wid,1,strike,tp);
end  
if Ud ~= 0
 [ux2,uy2,uz2] = ...
  calc_okada(1,Ud,x,y,nu,delta,d,len,wid,2,strike,tp);
end  
if Un ~= 0
 [ux3,uy3,uz3] = ...
  calc_okada(1,Un,x,y,nu,delta,d,len,wid,3,strike,tp);
end  
 ux = ux1 + ux2 + ux3;
 uy = uy1 + uy2 + uy3;
 uz = uz1 + uz2 + uz3;

end %switch source_type

delete(plate);
fprintf('press any key to continue\n\n');
