function [] = slipvects4(s,nve1,Lap,pm,pm1,pm2,pm3,pm4,sslip,dslip,nve,faults,G11,G21,G31,G41,G12,G22,G32,G42)
%
%	 slipvects(pm,slip,nve)
%
% plots slip on fault patches on 2-D grid
%
% 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
%
% Modification of plotpatchslip2 (P. Segall) by Kaj Johnson to plot slip vectors in the fault plane 


ns=sqrt(s(1:size(Lap,1)).^2+s(size(Lap,1)+1:end).^2);
%colormap(white)

depth=pm(:,3);


sslip1=sslip(1:size(pm1,1));
sslip2=sslip(size(pm1,1)+1:size(pm1,1)+size(pm2,1));
sslip3=sslip(size(pm1,1)+size(pm2,1)+1:size(pm1,1)+size(pm2,1)+size(pm3,1));
sslip4=sslip(size(pm1,1)+size(pm2,1)+size(pm3,1)+1:end);
dslip1=dslip(1:size(pm1,1));
dslip2=dslip(size(pm1,1)+1:size(pm1,1)+size(pm2,1));
dslip3=dslip(size(pm1,1)+size(pm2,1)+1:size(pm1,1)+size(pm2,1)+size(pm3,1));
dslip4=dslip(size(pm1,1)+size(pm2,1)+size(pm3,1)+1:end);

zslip1=zeros(size(sslip1));
zslip2=zeros(size(sslip2));
zslip3=zeros(size(sslip3));
zslip4=zeros(size(sslip4));


zslip=zeros(size(sslip));

dd1=zeros(size(depth));
dd2=dd1;dd3=dd1;dd4=dd1;

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));

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(pm1,1)/nve;
nhe2 = size(pm2,1)/nve;
nhe3 = size(pm3,1)/nve;
nhe4 = size(pm4,1)/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



%figure

%plots vectors at the center of the patch
xx1=(x1(:,1)+x1(:,2))/2;
xx2=(x2(:,1)+x2(:,2))/2;
xx3=(x3(:,1)+x3(:,2))/2;
xx4=(x4(:,1)+x4(:,2))/2;

yy1=(y1(:,1)+y1(:,2))/2;
yy2=(y2(:,1)+y2(:,2))/2;
yy3=(y3(:,1)+y3(:,2))/2;
yy4=(y4(:,1)+y4(:,2))/2;

hold on


% Shift the gridded planes so that the center of the top edge is at the origin
xcenter1=(xx1(1)+xx1(end))/2;
xcenter2=(xx2(1)+xx2(end))/2;
xcenter3=(xx3(1)+xx3(end))/2;
xcenter4=(xx4(1)+xx4(end))/2;


xx1=xx1-xcenter1;
xx2=xx2-xcenter2;
xx3=xx3-xcenter3;
xx4=xx4-xcenter4;


zz1=zeros(size(xx1));
zz2=zeros(size(xx2));
zz3=zeros(size(xx3));
zz4=zeros(size(xx4));

%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)];


   temp1=Ry1*Rx1*[xx1';yy1';zz1'];
   temp2=Ry2*Rx2*[xx2';yy2';zz2'];
   temp3=Ry3*Rx3*[xx3';yy3';zz3'];
   temp4=Ry4*Rx4*[xx4';yy4';zz4'];

   xx1=temp1(1,:);
   xx2=temp2(1,:);
   xx3=temp3(1,:);
   xx4=temp4(1,:);

   yy1=temp1(2,:);
   yy2=temp2(2,:);
   yy3=temp3(2,:);
   yy4=temp4(2,:);

   zz1=temp1(3,:);
   zz2=temp2(3,:);
   zz3=temp3(3,:);
   zz4=temp4(3,:);
   
   temp1=Ry1*Rx1*[sslip1';dslip1';zslip1'];
   temp2=Ry2*Rx2*[sslip2';dslip2';zslip2'];
   temp3=Ry3*Rx3*[sslip3';dslip3';zslip3'];
   temp4=Ry4*Rx4*[sslip4';dslip4';zslip4'];

   sslip1=temp1(1,:);
   sslip2=temp2(1,:);
   sslip3=temp3(1,:);
   sslip4=temp4(1,:);

   dslip1=temp1(2,:);
   dslip2=temp2(2,:);
   dslip3=temp3(2,:);
   dslip4=temp4(2,:);

   zslip1=temp1(3,:);
   zslip2=temp2(3,:);
   zslip3=temp3(3,:);
   zslip4=temp4(3,:);
   


% translate faults to actual position
xx1=xx1+faults(1,7);
xx2=xx2+faults(2,7);
xx3=xx3+faults(3,7);
xx4=xx4+faults(4,7);

zz1=zz1+faults(1,6);
zz2=zz2+faults(2,6);
zz3=zz3+faults(3,6);
zz4=zz4+faults(4,6);



%drop fault 4 vectors to actual position
yy4=yy4-faults(4,3);

S=.75; %arbitrary scaling of the slip vectors
quiver3(xx1,-zz1,yy1+.2,sslip1*S,-zslip1*S,dslip1*S,0,'k'), axis('equal')
hold on
quiver3(xx2,-zz2,yy2+.2,sslip2*S,-zslip2*S,dslip2*S,0,'k'), axis('equal')
quiver3(xx3,-zz3,yy3+.2,sslip3*S,-zslip3*S,dslip3*S,0,'k'), axis('equal')
quiver3(xx4,-zz4,yy4+.2,sslip4*S,-zslip4*S,dslip4*S,0,'k'), axis('equal')


%hold on


axis equal
grid on
view([-1 -1 1])
