function [fcc2bcc,bcc2fcc_odf] = determineor_triple(symmetries,ebsd,cond)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
% 
% symmetries{1} = ebsd.CSList{2};
% symmetries{2} = ebsd.CSList{3};

w = warning ('off','all');

symmetry1_string = symmetries{1}.mineral;

% Determine neighboring point pairs:
[V,F,I_FD] = spatialDecomposition([ebsd(symmetry1_string).prop.x(:), ebsd(symmetry1_string).prop.y(:)],ebsd(symmetry1_string).unitCell);
A_D = I_FD'*I_FD==1;
[Dl,Dr] = find(triu(A_D,1)); % lists of neighboring point pairs

[i,j] = find(A_D);  % list of all cells

% Reduce all orientations to quaternions:
O = quaternion(ebsd(symmetry1_string).orientations);

% Determine experimental misorientation set:
m_all = inv(O(Dl)).*O(Dr);

%%
% Determine the initial orientation relationship as Kurdjumov-Sachs:
fcc2bcc_KS = (orientation('map',Miller(1,1,1,symmetries{2}),Miller(0,1,1,symmetries{1}),Miller(-1,0,1,symmetries{2}),Miller(-1,-1,1,symmetries{1})));

%% Set up the necessary symmetries and the simulated intermartensitic misorientations:

%The assumption is that we're dealing with misorientations between cubic symmetry:
cs = crystalSymmetry('432');%432 is used instead of m3m because it only gives 24 symmetry operators, which is enough

% Set up symmetry operators in the order that is assumed to produce a
% misorientation set similar to that used by Morito et al. and others:
cs_symmetries = rotation(cs.properGroup);

symrots_true(1,1) = cs_symmetries.rot(1);
symrots_true(2,1) = cs_symmetries.rot(21);
symrots_true(3,1) = cs_symmetries.rot(11);
symrots_true(4,1) = cs_symmetries.rot(16);
symrots_true(5,1) = cs_symmetries.rot(24);
symrots_true(6,1) = cs_symmetries.rot(8);
symrots_true(7,1) = cs_symmetries.rot(18);
symrots_true(8,1) = cs_symmetries.rot(10);
symrots_true(9,1) = cs_symmetries.rot(20);
symrots_true(10,1) = cs_symmetries.rot(3);
symrots_true(11,1) = cs_symmetries.rot(7);
symrots_true(12,1) = cs_symmetries.rot(23);
symrots_true(13,1) = cs_symmetries.rot(19);
symrots_true(14,1) = cs_symmetries.rot(2);
symrots_true(15,1) = cs_symmetries.rot(9);
symrots_true(16,1) = cs_symmetries.rot(22);
symrots_true(17,1) = cs_symmetries.rot(17);
symrots_true(18,1) = cs_symmetries.rot(12);
symrots_true(19,1) = cs_symmetries.rot(5);
symrots_true(20,1) = cs_symmetries.rot(15);
symrots_true(21,1) = cs_symmetries.rot(4);
symrots_true(22,1) = cs_symmetries.rot(14);
symrots_true(23,1) = cs_symmetries.rot(6);
symrots_true(24,1) = cs_symmetries.rot(13);
symrots_true = rotation(symrots_true);

% Determine the full orientation relationship according to cubic symmetry
% as by Morito et al.:
fcc2bcc_true = orientation(symrots_true*rotation(fcc2bcc_KS),symmetries{1},1);

%%
l = 1;

if cond
    % Prepare a figure to display the evolution of the angular deviations between
    % close-packed planes and directions (see reconstruction articles by
    % Miyamoto et al.):
    figure
    h_1(1) = animatedline('Color','r');
    hold on
    h_1(2) = animatedline('Color','g');
    axis([0 20 0 5])
    z1 = legend(h_1,'Angle between {111}a and {011}m','Angle between <110>a and <111>m');
    set(z1,'interpreter','none');
end
%% Orientation relationship determination algorithm


while  l<3 || (angle_planes(end)-angle_planes(end-1))>0.01 % Determine end condition as reasonably diminished angular deviation between close-packed planes
    
    %Downsample the misorientation list for speed:
    if length(m_all)>100000
        g1 = m_all(randsample(length(m_all),100000));
    else
        g1 = m_all;
    end
    % Determine intermartensitic misorientations:
    bcc_trans =  inv(orientation(rotation(fcc2bcc_true),symmetries{1},1)) * (orientation(rotation(fcc2bcc_true(1)),symmetries{1},1));
    g2 = quaternion(bcc_trans);

    % Determine intermartensitic misorientations
    
    %Using a modified dot_outer_quat_cs, get the dot products, along with the index of the maximum dot product:
    [ind,d] = dot_outer_quat_cs_modified(g2,g1,cs,cs);
    
    %Calculate the angles corresponding to each dot product:
    omega_1 = [real(2*acos(abs(d))) / degree];
    
    % Get the minimum angle and index corresponding to the best-fitting
    % simulated intermartensitic misorientation:
    [omega_min_1, omega_min_1_I] = min(omega_1,[],2);
    
%     condition = omega_min_1<5 & omega_min_1_I~=1 & omega_min_1_I~=4;
        condition = omega_min_1<5 & omega_min_1_I~=1;

    % There are many misorientations corresponding to intragranular
    % misorientations (omega_min_1_I == 1). There are also some misorientations
    % that do not correspond to intermartensitic misorientations
    % (omega_min_1>arbitrary_cutoff_value). Let us remove these from consideration:
    cutoff_g1 = g1(condition);
    cutoff_ind = ind(find(condition),:);
    cutoff_omega_min_1 = omega_min_1(condition);
    cutoff_omega_min_1_I = omega_min_1_I(condition);
    
    % Determine all necessary symmetry operators
    
    %Using omega_min_1_I and its row number, we can pick the corresponding index from ind. First,
    %we need to alter omega_min_1_I to a corresponding linear index for ind:
    linearind_right_stuff = sub2ind(size(cutoff_ind),... %size of ind matrix
        [1:length(cutoff_omega_min_1)]',... % the row of each omega_min_1_I value
        [cutoff_omega_min_1_I]); % the index of the best-fitting intermartensitic misorientation on each row
    
    % dot_outer_quat_cs_modified works by picking the best option out of a
    % 24x24 matrix, from which we need the column value. All we get from ind is a linear index, however. %The rem command gives the division remainder for each linear index when divided by 24, giving the column value:
    ind_necessary_2 = rem(cutoff_ind(linearind_right_stuff),24);
    %The zero value given by the rem command actually corresponds to the 24th
    %symmetry operator:
    ind_necessary_2(ind_necessary_2 == 0) = 24;
    
    %To obtain the final necessary symmetry operator, we need to go back to the list of indexes given by dot_outer_quat_modified. This time we need the row,
    %which we can get by the command:
    ind_necessary_3 = ceil(cutoff_ind(linearind_right_stuff)/24);
    
    % Using the symmetry operators
    
    %We remove one symmetry problem using symmetry operator 2:
    primary_misos = symrots_true(ind_necessary_2).*orientation(cutoff_g1,cs,1);
        
    %The misorientations can be translated to the primary orientation relationship using symmetry operator 1:
    fcc2bcc_candidate = rotation(fcc2bcc_true(cutoff_omega_min_1_I)).*primary_misos;
    
    fcc2bcc_true = mean(fcc2bcc_candidate);
    fcc2bcc_true = symrots_true*project2FundamentalRegion(fcc2bcc_true,fcc2bcc_KS);
    
    % Calculate the angular deviations between close-packed planes and
    % angles:
    fcc2bcc_basis_miso = orientation(rotation(fcc2bcc_true(1)),symmetries{1},symmetries{2})
    angle_planes(l) = min(angle(fcc2bcc_basis_miso * (symrots_true*Miller(1,1,-1,symmetries{1})),Miller(0,1,1,symmetries{2}))) / degree;
    angle_directions(l) = min(angle(fcc2bcc_basis_miso * (symrots_true*Miller(1,0,1,symmetries{1},'uvw')),Miller(1,1,-1,symmetries{1},'uvw'))) / degree;

    % Update the graph:
    if cond
        addpoints(h_1(1),l,angle_planes(l));
        addpoints(h_1(2),l,angle_directions(l));
        drawnow
    end
    l = l+1;

end

fcc2bcc = fcc2bcc_true(1);
bcc2fcc_candidate = orientation(rotation(inv(fcc2bcc_candidate.project2FundamentalRegion(fcc2bcc_KS))),symmetries{2},1);
bcc2fcc_odf = calcODF(orientation(rotation(bcc2fcc_candidate),symmetries{2},1),'halfwidth',1*degree);

%%
    function [ind,d] = dot_outer_quat_cs_modified(g2,g1,cs,ss)
        % quick version that ignores inversion
        
        cs_symmetries = rotation(cs.properGroup);
        
        symrots_true(1,1) = cs_symmetries.rot(1);
        symrots_true(2,1) = cs_symmetries.rot(21);
        symrots_true(3,1) = cs_symmetries.rot(11);
        symrots_true(4,1) = cs_symmetries.rot(16);
        symrots_true(5,1) = cs_symmetries.rot(24);
        symrots_true(6,1) = cs_symmetries.rot(8);
        symrots_true(7,1) = cs_symmetries.rot(18);
        symrots_true(8,1) = cs_symmetries.rot(10);
        symrots_true(9,1) = cs_symmetries.rot(20);
        symrots_true(10,1) = cs_symmetries.rot(3);
        symrots_true(11,1) = cs_symmetries.rot(7);
        symrots_true(12,1) = cs_symmetries.rot(23);
        symrots_true(13,1) = cs_symmetries.rot(19);
        symrots_true(14,1) = cs_symmetries.rot(2);
        symrots_true(15,1) = cs_symmetries.rot(9);
        symrots_true(16,1) = cs_symmetries.rot(22);
        symrots_true(17,1) = cs_symmetries.rot(17);
        symrots_true(18,1) = cs_symmetries.rot(12);
        symrots_true(19,1) = cs_symmetries.rot(5);
        symrots_true(20,1) = cs_symmetries.rot(15);
        symrots_true(21,1) = cs_symmetries.rot(4);
        symrots_true(22,1) = cs_symmetries.rot(14);
        symrots_true(23,1) = cs_symmetries.rot(6);
        symrots_true(24,1) = cs_symmetries.rot(13);
        
        
        % g2rot = symmetrise_modified(quaternion(g2),cs,ss).'; % g2 x CS x SS
        g2rot = (inv(symrots_true)*((g2*symrots_true).'));
        % g2rot = (inv(symrots_true)*((g2*symrots_true).'));
        g2rot = reshape(g2rot,24 * 24,[]);
        g2rot = g2rot.';
        
        q1 = [g1.a(:) g1.b(:) g1.c(:) g1.d(:)];
        a2 = g2rot.a; b2 = g2rot.b; c2 = g2rot.c; d2 = g2rot.d;
        
        % this is implicite dot_outer
        d = abs(q1 * [a2(:,1).';b2(:,1).';c2(:,1).';d2(:,1).']);
        % allocate matrix for indices:
        zalathar = ones(length(d),24);
        
        for k=2:length(cs)*length(ss)
            d_compare = d;
            d = max(d,abs(q1 * [a2(:,k).';b2(:,k).';c2(:,k).';d2(:,k).'])); % g1 x g2 x CS * SS
            zalathar(d_compare<d) = k;
        end
        ind = zalathar;
    end

    function varargout = symmetrise_modified(q,CS,SS,varargin)
        % symmetrcially equivalent orientations
        %
        % Input
        %  q  - @quaternion
        %  CS - crystal @symmetry
        %  SS - specimen @symmetry
        %
        % Output
        %  q - symmetrically equivalent orientations CS x SS x q
        %
        % See also
        % CrystalSymmetries
        
        % q = (q * CS).'; % CS x M
        % if nargin>2 && length(SS)>1
        %   q = SS * q;     % SS x (CS X M)
        %   lSS = length(SS);
        % else
        %   lSS = 1;
        % end
        
        % Symmetry order correction according to Morito (this is not strictly necessary unless variant pair selection is of interest):
        cs_symmetries = rotation(CS.properGroup);
        
        symrots_true(1,1) = cs_symmetries.rot(1);
        symrots_true(2,1) = cs_symmetries.rot(21);
        symrots_true(3,1) = cs_symmetries.rot(11);
        symrots_true(4,1) = cs_symmetries.rot(16);
        symrots_true(5,1) = cs_symmetries.rot(24);
        symrots_true(6,1) = cs_symmetries.rot(8);
        symrots_true(7,1) = cs_symmetries.rot(18);
        symrots_true(8,1) = cs_symmetries.rot(10);
        symrots_true(9,1) = cs_symmetries.rot(20);
        symrots_true(10,1) = cs_symmetries.rot(3);
        symrots_true(11,1) = cs_symmetries.rot(7);
        symrots_true(12,1) = cs_symmetries.rot(23);
        symrots_true(13,1) = cs_symmetries.rot(19);
        symrots_true(14,1) = cs_symmetries.rot(2);
        symrots_true(15,1) = cs_symmetries.rot(9);
        symrots_true(16,1) = cs_symmetries.rot(22);
        symrots_true(17,1) = cs_symmetries.rot(17);
        symrots_true(18,1) = cs_symmetries.rot(12);
        symrots_true(19,1) = cs_symmetries.rot(5);
        symrots_true(20,1) = cs_symmetries.rot(15);
        symrots_true(21,1) = cs_symmetries.rot(4);
        symrots_true(22,1) = cs_symmetries.rot(14);
        symrots_true(23,1) = cs_symmetries.rot(6);
        symrots_true(24,1) = cs_symmetries.rot(13);
        
        q = (q * symrots_true).'; % CS x M
        q = symrots_true * q;     % SS x (CS X M)
        
        varargout{1} = reshape(q,24 * 24,[]); % (CSxSS) x M
    end

 end