samplingSphericalTriangles.jl

# "unitize" a vector
normalize(x) = x/norm(x)
# the [x|y] function
normcomp(x,y) = normalize(x - dot(x,y)*y)
# find the dihedral angle at a vertice `o`
dihedral(a,o,b) = acos(dot(normalize(cross(o,a)),normalize(cross(o,b))))
# the vertices of the triangle in the xy-plane
A2 = Float64[.1,.0,-2]
B2 = Float64[2,-.1,-2]
C2 = Float64[1.5,.9,-2]
# the vertices of the triangle on the sphere's surface
A = normalize(A2)
B = normalize(B2)
C = normalize(C2)
# the dihedral angles
α = dihedral(B,A,C)
β = dihedral(A,B,C)
γ = dihedral(A,C,B)
# the edge length
c = acos(dot(A,B))
# the triangle area
AA = α + β + γ - π
# some constants
cosα = cos(α)
sinα = sin(α)
cosc = cos(c)
normcompCA = normcomp(C,A)
# the main function
function sampleTriangle(ξ1::Float64,ξ2::Float64)
	# find a point on the sphere surface within a triangular window given two random numbers [0,1]
	AAhat = ξ1*AA
	s = sin(AAhat - α)
	t = cos(AAhat - α)
	u = t - cosα
	v = s + sinα*cosc
	q = ((v*t + u*s)*cosα - v)/((v*s + u*t)*sinα)
	Chat = q*A + sqrt(1 - q^2)*normcompCA
	z = 1 - ξ2*(1 - dot(Chat,B))
	return z*B + sqrt(1 - z^2)*normcomp(Chat,B)
end