%%%%%%%%%%%%%
%
% Kaj Johnson, Sept, 2000
%
% I altered plotpatchslip2(pm,slip,nve) to plot the patched fault surfaces in
% 3-D. This was designed to plot 4 faults surfaces. It must be slightly altered
% to plot more or less than 4 faults.
%
%
% Input:
%	pm = patchmodel (dislocation parameters)
%	pm(:,1) = length;
%	pm(:,2) = width;
%	pm(:,3) = depth;
%	pm(:,4) = dip;
%	pm(:,5) = strike;
%	pm(:,6) = East Offset;
%	pm(:,7) = North Offset;
%	slip = slip values
%	nve  = number of vertical elements
%

ns=s(1:end/2).^2+s(end/2+1
% grid three fault planes and then rotate into true (3D) position
nve=nve1;
slip1=ns(1:size(G11,2));
slip2=ns(1+size(G11,2):size(G11,2)+size(G21,2));
slip3=ns(1+size(G11,2)+size(G21,2):size(G11,2)+size(G21,2)+size(G31,2));
slip4=ns(1+size(G11,2)+size(G21,2)+size(G31,2):end);

depth1=pm1(:,3);
width1=pm1(:,2);
dip1=pm1(:,4);
dd1=zeros(size(depth1));

depth2=pm2(:,3);
width2=pm2(:,2);
dip2=pm2(:,4);
dd2=zeros(size(depth2));

depth3=pm3(:,3);
width3=pm3(:,2);
dip3=pm3(:,4);
dd3=zeros(size(depth3));

depth4=pm4(:,3);
width4=pm4(:,2);
dip4=pm4(:,4);
dd4=zeros(size(depth4));

% this needs to be changed if the depth of the fault planes are not all zero!!!
dd1(1:nve:end)=0*depth1(1:nve:end);
dd2(1:nve:end)=0*depth2(1:nve:end);
dd3(1:nve:end)=0*depth3(1:nve:end);
dd4(1:nve:end)=0*depth4(1:nve:end);

for k=2:nve
   dd1(k:nve:end)=dd1(k-1:nve:end) + pm1(1,2);
   dd2(k:nve:end)=dd2(k-1:nve:end) + pm2(1,2);
   dd3(k:nve:end)=dd3(k-1:nve:end) + pm3(1,2);
   dd4(k:nve:end)=dd4(k-1:nve:end) + pm4(1,2);
end


y1 = -[dd1 + pm1(:,2), dd1];
y2 = -[dd2 + pm2(:,2), dd2];
y3 = -[dd3 + pm3(:,2), dd3];
y4 = -[dd4 + pm4(:,2), dd4];


x1 = zeros(size(y1));
x2 = zeros(size(y2));
x3 = zeros(size(y3));
x4 = zeros(size(y4));

nhe1 = size(slip1)/nve;
nhe2 = size(slip2)/nve;
nhe3 = size(slip3)/nve;
nhe4 = size(slip4)/nve;


% for first nve elements
x1(1:nve,2) = pm1(1,1);
x2(1:nve,2) = pm2(1,1);
x3(1:nve,2) = pm3(1,1);
x4(1:nve,2) = pm4(1,1);

for k=2:nhe1
	index1 = nve*(k-1);
	index2 = nve*(k-2);
	x1( index1+1:index1+nve , :) = ...
	[x1(index2+1:index2+nve,2), ...
		 x1(index2+1:index2+nve,2) + pm1(index1+1,1)*ones(nve,1)];
end

for k=2:nhe2
	index1 = nve*(k-1);
	index2 = nve*(k-2);
	x2( index1+1:index1+nve , :) = ...
	[x2(index2+1:index2+nve,2), ...
		 x2(index2+1:index2+nve,2) + pm2(index1+1,1)*ones(nve,1)];
end

for k=2:nhe3
	index1 = nve*(k-1);
	index2 = nve*(k-2);
	x3( index1+1:index1+nve , :) = ...
	[x3(index2+1:index2+nve,2), ...
		 x3(index2+1:index2+nve,2) + pm3(index1+1,1)*ones(nve,1)];
end

for k=2:nhe4
	index1 = nve*(k-1);
	index2 = nve*(k-2);
	x4( index1+1:index1+nve , :) = ...
	[x4(index2+1:index2+nve,2), ...
		 x4(index2+1:index2+nve,2) + pm4(index1+1,1)*ones(nve,1)];
end



% Shift the gridded planes so that the center of the top edge is at the origin
xcenter1=(x1(1,1)+x1(end,2))/2;
xcenter2=(x2(1,1)+x2(end,2))/2;
xcenter3=(x3(1,1)+x3(end,2))/2;
xcenter4=(x4(1,1)+x4(end,2))/2;

X1 = [x1(:,1)'; x1(:,2)'; x1(:,2)'; x1(:,1)']-xcenter1;
Y1 = [y1(:,1)'; y1(:,1)'; y1(:,2)'; y1(:,2)'];
Z1=zeros(size(Y1));

X2 = [x2(:,1)'; x2(:,2)'; x2(:,2)'; x2(:,1)']-xcenter2;
Y2 = [y2(:,1)'; y2(:,1)'; y2(:,2)'; y2(:,2)'];
Z2=zeros(size(Y2));

X3 = [x3(:,1)'; x3(:,2)'; x3(:,2)'; x3(:,1)']-xcenter3;
Y3 = [y3(:,1)'; y3(:,1)'; y3(:,2)'; y3(:,2)'];
Z3=zeros(size(Y3));

X4 = [x4(:,1)'; x4(:,2)'; x4(:,2)'; x4(:,1)']-xcenter4;
Y4 = [y4(:,1)'; y4(:,1)'; y4(:,2)'; y4(:,2)'];
Z4=zeros(size(Y4));


%rotate and translate gridded faults to actual position (3D)

%rotate about x-axis (90-dip) degrees

theta1=(90-pm1(1,4))*pi/180;
theta2=(90-pm2(1,4))*pi/180;
theta3=(90-pm3(1,4))*pi/180;
theta4=(90-pm4(1,4))*pi/180;
strike1=-faults(1,5)*pi/180;
strike2=-faults(2,5)*pi/180;
strike3=-faults(3,5)*pi/180;
strike4=-faults(4,5)*pi/180;

Rx1=[1 0 0;0 cos(theta1) sin(theta1);0 -sin(theta1) cos(theta1)];
Rx2=[1 0 0;0 cos(theta2) sin(theta2);0 -sin(theta2) cos(theta2)];
Rx3=[1 0 0;0 cos(theta3) sin(theta3);0 -sin(theta3) cos(theta3)];
Rx4=[1 0 0;0 cos(theta4) sin(theta4);0 -sin(theta4) cos(theta4)];

Ry1=[cos(strike1) 0 sin(strike1);0 1 0;-sin(strike1) 0 cos(strike1)];
Ry2=[cos(strike2) 0 sin(strike2);0 1 0;-sin(strike2) 0 cos(strike2)];
Ry3=[cos(strike3) 0 sin(strike3);0 1 0;-sin(strike3) 0 cos(strike3)];
Ry4=[cos(strike4) 0 sin(strike4);0 1 0;-sin(strike4) 0 cos(strike4)];


for m=1:4
   temp1=Ry1*Rx1*[X1(m,:);Y1(m,:);Z1(m,:)];
   temp2=Ry2*Rx2*[X2(m,:);Y2(m,:);Z2(m,:)];
   temp3=Ry3*Rx3*[X3(m,:);Y3(m,:);Z3(m,:)];
   temp4=Ry4*Rx4*[X4(m,:);Y4(m,:);Z4(m,:)];

   X1(m,:)=temp1(1,:);
   X2(m,:)=temp2(1,:);
   X3(m,:)=temp3(1,:);
   X4(m,:)=temp4(1,:);
   Y1(m,:)=temp1(2,:);
   Y2(m,:)=temp2(2,:);
   Y3(m,:)=temp3(2,:);
   Y4(m,:)=temp4(2,:);
   Z1(m,:)=temp1(3,:);
   Z2(m,:)=temp2(3,:);
   Z3(m,:)=temp3(3,:);
   Z4(m,:)=temp4(3,:);   
end


% translate faults to actual position
X1=X1+faults(1,7);
X2=X2+faults(2,7);
X3=X3+faults(3,7);
X4=X4+faults(4,7);
Z1=Z1+faults(1,6);
Z2=Z2+faults(2,6);
Z3=Z3+faults(3,6);
Z4=Z4+faults(4,6);
Y4=Y4-faults(4,3);

figure
patch(X1,-Z1,Y1,slip1'), axis('equal')
hold on
patch(X2,-Z2,Y2,slip2')
patch(X3,-Z3,Y3,slip3')
patch(X4,-Z4,Y4,slip4')

colormap(flipud(hot))
colorbar
grid on
view([-1 -1 1])