From 4c368ac903909aa98f07c9e7eb9131b3c30422d8 Mon Sep 17 00:00:00 2001 From: Nuwan Yapa Date: Wed, 17 Apr 2024 16:25:06 -0400 Subject: [PATCH] Proper mu and omega constants --- XZ_technique.jl | 6 +++--- ho_basis.jl | 8 ++++---- ho_basis_2body.jl | 6 +++--- ho_basis_3body.jl | 10 +++++----- ho_basis_benchmark.jl | 6 +++--- 5 files changed, 18 insertions(+), 18 deletions(-) diff --git a/XZ_technique.jl b/XZ_technique.jl index a4e5158..f4a13f8 100644 --- a/XZ_technique.jl +++ b/XZ_technique.jl @@ -2,7 +2,7 @@ using LinearAlgebra, Plots include("ho_basis.jl") include("p_space.jl") -ω = 0.5 * exp(-1im * 0.47 * pi) +μω_gen = 0.5 * exp(-1im * 0.47 * pi) μ = 0.5 l = 0 V1 = -5 @@ -14,8 +14,8 @@ n_max = 15 ns = collect(0:n_max) ls = fill(l, n_max + 1) -T = sp_T_matrix(ns, ls; ω=ω, μ=μ) -V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); ω=ω) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); ω=ω) +T = sp_T_matrix(ns, ls; μω_gen=μω_gen, μ=μ) +V = V1 .* V_Gaussian.(R1, l, ns, transpose(ns); μω_gen=μω_gen) + V2 .* V_Gaussian.(R2, l, ns, transpose(ns); μω_gen=μω_gen) n_EC = 8 train_cs = (0.7 .+ 0.05 * randn(n_EC)) - 1im * (0.2 .+ 0.05 * randn(n_EC)) diff --git a/ho_basis.jl b/ho_basis.jl index e51ebd0..97a3ca8 100644 --- a/ho_basis.jl +++ b/ho_basis.jl @@ -7,8 +7,8 @@ include("helper.jl") inv_factorial(n) = Iterators.prod(inv.(1:n)) sqrt_factorial(n) = Iterators.prod(sqrt.(n:-1:1)) N_lnk(l, n, k) = 1/sqrt_factorial(n) * binomial(n, k) * sqrt(gamma(n + l + 3/2)) / gamma(k + l + 3/2) -Talmi(l, R, k1, k2; ω=1.0) = (-1)^(k1 + k2) * (1 + 1/(ω * R^2))^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2) -V_Gaussian(R, l, n1, n2; ω=1.0) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2; ω=ω) for (k1, k2) in Iterators.product(0:n1, 0:n2)]) +Talmi(l, R, k1, k2; μω_gen=1.0) = (-1)^(k1 + k2) * (1 + 1/(μω_gen * R^2))^-(3/2 + l + k1 + k2) * gamma(3/2 + l + k1 + k2) +V_Gaussian(R, l, n1, n2; μω_gen=1.0) = (-1)^(n1 + n2) * better_sum([N_lnk(l, n1, k1) * N_lnk(l, n2, k2) * Talmi(l, R, k1, k2; μω_gen=μω_gen) for (k1, k2) in Iterators.product(0:n1, 0:n2)]) function get_sp_basis(E_max) Es = Int[] @@ -54,7 +54,7 @@ function get_2p_basis(E_max) return (Es, n1s, l1s, n2s, l2s) end -function sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), ω=1.0, μ=1.0) +function sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), μω_gen=1.0, μ=1.0) mat = spzeros(length(ns), length(ns)) for idx in CartesianIndices(mat) if !mask[idx]; continue; end @@ -68,7 +68,7 @@ function sp_T_matrix(ns, ls; mask=trues(length(ns),length(ns)), ω=1.0, μ=1.0) end end end - return (ω / μ) .* mat + return (μω_gen / μ) .* mat end function sp_V_matrix(V_l, ns, ls; mask=trues(length(ns),length(ns)), dtype=Float64) diff --git a/ho_basis_2body.jl b/ho_basis_2body.jl index f1ed433..2a2cc04 100644 --- a/ho_basis_2body.jl +++ b/ho_basis_2body.jl @@ -3,7 +3,7 @@ include("ho_basis.jl") include("p_space.jl") E_max = 30 -ω = 0.2 +μω_gen = 0.2 Λ = 0 m = 1.0 @@ -18,10 +18,10 @@ Es, n1s, l1s = get_sp_basis(E_max) println("Basis size = ", length(Es)) println("Constructing KE matrices") -@time "T1" T1 = sp_T_matrix(n1s, l1s; ω=ω, μ=μ1) +@time "T1" T1 = sp_T_matrix(n1s, l1s; μω_gen=μω_gen, μ=μ1) println("Constructing PE matrices") -V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; ω=ω) +V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μω_gen) @time "V1" V1 = sp_V_matrix(V1_elem, n1s, l1s) println("Calculating spectrum") diff --git a/ho_basis_3body.jl b/ho_basis_3body.jl index 8e20791..71a5c20 100644 --- a/ho_basis_3body.jl +++ b/ho_basis_3body.jl @@ -3,7 +3,7 @@ include("ho_basis.jl") include("p_space.jl") E_max = 20 -ω = 0.2 +μω_gen = 0.2 Λ = 0 m = 1.0 @@ -25,12 +25,12 @@ end println("Basis size = ", length(Es)) println("Constructing KE matrices") -@time "T1" T1 = sp_T_matrix(n1s, l1s; mask=mask1, ω=ω, μ=μ1) -@time "T2" T2 = sp_T_matrix(n2s, l2s; mask=mask2, ω=ω, μ=μ2) +@time "T1" T1 = sp_T_matrix(n1s, l1s; mask=mask1, μω_gen=μω_gen, μ=μ1) +@time "T2" T2 = sp_T_matrix(n2s, l2s; mask=mask2, μω_gen=μω_gen, μ=μ2) println("Constructing PE matrices") -V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; ω=ω) -V_relative_elem(l, n1, n2) = Va * V_Gaussian(Ra / c, l, n1, n2; ω=ω) +V1_elem(l, n1, n2) = Va * V_Gaussian(Ra, l, n1, n2; μω_gen=μω_gen) +V_relative_elem(l, n1, n2) = Va * V_Gaussian(Ra / c, l, n1, n2; μω_gen=μω_gen) @time "V1" V1 = sp_V_matrix(V1_elem, n1s, l1s; mask=mask1) @time "V relative" V_relative = sp_V_matrix(V_relative_elem, n1s, l1s; mask=mask1) + sp_V_matrix(V_relative_elem, n2s, l2s; mask=mask2) @time "Moshinsky brackets" U = Moshinsky_transform(Es, n1s, l1s, n2s, l2s, Λ) diff --git a/ho_basis_benchmark.jl b/ho_basis_benchmark.jl index 58c3309..15fad92 100644 --- a/ho_basis_benchmark.jl +++ b/ho_basis_benchmark.jl @@ -2,7 +2,7 @@ using LinearAlgebra, Plots include("ho_basis.jl") include("p_space.jl") -ω = 0.5 * exp(-1im * 0.45 * pi) +μω_gen = 0.5 * exp(-1im * 0.45 * pi) μ = 0.5 l = 0 V1 = -5 @@ -14,9 +14,9 @@ n_max = 15 ns = collect(0:n_max) ls = fill(l, n_max + 1) -T = sp_T_matrix(ns, ls; ω=ω, μ=μ) +T = sp_T_matrix(ns, ls; μω_gen=μω_gen, μ=μ) -V_l(l, n1, n2) = V1 * V_Gaussian(R1, l, n1, n2; ω=ω) + V2 * V_Gaussian(R2, l, n1, n2; ω=ω) +V_l(l, n1, n2) = V1 * V_Gaussian(R1, l, n1, n2; μω_gen=μω_gen) + V2 * V_Gaussian(R2, l, n1, n2; μω_gen=μω_gen) V = sp_V_matrix(V_l, ns, ls; dtype=ComplexF64) cs = range(1.25, 0.25, 10)