function [positionMatrix probabilityVector] = sampleRay(...
    nSample, ray, radius)
% sample from n-dim cone (disc base)
%
% Syntax: [pointMatrix, probabilityVector] = sampleNormal(nSample, ray, radius)
%
% THIS IS NO USER FUNCTION

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2013 Alexander Reinhold
%
% 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 (at your 
%   option) 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 <http://www.gnu.org/licenses/>.

%% Algorithm
% A point of the cone can be sampled by first determine the level at
%   the ray, next the distance from the ray and finally a point at the
%   boundary for the direction. This topic is handled for n-dimensional
%   cases.
% There must be linearly increasing more points at the base of the cone
%   than at the tip. The probability distribution for points at the ray
%   (l=0 -> tip, l=1 -> base) can be given by:
%      P(l) = 2*l
%   so that the cumulative density function is given by:
%      PC(l) = l^2
%   from which can be sampled using a uniform distribution for PC to
%   get l:
%      l = sqrt(PC)
% The distance from the ray shall also be chosen in consideration of
%   uniformly distributed points in the cone. The sphere at the
%   corresponding radius increases its surface polynomial with r^(n-2)
%   so that the radius is chosen from the uniformly distributed
%   variable PC:
%     r = PC^(1/(n-2))
%   Special considerations are required for n=1.
% The direction that is given by a point of the boundary of the
%   crosssectional hypersphere is detemined by a formula that can be
%   found in D. E. Knuth "The art of computer Programming, vol. 2
%   Seminumerical algorithms". Therefore, (n-1) normally distributed
%   variables x are used to obtain the points of a (n-1)-sphere
%   embedded in the n-dimensional space:
%     y = [x / norm(x); 0];
% We that need a rotation matrix that transforms the n-th axis to the
%   ray that will be U and obtain the final random points:
%     pivotPosition = U*(r*l)*y + l*ray;

%% input handling
% adjust to column vector
ray = ray(:);
% problem dimension
nDim = length(ray);

%% required random numbers
% radial distance from cone
r = rand(1, nSample).^(1/(nDim-1));
% position on ray (length)
l = rand(1, nSample).^(1/nDim);
% normal distributed points (equally distributed over angle of cone if
%   projected to unit length)
x = randn(nDim-1, nSample);

%% building a cone (direction to x-axis
% ray along x-axis, radius 1
positionMatrix = [l*norm(ray); bsxfun(@times, radius*(r)./sqrt(sum(x.*x, 1)).*l, x)];
rayUnit = zeros(nDim, 1);
rayUnit(1) = 1;

%% rotation matrix (U)
[U, S, V] = svd([ray, zeros(nDim, nDim-1)]);
% correct direction of rotation - svd does not care about the sign of new
%   coordinates, but here we do. Tolerance is not essential, however.
if ~all(isZero(ray - U*rayUnit*norm(ray), 1e-10*norm(ray)))
    U = -1*U;
end
positionMatrix = U*positionMatrix;

%% density
volume = computeConeVolume(ray, radius);
probabilityVector = positionMatrix(1,:)*0 + 1/volume;