%
% comparison of random distributions -- for rotations of double-couples (DC)
% Tests equal-area projections
%
% Output ROTDCE1 program, 10000 points, 
%      or ROTDCG_TEST3.FOR -- to see influence of 
%      angle discretization in T- and P-axes
%
%  Aug 24, 2004
%
clear all
%
fs1 = 14
fs2 = 14
%
rad = 180/pi
%
rr = 2*sin(0.5*acos(sqrt(2/3)))
x1 = [0 0]
y1 = [0 -rr]
x2 = [0 sqrt(3)*rr/2]
y2 = [0 rr/2]
x3 = [0 -sqrt(3)*rr/2]
y3 = [0 rr/2]
%
nx = 500;
nx = 100;
%
z3l1 = sin(0/rad)
z2l1 = linspace(0,sqrt(1-z3l1^2),nx);
for i=1:nx;
z1l1(i) = sqrt(1-z3l1^2-z2l1(i)^2);
rrl1(i) = 2*sin(0.5*acos((z3l1+z2l1(i)+z1l1(i))*sqrt(1/3)));
nl1(i) = sqrt(2*((z3l1-z2l1(i))^2+(z3l1-z1l1(i))^2+(z1l1(i)-z2l1(i))^2));
xl1(i) = sqrt(3)*(rrl1(i)/nl1(i))*(z1l1(i)-z2l1(i));
yl1(i) = (rrl1(i)/nl1(i))*(2*z3l1-z2l1(i)-z1l1(i));
end
%
z3l2 = sin(30/rad)
z2l2 = linspace(0,sqrt(1-z3l2^2),nx);
for i=1:nx;
z1l2(i) = sqrt(1-z3l2^2-z2l2(i)^2);
rrl2(i) = 2*sin(0.5*acos((z3l2+z2l2(i)+z1l2(i))*sqrt(1/3)));
nl2(i) = sqrt(2*((z3l2-z2l2(i))^2+(z3l2-z1l2(i))^2+(z1l2(i)-z2l2(i))^2));
xl2(i) = sqrt(3)*(rrl2(i)/nl2(i))*(z1l2(i)-z2l2(i));
yl2(i) = (rrl2(i)/nl2(i))*(2*z3l2-z2l2(i)-z1l2(i));
end
%
z3l3 = sin(60/rad)
z2l3 = linspace(0,sqrt(1-z3l3^2),nx);
for i=1:nx;
z1l3(i) = sqrt(1-z3l3^2-z2l3(i)^2);
rrl3(i) = 2*sin(0.5*acos((z3l3+z2l3(i)+z1l3(i))*sqrt(1/3)));
nl3(i) = sqrt(2*((z3l3-z2l3(i))^2+(z3l3-z1l3(i))^2+(z1l3(i)-z2l3(i))^2));
xl3(i) = sqrt(3)*(rrl3(i)/nl3(i))*(z1l3(i)-z2l3(i));
yl3(i) = (rrl3(i)/nl3(i))*(2*z3l3-z2l3(i)-z1l3(i));
end
%
%
z2n1 = sin(0/rad)
z3n1 = linspace(0,sqrt(1-z2n1^2),nx);
for i=1:nx;
z1n1(i) = sqrt(1-z3n1(i)^2-z2n1^2);
rrn1(i) = 2*sin(0.5*acos((z3n1(i)+z2n1+z1n1(i))*sqrt(1/3)));
nn1(i) = sqrt(2*((z3n1(i)-z2n1)^2+(z3n1(i)-z1n1(i))^2+(z1n1(i)-z2n1)^2));
xn1(i) = sqrt(3)*(rrn1(i)/nn1(i))*(z1n1(i)-z2n1);
yn1(i) = (rrn1(i)/nn1(i))*(2*z3n1(i)-z2n1-z1n1(i));
end
%
z2n2 = sin(30/rad)
z3n2 = linspace(0,sqrt(1-z2n2^2),nx);
for i=1:nx;
z1n2(i) = sqrt(1-z3n2(i)^2-z2n2^2);
rrn2(i) = 2*sin(0.5*acos((z3n2(i)+z2n2+z1n2(i))*sqrt(1/3)));
nn2(i) = sqrt(2*((z3n2(i)-z2n2)^2+(z3n2(i)-z1n2(i))^2+(z1n2(i)-z2n2)^2));
xn2(i) = sqrt(3)*(rrn2(i)/nn2(i))*(z1n2(i)-z2n2);
yn2(i) = (rrn2(i)/nn2(i))*(2*z3n2(i)-z2n2-z1n2(i));
end
%
z2n3 = sin(60/rad)
z3n3 = linspace(0,sqrt(1-z2n3^2),nx);
for i=1:nx;
z1n3(i) = sqrt(1-z3n3(i)^2-z2n3^2);
rrn3(i) = 2*sin(0.5*acos((z3n3(i)+z2n3+z1n3(i))*sqrt(1/3)));
nn3(i) = sqrt(2*((z3n3(i)-z2n3)^2+(z3n3(i)-z1n3(i))^2+(z1n3(i)-z2n3)^2));
xn3(i) = sqrt(3)*(rrn3(i)/nn3(i))*(z1n3(i)-z2n3);
yn3(i) = (rrn3(i)/nn3(i))*(2*z3n3(i)-z2n3-z1n3(i));
end
%
%
z1t1 = sin(0/rad)
z3t1 = linspace(0,sqrt(1-z1t1^2),nx);
for i=1:nx;
z2t1(i) = sqrt(1-z3t1(i)^2-z1t1^2);
rrt1(i) = 2*sin(0.5*acos((z3t1(i)+z2t1(i)+z1t1)*sqrt(1/3)));
nt1(i) = sqrt(2*((z3t1(i)-z2t1(i))^2+(z3t1(i)-z1t1)^2+(z1t1-z2t1(i))^2));
xt1(i) = sqrt(3)*(rrt1(i)/nt1(i))*(z1t1-z2t1(i));
yt1(i) = (rrt1(i)/nt1(i))*(2*z3t1(i)-z2t1(i)-z1t1);
end
%
z1t2 = sin(30/rad)
z3t2 = linspace(0,sqrt(1-z1t2^2),nx);
for i=1:nx;
z2t2(i) = sqrt(1-z3t2(i)^2-z1t2^2);
rrt2(i) = 2*sin(0.5*acos((z3t2(i)+z2t2(i)+z1t2)*sqrt(1/3)));
nt2(i) = sqrt(2*((z3t2(i)-z2t2(i))^2+(z3t2(i)-z1t2)^2+(z1t2-z2t2(i))^2));
xt2(i) = sqrt(3)*(rrt2(i)/nt2(i))*(z1t2-z2t2(i));
yt2(i) = (rrt2(i)/nt2(i))*(2*z3t2(i)-z2t2(i)-z1t2);
end
%
z1t3 = sin(60/rad)
z3t3 = linspace(0,sqrt(1-z1t3^2),nx);
for i=1:nx;
z2t3(i) = sqrt(1-z3t3(i)^2-z1t3^2);
rrt3(i) = 2*sin(0.5*acos((z3t3(i)+z2t3(i)+z1t3)*sqrt(1/3)));
nt3(i) = sqrt(2*((z3t3(i)-z2t3(i))^2+(z3t3(i)-z1t3)^2+(z1t3-z2t3(i))^2));
xt3(i) = sqrt(3)*(rrt3(i)/nt3(i))*(z1t3-z2t3(i));
yt3(i) = (rrt3(i)/nt3(i))*(2*z3t3(i)-z2t3(i)-z1t3);
end
%
plot(x1,y1,'k--', x2,y2,'k--', x3,y3,'k--', 'LineWidth',2.0)
hold on
plot(xl1,yl1,'k-', xl2,yl2,'k-', xl3,yl3,'k-', ...
xn1,yn1,'k-', xn2,yn2,'k-', xn3,yn3,'k-', ...
xt1,yt1,'k-', xt2,yt2,'k-', xt3,yt3,'k-', 'LineWidth', 0.75)
%
% title('Random rotation (eigenvectors) of DC source, ROTDCE2, Octovue projection', 'FontSize', fs1)
title('Rotation of DC source, Octovue projection', 'FontSize', fs1)
xlabel('X', 'FontSize', fs1)
ylabel('Y', 'FontSize', fs1)
%
grid
axis square
%
load for018_2500.dat
for018 = for018_2500;
%
z = for018(:,6:6);
x = for018(:,7:7);
y = for018(:,8:8);
%
%
[narrf, narrf1] = size(z)
hmin = min(x)
hmax = max(x)
vmax = max(y)
vmin = min(y)
steph = (hmax - hmin)/(narrf - 1)
stepv = (vmax - vmin)/(narrf - 1)
%
[xi, yi] = meshgrid(hmin:steph:hmax,vmin:stepv:vmax);      %(b)
zi = griddata(x,y,z,xi,yi,'cubic');
%
% [c,h] = contour(xi,yi,zi,v); 
[c,h] = contour(xi,yi,zi,'-r');
clabel(c,h)
hold off