diff --git a/berggren.jl b/berggren.jl index 2b2c7b6..d41ee58 100644 --- a/berggren.jl +++ b/berggren.jl @@ -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) diff --git a/berggren_3body.jl b/berggren_3body.jl index 5c74246..010c5a4 100644 --- a/berggren_3body.jl +++ b/berggren_3body.jl @@ -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") diff --git a/helper.jl b/helper.jl index ae1bb0e..dbd70a5 100644 --- a/helper.jl +++ b/helper.jl @@ -26,4 +26,7 @@ nearest(list::Array, ref) = list[nearestIndex(list, ref)] 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 \ No newline at end of file +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)...))[:]) \ No newline at end of file