三维方向统计：TomoFab

---

1 方向的定义与选取

参考文献 Tauxe, L., 2003. Paleomagnetic Principles and Practice, Modern Approaches in Geophysics. Kluwer Academic Publishers, Dordrecht.

1.1 坐标系

地理坐标系如下图

1.2 倾向(Dip direction)和倾角(Dip angle)

XY平面（North-East平面）为水平面。

1. 走向：断层面与水平面的交线称为走向线，走向线指示的地理方位（与地理北极沿顺时针方向的夹角）叫走向。走向线有无数条平行线，但走向只有两个，且相差180°。
2. 倾向：断层面上与走向线相垂直的线叫倾斜线或真倾斜线,它在水平面上的投影所指的沿平面向下倾斜的方位即为断层的倾向。倾向与走向相差90°或270°，但倾向确定后，走向就可以确定，走向确定后，倾向不一定确定。
3. 倾角：岩层层面上的真倾斜线与其在水平面上投影线的夹角叫倾角。它表示在垂直断层面走向的直立剖面上该层面与水平面间的夹角。

由定义可知，地球上的经线所在的大圆平面（与水平面的交线就是北极轴）的走向为0倾向为90度。纬线的倾向为0度。

1.3 Lambert投影(Schmidt投影)

为了保持不同区域的面积等比例，要求圆球的半径为根号2，因而投影后圆的半径为2。并且是z轴大于0的半球叫做下半球(lower hemisphere)。

2 Fisher置信区间

3D rock fabric analysis using micro-tomography: An introduction to the open-source TomoFab MATLAB code

2.1数据预处理 ImportFisherVector3d.m

输入为n个坐标系，v1v2v3对应每个坐标系的三个坐标轴的笛卡尔方向，n*3
输出为对接tomofab的格式数据
orentation tensor linear fibric的定义参照文献中的公式(3)和(4)
对接tomofab的结果如下

function [V1,V2,V3,Vtot]=ImportFisherVector3d(v1,v2,v3)
%输入为n个坐标系，v1v2v3对应每个坐标系的三个坐标轴的笛卡尔方向，n*3
%输出为对接tomofab的格式数据
%数据预先清理
%将向量转化为单位长度的向量
V1.cosine=v1./ sqrt(sum(v1.^2, 2));
V2.cosine=v2./ sqrt(sum(v2.^2, 2));
V3.cosine=v3./ sqrt(sum(v3.^2, 2));
% replace all directional cosines to the lower hemisphere
for i=1:size(V1.cosine,1) 
    if V1.cosine(i,3)>0
        V1.cosine(i,1)=-V1.cosine(i,1);
        V1.cosine(i,2)=-V1.cosine(i,2);
        V1.cosine(i,3)=-V1.cosine(i,3);
    end
    if V2.cosine(i,3)>0
        V2.cosine(i,1)=-V2.cosine(i,1);
        V2.cosine(i,2)=-V2.cosine(i,2);
        V2.cosine(i,3)=-V2.cosine(i,3);
    end
    if V3.cosine(i,3)>0
        V3.cosine(i,1)=-V3.cosine(i,1);
        V3.cosine(i,2)=-V3.cosine(i,2);
        V3.cosine(i,3)=-V3.cosine(i,3);
    end
end
[V1,V2,V3,Vtot]=OrientationStatistics(V1,V2,V3);
end

function [V1,V2,V3,Vtot]=OrientationStatistics(V1,V2,V3)
%
% function OrientationStatistics: Calculates orientation statistics in various ways. See Tauxe (2002) for details.
% input:   - V1, V2, V3: structures containing the selected set of 
%           directional data.
%           - params: structure containing calculated parameters and
%           methods of calculation
%          - UIcontrol recording projection type
% output: - V1, V2, V3: structures containing the selected set of 
%           directional data.
%           - Vtot: structure containing total orientation tensor and
%           associated parameters
%           - params: structure containing calculated parameters and
%           methods of calculation
% references:
% - Fisher, R., 1953. Dispersion on a sphere. Proc. R. Soc. London. Ser.
%   A. Math. Phys. Sci. 217, 295 LP-305.
% - Hext, G.R., 1963. The Estimation of Second-Order Tensors, with 
%   Related Tests and Designs. Biometrika 50, 353-373.
%   https://doi.org/10.2307/2333905
% - Tauxe, L., 2002. Paleomagnetic principles and practice.
%   Kluwer, Dodrecht.
% - Petri, B., Almqvist, B.S.G., Pistone, M., 2020. 3D rock fabric analysis
%   using micro-tomography: An introduction to the open-source TomoFab
%   MATLAB code. Comput. Geosci. 138, 104444.
%   https://doi.org/10.1016/j.cageo.2020.104444
% - Scheidegger, A.E., 1965. On the statistics of the orientation of 
%   bedding planes, grain axes, and similar sedimentological data. US
%   Geol. Surv. Prof. Pap. 525, 164-167.
% - Jelinek, V., 1978. Statistical processing of anisotropy of magnetic
%   susceptibility measured on groups of specimens. Stud. Geophys. 
%   Geod. 22, 50-62. https://doi.org/10.1007/BF01613632
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     This file is part of TOMOFAB. Copyright (C) 2018-2021  Benoit Petri
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.
%      
%     Please report any bug, error or suggestion to bpetri@unistra.fr
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%% Calculation of orientation tensors
% Compiling orientation tensor for V1 directions
V1.ortens=zeros(3,3);

for i=1:size(V1.cosine,1)
    V1.ortens(1,1) = V1.ortens(1,1)+V1.cosine(i,1)^2;
    V1.ortens(2,2) = V1.ortens(2,2)+V1.cosine(i,2)^2;
    V1.ortens(3,3) = V1.ortens(3,3)+V1.cosine(i,3)^2;
    V1.ortens(1,2) = V1.ortens(1,2)+V1.cosine(i,1)*V1.cosine(i,2);
    V1.ortens(1,3) = V1.ortens(1,3)+V1.cosine(i,1)*V1.cosine(i,3);
    V1.ortens(2,3) = V1.ortens(2,3)+V1.cosine(i,2)*V1.cosine(i,3);
end
V1.ortens(2,1) = V1.ortens(1,2); % The orientation tensor is symmetric, hence can be constructed after the loop
V1.ortens(3,1) = V1.ortens(1,3);
V1.ortens(3,2) = V1.ortens(2,3);
V1.ortens = V1.ortens/(size(V1.cosine,1));
V1.oteigval=abs(sort(eig(V1.ortens)));

% Compiling orientation tensor for V2 directions
V2.ortens=zeros(3,3);
for i=1:size(V2.cosine,1)
    V2.ortens(1,1) = V2.ortens(1,1)+V2.cosine(i,1)^2;
    V2.ortens(2,2) = V2.ortens(2,2)+V2.cosine(i,2)^2;
    V2.ortens(3,3) = V2.ortens(3,3)+V2.cosine(i,3)^2;
    V2.ortens(1,2) = V2.ortens(1,2)+V2.cosine(i,1)*V2.cosine(i,2);
    V2.ortens(1,3) = V2.ortens(1,3)+V2.cosine(i,1)*V2.cosine(i,3);
    V2.ortens(2,3) = V2.ortens(2,3)+V2.cosine(i,2)*V2.cosine(i,3);
end
V2.ortens(2,1) = V2.ortens(1,2);
V2.ortens(3,1) = V2.ortens(1,3);
V2.ortens(3,2) = V2.ortens(2,3);
V2.ortens = V2.ortens/(size(V2.cosine,1));
V2.oteigval=abs(sort(eig(V2.ortens)));

% Compiling orientation tensor for V3 directions
V3.ortens=zeros(3,3);
for i=1:size(V3.cosine,1)
    V3.ortens(1,1) = V3.ortens(1,1)+V3.cosine(i,1)^2;
    V3.ortens(2,2) = V3.ortens(2,2)+V3.cosine(i,2)^2;
    V3.ortens(3,3) = V3.ortens(3,3)+V3.cosine(i,3)^2;
    V3.ortens(1,2) = V3.ortens(1,2)+V3.cosine(i,1)*V3.cosine(i,2);
    V3.ortens(1,3) = V3.ortens(1,3)+V3.cosine(i,1)*V3.cosine(i,3);
    V3.ortens(2,3) = V3.ortens(2,3)+V3.cosine(i,2)*V3.cosine(i,3);
end
V3.ortens(2,1) = V3.ortens(1,2);
V3.ortens(3,1) = V3.ortens(1,3);
V3.ortens(3,2) = V3.ortens(2,3);
V3.ortens = V3.ortens/(size(V3.cosine,1));
V3.oteigval=abs(sort(eig(V3.ortens)));

% Compiling mean fabric tensor
Vtot.fabtens=zeros(3,3); % Mean fabric tensor
for i=1:size(V1.cosine,1)
    Vtot.fabtens(1,1) = Vtot.fabtens(1,1)+V1.cosine(i,1)^2+V2.cosine(i,1)^2+V3.cosine(i,1)^2;
    Vtot.fabtens(2,2) = Vtot.fabtens(2,2)+V1.cosine(i,2)^2+V2.cosine(i,2)^2+V3.cosine(i,2)^2;
    Vtot.fabtens(3,3) = Vtot.fabtens(3,3)+V1.cosine(i,3)^2+V2.cosine(i,3)^2+V3.cosine(i,3)^2;
    Vtot.fabtens(1,2) = Vtot.fabtens(1,2)+V1.cosine(i,1)*V1.cosine(i,2)+V2.cosine(i,1)*V2.cosine(i,2)+V3.cosine(i,1)*V3.cosine(i,2);
    Vtot.fabtens(1,3) = Vtot.fabtens(1,3)+V1.cosine(i,1)*V1.cosine(i,3)+V2.cosine(i,1)*V2.cosine(i,3)+V3.cosine(i,1)*V3.cosine(i,3);
    Vtot.fabtens(2,3) = Vtot.fabtens(2,3)+V1.cosine(i,2)*V1.cosine(i,3)+V2.cosine(i,2)*V2.cosine(i,3)+V3.cosine(i,2)*V3.cosine(i,3);
end

Vtot.fabtens(2,1) = Vtot.fabtens(1,2);
Vtot.fabtens(3,1) = Vtot.fabtens(1,3);
Vtot.fabtens(3,2) = Vtot.fabtens(2,3);
Vtot.fabtens=Vtot.fabtens/sum(eig(Vtot.fabtens))*3; % Fabric tensor normalization by its eigenvalues

%Fabric tensor approach. All directions distribution are dependent. Fabric tensor calculated following the approach of Petri et al. (2020).
% Summ of the fabric tensor of each sample. Allows Hext (1953) and Jelinek (1978)
% confidence cones calculation.

    
    [eigvec, eigval] = eig(Vtot.fabtens);
    V1.meancart = [eigvec(1,3) eigvec(2,3) eigvec(3,3)];
    if V1.meancart(3)>0
        V1.meancart=-V1.meancart;
    end
    [V1.meandeg]=Cart2Deg(V1.meancart);
    
    V2.meancart = [eigvec(1,2) eigvec(2,2) eigvec(3,2)];
    if V2.meancart(3)>0
        V2.meancart=-V2.meancart;
    end
    [V2.meandeg]=Cart2Deg(V2.meancart);
    
    V3.meancart = [eigvec(1,1) eigvec(2,1) eigvec(3,1)];
    if V3.meancart(3)>0
        V3.meancart=-V3.meancart;
    end
    [V3.meandeg]=Cart2Deg(V3.meancart);
    
    Vtot.fteigval=[eigval(3,3) eigval(2,2) eigval(1,1)];
    
    [Vtot.plane12, Vtot.plane13, Vtot.plane23] = GCSPlane(V1.meancart, V2.meancart, V3.meancart);

end

2.2 ConfidenceFisher.m

function [a95,Ell] = ConfidenceFisher(Points)
%
% function ConfidenceFisher: Calulates the opening angle of a confidence
% cone around mean direction calculated by vector analysis following Fisher
% (1953). See Tauxe (2002) for synthesis.
% input:    - Points: set of points to analyse, in cartesian coordinates
% output:   - a95: half opening angle of the confidence circle at a 95%
%              level of conficense
%           - ak95: k-parameter of the Fisher statistics, indicates
%           clustered or scattered distributions
% references:
% - Fisher, R., 1953. Dispersion on a sphere. Proc. R. Soc.
% London. Ser. A. Math. Phys. Sci. 217, 295 LP-305.
% - Tauxe, L., 2002. Paleomagnetic principles and practice. Kluwer, Dodrecht.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     This file is part of TOMOFAB. Copyright (C) 2018-2021  Benoit Petri
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.
%      
%     Please report any bug, error or suggestion to bpetri@unistra.fr
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sumx=0;
sumy=0;
sumz=0;
for i=1:size(Points,1)
    sumx=sumx+Points(i,1);
    sumy=sumy+Points(i,2);
    sumz=sumz+Points(i,3);
end
R=sqrt(sumx^2+sumy^2+sumz^2); % Calculate the R parameter, eq 4.2 of Tauxe (2002)
N=size(Points,1); % Number of points
p=0.05; % 95% level of confidence
a95=acosd(1-(N-R)/R*((1/p)^(1/(N-1))-1)); % Calculate the opening half angle, eq 4.4 of Tauxe (2002)
% ak95=(N-1)/(N-R); %% Estimate the degree of scatter, 1 = high scatter, inf = tight cluster

r=sind(a95); % Calculates the circle radius
th = linspace(0,2*pi,100);
xunit = r .* cos(th);
yunit = r .* sin(th);
zpos=cosd(a95); % Calculates the z position of the circle
zunit=-zpos*ones(size(th)); 
Ell=[xunit' yunit' zunit']; % Put everything together

end

2.3 画图展示

导入数据

[V1,V2,V3,Vtot]=ImportFisherVector3d(v1,v2,v3);

步骤1：先绘制lambert投影图

figure
plotlambertbackground()
function plotlambertbackground()
    % Spherical projection plot section
    N = 50;
    cx = cos(0:pi/N:2*pi);
    cy = sin(0:pi/N:2*pi);
    xh = [-1 1];
    yh = [0 0];
    xv = [0 0];
    yv = [-1 1];
    axis([-1 1 -1 1]);
    axis('square');
    fill(cx,cy,'w'); % Plot a white filled circle
    hold on
    plot(xh,yh,xv,yv,'Color', [0.7 0.7 0.7]); % Plot axes
    axis off;
    hold on;
    for Dip = 10:10:80
        DipDir=90;
        [x,y]=GCSchmidt(DipDir,Dip);
        plot(x,y,-x,y,'Color', [0.7 0.7 0.7]); % Plot great circles
    end
    for Dip = 10:10:80
        DipDir=0;
        [x,y]=GCSchmidt(DipDir,Dip);
        plot(x,y,x,-y,'Color', [0.7 0.7 0.7]); % Plot small circles
    end
    plot(cx,cy,'-k'); % Plot black contour on the spherical projection plot
    axis('square');
end

步骤2：绘制数据点

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%绘制数据点
[Points]=Cart2Deg(V1.cosine);
% Equal area case
[xp1,yp1]=SchmidtParams(Points);
plot(xp1,yp1,'or','MarkerSize',2,'MarkerFaceColor','r');
hold on
[Points]=Cart2Deg(V2.cosine);
% Equal area case
[xp2,yp2]=SchmidtParams(Points);
plot(xp2,yp2,'o','MarkerSize',2,'MarkerFaceColor',[0.127 0.746 0.127],'MarkerEdgeColor',[0.127 0.746 0.127]);
hold on
[Points]=Cart2Deg(V3.cosine);
% Equal area case
[xp3,yp3]=SchmidtParams(Points);
plot(xp3,yp3,'ob','MarkerSize',2,'MarkerFaceColor','b');
hold on

步骤3：绘制置信区间轮廓

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%绘制置信区间
% [V1,V2,V3,Vtot]=ConfidenceHext(V1, V2, V3, Vtot);
%fisher method
V1.confang=[]; V2.confang=[]; V3.confang=[];
V1.confell=[]; V2.confell=[]; V3.confell=[];
[V1.confang,V1.confell]=ConfidenceFisher(V1.cosine); % Call Fisher stats calculation for V1
[V2.confang,V2.confell]=ConfidenceFisher(V2.cosine); % Call Fisher stats calculation for V2
[V3.confang,V3.confell]=ConfidenceFisher(V3.cosine); % Call Fisher stats calculation for V3
[V1.confell]=ConfidencePlot(V1.confell,V1.meandeg,'r'); % Send ellipses coordinates to plotting function
[V2.confell]=ConfidencePlot(V2.confell,V2.meandeg,[0.127 0.746 0.127]); % Send ellipses coordinates to plotting function
[V3.confell]=ConfidencePlot(V3.confell,V3.meandeg,'b'); % Send ellipses coordinates to plotting function

步骤4：绘制平均方向及平面

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%绘制平均方向及平面
% plot the great circles related to opening angles
[x,y]=GCSchmidt(Vtot.plane12(1),Vtot.plane12(2));
plot(x,y,'--k');
hold on
clearvars x y
[x,y]=GCSchmidt(Vtot.plane13(1),Vtot.plane13(2));
plot(x,y,'--k');
hold on
clearvars x y
[x,y]=GCSchmidt(Vtot.plane23(1),Vtot.plane23(2));
plot(x,y,'--k');
clearvars x y
hold on
%%%
Points=V1.meandeg; % For V1 mean direction
[xp,yp]=SchmidtParams(Points);    
plot(xp,yp,'sk','MarkerFaceColor',[1 0 0],'MarkerSize',9);
clearvars Points xp yp
Points=V2.meandeg; % For V2 mean direction
[xp,yp]=SchmidtParams(Points);
plot(xp,yp,'^k','MarkerFaceColor',[0.127 0.746 0.127],'MarkerSize',8);
clearvars Points xp yp
Points=V3.meandeg; % For V3 mean direction
[xp,yp]=SchmidtParams(Points);
plot(xp,yp,'ok','MarkerFaceColor',[0 0 1],'MarkerSize',8);
clearvars Points xp yp

附录

几个公用的文件

1. Cart2Deg.m [PointsDeg] = Cart2Deg(PointsCart) 将笛卡尔坐标n×3的输入，转为n×2的地理角度输出（DipDir / Dip）
2. ConfidencePlot [Ell]=ConfidencePlot(Ell,AxisCoord,Color) Plots confidence ellipses on a spherical projection plot
3. GCSchmidt.m [x,y] = GCSchmidt(DipDir,Dip) 该函数针对平面，知道平面的DipDir,Dip，将其转为lambert投影平面上的坐标，从而绘制大圆。
4. GCSPlane [V1V2Deg, V1V3Deg, V2V3Deg] = GCSPlane(V1Cart, V2Cart, V3Cart)知道三个坐标轴的直角坐标，拿到三个平面（xy、xz、yz）的DipDir,Dip
   4.SchmidtParams [xp,yp] = SchmidtParams(Points) 针对散点的绘制，将角度坐标的散点改成lambert投影平面下的坐标

Cart2Deg.m

function [PointsDeg] = Cart2Deg(PointsCart)
%
% function Cart2Deg: Transforms directional cosines (cartesian coordinates)
% to Dip direction / Dip angles coordinates
% input:    - PointsCart: Coordinates of points in cartesian coordinates
% output:   - PointsDeg: Coordinates of points in DipDir / Dip
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     This file is part of TOMOFAB. Copyright (C) 2018-2021  Benoit Petri
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.
%      
%     Please report any bug, error or suggestion to bpetri@unistra.fr
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PointsSph=zeros(size(PointsCart,1),3);
PointsDeg=zeros(size(PointsCart,1),2);

for i=1:size(PointsCart,1)
    [PointsSph(i,1), PointsSph(i,2), PointsSph(i,3)]=cart2sph(PointsCart(i,1),PointsCart(i,2),PointsCart(i,3));
    PointsDeg(i,1)=round(-PointsSph(i,1)*180/pi+90,1);
    PointsDeg(i,2)=round(-PointsSph(i,2)*180/pi,1);
    if PointsDeg(i,1)<0
        PointsDeg(i,1)=PointsDeg(i,1)+360;
    end
end

ConfidencePlot.m

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

function [Ell]=ConfidencePlot(Ell,AxisCoord,Color)
%
% function ConfidencePlot: Plots confidence ellipses on a spherical
% projection plot
% input:    - Ell: coordinatees of points constituting the confidence
%           ellipses
%           - AxisCoord: coordinates of the axis anchoring the ellipse
%           - Color: color of the ellipse
%           - chb1: requested projection
% output:   - Ell: recalculated ellipse coordinates
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     This file is part of TOMOFAB. Copyright (C) 2018-2021  Benoit Petri
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.
%      
%     Please report any bug, error or suggestion to bpetri@unistra.fr
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Ell, revrec]=Rotator([AxisCoord 0 0],Ell,1); % anchor the ellipse to its right point
[EllDeg]=Cart2Deg(Ell);
[xel,yel]=SchmidtParams(EllDeg);% Calculates the coordinates on the spherical projection plot


fsize=1; 
for k=2:size(xel) % loop to quantify the number of rows requested to draw nice ellipses
    if revrec(k-1) ~= revrec(k)
        fsize=fsize+1;
    end
end

cutxel=cell(fsize,1);
cutyel=cell(fsize,1);

cutxel{1}=xel(1);
cutyel{1}=yel(1);
nameinc=1;

for k=2:size(xel) % Cut the ellipse coordinates into indepentent curves to avoid connecting lines if they are separated in different segments
    if revrec(k-1) == revrec(k)
        A = cutxel{nameinc};
        B = cutyel{nameinc};
        cutxel{nameinc}=[A xel(k)];
        cutyel{nameinc}=[B yel(k)];
    else
        nameinc=nameinc+1;
        A = cutxel{nameinc};
        B = cutyel{nameinc};
        cutxel{nameinc}=[A xel(k)];
        cutyel{nameinc}=[B yel(k)];
    end
end

for i=1:fsize % Plot all segments
    plot(cutxel{i},cutyel{i},'Color',Color) 
end


end



function [Data,Revrec]=Rotator(CoreOr, Data, IsCos)
%
% function Rotator: Apply several rotations to the dataset to calculate cartesian
% coordinates in the geographic reference frame.
% input:    - CoreOr: core orientation (DipDir, Dip, reverse position, misorientation angle)
%           - Data: data onto which rotations should be applied.
%           - IsCos: 1 if input data are directional cosines, 0 if not.
% output:   - Data: set of rotated points in cartesian coordinates
%           - Revrec: array recording if points had z > 0 after rotation.
%           Important for ellipses segmentation during plotting.
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     This file is part of TOMOFAB. Copyright (C) 2018-2021  Benoit Petri
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.
%      
%     Please report any bug, error or suggestion to bpetri@unistra.fr
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Revrec=zeros(size(Data,1),1);

CoreDeg = [CoreOr(1) CoreOr(2)];
RevPosition = CoreOr(3);
MisAngle = CoreOr(4);


%% Change of coordinates system
CoreDeg = -CoreDeg;
MisAngle = -MisAngle;

for i=1:size(Data,1)
    DataCart=[Data(i,1) Data(i,2) Data(i,3)];
    
    %%
    %% Rotation phase 1: correction for acquisition misorientation (core centering and bottom down reorientation)
    %%
    
    %%
    %% Rotation 1: correction for sample misorientation: around Z, amount = misorientation angle
    %%
    
    thetaR1 = -MisAngle; % Rotation angle
    thetaR1 = -pi*(thetaR1)/180; % transformation into rad
    R1= ([cos(thetaR1) -sin(thetaR1) 0 ; sin(thetaR1) cos(thetaR1) 0 ; 0 0 1]); % Rotation matrix
    
    DataCartT1 = R1*DataCart';

    
    
    %%
    %% Rotation 2: correction for bottom up acquisition: for reversed samples, around X, amount = 180 
    %%
    
    if RevPosition == 1
        u2 = [1 0 0]; % Unitary vector of the rotation axis
        thetaR2 = -180; % Rotation angle
        thetaR2 = -pi*(thetaR2)/180; % transformation into rad
        c2 = cos(thetaR2);
        s2 = sin(thetaR2);
        R2 = ([u2(1)^2*(1-c2)+c2 u2(1)*u2(2)*(1-c2)+u2(3)*s2 u2(1)*u2(3)*(1-c2)-u2(2)*s2 ; u2(1)*u2(2)*(1-c2)-u2(3)*s2 u2(2)^2*(1-c2)+c2 u2(2)*u2(3)*(1-c2)+u2(1)*s2 ; u2(1)*u2(3)*(1-c2)+u2(2)*s2 u2(2)*u2(3)*(1-c2)-u2(1)*s2 u2(3)^2*(1-c2)+c2]); % Rotation matrix
        
        DataCartT2 = R2*DataCartT1;
        DataCartT2 = -DataCartT2;
    else
        DataCartT2 = DataCartT1;
    end
    
    
    %%
    %% Rotation phase 2: core and axes reorientation in geographic reference frame
    %%
    
    %%
    %% Rotation 3: Strike correction: along the z axis, amount = dip direction
    %%
    
    
    thetaR3 = -CoreDeg(1); % Rotation angle
    thetaR3 = -pi*(thetaR3)/180; % transformation into rad
    R3= ([cos(thetaR3) -sin(thetaR3) 0 ; sin(thetaR3) cos(thetaR3) 0 ; 0 0 1]); % Rotation matrix

    DataCartT3 = R3*DataCartT2;
    
    
    %%
    %% Rotation 4: Dip correction: along the azimuth, amount = dip
    %%
    
    [u4(1),u4(2),u4(3)] = sph2cart(pi*(CoreDeg(1))/180, 0, 1); % Unitary vector of the rotation axis
    thetaR4 = -CoreDeg(2); % Rotation angle
    thetaR4 = -pi*(90-thetaR4)/180; % transformation into rad
    c4 = cos(thetaR4);
    s4 = sin(thetaR4);
    R4 = ([u4(1)^2*(1-c4)+c4 u4(1)*u4(2)*(1-c4)+u4(3)*s4 u4(1)*u4(3)*(1-c4)-u4(2)*s4 ; u4(1)*u4(2)*(1-c4)-u4(3)*s4 u4(2)^2*(1-c4)+c4 u4(2)*u4(3)*(1-c4)+u4(1)*s4 ; u4(1)*u4(3)*(1-c4)+u4(2)*s4 u4(2)*u4(3)*(1-c4)-u4(1)*s4 u4(3)^2*(1-c4)+c4]); % Rotation matrix
    DataCartT4 = R4*DataCartT3;
    
       
    if IsCos == 1 && DataCartT4(3)>0 % Correction for upper hemisphere data, noted -1 in Revrec
        DataCartT4 = -DataCartT4;
        Revrec(i)=-1;
    end
    
    Data(i,1)=round(DataCartT4(1),6);
    Data(i,2)=round(DataCartT4(2),6);
    Data(i,3)=round(DataCartT4(3),6);

end

return
end

GCSchmidt.m

function [x,y] = GCSchmidt(DipDir,Dip)
%
% function GCWulff: Calculate coordinates of a great circle on an equal-area spherical projection plot (Schmidt net)
% input: - DipDir, Dip: Coordinates of the plane
% output: - x, y: coordinates of the great circle that will be plot
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     This file is part of TOMOFAB. Copyright (C) 2018-2021  Benoit Petri
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.
%      
%     Please report any bug, error or suggestion to bpetri@unistra.fr
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

psi=0:pi/50:pi;
rdip = Dip*(pi/180);
radip = atan(tan(rdip)*sin(psi));
rproj = sqrt(2)*sin((pi/2 - radip)/2);
x = rproj .* sin(psi);
y = rproj .* cos(psi);

thetaR1 = -90+DipDir; % Rotation angle
thetaR1 = -pi*(thetaR1)/180; % Transformation into rad
rot= ([cos(thetaR1) -sin(thetaR1) 0 ; sin(thetaR1) cos(thetaR1) 0 ; 0 0 1]); % Rotation matrix

for i=1:size(x,2)
    tmp= rot*[x(i) y(i) 0]';
    x(i)=tmp(1);
    y(i)=tmp(2);
end

end

GCSPlane.m

function [V1V2Deg, V1V3Deg, V2V3Deg] = GCSPlane(V1Cart, V2Cart, V3Cart)
%
%根据三个主方向的笛卡尔坐标系，拿到三个平面在GCS下的表示
% function Fabrics: Calculates the V1-V2, V1-V3, V2-V3 planes
% input:    - V1, V2, V3 in cartesian coordinates
% output:   - planes dip direction / dip angle in degrees
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     This file is part of TOMOFAB. Copyright (C) 2018-2021  Benoit Petri
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.
%      
%     Please report any bug, error or suggestion to bpetri@unistra.fr
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%% Coordinates transformation from cartesian to Dip Dir and Dip
%%
V1Sph = [0 0 0];
[V1Sph(1), V1Sph(2), V1Sph(3)]=cart2sph(V1Cart(1),V1Cart(2),V1Cart(3));
V2Sph = [0 0 0];
[V2Sph(1), V2Sph(2), V2Sph(3)]=cart2sph(V2Cart(1),V2Cart(2),V2Cart(3));
V3Sph = [0 0 0];
[V3Sph(1), V3Sph(2), V3Sph(3)]=cart2sph(V3Cart(1),V3Cart(2),V3Cart(3));
V1Deg = [0 0];
V1Deg(1)=-V1Sph(1)*180/pi+90;
V1Deg(2)=-V1Sph(2)*180/pi;
V2Deg = [0 0];
V2Deg(1)=-V2Sph(1)*180/pi+90;
V2Deg(2)=-V2Sph(2)*180/pi;
V3Deg = [0 0];
V3Deg(1)=-V3Sph(1)*180/pi+90;
V3Deg(2)=-V3Sph(2)*180/pi;
if V1Deg(1)<0
    V1Deg(1)=V1Deg(1)+360;
end
if V2Deg(1)<0
    V2Deg(1)=V2Deg(1)+360;
end
if V3Deg(1)<0
    V3Deg(1)=V3Deg(1)+360;
end

%% V1V2 plane calculation
V1V2Deg = zeros(size(V3Deg));

NormVec=cross(V1Cart,V2Cart);
[FolDipDir, FolDip, ~]=cart2sph(NormVec(1),NormVec(2),NormVec(3));

V1V2Deg(1) = FolDipDir;
V1V2Deg(2) = FolDip;

if V1V2Deg(2) <= 0
    V1V2Deg(1)=-V1V2Deg(1)*180/pi+90+180;
    V1V2Deg(2)=90+V1V2Deg(2)*180/pi;
else
    V1V2Deg(1)=V1V2Deg(1)-pi;
    V1V2Deg(2)=-V1V2Deg(2);
    V1V2Deg(1)=-V1V2Deg(1)*180/pi+90+180;
    V1V2Deg(2)=90+V1V2Deg(2)*180/pi;
end

if V1V2Deg(1) > 360
    V1V2Deg(1) = V1V2Deg(1)-360;
end

V1V2Deg = round(V1V2Deg,2);


%% V1V3 plane calculation
V1V3Deg = zeros(size(V2Deg));

NormVec=cross(V1Cart,V3Cart);
[FolDipDir, FolDip, ~]=cart2sph(NormVec(1),NormVec(2),NormVec(3));


V1V3Deg(1) = FolDipDir;
V1V3Deg(2) = FolDip;

if V1V3Deg(2) <= 0
    V1V3Deg(1)=-V1V3Deg(1)*180/pi+90+180;
    V1V3Deg(2)=90+V1V3Deg(2)*180/pi;
else
    V1V3Deg(1)=V1V3Deg(1)-pi;
    V1V3Deg(2)=-V1V3Deg(2);
    V1V3Deg(1)=-V1V3Deg(1)*180/pi+90+180;
    V1V3Deg(2)=90+V1V3Deg(2)*180/pi;
end

if V1V3Deg(1) > 360
    V1V3Deg(1) = V1V3Deg(1)-360;
end

V1V3Deg = round(V1V3Deg,2);


%% V2V3 plane calculation
V2V3Deg = zeros(size(V1Deg));

NormVec=cross(V2Cart,V3Cart);
[FolDipDir, FolDip, ~]=cart2sph(NormVec(1),NormVec(2),NormVec(3));

V2V3Deg(1) = FolDipDir;
V2V3Deg(2) = FolDip;

if V2V3Deg(2) <= 0
    V2V3Deg(1)=-V2V3Deg(1)*180/pi+90+180;
    V2V3Deg(2)=90+V2V3Deg(2)*180/pi;
else
    V2V3Deg(1)=V2V3Deg(1)-pi;
    V2V3Deg(2)=-V2V3Deg(2);
    V2V3Deg(1)=-V2V3Deg(1)*180/pi+90+180;
    V2V3Deg(2)=90+V2V3Deg(2)*180/pi;
end

if V2V3Deg(1) > 360
    V2V3Deg(1) = V2V3Deg(1)-360;
end

V2V3Deg = round(V2V3Deg,2);
end

SchmidtParams.m

function [xp,yp] = SchmidtParams(Points)
%
% function SchmidtParams: Calculates coordinates of points on an equal-area spherical projection plot (Schmidt net)
% input: - Points directions in degrees
% output: - xp, yp: coordinates of points on the stereonet
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     This file is part of TOMOFAB. Copyright (C) 2018-2021  Benoit Petri
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.
%      
%     Please report any bug, error or suggestion to bpetri@unistra.fr
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

theta = pi*(Points(:,1))/180;
rho = pi*(Points(:,2))/180;
xp=zeros(size(theta));
yp=zeros(size(theta));
for i=1:size(theta)
    trd=theta(i);
    plg=rho(i);
    if plg < 0.0
        trd = ZeroTwoPi(trd+pi);
        plg = -plg;
    end
    piS4 = pi/4.0;
    s2 = sqrt(2.0);
    plgS2 = plg/2.0;
    xpt = s2*sin(piS4 - plgS2)*sin(trd);
    ypt = s2*sin(piS4 - plgS2)*cos(trd);
    xp(i)=xpt;
    yp(i)=ypt;
end

end