Fixed major bug related to index ordering

This commit is contained in:
ysyapa 2024-06-05 23:22:40 +00:00
parent 03f9ae6789
commit 25150cbf87
3 changed files with 8 additions and 5 deletions

View File

@ -7,7 +7,7 @@ function get_W_matrix(basis_p, basis_HO, μ1ω1, μ2ω2=μ1ω1; weights=true)
W = zeros(ComplexF64, length(basis_p), length(Es))
Threads.@threads for idx in CartesianIndices(W)
(i1, i2) = Tuple(idx)
((k1, w1), (k2, w2), (j1, j2)) = basis_p[i1]
((j1, j2), (k1, w1), (k2, w2)) = basis_p[i1]
if j1 == l1s[i2] && j2 == l2s[i2]
elem1 = 1/sqrt(sqrt(μ1ω1)) * (-1)^n1s[i2] * ho_basis(j1, n1s[i2], 1/sqrt(μ1ω1) * k1)
elem2 = 1/sqrt(sqrt(μ2ω2)) * (-1)^n2s[i2] * ho_basis(j2, n2s[i2], 1/sqrt(μ2ω2) * k2)

View File

@ -24,10 +24,10 @@ ks, ws = get_mesh(vertices, subdivisions)
jmax = 4
tri((j1, j2)) = triangle_ineq(j1, j2, Λ)
js = collect(Iterators.filter(tri, Iterators.product(0:jmax, 0:jmax)))
js = collect(Iterators.filter(tri, iter_prod(0:jmax, 0:jmax)))
basis = collect(Iterators.product(zip(ks, ws), zip(ks, ws), js))[:]
basis_size = length(ks)^2 * length(js)
basis = iter_prod(js, zip(ks, ws), zip(ks, ws)) # basis = ((j1, j2), (k1, w1), (k2, w2))
basis_size = length(js) * length(ks)^2
@assert length(basis) == basis_size "Something wrong with the basis"
println("Basis size = $basis_size")

View File

@ -27,3 +27,6 @@ function kron_sum(A::AbstractMatrix, B::AbstractMatrix)
@assert size(A, 1) == size(A, 2) && size(B, 1) == size(B, 2) "Matrices should be square"
return kron(A, I(size(B, 1))) + kron(I(size(A, 1)), B)
end
"Flattened vector version of Iterators.product(...) with index hierachy reversed -- leftmost index has the highest hierachy"
iter_prod(args...) = reverse.(collect(Iterators.product(reverse(args)...))[:])