function [q1,q2]=dpfg5design(e1,e2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Written by 
% Sune Nørhøj Jespersen
% CFIN/MindLab and Dept. of Physics and Astronomy
% Aarhus Universitet
% Nørrebrogade 44, bygn 10G, 5. sal
% 8000 Århus C
% Denmark
% Phone: +45 78463334
% Cell: +45 60896642
% E-mail: sune@cfin.au.dk
% Web: http://www.cfin.au.dk/~sune
%
% This Matlab function computes 60 pairs of diffusion directions for
% rotationally invariant sampling of double PFG MR diffusion experiments,
% as described in ref. 2. This system is based on the theory of
% quadrature on the sphere, more specifically the framework proposed in
% Ref. 1. It uses 60 rotations to perform 60 rotatiosn of the input
% vectors. The input vectors can be arbitrary unit vectors. Free usage is
% permitted, but please acknowledge Refs 1 and 2.
%
% % References:
% 1) Graf, M., Potts, D., 2009. Sampling Sets and Quadrature Formulae on the Rotation Group. Numerical
% Functional Analysis and Optimization 30, 665-688
% 2) Jespersen, S., Lundell, H., Sønderby C. K., Dyrby, T.B. NMR in
% Biomedicine (2013).
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Input arguments: 2 unit vectors upon which the the 60 rotations will
%  act. If no input is supplied, these 2 vectors will default to (1,0,0)
%  and (0,0,1).
% Output: 60 pairs of diffusion directions.

if nargin==0
    e1=[1,0,0]';
    e2=[0,0,1]';
elseif nargin==2
    e1=e1(:);
    e2=e2(:);
else    error('Syntax error. See help');
    return;
end

phi=[0 2*pi/5*[1:5] (pi/5+2*pi/5*[1:5]) 0];
theta=[0 ones(1,5)*atan(2) pi-ones(1,5)*atan(2) pi];
psi=zeros(5,12);
for i=1:12
    for k=1:5
        psi(k,i)=pi/5*not(i==1)*not(i==12)+2/5*pi*(k-1)-phi(i);
    end
end
phi=repmat(phi,5,1);
theta=repmat(theta,5,1);
phi=reshape(phi,[1,60]);
theta=reshape(theta,[1,60]);
psi=reshape(psi,[1,60]);

q1=zeros(3,60);
q2=q1;
R=zeros(3,3,60);
for i=1:60;
R1(:,:,i)=[cos(phi(i)) -sin(phi(i)) 0; sin(phi(i)) cos(phi(i)) 0; 0 0 1];
R2(:,:,i)=[cos(theta(i)) 0 sin(theta(i));0 1 0; -sin(theta(i)) 0 cos(theta(i))];
R3(:,:,i)=[cos(psi(i)) -sin(psi(i)) 0; sin(psi(i)) cos(psi(i)) 0; 0 0 1];
R(:,:,i)=R1(:,:,i)*R2(:,:,i)*R3(:,:,i);
        q1(:,i)=R(:,:,i)*e1;
        q2(:,i)=R(:,:,i)*e2;
end

    
