Skip to content

Instantly share code, notes, and snippets.

@yakir12
Last active August 29, 2015 14:26
Show Gist options
  • Save yakir12/d4e557b35e927b38206c to your computer and use it in GitHub Desktop.
Save yakir12/d4e557b35e927b38206c to your computer and use it in GitHub Desktop.
# "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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment