diff --git a/Project.toml b/Project.toml index d151dc7..9735fec 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Antique" uuid = "be6e5d0e-34a5-4c8f-af83-e1b5389203d8" authors = ["Shuhei Ohno"] -version = "0.1.2" +version = "0.2.0" [deps] SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" diff --git a/README.md b/README.md index 66b9ea0..fc1ab5e 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ using Pkg; Pkg.add("Antique") ## Usage & Examples -Install Antique.jl for the first use and run `using Antique` before each use. The function `antique(model, parameters...)` returns a module that has `E`, `ψ`, `V` and some other functions. Here are examples in hydrogen-like atom. The analytical notation of energy (eigen value of the Hamiltonian) is written as +Install Antique.jl for the first use and run `using Antique` before each use. The energy `E()`, wavefunction `ψ()`, potential `V()` and some other functions are suppoted. Here are examples in hydrogen-like atom. The analytical notation of energy (eigen value of the Hamiltonian) is written as ```math E_n = -\frac{Z^2}{2n^2} E_\mathrm{h}. @@ -26,8 +26,8 @@ Hydrogen atom has symbol $\mathrm{H}$ and atomic number 1 ($Z=1$). Therefore the ```julia using Antique -H = antique(:HydrogenAtom, Z=1) -H.E(n=1) +H = HydrogenAtom(Z=1) +E(H) # output> -0.5 ``` @@ -35,8 +35,8 @@ Helium cation has symbol $\mathrm{He}^+$ and atomic number 2 ($Z=2$). Therefore ```julia using Antique -He⁺ = antique(:HydrogenAtom, Z=2) -He⁺.E(n=1) +He⁺ = HydrogenAtom(Z=2) +E(He⁺) # output> -2.0 ``` @@ -44,10 +44,10 @@ There are more examples on each model page. ## Supported Models -- [Infinite Potential Well](https://ohno.github.io/Antique.jl/dev/InfinitePotentialWell/) `:InfinitePotentialWell` -- [Harmonic Oscillator](https://ohno.github.io/Antique.jl/dev/HarmonicOscillator/) `:HarmonicOscillator` -- [Morse Potential](https://ohno.github.io/Antique.jl/dev/MorsePotential/) `:MorsePotential` -- [Hydrogen Atom](https://ohno.github.io/Antique.jl/dev/HydrogenAtom/) `:HydrogenAtom` +- [Infinite Potential Well](https://ohno.github.io/Antique.jl/dev/InfinitePotentialWell/) `InfinitePotentialWell` +- [Harmonic Oscillator](https://ohno.github.io/Antique.jl/dev/HarmonicOscillator/) `HarmonicOscillator` +- [Morse Potential](https://ohno.github.io/Antique.jl/dev/MorsePotential/) `MorsePotential` +- [Hydrogen Atom](https://ohno.github.io/Antique.jl/dev/HydrogenAtom/) `HydrogenAtom` ## Future Works diff --git a/developer/dev.ipynb b/developer/dev.ipynb index 0ac587b..e025792 100644 --- a/developer/dev.ipynb +++ b/developer/dev.ipynb @@ -10,21 +10,66 @@ "name": "stderr", "output_type": "stream", "text": [ - "\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m new project at `C:\\Users\\user\\Desktop\\GitHub\\Antique.jl\\developer`\n" + "\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m project at `C:\\Users\\user\\Desktop\\Antique.jl`\n", + "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mPrecompiling Antique [be6e5d0e-34a5-4c8f-af83-e1b5389203d8]\n" ] } ], "source": [ - "include(\"revice.jl\")" + "# run `include(\"./developer/revice.jl\")`\n", + "using Pkg\n", + "Pkg.instantiate()\n", + "# Pkg.add(\"Revise\")\n", + "using Revise\n", + "Pkg.activate(\"../\")\n", + "using Antique" ] }, { "cell_type": "code", - "execution_count": null, - "id": "97a634d3", + "execution_count": 2, + "id": "526ca10f", "metadata": {}, - "outputs": [], - "source": [] + "outputs": [ + { + "data": { + "text/plain": [ + "-0.5" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using Antique\n", + "H = HydrogenAtom(Z=1)\n", + "E(H)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "4e1782eb", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "-2.0" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using Antique\n", + "He⁺ = HydrogenAtom(Z=2)\n", + "E(He⁺)" + ] } ], "metadata": { diff --git a/docs/jmd2md.jl b/docs/jmd2md.jl index f6929f6..0f2df33 100644 --- a/docs/jmd2md.jl +++ b/docs/jmd2md.jl @@ -1,7 +1,7 @@ using Antique using Weave -for file in Antique.models # [:InfinitePotentialWell] +for file in Antique.models # [:InfinitePotentialWell :HarmonicOscillator :MorsePotential :HydrogenAtom] weave("./src/jmd/$file.jmd", doctype="github", out_path="./src/", fig_path="./assets/fig/") text = Antique.load("./src/$file.md") # remove ``` after include(...s) diff --git a/docs/src/HarmonicOscillator.md b/docs/src/HarmonicOscillator.md index 43549f0..cf90e5e 100644 --- a/docs/src/HarmonicOscillator.md +++ b/docs/src/HarmonicOscillator.md @@ -25,7 +25,7 @@ The harmonic oscillator is the most frequently used model in quantum physics. ``` #### Potential -`V(x; k=k, m=m, ω=sqrt(k/m))` +`V(model::HarmonicOscillator; x)` ```math V(x) = \frac{1}{2} k x^2 @@ -34,19 +34,19 @@ The harmonic oscillator is the most frequently used model in quantum physics. ``` #### Eigen Values -`E(n; k, m, ω=sqrt(k/m), ℏ=ℏ)` +`E(model::HarmonicOscillator; n=0)` ```math E_n = \hbar \omega \left( n + \frac{1}{2} \right) ``` #### Eigen Functions -`ψ(n, x; k=k, m=m, ω=sqrt(k/m), ℏ=ℏ)` +`ψ(model::HarmonicOscillator, x; n=0)` ```math \psi_n(x) = A_n H_n(\xi) \exp{\left( -\frac{\xi^2}{2} \right)} ``` #### Hermite Polynomials -`H(x; n=0)` +`H(model::HarmonicOscillator, x; n=0)` Rodrigues' formula & closed-form: ```math @@ -88,11 +88,11 @@ Examples: ## Usage & Examples -[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The function `antique(model, parameters...)` returns a module that has `E()`, `ψ(x)`, `V(x)` and some other functions. In this system, the model name is specified by `:HarmonicOscillator` and several parameters `k`, `m` and `ℏ` are set as optional arguments. +[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The energy `E()`, wavefunction `ψ()`, potential `V()` and some other functions are suppoted. In this system, the model is generated by `HarmonicOscillator` and several parameters `k`, `m` and `ℏ` are set as optional arguments. ```julia using Antique -HO = antique(:HarmonicOscillator, k=1.0, m=1.0, ℏ=1.0) +HO = HarmonicOscillator(k=1.0, m=1.0, ℏ=1.0) ``` @@ -116,10 +116,10 @@ julia> HO.ℏ Eigen values: ```julia -julia> HO.E(n=0) +julia> E(HO, n=0) 0.5 -julia> HO.E(n=1) +julia> E(HO, n=1) 1.5 ``` @@ -129,7 +129,7 @@ Potential energy curve: ```julia using Plots -plot(-5:0.1:5, x -> HO.V(x), lw=2, label="", xlabel="x", ylabel="V(x)") +plot(-5:0.1:5, x -> V(HO, x), lw=2, label="", xlabel="x", ylabel="V(x)") ``` ![](./assets/fig//HarmonicOscillator_4_1.png) @@ -141,11 +141,11 @@ Wave functions: ```julia using Plots plot(xlim=(-5,5), xlabel="x", ylabel="ψ(x)") -plot!(x -> HO.ψ(x, n=0), label="n=0", lw=2) -plot!(x -> HO.ψ(x, n=1), label="n=1", lw=2) -plot!(x -> HO.ψ(x, n=2), label="n=2", lw=2) -plot!(x -> HO.ψ(x, n=3), label="n=3", lw=2) -plot!(x -> HO.ψ(x, n=4), label="n=4", lw=2) +plot!(x -> ψ(HO, x, n=0), label="n=0", lw=2) +plot!(x -> ψ(HO, x, n=1), label="n=1", lw=2) +plot!(x -> ψ(HO, x, n=2), label="n=2", lw=2) +plot!(x -> ψ(HO, x, n=3), label="n=3", lw=2) +plot!(x -> ψ(HO, x, n=4), label="n=4", lw=2) ``` ![](./assets/fig//HarmonicOscillator_5_1.png) @@ -159,13 +159,13 @@ using Plots plot(xlim=(-5.5,5.5), ylim=(-0.2,5.4), xlabel="\$x\$", ylabel="\$V(x),~E_n,~\\psi_n(x)\\times0.5+E_n\$", size=(480,400), dpi=300) for n in 0:4 # energy - hline!([HO.E(n=n)], lc=:black, ls=:dash, label="") - plot!([-sqrt(2*HO.k*HO.E(n=n)),sqrt(2*HO.k*HO.E(n=n))], fill(HO.E(n=n),2), lc=:black, lw=2, label="") + hline!([E(HO, n=n)], lc=:black, ls=:dash, label="") + plot!([-sqrt(2*HO.k*E(HO, n=n)),sqrt(2*HO.k*E(HO, n=n))], fill(E(HO, n=n),2), lc=:black, lw=2, label="") # wave function - plot!(x -> HO.E(n=n) + 0.5*HO.ψ(x,n=n), lc=n+1, lw=2, label="") + plot!(x -> E(HO, n=n) + 0.5*ψ(HO, x,n=n), lc=n+1, lw=2, label="") end # potential -plot!(x -> HO.V(x), lc=:black, lw=2, label="") +plot!(x -> V(HO, x), lc=:black, lw=2, label="") ``` ![](./assets/fig//HarmonicOscillator_6_1.png) diff --git a/docs/src/HydrogenAtom.md b/docs/src/HydrogenAtom.md index fe64521..cc85334 100644 --- a/docs/src/HydrogenAtom.md +++ b/docs/src/HydrogenAtom.md @@ -22,7 +22,7 @@ The hydrogen atom is the simplest 2-body Coulomb system. Where, $\mu=\left(\frac{1}{m_\mathrm{e}}+\frac{1}{m_\mathrm{p}}\right)^{-1}$ is the reduced mass of electron $\mathrm{e}$ and proton $\mathrm{p}$. $\mu = m_\mathrm{e}$ holds in the limit $m_\mathrm{p}\rightarrow\infty$. #### Potential -`V(r; Z=Z, a₀=a₀)` +`V(model::HydrogenAtom, r)` ```math \begin{aligned} V(r) @@ -33,27 +33,27 @@ Where, $\mu=\left(\frac{1}{m_\mathrm{e}}+\frac{1}{m_\mathrm{p}}\right)^{-1}$ is ``` #### Eigen Values -`E(; n=1, Z=Z, Eₕ=Eₕ)` +`E(model::HydrogenAtom; n=1)` ```math E_n = -\frac{m_\mathrm{e} e^4 Z^2}{2n^2(4\pi\varepsilon_0)^2\hbar^2} = -\frac{Z^2}{2n^2} E_\mathrm{h} ``` Where, ``E_\mathrm{h}`` is the Hartree energy, one of atomic unit. About atomic units, see section 3.9.2 of the [IUPAC GreenBook](https://iupac.org/what-we-do/books/greenbook/). In other units, ``E_\mathrm{h} = 27.211~386~245~988(53)~\mathrm{eV}`` from [here](https://physics.nist.gov/cgi-bin/cuu/Value?hrev). #### Eigen Functions -`ψ(r, θ, φ; n=1, l=0, m=0, Z=Z, a₀=a₀)` +`ψ(model::HydrogenAtom, r, θ, φ; n=1, l=0, m=0)` ```math \psi_{nlm}(\pmb{r}) = R_{nl}(r) Y_{lm}(\theta,\varphi) ``` #### Radial Functions -`R(r; n=1, l=0, Z=Z, a₀=a₀)` +`R(model::HydrogenAtom, r; n=1, l=0)` ```math R_{nl}(r) = -\sqrt{\frac{(n-l-1)!}{2n(n+l)!} \left(\frac{2Z}{n a_0}\right)^3} \left(\frac{2Zr}{n a_0}\right)^l \exp \left(-\frac{Zr}{n a_0}\right) L_{n+l}^{2l+1} \left(\frac{2Zr}{n a_0}\right) ``` Where, Laguerre polynomials are defined as ``L_n(x) = \frac{1}{n!} \mathrm{e}^x \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right)``, and associated Laguerre polynomials are defined as ``L_n^{k}(x) = \frac{\mathrm{d}^k}{\mathrm{d}x^k} L_n(x)``. Note that, replace ``2n(n+l)!`` with ``2n[(n+l)!]^3`` if Laguerre polynomials are defined as ``L_n(x) = \mathrm{e}^x \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right)``. #### Associated Laguerre Polynomials -`L(x; n=0, k=0)` +`L(model::HydrogenAtom, x; n=0, k=0)` Rodrigues' formula & closed-form: ```math @@ -91,7 +91,7 @@ Examples: ``` #### Spherical Harmonics -`Y(θ, φ; l=0, m=0)` +`Y(model::HydrogenAtom, θ, φ; l=0, m=0)` ```math Y_{lm}(\theta,\varphi) = (-1)^{\frac{|m|+m}{2}} \sqrt{\frac{2l+1}{4\pi} \frac{(l-|m|)!}{(l+|m|)!}} P_l^{|m|} (\cos\theta) \mathrm{e}^{im\varphi} ``` @@ -101,7 +101,7 @@ i^{|m|+m} \sqrt{\frac{(l-|m|)!}{(l+|m|)!}} P_l^{|m|} = (-1)^{\frac{|m|+m}{2}} \s ``` #### Associated Legendre Polynomials -`P(x; n=0, m=0)` +`P(model::HydrogenAtom, x; n=0, m=0)` Rodrigues' formula & closed-form: ```math @@ -151,11 +151,11 @@ Examples: ## Usage & Examples -[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The function `antique(model, parameters...)` returns a module that has `E()`, `ψ(r)`, `V(r)` and some other functions. In this system, the model name is specified by `:HydrogenAtom` and several parameters `Z`, `Eₕ`, `mₑ`, `a₀` and `ℏ` are set as optional arguments. +[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The energy `E()`, wavefunction `ψ()`, potential `V()` and some other functions are suppoted. In this system, the model is generated by `HydrogenAtom` and several parameters `Z`, `Eₕ`, `mₑ`, `a₀` and `ℏ` are set as optional arguments. ```julia using Antique -H = antique(:HydrogenAtom, Z=1, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0) +H = HydrogenAtom(Z=1, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0) ``` @@ -185,10 +185,10 @@ julia> H.ℏ Eigen values: ```julia -julia> H.E(n=1) +julia> E(H, n=1) -0.5 -julia> H.E(n=2) +julia> E(H, n=2) -0.125 ``` @@ -198,8 +198,8 @@ Wave length ($n=2\rightarrow1$, the first line of the Lyman series): ```julia Eₕ2nm⁻¹ = 2.1947463136320e-2 # https://physics.nist.gov/cgi-bin/cuu/CCValue?hrminv -println("ΔE = ", H.E(n=2) - H.E(n=1), " Eₕ") -println("λ = ", ((H.E(n=2)-H.E(n=1))*Eₕ2nm⁻¹)^-1, " nm") +println("ΔE = ", E(H,n=2) - E(H,n=1), " Eₕ") +println("λ = ", ((E(H,n=2)-E(H,n=1))*Eₕ2nm⁻¹)^-1, " nm") ``` ``` @@ -226,7 +226,7 @@ ge = 2.00231930436256 # https://physics.nist.gov/cgi-bin/cuu/Value?gem gp = 5.5856946893 # https://physics.nist.gov/cgi-bin/cuu/Value?gp # calculation: https://doi.org/10.1119/1.12733 -δ = abs(H.ψ(0,0,0))^2 +δ = abs(ψ(H,0,0,0))^2 ΔE = 2 / 3 * µ0 * µN * µB * gp * ge * δ * a0^(-3) println("1/π = ", 1/π) println("<δ(r)> = ", δ, " a₀⁻³") @@ -254,7 +254,7 @@ Potential energy curve: ```julia using Plots plot(xlims=(0.0,15.0), ylims=(-0.6,0.05), xlabel="\$r~/~a_0\$", ylabel="\$V(r)/E_\\mathrm{h},~E_n/E_\\mathrm{h}\$", legend=:bottomright, size=(480,400)) -plot!(0.1:0.01:15, r -> H.V(r), lc=:black, lw=2, label="") # potential +plot!(0.1:0.01:15, r -> V(H,r), lc=:black, lw=2, label="") # potential ``` ![](./assets/fig//HydrogenAtom_6_1.png) @@ -267,9 +267,9 @@ Potential energy curve, Energy levels: using Plots plot(xlims=(0.0,15.0), ylims=(-0.6,0.05), xlabel="\$r~/~a_0\$", ylabel="\$V(r)/E_\\mathrm{h}\$", legend=:bottomright, size=(480,400)) for n in 0:10 - plot!(0.0:0.01:15, r -> H.E(n=n) > H.V(r) ? H.E(n=n) : NaN, lc=n, lw=1, label="") # energy level + plot!(0.0:0.01:15, r -> E(H,n=n) > V(H,r) ? E(H,n=n) : NaN, lc=n, lw=1, label="") # energy level end -plot!(0.1:0.01:15, r -> H.V(r), lc=:black, lw=2, label="") # potential +plot!(0.1:0.01:15, r -> V(H,r), lc=:black, lw=2, label="") # potential ``` ![](./assets/fig//HydrogenAtom_7_1.png) @@ -283,7 +283,7 @@ using Plots plot(xlabel="\$r~/~a_0\$", ylabel="\$r^2|R_{nl}(r)|^2~/~a_0^{-1}\$", ylims=(-0.01,0.55), xticks=0:1:20, size=(480,400), dpi=300) for n in 1:3 for l in 0:n-1 - plot!(0:0.01:20, r->r^2*H.R(r,n=n,l=l)^2, lc=n, lw=2, ls=[:solid,:dash,:dot,:dashdot,:dashdotdot][l+1], label="\$n = $n, l=$l\$") + plot!(0:0.01:20, r->r^2*R(H,r,n=n,l=l)^2, lc=n, lw=2, ls=[:solid,:dash,:dot,:dashdot,:dashdotdot][l+1], label="\$n = $n, l=$l\$") end end plot!() diff --git a/docs/src/InfinitePotentialWell.md b/docs/src/InfinitePotentialWell.md index f0e7848..72f230b 100644 --- a/docs/src/InfinitePotentialWell.md +++ b/docs/src/InfinitePotentialWell.md @@ -21,7 +21,7 @@ The infinite potential well (particle in a box) is the simplest model for quantu ``` #### Potential -`V(x; L=L)` +`V(model::InfinitePotentialWell; x)` ```math V(x) = \left\{ @@ -33,13 +33,13 @@ The infinite potential well (particle in a box) is the simplest model for quantu ``` #### Eigen Values -`E(; n=0, L=L, m=m, ℏ=ℏ)` +`E(model::InfinitePotentialWell; n=0)` ```math E_n = \frac{\hbar^2 n^2 \pi^2}{2 m L^2} ``` #### Eigen Functions -`ψ(x; n=0, L=L)` +`ψ(model::InfinitePotentialWell, x; n=0)` ```math \psi_n(x) = \sqrt{\frac{2}{L}} \sin \frac{n\pi x}{L} ``` @@ -50,11 +50,11 @@ The infinite potential well (particle in a box) is the simplest model for quantu ## Usage & Examples -[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The function `antique(model, parameters...)` returns a module that has `E()`, `ψ(x)`, `V(x)` and some other functions. In this system, the model name is specified by `:InfinitePotentialWell` and several parameters `L`, `m` and `ℏ` are set as optional arguments. +[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The energy `E()`, wavefunction `ψ()`, potential `V()` and some other functions are suppoted. In this system, the model is generated by `InfinitePotentialWell` and several parameters `L`, `m` and `ℏ` are set as optional arguments. ```julia using Antique -IPW = antique(:InfinitePotentialWell, L=1.0, m=1.0, ℏ=1.0) +IPW = InfinitePotentialWell(L=1.0, m=1.0, ℏ=1.0) ``` @@ -78,10 +78,10 @@ julia> IPW.ℏ Eigen values: ```julia -julia> IPW.E(n=1) +julia> E(IPW, n=1) 4.934802200544679 -julia> IPW.E(n=2) +julia> E(IPW, n=2) 19.739208802178716 ``` @@ -92,11 +92,11 @@ Wave functions: ```julia using Plots plot(xlim=(0,1), xlabel="x", ylabel="ψ(x)") -plot!(x -> IPW.ψ(x, n=1), label="n=1", lw=2) -plot!(x -> IPW.ψ(x, n=2), label="n=2", lw=2) -plot!(x -> IPW.ψ(x, n=3), label="n=3", lw=2) -plot!(x -> IPW.ψ(x, n=4), label="n=4", lw=2) -plot!(x -> IPW.ψ(x, n=5), label="n=5", lw=2) +plot!(x -> ψ(IPW, x, n=1), label="n=1", lw=2) +plot!(x -> ψ(IPW, x, n=2), label="n=2", lw=2) +plot!(x -> ψ(IPW, x, n=3), label="n=3", lw=2) +plot!(x -> ψ(IPW, x, n=4), label="n=4", lw=2) +plot!(x -> ψ(IPW, x, n=5), label="n=5", lw=2) ``` ![](./assets/fig//InfinitePotentialWell_4_1.png) @@ -111,9 +111,9 @@ using Plots plot(xlim=(-0.5,1.5), ylim=(-5,140), xlabel="\$x\$", ylabel="\$V(x),~E_n,~\\psi_n(x)\\times5+E_n\$", size=(480,400), dpi=300) for n in 1:5 # energy - plot!([0,L], fill(IPW.E(n=n,L=L),2), lc=:black, lw=2, label="") + plot!([0,L], fill(E(IPW,n=n),2), lc=:black, lw=2, label="") # wave function - plot!(0:0.01:L, x->IPW.E(n=n,L=L)+5*IPW.ψ(x,n=n,L=L), lc=n, lw=2, label="") + plot!(0:0.01:L, x->E(IPW,n=n) + 5*ψ(IPW,x,n=n), lc=n, lw=2, label="") end # potential plot!([0,0,L,L], [140,0,0,140], lc=:black, lw=2, label="") diff --git a/docs/src/MorsePotential.md b/docs/src/MorsePotential.md index 16ad77c..f748075 100644 --- a/docs/src/MorsePotential.md +++ b/docs/src/MorsePotential.md @@ -22,25 +22,25 @@ The Morse potential is a model for inter-nuclear anharmonic vibration in a diato ``` #### Potential -`V(r; rₑ=rₑ, Dₑ=Dₑ, k=k, a=sqrt(k/(2*Dₑ)))` +`V(model::MorsePotential, r)` ```math V(r) = D_\mathrm{e} \left( \mathrm{e}^{-2a(r-r_e)} - 2\mathrm{e}^{-a(r-r_e)} \right) ``` #### Eigen Values -`E(n; rₑ=rₑ, Dₑ=Dₑ, k=k, a=sqrt(k/(2*Dₑ)), µ=µ, ω=sqrt(k/µ), χ=ℏ*ω/(4*Dₑ), ℏ=ℏ)` +`E(model::MorsePotential; n=0)` ```math E_n = - D_\mathrm{e} + \hbar \omega \left( n + \frac{1}{2} \right) - \chi \hbar \omega \left( n + \frac{1}{2} \right)^2 ``` #### Eigen Functions -`ψ(n, r; rₑ=rₑ, Dₑ=Dₑ, k=k, a=sqrt(k/(2*Dₑ)), µ=µ, ω=sqrt(k/µ), χ=ℏ*ω/(4*Dₑ), ℏ=ℏ)` +`ψ(model::MorsePotential, r; n=0)` ```math \psi_n(r) = N_n z^{\lambda-n-1/2} \mathrm{e}^{-z/2} L_n^{(2\lambda-2n-1)}(\xi) ``` #### Generalized Laguerre Polynomials -`L(x; n=0, α=0)` +`L(model::MorsePotential, x; n=0, α=0)` Rodrigues' formula & closed-form: ```math @@ -81,11 +81,7 @@ Examples: ## Usage & Examples -[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The function `antique(model, parameters...)` returns a module that has `E()`, `ψ(r)`, `V(r)` and some other functions. In this system, the model name is specified by `:MorsePotential` and several parameters `rₑ`, `Dₑ`, `k`, `µ` and `ℏ` are set as optional arguments. - -!!! warning - - Run `using Pkg; Pkg.add("SpecialFunctions")` if the following returns an error. +[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The energy `E()`, wavefunction `ψ()`, potential `V()` and some other functions are suppoted. In this system, the model is generated by `MorsePotential` and several parameters `rₑ`, `Dₑ`, `k`, `µ` and `ℏ` are set as optional arguments. ```julia # Parameters for H₂⁺ @@ -93,14 +89,13 @@ Examples: # https://doi.org/10.5281/zenodo.5047817 # https://physics.nist.gov/cgi-bin/cuu/Value?mpsme rₑ = 1.997193319969992120068298141276 -Vₑ = -0.602634619106539878727562156289 -Dₑ = - 0.5 - Vₑ -k = 2*((-1.1026342144949464615+1/2.00) - Vₑ) / (2.00 - rₑ)^2 +Dₑ = - 0.5 - (-0.602634619106539878727562156289) +k = 2*((-1.1026342144949464615+1/2.00) - (-0.602634619106539878727562156289)) / (2.00 - rₑ)^2 µ = 1/(1/1836.15267343 + 1/1836.15267343) ℏ = 1.0 using Antique -MP = antique(:MorsePotential, rₑ=rₑ, Vₑ=Vₑ, Dₑ=Dₑ, k=k, µ =µ, ℏ =ℏ) +MP = MorsePotential(rₑ=rₑ, Dₑ=Dₑ, k=k, µ=µ, ℏ=ℏ) ``` @@ -130,10 +125,10 @@ julia> MP.ℏ Eigen values: ```julia -julia> MP.E(n=0) +julia> E(MP, n=0) -0.09741377794418261 -julia> MP.E(n=1) +julia> E(MP, n=1) -0.08738092406760907 ``` @@ -143,7 +138,7 @@ Potential energy curve: ```julia using Plots -plot(0.1:0.01:15, r -> MP.V(r), lw=2, label="", xlims=(0.1,9.1), ylims=(-0.11,0.01), xlabel="r", ylabel="V(r)") +plot(0.1:0.01:15, r -> V(MP, r), lw=2, label="", xlims=(0.1,9.1), ylims=(-0.11,0.01), xlabel="r", ylabel="V(r)") ``` ![](./assets/fig//MorsePotential_4_1.png) @@ -155,12 +150,12 @@ Wave functions: ```julia using Plots plot(xlim=(0,5), xlabel="x", ylabel="ψ(x)") -plot!(x -> MP.ψ(x, n=0), label="n=0", lw=2) -plot!(x -> MP.ψ(x, n=1), label="n=1", lw=2) -plot!(x -> MP.ψ(x, n=2), label="n=2", lw=2) -plot!(x -> MP.ψ(x, n=3), label="n=3", lw=2) -plot!(x -> MP.ψ(x, n=4), label="n=4", lw=2) -plot!(x -> MP.ψ(x, n=5), label="n=5", lw=2) +plot!(x -> ψ(MP, x, n=0), label="n=0", lw=2) +plot!(x -> ψ(MP, x, n=1), label="n=1", lw=2) +plot!(x -> ψ(MP, x, n=2), label="n=2", lw=2) +plot!(x -> ψ(MP, x, n=3), label="n=3", lw=2) +plot!(x -> ψ(MP, x, n=4), label="n=4", lw=2) +plot!(x -> ψ(MP, x, n=5), label="n=5", lw=2) ``` ![](./assets/fig//MorsePotential_5_1.png) @@ -170,20 +165,20 @@ plot!(x -> MP.ψ(x, n=5), label="n=5", lw=2) Potential energy curve, Energy levels, Comparison with harmonic oscillator: ```julia -MP = antique(:MorsePotential) -HO = antique(:HarmonicOscillator, k=MP.k, m=MP.μ) +MP = MorsePotential() +HO = HarmonicOscillator(k=MP.k, m=MP.μ) using Plots plot(xlims=(0.1,9.1), ylims=(-0.11,0.01), xlabel="\$r\$", ylabel="\$V(r), E_n\$", legend=:bottomright, size=(480,400), dpi=300) -for n in 0:MP.nₘₐₓ() +for n in 0:nₘₐₓ(MP) # energy - EM = MP.E(n=n) - EH = HO.E(n=n) - MP.Dₑ - plot!(0.1:0.01:15, r -> EH > HO.V(r-MP.rₑ) - MP.Dₑ ? EH : NaN, lc="#BC1C5F", lw=1, label="") - plot!(0.1:0.01:15, r -> EM > MP.V(r) ? EM : NaN, lc="#578FC7", lw=1, label="") + EM = E(MP, n=n) + EH = E(HO, n=n) - MP.Dₑ + plot!(0.1:0.01:15, r -> EH > V(HO, r-MP.rₑ) - MP.Dₑ ? EH : NaN, lc="#BC1C5F", lw=1, label="") + plot!(0.1:0.01:15, r -> EM > V(MP, r) ? EM : NaN, lc="#578FC7", lw=1, label="") end # potential -plot!(0.1:0.01:15, r -> HO.V(r-MP.rₑ) - MP.Dₑ, lc="#BC1C5F", lw=2, label="Harmonic Oscillator") -plot!(0.1:0.01:15, r -> MP.V(r), lc="#578FC7", lw=2, label="Morse Potential") +plot!(0.1:0.01:15, r -> V(HO, r-MP.rₑ) - MP.Dₑ, lc="#BC1C5F", lw=2, label="Harmonic Oscillator") +plot!(0.1:0.01:15, r -> V(MP, r), lc="#578FC7", lw=2, label="Morse Potential") ``` ![](./assets/fig//MorsePotential_6_1.png) @@ -773,55 +768,55 @@ Lₙ⁽ᵅ⁾(x) = x⁻ᵅeˣ/n! dⁿ/dxⁿ xⁿ⁺ᵅe⁻ˣ | 15 15 1.00 | 0 | 0 | 1.000000000000 | 1.000000000000 ✔ 1.00 | 0 | 1 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 0 | 2 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 0 | 3 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 0 | 4 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 0 | 5 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 0 | 3 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 0 | 4 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 0 | 5 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 0 | 6 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 0 | 7 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 0 | 7 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 0 | 8 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 0 | 9 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 0 | 9 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 1 | 0 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 1 | 1 | 2.000000000000 | 2.000000000000 ✔ 1.00 | 1 | 2 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 1 | 3 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 1 | 4 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 1 | 5 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 1 | 3 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 1 | 4 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 1 | 5 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 1 | 6 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 1 | 7 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 1 | 8 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 1 | 9 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 1 | 7 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 1 | 8 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 1 | 9 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 2 | 0 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 2 | 1 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 2 | 2 | 3.000000000000 | 3.000000000000 ✔ -1.00 | 2 | 3 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 2 | 3 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 2 | 4 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 2 | 5 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 2 | 6 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 2 | 7 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 2 | 6 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 2 | 7 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 2 | 8 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 2 | 9 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 3 | 0 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 3 | 1 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 3 | 2 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 3 | 0 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 3 | 1 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 3 | 2 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 3 | 3 | 4.000000000000 | 4.000000000000 ✔ -1.00 | 3 | 4 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 3 | 4 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 3 | 5 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 3 | 6 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 3 | 7 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 3 | 8 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 3 | 9 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 4 | 0 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 4 | 1 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 3 | 9 | 0.000000000000 | -0.000000000001 ✔ +1.00 | 4 | 0 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 4 | 1 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 4 | 2 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 4 | 3 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 4 | 3 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 4 | 4 | 5.000000000000 | 4.999999999999 ✔ 1.00 | 4 | 5 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 4 | 6 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 4 | 7 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 4 | 8 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 4 | 8 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 4 | 9 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 5 | 0 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 5 | 1 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 5 | 0 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 5 | 1 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 5 | 2 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 5 | 3 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 5 | 4 | 0.000000000000 | -0.000000000000 ✔ @@ -829,10 +824,10 @@ Lₙ⁽ᵅ⁾(x) = x⁻ᵅeˣ/n! dⁿ/dxⁿ xⁿ⁺ᵅe⁻ˣ | 15 15 1.00 | 5 | 6 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 5 | 7 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 5 | 8 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 5 | 9 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 5 | 9 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 6 | 0 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 6 | 1 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 6 | 2 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 6 | 2 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 6 | 3 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 6 | 4 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 6 | 5 | 0.000000000000 | -0.000000000000 ✔ @@ -840,9 +835,9 @@ Lₙ⁽ᵅ⁾(x) = x⁻ᵅeˣ/n! dⁿ/dxⁿ xⁿ⁺ᵅe⁻ˣ | 15 15 1.00 | 6 | 7 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 6 | 8 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 6 | 9 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 7 | 0 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 7 | 1 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 7 | 2 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 7 | 0 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 7 | 1 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 7 | 2 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 7 | 3 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 7 | 4 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 7 | 5 | 0.000000000000 | -0.000000000000 ✔ @@ -851,25 +846,25 @@ Lₙ⁽ᵅ⁾(x) = x⁻ᵅeˣ/n! dⁿ/dxⁿ xⁿ⁺ᵅe⁻ˣ | 15 15 1.00 | 7 | 8 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 7 | 9 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 8 | 0 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 8 | 1 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 8 | 1 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 8 | 2 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 8 | 3 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 8 | 4 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 8 | 4 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 8 | 5 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 8 | 6 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 8 | 7 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 8 | 8 | 9.000000000000 | 9.000000000001 ✔ +1.00 | 8 | 8 | 9.000000000000 | 9.000000000000 ✔ 1.00 | 8 | 9 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 9 | 0 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 9 | 1 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 9 | 0 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 9 | 1 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 9 | 2 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 9 | 3 | 0.000000000000 | -0.000000000000 ✔ +1.00 | 9 | 3 | 0.000000000000 | -0.000000000001 ✔ 1.00 | 9 | 4 | 0.000000000000 | 0.000000000000 ✔ -1.00 | 9 | 5 | 0.000000000000 | 0.000000000000 ✔ +1.00 | 9 | 5 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 9 | 6 | 0.000000000000 | 0.000000000000 ✔ 1.00 | 9 | 7 | 0.000000000000 | -0.000000000000 ✔ 1.00 | 9 | 8 | 0.000000000000 | -0.000000000000 ✔ -1.00 | 9 | 9 | 10.000000000000 | 10.000000000001 ✔ +1.00 | 9 | 9 | 10.000000000000 | 10.000000000002 ✔ Test Summary: | Pass Total ∫Lᵢ⁽ᵅ⁾(x)Lⱼ⁽ᵅ⁾(x)xᵅexp(-x)dx = Γ(i+α+1)/i! δᵢⱼ | 500 500 ``` diff --git a/docs/src/index.md b/docs/src/index.md index 614ec63..4a49129 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -16,7 +16,7 @@ using Pkg; Pkg.add("Antique") ## Usage & Examples -[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The function `antique(model, parameters...)` returns a module that has `E`, `ψ`, `V` and some other functions. Here are examples in hydrogen-like atom. The analytical notation of energy (eigen value of the Hamiltonian) is written as +[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The energy `E()`, wavefunction `ψ()`, potential `V()` and some other functions are suppoted. Here are examples in hydrogen-like atom. The analytical notation of energy (eigen value of the Hamiltonian) is written as ```math E_n = -\frac{Z^2}{2n^2} E_\mathrm{h}. @@ -26,8 +26,8 @@ Hydrogen atom has symbol $\mathrm{H}$ and atomic number 1 ($Z=1$). Therefore the ```julia using Antique -H = antique(:HydrogenAtom, Z=1) -H.E(n=1) +H = HydrogenAtom(Z=1) +E(H) # output> -0.5 ``` @@ -35,8 +35,8 @@ Helium cation has symbol $\mathrm{He}^+$ and atomic number 2 ($Z=2$). Therefore ```julia using Antique -He⁺ = antique(:HydrogenAtom, Z=2) -He⁺.E(n=1) +He⁺ = HydrogenAtom(Z=2) +E(He⁺) # output> -2.0 ``` @@ -50,25 +50,25 @@ There are more examples on each model page. InfinitePotentialWell - :InfinitePotentialWell + InfinitePotentialWell
HarmonicOscillator - :HarmonicOscillator + HarmonicOscillator
MorsePotential - :MorsePotential + MorsePotential
HydrogenAtom - :HydrogenAtom + HydrogenAtom
``` diff --git a/docs/src/jmd/HarmonicOscillator.jmd b/docs/src/jmd/HarmonicOscillator.jmd index 65da077..9faa8c6 100644 --- a/docs/src/jmd/HarmonicOscillator.jmd +++ b/docs/src/jmd/HarmonicOscillator.jmd @@ -25,7 +25,7 @@ The harmonic oscillator is the most frequently used model in quantum physics. ``` #### Potential -`V(x; k=k, m=m, ω=sqrt(k/m))` +`V(model::HarmonicOscillator; x)` ```math V(x) = \frac{1}{2} k x^2 @@ -34,19 +34,19 @@ The harmonic oscillator is the most frequently used model in quantum physics. ``` #### Eigen Values -`E(n; k, m, ω=sqrt(k/m), ℏ=ℏ)` +`E(model::HarmonicOscillator; n=0)` ```math E_n = \hbar \omega \left( n + \frac{1}{2} \right) ``` #### Eigen Functions -`ψ(n, x; k=k, m=m, ω=sqrt(k/m), ℏ=ℏ)` +`ψ(model::HarmonicOscillator, x; n=0)` ```math \psi_n(x) = A_n H_n(\xi) \exp{\left( -\frac{\xi^2}{2} \right)} ``` #### Hermite Polynomials -`H(x; n=0)` +`H(model::HarmonicOscillator, x; n=0)` Rodrigues' formula & closed-form: ```math @@ -88,11 +88,11 @@ Examples: ## Usage & Examples -[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The function `antique(model, parameters...)` returns a module that has `E()`, `ψ(x)`, `V(x)` and some other functions. In this system, the model name is specified by `:HarmonicOscillator` and several parameters `k`, `m` and `ℏ` are set as optional arguments. +[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The energy `E()`, wavefunction `ψ()`, potential `V()` and some other functions are suppoted. In this system, the model is generated by `HarmonicOscillator` and several parameters `k`, `m` and `ℏ` are set as optional arguments. ```julia; cache = :all; results = "hidden" using Antique -HO = antique(:HarmonicOscillator, k=1.0, m=1.0, ℏ=1.0) +HO = HarmonicOscillator(k=1.0, m=1.0, ℏ=1.0) ``` Parameters: @@ -106,15 +106,15 @@ HO.ℏ Eigen values: ```julia; term = true -HO.E(n=0) -HO.E(n=1) +E(HO, n=0) +E(HO, n=1) ``` Potential energy curve: ```julia using Plots -plot(-5:0.1:5, x -> HO.V(x), lw=2, label="", xlabel="x", ylabel="V(x)") +plot(-5:0.1:5, x -> V(HO, x), lw=2, label="", xlabel="x", ylabel="V(x)") ``` Wave functions: @@ -122,11 +122,11 @@ Wave functions: ```julia using Plots plot(xlim=(-5,5), xlabel="x", ylabel="ψ(x)") -plot!(x -> HO.ψ(x, n=0), label="n=0", lw=2) -plot!(x -> HO.ψ(x, n=1), label="n=1", lw=2) -plot!(x -> HO.ψ(x, n=2), label="n=2", lw=2) -plot!(x -> HO.ψ(x, n=3), label="n=3", lw=2) -plot!(x -> HO.ψ(x, n=4), label="n=4", lw=2) +plot!(x -> ψ(HO, x, n=0), label="n=0", lw=2) +plot!(x -> ψ(HO, x, n=1), label="n=1", lw=2) +plot!(x -> ψ(HO, x, n=2), label="n=2", lw=2) +plot!(x -> ψ(HO, x, n=3), label="n=3", lw=2) +plot!(x -> ψ(HO, x, n=4), label="n=4", lw=2) ``` Potential energy curve, Energy levels, Wave functions: @@ -136,13 +136,13 @@ using Plots plot(xlim=(-5.5,5.5), ylim=(-0.2,5.4), xlabel="\$x\$", ylabel="\$V(x),~E_n,~\\psi_n(x)\\times0.5+E_n\$", size=(480,400), dpi=300) for n in 0:4 # energy - hline!([HO.E(n=n)], lc=:black, ls=:dash, label="") - plot!([-sqrt(2*HO.k*HO.E(n=n)),sqrt(2*HO.k*HO.E(n=n))], fill(HO.E(n=n),2), lc=:black, lw=2, label="") + hline!([E(HO, n=n)], lc=:black, ls=:dash, label="") + plot!([-sqrt(2*HO.k*E(HO, n=n)),sqrt(2*HO.k*E(HO, n=n))], fill(E(HO, n=n),2), lc=:black, lw=2, label="") # wave function - plot!(x -> HO.E(n=n) + 0.5*HO.ψ(x,n=n), lc=n+1, lw=2, label="") + plot!(x -> E(HO, n=n) + 0.5*ψ(HO, x,n=n), lc=n+1, lw=2, label="") end # potential -plot!(x -> HO.V(x), lc=:black, lw=2, label="") +plot!(x -> V(HO, x), lc=:black, lw=2, label="") ``` ## Testing diff --git a/docs/src/jmd/HydrogenAtom.jmd b/docs/src/jmd/HydrogenAtom.jmd index d78510d..af691d0 100644 --- a/docs/src/jmd/HydrogenAtom.jmd +++ b/docs/src/jmd/HydrogenAtom.jmd @@ -22,7 +22,7 @@ The hydrogen atom is the simplest 2-body Coulomb system. Where, $\mu=\left(\frac{1}{m_\mathrm{e}}+\frac{1}{m_\mathrm{p}}\right)^{-1}$ is the reduced mass of electron $\mathrm{e}$ and proton $\mathrm{p}$. $\mu = m_\mathrm{e}$ holds in the limit $m_\mathrm{p}\rightarrow\infty$. #### Potential -`V(r; Z=Z, a₀=a₀)` +`V(model::HydrogenAtom, r)` ```math \begin{aligned} V(r) @@ -33,27 +33,27 @@ Where, $\mu=\left(\frac{1}{m_\mathrm{e}}+\frac{1}{m_\mathrm{p}}\right)^{-1}$ is ``` #### Eigen Values -`E(; n=1, Z=Z, Eₕ=Eₕ)` +`E(model::HydrogenAtom; n=1)` ```math E_n = -\frac{m_\mathrm{e} e^4 Z^2}{2n^2(4\pi\varepsilon_0)^2\hbar^2} = -\frac{Z^2}{2n^2} E_\mathrm{h} ``` Where, ``E_\mathrm{h}`` is the Hartree energy, one of atomic unit. About atomic units, see section 3.9.2 of the [IUPAC GreenBook](https://iupac.org/what-we-do/books/greenbook/). In other units, ``E_\mathrm{h} = 27.211~386~245~988(53)~\mathrm{eV}`` from [here](https://physics.nist.gov/cgi-bin/cuu/Value?hrev). #### Eigen Functions -`ψ(r, θ, φ; n=1, l=0, m=0, Z=Z, a₀=a₀)` +`ψ(model::HydrogenAtom, r, θ, φ; n=1, l=0, m=0)` ```math \psi_{nlm}(\pmb{r}) = R_{nl}(r) Y_{lm}(\theta,\varphi) ``` #### Radial Functions -`R(r; n=1, l=0, Z=Z, a₀=a₀)` +`R(model::HydrogenAtom, r; n=1, l=0)` ```math R_{nl}(r) = -\sqrt{\frac{(n-l-1)!}{2n(n+l)!} \left(\frac{2Z}{n a_0}\right)^3} \left(\frac{2Zr}{n a_0}\right)^l \exp \left(-\frac{Zr}{n a_0}\right) L_{n+l}^{2l+1} \left(\frac{2Zr}{n a_0}\right) ``` Where, Laguerre polynomials are defined as ``L_n(x) = \frac{1}{n!} \mathrm{e}^x \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right)``, and associated Laguerre polynomials are defined as ``L_n^{k}(x) = \frac{\mathrm{d}^k}{\mathrm{d}x^k} L_n(x)``. Note that, replace ``2n(n+l)!`` with ``2n[(n+l)!]^3`` if Laguerre polynomials are defined as ``L_n(x) = \mathrm{e}^x \frac{\mathrm{d}^n}{\mathrm{d}x ^n} \left( \mathrm{e}^{-x} x^n \right)``. #### Associated Laguerre Polynomials -`L(x; n=0, k=0)` +`L(model::HydrogenAtom, x; n=0, k=0)` Rodrigues' formula & closed-form: ```math @@ -91,7 +91,7 @@ Examples: ``` #### Spherical Harmonics -`Y(θ, φ; l=0, m=0)` +`Y(model::HydrogenAtom, θ, φ; l=0, m=0)` ```math Y_{lm}(\theta,\varphi) = (-1)^{\frac{|m|+m}{2}} \sqrt{\frac{2l+1}{4\pi} \frac{(l-|m|)!}{(l+|m|)!}} P_l^{|m|} (\cos\theta) \mathrm{e}^{im\varphi} ``` @@ -101,7 +101,7 @@ i^{|m|+m} \sqrt{\frac{(l-|m|)!}{(l+|m|)!}} P_l^{|m|} = (-1)^{\frac{|m|+m}{2}} \s ``` #### Associated Legendre Polynomials -`P(x; n=0, m=0)` +`P(model::HydrogenAtom, x; n=0, m=0)` Rodrigues' formula & closed-form: ```math @@ -151,11 +151,11 @@ Examples: ## Usage & Examples -[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The function `antique(model, parameters...)` returns a module that has `E()`, `ψ(r)`, `V(r)` and some other functions. In this system, the model name is specified by `:HydrogenAtom` and several parameters `Z`, `Eₕ`, `mₑ`, `a₀` and `ℏ` are set as optional arguments. +[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The energy `E()`, wavefunction `ψ()`, potential `V()` and some other functions are suppoted. In this system, the model is generated by `HydrogenAtom` and several parameters `Z`, `Eₕ`, `mₑ`, `a₀` and `ℏ` are set as optional arguments. ```julia; cache = :all; results = "hidden" using Antique -H = antique(:HydrogenAtom, Z=1, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0) +H = HydrogenAtom(Z=1, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0) ``` Parameters: @@ -171,16 +171,16 @@ H.ℏ Eigen values: ```julia; term = true -H.E(n=1) -H.E(n=2) +E(H, n=1) +E(H, n=2) ``` Wave length ($n=2\rightarrow1$, the first line of the Lyman series): ```julia Eₕ2nm⁻¹ = 2.1947463136320e-2 # https://physics.nist.gov/cgi-bin/cuu/CCValue?hrminv -println("ΔE = ", H.E(n=2) - H.E(n=1), " Eₕ") -println("λ = ", ((H.E(n=2)-H.E(n=1))*Eₕ2nm⁻¹)^-1, " nm") +println("ΔE = ", E(H,n=2) - E(H,n=1), " Eₕ") +println("λ = ", ((E(H,n=2)-E(H,n=1))*Eₕ2nm⁻¹)^-1, " nm") ``` Hyperfine Splitting: @@ -198,7 +198,7 @@ ge = 2.00231930436256 # https://physics.nist.gov/cgi-bin/cuu/Value?gem gp = 5.5856946893 # https://physics.nist.gov/cgi-bin/cuu/Value?gp # calculation: https://doi.org/10.1119/1.12733 -δ = abs(H.ψ(0,0,0))^2 +δ = abs(ψ(H,0,0,0))^2 ΔE = 2 / 3 * µ0 * µN * µB * gp * ge * δ * a0^(-3) println("1/π = ", 1/π) println("<δ(r)> = ", δ, " a₀⁻³") @@ -213,7 +213,7 @@ Potential energy curve: ```julia using Plots plot(xlims=(0.0,15.0), ylims=(-0.6,0.05), xlabel="\$r~/~a_0\$", ylabel="\$V(r)/E_\\mathrm{h},~E_n/E_\\mathrm{h}\$", legend=:bottomright, size=(480,400)) -plot!(0.1:0.01:15, r -> H.V(r), lc=:black, lw=2, label="") # potential +plot!(0.1:0.01:15, r -> V(H,r), lc=:black, lw=2, label="") # potential ``` Potential energy curve, Energy levels: @@ -222,9 +222,9 @@ Potential energy curve, Energy levels: using Plots plot(xlims=(0.0,15.0), ylims=(-0.6,0.05), xlabel="\$r~/~a_0\$", ylabel="\$V(r)/E_\\mathrm{h}\$", legend=:bottomright, size=(480,400)) for n in 0:10 - plot!(0.0:0.01:15, r -> H.E(n=n) > H.V(r) ? H.E(n=n) : NaN, lc=n, lw=1, label="") # energy level + plot!(0.0:0.01:15, r -> E(H,n=n) > V(H,r) ? E(H,n=n) : NaN, lc=n, lw=1, label="") # energy level end -plot!(0.1:0.01:15, r -> H.V(r), lc=:black, lw=2, label="") # potential +plot!(0.1:0.01:15, r -> V(H,r), lc=:black, lw=2, label="") # potential ``` Radial functions: @@ -234,7 +234,7 @@ using Plots plot(xlabel="\$r~/~a_0\$", ylabel="\$r^2|R_{nl}(r)|^2~/~a_0^{-1}\$", ylims=(-0.01,0.55), xticks=0:1:20, size=(480,400), dpi=300) for n in 1:3 for l in 0:n-1 - plot!(0:0.01:20, r->r^2*H.R(r,n=n,l=l)^2, lc=n, lw=2, ls=[:solid,:dash,:dot,:dashdot,:dashdotdot][l+1], label="\$n = $n, l=$l\$") + plot!(0:0.01:20, r->r^2*R(H,r,n=n,l=l)^2, lc=n, lw=2, ls=[:solid,:dash,:dot,:dashdot,:dashdotdot][l+1], label="\$n = $n, l=$l\$") end end plot!() diff --git a/docs/src/jmd/InfinitePotentialWell.jmd b/docs/src/jmd/InfinitePotentialWell.jmd index d23a5f6..d548ce1 100644 --- a/docs/src/jmd/InfinitePotentialWell.jmd +++ b/docs/src/jmd/InfinitePotentialWell.jmd @@ -21,7 +21,7 @@ The infinite potential well (particle in a box) is the simplest model for quantu ``` #### Potential -`V(x; L=L)` +`V(model::InfinitePotentialWell; x)` ```math V(x) = \left\{ @@ -33,13 +33,13 @@ The infinite potential well (particle in a box) is the simplest model for quantu ``` #### Eigen Values -`E(; n=0, L=L, m=m, ℏ=ℏ)` +`E(model::InfinitePotentialWell; n=0)` ```math E_n = \frac{\hbar^2 n^2 \pi^2}{2 m L^2} ``` #### Eigen Functions -`ψ(x; n=0, L=L)` +`ψ(model::InfinitePotentialWell, x; n=0)` ```math \psi_n(x) = \sqrt{\frac{2}{L}} \sin \frac{n\pi x}{L} ``` @@ -50,11 +50,11 @@ The infinite potential well (particle in a box) is the simplest model for quantu ## Usage & Examples -[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The function `antique(model, parameters...)` returns a module that has `E()`, `ψ(x)`, `V(x)` and some other functions. In this system, the model name is specified by `:InfinitePotentialWell` and several parameters `L`, `m` and `ℏ` are set as optional arguments. +[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The energy `E()`, wavefunction `ψ()`, potential `V()` and some other functions are suppoted. In this system, the model is generated by `InfinitePotentialWell` and several parameters `L`, `m` and `ℏ` are set as optional arguments. ```julia; cache = :all; results = "hidden" using Antique -IPW = antique(:InfinitePotentialWell, L=1.0, m=1.0, ℏ=1.0) +IPW = InfinitePotentialWell(L=1.0, m=1.0, ℏ=1.0) ``` Parameters: @@ -68,8 +68,8 @@ IPW.ℏ Eigen values: ```julia; term = true -IPW.E(n=1) -IPW.E(n=2) +E(IPW, n=1) +E(IPW, n=2) ``` Wave functions: @@ -77,11 +77,11 @@ Wave functions: ```julia using Plots plot(xlim=(0,1), xlabel="x", ylabel="ψ(x)") -plot!(x -> IPW.ψ(x, n=1), label="n=1", lw=2) -plot!(x -> IPW.ψ(x, n=2), label="n=2", lw=2) -plot!(x -> IPW.ψ(x, n=3), label="n=3", lw=2) -plot!(x -> IPW.ψ(x, n=4), label="n=4", lw=2) -plot!(x -> IPW.ψ(x, n=5), label="n=5", lw=2) +plot!(x -> ψ(IPW, x, n=1), label="n=1", lw=2) +plot!(x -> ψ(IPW, x, n=2), label="n=2", lw=2) +plot!(x -> ψ(IPW, x, n=3), label="n=3", lw=2) +plot!(x -> ψ(IPW, x, n=4), label="n=4", lw=2) +plot!(x -> ψ(IPW, x, n=5), label="n=5", lw=2) ``` Potential energy curve, Energy levels, Wave functions: @@ -92,9 +92,9 @@ using Plots plot(xlim=(-0.5,1.5), ylim=(-5,140), xlabel="\$x\$", ylabel="\$V(x),~E_n,~\\psi_n(x)\\times5+E_n\$", size=(480,400), dpi=300) for n in 1:5 # energy - plot!([0,L], fill(IPW.E(n=n,L=L),2), lc=:black, lw=2, label="") + plot!([0,L], fill(E(IPW,n=n),2), lc=:black, lw=2, label="") # wave function - plot!(0:0.01:L, x->IPW.E(n=n,L=L)+5*IPW.ψ(x,n=n,L=L), lc=n, lw=2, label="") + plot!(0:0.01:L, x->E(IPW,n=n) + 5*ψ(IPW,x,n=n), lc=n, lw=2, label="") end # potential plot!([0,0,L,L], [140,0,0,140], lc=:black, lw=2, label="") diff --git a/docs/src/jmd/MorsePotential.jmd b/docs/src/jmd/MorsePotential.jmd index a09d10b..97b7740 100644 --- a/docs/src/jmd/MorsePotential.jmd +++ b/docs/src/jmd/MorsePotential.jmd @@ -22,25 +22,25 @@ The Morse potential is a model for inter-nuclear anharmonic vibration in a diato ``` #### Potential -`V(r; rₑ=rₑ, Dₑ=Dₑ, k=k, a=sqrt(k/(2*Dₑ)))` +`V(model::MorsePotential, r)` ```math V(r) = D_\mathrm{e} \left( \mathrm{e}^{-2a(r-r_e)} - 2\mathrm{e}^{-a(r-r_e)} \right) ``` #### Eigen Values -`E(n; rₑ=rₑ, Dₑ=Dₑ, k=k, a=sqrt(k/(2*Dₑ)), µ=µ, ω=sqrt(k/µ), χ=ℏ*ω/(4*Dₑ), ℏ=ℏ)` +`E(model::MorsePotential; n=0)` ```math E_n = - D_\mathrm{e} + \hbar \omega \left( n + \frac{1}{2} \right) - \chi \hbar \omega \left( n + \frac{1}{2} \right)^2 ``` #### Eigen Functions -`ψ(n, r; rₑ=rₑ, Dₑ=Dₑ, k=k, a=sqrt(k/(2*Dₑ)), µ=µ, ω=sqrt(k/µ), χ=ℏ*ω/(4*Dₑ), ℏ=ℏ)` +`ψ(model::MorsePotential, r; n=0)` ```math \psi_n(r) = N_n z^{\lambda-n-1/2} \mathrm{e}^{-z/2} L_n^{(2\lambda-2n-1)}(\xi) ``` #### Generalized Laguerre Polynomials -`L(x; n=0, α=0)` +`L(model::MorsePotential, x; n=0, α=0)` Rodrigues' formula & closed-form: ```math @@ -81,11 +81,7 @@ Examples: ## Usage & Examples -[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The function `antique(model, parameters...)` returns a module that has `E()`, `ψ(r)`, `V(r)` and some other functions. In this system, the model name is specified by `:MorsePotential` and several parameters `rₑ`, `Dₑ`, `k`, `µ` and `ℏ` are set as optional arguments. - -!!! warning - - Run `using Pkg; Pkg.add("SpecialFunctions")` if the following returns an error. +[Install Antique.jl](@ref Install) for the first use and run `using Antique` before each use. The energy `E()`, wavefunction `ψ()`, potential `V()` and some other functions are suppoted. In this system, the model is generated by `MorsePotential` and several parameters `rₑ`, `Dₑ`, `k`, `µ` and `ℏ` are set as optional arguments. ```julia; cache = :all; results = "hidden" # Parameters for H₂⁺ @@ -93,14 +89,13 @@ Examples: # https://doi.org/10.5281/zenodo.5047817 # https://physics.nist.gov/cgi-bin/cuu/Value?mpsme rₑ = 1.997193319969992120068298141276 -Vₑ = -0.602634619106539878727562156289 -Dₑ = - 0.5 - Vₑ -k = 2*((-1.1026342144949464615+1/2.00) - Vₑ) / (2.00 - rₑ)^2 +Dₑ = - 0.5 - (-0.602634619106539878727562156289) +k = 2*((-1.1026342144949464615+1/2.00) - (-0.602634619106539878727562156289)) / (2.00 - rₑ)^2 µ = 1/(1/1836.15267343 + 1/1836.15267343) ℏ = 1.0 using Antique -MP = antique(:MorsePotential, rₑ=rₑ, Vₑ=Vₑ, Dₑ=Dₑ, k=k, µ =µ, ℏ =ℏ) +MP = MorsePotential(rₑ=rₑ, Dₑ=Dₑ, k=k, µ=µ, ℏ=ℏ) ``` Parameters: @@ -116,15 +111,15 @@ MP.ℏ Eigen values: ```julia; term = true -MP.E(n=0) -MP.E(n=1) +E(MP, n=0) +E(MP, n=1) ``` Potential energy curve: ```julia using Plots -plot(0.1:0.01:15, r -> MP.V(r), lw=2, label="", xlims=(0.1,9.1), ylims=(-0.11,0.01), xlabel="r", ylabel="V(r)") +plot(0.1:0.01:15, r -> V(MP, r), lw=2, label="", xlims=(0.1,9.1), ylims=(-0.11,0.01), xlabel="r", ylabel="V(r)") ``` Wave functions: @@ -132,31 +127,31 @@ Wave functions: ```julia using Plots plot(xlim=(0,5), xlabel="x", ylabel="ψ(x)") -plot!(x -> MP.ψ(x, n=0), label="n=0", lw=2) -plot!(x -> MP.ψ(x, n=1), label="n=1", lw=2) -plot!(x -> MP.ψ(x, n=2), label="n=2", lw=2) -plot!(x -> MP.ψ(x, n=3), label="n=3", lw=2) -plot!(x -> MP.ψ(x, n=4), label="n=4", lw=2) -plot!(x -> MP.ψ(x, n=5), label="n=5", lw=2) +plot!(x -> ψ(MP, x, n=0), label="n=0", lw=2) +plot!(x -> ψ(MP, x, n=1), label="n=1", lw=2) +plot!(x -> ψ(MP, x, n=2), label="n=2", lw=2) +plot!(x -> ψ(MP, x, n=3), label="n=3", lw=2) +plot!(x -> ψ(MP, x, n=4), label="n=4", lw=2) +plot!(x -> ψ(MP, x, n=5), label="n=5", lw=2) ``` Potential energy curve, Energy levels, Comparison with harmonic oscillator: ```julia -MP = antique(:MorsePotential) -HO = antique(:HarmonicOscillator, k=MP.k, m=MP.μ) +MP = MorsePotential() +HO = HarmonicOscillator(k=MP.k, m=MP.μ) using Plots plot(xlims=(0.1,9.1), ylims=(-0.11,0.01), xlabel="\$r\$", ylabel="\$V(r), E_n\$", legend=:bottomright, size=(480,400), dpi=300) -for n in 0:MP.nₘₐₓ() +for n in 0:nₘₐₓ(MP) # energy - EM = MP.E(n=n) - EH = HO.E(n=n) - MP.Dₑ - plot!(0.1:0.01:15, r -> EH > HO.V(r-MP.rₑ) - MP.Dₑ ? EH : NaN, lc="#BC1C5F", lw=1, label="") - plot!(0.1:0.01:15, r -> EM > MP.V(r) ? EM : NaN, lc="#578FC7", lw=1, label="") + EM = E(MP, n=n) + EH = E(HO, n=n) - MP.Dₑ + plot!(0.1:0.01:15, r -> EH > V(HO, r-MP.rₑ) - MP.Dₑ ? EH : NaN, lc="#BC1C5F", lw=1, label="") + plot!(0.1:0.01:15, r -> EM > V(MP, r) ? EM : NaN, lc="#578FC7", lw=1, label="") end # potential -plot!(0.1:0.01:15, r -> HO.V(r-MP.rₑ) - MP.Dₑ, lc="#BC1C5F", lw=2, label="Harmonic Oscillator") -plot!(0.1:0.01:15, r -> MP.V(r), lc="#578FC7", lw=2, label="Morse Potential") +plot!(0.1:0.01:15, r -> V(HO, r-MP.rₑ) - MP.Dₑ, lc="#BC1C5F", lw=2, label="Harmonic Oscillator") +plot!(0.1:0.01:15, r -> V(MP, r), lc="#578FC7", lw=2, label="Morse Potential") ``` where, the potential of harmonic oscillator is defined as $V(r) \simeq \frac{1}{2} k (r - r_\mathrm{e})^2 + V_0$. diff --git a/src/Antique.jl b/src/Antique.jl index 0a3006c..2826d60 100644 --- a/src/Antique.jl +++ b/src/Antique.jl @@ -1,7 +1,5 @@ module Antique - export antique - # Update this list when you add a model. models = [ :InfinitePotentialWell, @@ -15,6 +13,23 @@ module Antique include("./$(model).jl") end + # --------------- for old version --------------- + + export antique + + # Update this list when you add a model. + oldmodels = [ + :InfinitePotentialWell, + :HarmonicOscillator, + :MorsePotential, + :HydrogenAtom, + ] + + # include statements + for model in oldmodels + include("./old/$(model).jl") + end + # I/O function function save(path, text) mkpath(dirname(path)) @@ -33,15 +48,17 @@ module Antique # main functions function antique(model; parameters...) + # alart + @warn "The function `antique(model; parameters...)` is deprecated. Please check the latest document: https://ohno.github.io/Antique.jl/stable/" # check existence of model if model ∉ models throw(ErrorException("\`:$(model)\` is not in the list of supported models $(models).")) end # load source code file if Sys.iswindows() - source = load("$(dirname(@__FILE__()))\\\\$(model).jl") + source = load("$(dirname(@__FILE__()))\\\\old\\\\$(model).jl") else - source = load("$(dirname(@__FILE__()))/$(model).jl") + source = load("$(dirname(@__FILE__()))/old/$(model).jl") end # replace name of module for m in eachmatch(Regex("module $(model)"), source) diff --git a/src/HarmonicOscillator.jl b/src/HarmonicOscillator.jl index bfe3e32..ef7ac47 100644 --- a/src/HarmonicOscillator.jl +++ b/src/HarmonicOscillator.jl @@ -1,24 +1,39 @@ -module HarmonicOscillator +export HarmonicOscillator, V, E, ψ, H - # Default - k = 1.0 # change here! - m = 1.0 # change here! - ℏ = 1.0 # change here! - - # Potential - V(x; k=k, m=m, ω=sqrt(k/m)) = 1/2*k*x^2 +# Parameters +Base.@kwdef struct HarmonicOscillator + k = 1.0 + m = 1.0 + ℏ = 1.0 +end - # Energy - E(; n=0, k=k, m=m, ω=sqrt(k/m), ℏ=ℏ) = ℏ*ω*(n+1/2) +# Potential +function V(model::HarmonicOscillator, x) + k = model.k + return 1//2 * k * x^2 +end - # Wave Function - function ψ(x; n=0, k=k, m=m, ω=sqrt(k/m), ℏ=ℏ) - A = sqrt(1/(factorial(n)*2^n)*sqrt(m*ω/(π*ℏ))) - ξ = sqrt(m*ω/ℏ) * x - return A*H(ξ,n=n)*exp(-ξ^2/2) - end +# Energy +function E(model::HarmonicOscillator; n=0) + k = model.k + m = model.m + ℏ = model.ℏ + ω = sqrt(k/m) + return ℏ * ω * (n+1//2) +end - # Hermite polynomials - H(x; n=0) = factorial(n) * sum(i -> (-1)^i // (factorial(i) * factorial(n-2*i)) * (2*x)^(n-2*i), 0:Int(floor(n/2))) +# Wave Function +function ψ(model::HarmonicOscillator, x; n=0) + k = model.k + m = model.m + ℏ = model.ℏ + ω = sqrt(k/m) + A = sqrt(1//(factorial(n)*2^n)*sqrt(m*ω/(π*ℏ))) + ξ = sqrt(m*ω/ℏ) * x + return A * H(model,ξ,n=n) * exp(-ξ^2/2) +end +# Hermite Polynomials +function H(model::HarmonicOscillator, x; n=0) + return factorial(n) * sum(i -> (-1)^i // (factorial(i) * factorial(n-2*i)) * (2*x)^(n-2*i), 0:Int(floor(n/2))) end diff --git a/src/HydrogenAtom.jl b/src/HydrogenAtom.jl index 7fb2d03..6cc0a3f 100644 --- a/src/HydrogenAtom.jl +++ b/src/HydrogenAtom.jl @@ -1,41 +1,65 @@ -module HydrogenAtom - - # Default - Z = 1 # change here! - ℏ = 1.0 # change here! - Eₕ = 1.0 # change here! - a₀ = 1.0 # change here! - mₑ = 1.0 # change here! - - # Potential - V(r; Z=Z, a₀=a₀) = r<0 ? throw(DomainError(r, "r=$r is out of the domain (0≦r)")) : -Z/abs(r/a₀) - - # Energy - E(; n=1, Z=Z, Eₕ=Eₕ) = -Z^2/(2*n^2) * Eₕ - - # Wave Function - ψ(r, θ, φ; n=1, l=0, m=0, Z=Z, a₀=a₀) = R(r, n=n, l=l, Z=Z, a₀=a₀) * Y(θ, φ, l=l, m=m) - - # Radial Function - function R(r; n=1, l=0, Z=Z, a₀=a₀) - if r<0 - throw(DomainError(r, "r=$r is out of the domain (0≦r)")) - end - ρ = 2*Z*abs(r)/(n*a₀) - N = -sqrt( factorial(n-l-1)/(2*n*factorial(n+l)) * (2*Z/(n*a₀))^3 ) - return N*ρ^l * exp(-ρ/2) * L(ρ, n=n+l, k=2*l+1) - end - - # Associated Laguerre Polynomials - L(x; n=0, k=0) = sum(m -> (-1)^(m+k) * factorial(n) // (factorial(m) * factorial(m+k) * factorial(n-m-k)) * x^m, 0:n-k) - - # Spherical Harmonics - function Y(θ, φ; l=0, m=0) - N = (im)^(m+abs(m)) * sqrt( (2*l+1)*factorial(l-Int(abs(m))) / (2*factorial(l+Int(abs(m)))) ) - return N * P(cos(θ),n=l,m=Int(abs(m))) * exp(im*m*φ) / sqrt(2*π) - end - - # Associated Legendre Polynomials - P(x; n=0, m=0) = (1//2)^n * (1-x^2)^(m//2) * sum(j -> (-1)^j * factorial(2*n-2*j) // (factorial(j) * factorial(n-j) * factorial(n-2*j-m)) * x^(n-2*j-m), 0:Int(floor((n-m)/2))) +export HydrogenAtom, V, E, ψ, R, L, Y, P +# Parameters +Base.@kwdef struct HydrogenAtom + Z = 1 + ℏ = 1.0 + Eₕ = 1.0 + a₀ = 1.0 + mₑ = 1.0 end + +# Potential +function V(model::HydrogenAtom, r) + # if r<0 + # throw(DomainError(r, "r=$r is out of the domain (0≦r)")) + # end + Z = model.Z + a₀ = model.a₀ + return -Z/abs(r/a₀) +end + +# Energy +function E(model::HydrogenAtom; n=1) + Z = model.Z + Eₕ = model.Eₕ + return -Z^2/(2*n^2) * Eₕ +end + +# Wave Function +function ψ(model::HydrogenAtom, r, θ, φ; n=1, l=0, m=0) + # if r<0 + # throw(DomainError(r, "r=$r is out of the domain (0≦r)")) + # end + Z = model.Z + a₀ = model.a₀ + return R(model, r, n=n, l=l) * Y(model, θ, φ, l=l, m=m) +end + +# Radial Function +function R(model::HydrogenAtom, r; n=1, l=0) + # if r<0 + # throw(DomainError(r, "r=$r is out of the domain (0≦r)")) + # end + Z = model.Z + a₀ = model.a₀ + ρ = 2*Z*abs(r)/(n*a₀) + N = -sqrt( factorial(n-l-1)/(2*n*factorial(n+l)) * (2*Z/(n*a₀))^3 ) + return N*ρ^l * exp(-ρ/2) * L(model, ρ, n=n+l, k=2*l+1) +end + +# Associated Laguerre Polynomials +function L(model::HydrogenAtom, x; n=0, k=0) + return sum(m -> (-1)^(m+k) * factorial(n) // (factorial(m) * factorial(m+k) * factorial(n-m-k)) * x^m, 0:n-k) +end + +# Spherical Harmonics +function Y(model::HydrogenAtom, θ, φ; l=0, m=0) + N = (im)^(m+abs(m)) * sqrt( (2*l+1)*factorial(l-Int(abs(m))) / (2*factorial(l+Int(abs(m)))) ) + return N * P(model,cos(θ), n=l, m=Int(abs(m))) * exp(im*m*φ) / sqrt(2*π) +end + +# Associated Legendre Polynomials +function P(model::HydrogenAtom, x; n=0, m=0) + return (1//2)^n * (1-x^2)^(m//2) * sum(j -> (-1)^j * factorial(2*n-2*j) // (factorial(j) * factorial(n-j) * factorial(n-2*j-m)) * x^(n-2*j-m), 0:Int(floor((n-m)/2))) +end \ No newline at end of file diff --git a/src/InfinitePotentialWell.jl b/src/InfinitePotentialWell.jl index 8ac2b14..27e6a6f 100644 --- a/src/InfinitePotentialWell.jl +++ b/src/InfinitePotentialWell.jl @@ -1,17 +1,28 @@ -module InfinitePotentialWell +export InfinitePotentialWell, V, E, ψ - # Default - L = 1.0 # change here! - m = 1.0 # change here! - ℏ = 1.0 # change here! +# Parameters +Base.@kwdef struct InfinitePotentialWell + L = 1.0 + m = 1.0 + ℏ = 1.0 +end - # Potential - V(x; L=L) = 0 (-1)^(k) * (gamma(α+n+1) / (gamma(α+1+k)*gamma(n-k+1))) * x^k / factorial(k), 0:n) - Lαint(x; n=0, α::Int=0) = sum(k -> (-1)^(k) * (Int(gamma(α+n+1)) // Int((gamma(α+1+k)*gamma(n-k+1)))) * x^k // factorial(k), 0:n) # α::Int +# Wave Function +function ψ(model::MorsePotential, r; n=0) + if r<0 + throw(DomainError(r, "r=$r is out of the domain (0≦r)")) + end + rₑ = model.rₑ + Dₑ = model.Dₑ + k = model.k + µ = model.µ + ℏ = model.ℏ + a = sqrt(k/(2*Dₑ)) + λ = sqrt(2*µ*Dₑ) / (a*ℏ) + ξ = 2*λ*exp(-a*(r-rₑ)) + s = 2*λ - 2*n - 1 + N = sqrt(factorial(n) * s * a / gamma(s+n+1)) + return N * ξ^(s/2) * exp(-ξ/2) * L(model, ξ, n=n, α=s) +end +# Generalized Laguerre Polynomials +function L(model::MorsePotential, x; n=0, α=0) + if isinteger(α) + return sum(k -> (-1)^(k) * (Int(gamma(α+n+1)) // Int((gamma(α+1+k)*gamma(n-k+1)))) * x^k * 1 // factorial(k), 0:n) + else + return sum(k -> (-1)^(k) * (gamma(α+n+1) / (gamma(α+1+k)*gamma(n-k+1))) * x^k / factorial(k), 0:n) + end end \ No newline at end of file diff --git a/src/old/HarmonicOscillator.jl b/src/old/HarmonicOscillator.jl new file mode 100644 index 0000000..1c5891b --- /dev/null +++ b/src/old/HarmonicOscillator.jl @@ -0,0 +1,24 @@ +module OldHarmonicOscillator + + # Default + k = 1.0 # change here! + m = 1.0 # change here! + ℏ = 1.0 # change here! + + # Potential + V(x; k=k, m=m, ω=sqrt(k/m)) = 1/2*k*x^2 + + # Energy + E(; n=0, k=k, m=m, ω=sqrt(k/m), ℏ=ℏ) = ℏ*ω*(n+1/2) + + # Wave Function + function ψ(x; n=0, k=k, m=m, ω=sqrt(k/m), ℏ=ℏ) + A = sqrt(1/(factorial(n)*2^n)*sqrt(m*ω/(π*ℏ))) + ξ = sqrt(m*ω/ℏ) * x + return A*H(ξ,n=n)*exp(-ξ^2/2) + end + + # Hermite polynomials + H(x; n=0) = factorial(n) * sum(i -> (-1)^i // (factorial(i) * factorial(n-2*i)) * (2*x)^(n-2*i), 0:Int(floor(n/2))) + +end diff --git a/src/old/HydrogenAtom.jl b/src/old/HydrogenAtom.jl new file mode 100644 index 0000000..294245d --- /dev/null +++ b/src/old/HydrogenAtom.jl @@ -0,0 +1,41 @@ +module OldHydrogenAtom + + # Default + Z = 1 # change here! + ℏ = 1.0 # change here! + Eₕ = 1.0 # change here! + a₀ = 1.0 # change here! + mₑ = 1.0 # change here! + + # Potential + V(r; Z=Z, a₀=a₀) = r<0 ? throw(DomainError(r, "r=$r is out of the domain (0≦r)")) : -Z/abs(r/a₀) + + # Energy + E(; n=1, Z=Z, Eₕ=Eₕ) = -Z^2/(2*n^2) * Eₕ + + # Wave Function + ψ(r, θ, φ; n=1, l=0, m=0, Z=Z, a₀=a₀) = R(r, n=n, l=l, Z=Z, a₀=a₀) * Y(θ, φ, l=l, m=m) + + # Radial Function + function R(r; n=1, l=0, Z=Z, a₀=a₀) + if r<0 + throw(DomainError(r, "r=$r is out of the domain (0≦r)")) + end + ρ = 2*Z*abs(r)/(n*a₀) + N = -sqrt( factorial(n-l-1)/(2*n*factorial(n+l)) * (2*Z/(n*a₀))^3 ) + return N*ρ^l * exp(-ρ/2) * L(ρ, n=n+l, k=2*l+1) + end + + # Associated Laguerre Polynomials + L(x; n=0, k=0) = sum(m -> (-1)^(m+k) * factorial(n) // (factorial(m) * factorial(m+k) * factorial(n-m-k)) * x^m, 0:n-k) + + # Spherical Harmonics + function Y(θ, φ; l=0, m=0) + N = (im)^(m+abs(m)) * sqrt( (2*l+1)*factorial(l-Int(abs(m))) / (2*factorial(l+Int(abs(m)))) ) + return N * P(cos(θ),n=l,m=Int(abs(m))) * exp(im*m*φ) / sqrt(2*π) + end + + # Associated Legendre Polynomials + P(x; n=0, m=0) = (1//2)^n * (1-x^2)^(m//2) * sum(j -> (-1)^j * factorial(2*n-2*j) // (factorial(j) * factorial(n-j) * factorial(n-2*j-m)) * x^(n-2*j-m), 0:Int(floor((n-m)/2))) + +end diff --git a/src/old/InfinitePotentialWell.jl b/src/old/InfinitePotentialWell.jl new file mode 100644 index 0000000..3449f9d --- /dev/null +++ b/src/old/InfinitePotentialWell.jl @@ -0,0 +1,17 @@ +module OldInfinitePotentialWell + + # Default + L = 1.0 # change here! + m = 1.0 # change here! + ℏ = 1.0 # change here! + + # Potential + V(x; L=L) = 0 (-1)^(k) * (gamma(α+n+1) / (gamma(α+1+k)*gamma(n-k+1))) * x^k / factorial(k), 0:n) + Lαint(x; n=0, α::Int=0) = sum(k -> (-1)^(k) * (Int(gamma(α+n+1)) // Int((gamma(α+1+k)*gamma(n-k+1)))) * x^k // factorial(k), 0:n) # α::Int + +end \ No newline at end of file diff --git a/test/HarmonicOscillator.jl b/test/HarmonicOscillator.jl index d764ab0..17be73f 100644 --- a/test/HarmonicOscillator.jl +++ b/test/HarmonicOscillator.jl @@ -6,7 +6,7 @@ using QuadGK using Symbolics using Latexify using LaTeXStrings -HO = antique(:HarmonicOscillator, k=1.0, m=1.0, ℏ=1.0) +HO = HarmonicOscillator(k=1.0, m=1.0, ℏ=1.0) # Hₙ(x) = (-1)ⁿ exp(x²) dⁿ/dxⁿ　exp(-x²) = ... @@ -34,7 +34,7 @@ println(raw""" c = a * D(b) # Rodrigues' formula d = expand_derivatives(c) # expand dⁿ/dxⁿ e = simplify(d, expand=true) # simplify - f = simplify(HO.H(x, n=n), expand=true) # closed-form + f = simplify(H(HO, x, n=n), expand=true) # closed-form # latexify eq1 = latexify(e, env=:raw) eq2 = latexify(f, env=:raw) @@ -79,7 +79,7 @@ println(raw""" for i in 0:9 for j in 0:9 analytical = sqrt(π)*2^j*factorial(j)*(i == j ? 1 : 0) - numerical = quadgk(x -> HO.H(x, n=j) * HO.H(x, n=i)* exp(-x^2), -Inf, Inf, maxevals=10^3)[1] + numerical = quadgk(x -> H(HO, x, n=j) * H(HO, x, n=i)* exp(-x^2), -Inf, Inf, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%2d | %2d | %17.12f | %17.12f %s\n", i, j, analytical, numerical, acceptance ? "✔" : "✗") @@ -109,7 +109,7 @@ println(raw""" for i in 0:9 for j in 0:9 analytical = (i == j ? 1 : 0) - numerical = quadgk(x -> conj(HO.ψ(x, n=i)) * HO.ψ(x, n=j), -Inf, Inf, maxevals=10^3)[1] + numerical = quadgk(x -> conj(ψ(HO, x, n=i)) * ψ(HO, x, n=j), -Inf, Inf, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%2d | %2d | %17.12f | %17.12f %s\n", i, j, analytical, numerical, acceptance ? "✔" : "✗") @@ -140,8 +140,8 @@ The virial theorem $\langle T \rangle = \langle V \rangle$ and the definition of println("--- | -- | ----------------- | ----------------- ") for k in [0.1,0.5,1.0,5.0] for n in 0:9 - analytical = HO.E(n=n) - numerical = 2 * quadgk(x -> conj(HO.ψ(x, n=n)) * HO.V(x) * HO.ψ(x, n=n), -Inf, Inf, maxevals=10^3)[1] + analytical = E(HO, n=n) + numerical = 2 * quadgk(x -> conj(ψ(HO, x, n=n)) * V(HO, x) * ψ(HO, x, n=n), -Inf, Inf, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%.1f | %2d | %17.12f | %17.12f %s\n", k, n, analytical, numerical, acceptance ? "✔" : "✗") @@ -222,19 +222,21 @@ are given by the sum of 2 Taylor series: ```""") @testset "∫ψₙ*Hψₙdx = <ψₙ|H|ψₙ> = Eₙ" begin - ψHψ(x; n=0, k=HO.k, m=HO.m, ℏ=HO.ℏ, Δx=0.005) = HO.V(x,k=k,m=m)*HO.ψ(x,n=n,k=k,m=m,ℏ=ℏ)^2 - ℏ^2/(2*m)*conj(HO.ψ(x,n=n,k=k,m=m,ℏ=ℏ))*(HO.ψ(x+Δx,n=n,k=k,m=m,ℏ=ℏ)-2*HO.ψ(x,n=n,k=k,m=m,ℏ=ℏ)+HO.ψ(x-Δx,n=n,k=k,m=m,ℏ=ℏ))/Δx^2 + ψHψ(HO, x; n=0, Δx=0.005) = V(HO,x)*ψ(HO,x,n=n)^2 - HO.ℏ^2/(2*HO.m)*conj(ψ(HO,x,n=n))*(ψ(HO,x+Δx,n=n)-2*ψ(HO,x,n=n)+ψ(HO,x-Δx,n=n))/Δx^2 println(" k | n | analytical | numerical ") println("--- | -- | ----------------- | ----------------- ") for k in [0.1,0.5,1.0,5.0] for n in 0:9 - analytical = HO.E(n=n, k=k) - numerical = quadgk(x -> ψHψ(x, n=n, k=k, Δx=0.001), -Inf, Inf, maxevals=10^3)[1] + HO = HarmonicOscillator(k=k) + analytical = E(HO, n=n) + numerical = quadgk(x -> ψHψ(HO, x, n=n, Δx=0.001), -Inf, Inf, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%.1f | %2d | %17.12f | %17.12f %s\n", k, n, analytical, numerical, acceptance ? "✔" : "✗") end end end +HO = HarmonicOscillator(k=1.0, m=1.0, ℏ=1.0) println("""``` """) diff --git a/test/HydrogenAtom.jl b/test/HydrogenAtom.jl index 4f623d0..f150e123a 100644 --- a/test/HydrogenAtom.jl +++ b/test/HydrogenAtom.jl @@ -6,8 +6,8 @@ using QuadGK using Symbolics using Latexify using LaTeXStrings -HA = antique(:HydrogenAtom, Z=1, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0) -MP = antique(:MorsePotential) +HA = HydrogenAtom(Z=1, Eₕ=1.0, a₀=1.0, mₑ=1.0, ℏ=1.0) +MP = MorsePotential() # Pₙᵐ(x) = √(1-x²)ᵐ dᵐ/dxᵐ Pₙ(x); Pₙ(x) = 1/(2ⁿn!) dⁿ/dxⁿ (x²-1)ⁿ @@ -38,7 +38,7 @@ println(raw""" c = (1 - x^2)^(m//2) * Dm(a * Dn(b)) # Rodrigues' formula d = expand_derivatives(c) # expand dⁿ/dxⁿ and dᵐ/dxᵐ e = simplify(d, expand=true) # simplify - f = simplify(HA.P(x, n=n, m=m), expand=true) # closed-form + f = simplify(P(HA, x, n=n, m=m), expand=true) # closed-form # latexify eq1 = latexify(e, env=:raw) eq2 = latexify(f, env=:raw) @@ -84,7 +84,7 @@ println(raw""" for i in m:9 for j in m:9 analytical = 2*factorial(j+m)/(2*j+1)/factorial(j-m)*(i == j ? 1 : 0) - numerical = quadgk(x -> HA.P(x, n=i, m=m) * HA.P(x, n=j, m=m), -1, 1, maxevals=10^3)[1] + numerical = quadgk(x -> P(HA, x, n=i, m=m) * P(HA, x, n=j, m=m), -1, 1, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%2d | %2d | %2d | %17.12f | %17.12f %s\n", m, i, j, analytical, numerical, acceptance ? "✔" : "✗") @@ -123,7 +123,7 @@ Y_{lm}(\theta,\varphi)^* Y_{l'm'}(\theta,\varphi) \sin(\theta) numerical = real( quadgk(φ -> quadgk(θ -> - conj(HA.Y(θ,φ,l=l1,m=m1)) * HA.Y(θ,φ,l=l2,m=m2) * sin(θ) + conj(Y(HA,θ,φ,l=l1,m=m1)) * Y(HA,θ,φ,l=l2,m=m2) * sin(θ) , 0, π, maxevals=50)[1] , 0, 2π, maxevals=100)[1] ) @@ -162,15 +162,15 @@ println(raw""" for k in 0:n # Rodriguesの公式の展開 @variables x - Dn = n==0 ? x->x : Differential(x)^n # dⁿ/dxⁿ - Dk = k==0 ? x->x : Differential(x)^k # dᵐ/dxᵐ - a = exp(x) / factorial(n) # left - b = exp(-x) * x^n # right - c = Dk(a * Dn(b)) # Rodrigues' formula - d = expand_derivatives(c) # expand dⁿ/dxⁿ and dᵐ/dxᵐ - e = simplify(d, expand=true) # simplify - f = simplify(HA.L(x, n=n, k=k), expand=true) # closed-form - g = simplify((-1)^k * MP.Lαint(x, n=n-k, α=k), expand=true) # closed-form + Dn = n==0 ? x->x : Differential(x)^n # dⁿ/dxⁿ + Dk = k==0 ? x->x : Differential(x)^k # dᵐ/dxᵐ + a = exp(x) / factorial(n) # left + b = exp(-x) * x^n # right + c = Dk(a * Dn(b)) # Rodrigues' formula + d = expand_derivatives(c) # expand dⁿ/dxⁿ and dᵐ/dxᵐ + e = simplify(d, expand=true) # simplify + f = simplify(L(HA, x, n=n, k=k), expand=true) # closed-form + g = simplify((-1)^k * L(MP, x, n=n-k, α=k), expand=true) # closed-form # latexify eq1 = latexify(e, env=:raw) eq2 = latexify(f, env=:raw) @@ -220,7 +220,7 @@ Replace $n+k$ with $n$ for [the definition of Wolfram MathWorld](https://mathwor for j in 0:7 for k in 0:min(i,j) analytical = factorial(i) / factorial(i-k) * (i == j ? 1 : 0) - numerical = quadgk(x -> exp(-x) * x^k * HA.L(x, n=i, k=k) * HA.L(x, n=j, k=k), 0, Inf, maxevals=10^3)[1] + numerical = quadgk(x -> exp(-x) * x^k * L(HA, x, n=i, k=k) * L(HA, x, n=j, k=k), 0, Inf, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%2d | %2d | %2d | %17.12f | %17.12f %s\n", i, j, k, analytical, numerical, acceptance ? "✔" : "✗") @@ -250,7 +250,7 @@ println(raw""" for n in 1:9 for l in 0:n-1 analytical = 1 - numerical = quadgk(r -> r^2 * HA.R(r,n=n,l=l)^2, 0, Inf, maxevals=10^3)[1] + numerical = quadgk(r -> r^2 * R(HA,r,n=n,l=l)^2, 0, Inf, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%2d | %2d | %17.12f | %17.12f %s\n", n, l, analytical, numerical, acceptance ? "✔" : "✗") @@ -288,7 +288,7 @@ Reference: for n in 1:9 for l in 0:n-1 analytical = HA.a₀/2/HA.Z * (3*n^2-l*(l+1)) - numerical = quadgk(r -> r^3 * HA.R(r,n=n,l=l)^2, 0, Inf, maxevals=10^3)[1] + numerical = quadgk(r -> r^3 * R(HA,r,n=n,l=l)^2, 0, Inf, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%2d | %2d | %17.12f | %17.12f %s\n", n, l, analytical, numerical, acceptance ? "✔" : "✗") @@ -325,7 +325,7 @@ Reference: for n in 1:9 for l in 0:n-1 analytical = HA.a₀^2/2/HA.Z^2 * n^2*(5*n^2+1-3*l*(l+1)) - numerical = quadgk(r -> r^4 * HA.R(r,n=n,l=l)^2, 0, Inf, maxevals=10^3)[1] + numerical = quadgk(r -> r^4 * R(HA,r,n=n,l=l)^2, 0, Inf, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%2d | %2d | %17.12f | %17.12f %s\n", n, l, analytical, numerical, acceptance ? "✔" : "✗") @@ -354,8 +354,8 @@ The virial theorem $2\langle T \rangle + \langle V \rangle = 0$ and the definiti println(" n | analytical | numerical ") println("-- | ----------------- | ----------------- ") for n in 1:10 - analytical = HA.E(n=n) * 2 - numerical = quadgk(r -> 4*π*r^2 * conj(HA.ψ(r,0,0, n=n)) * HA.V(r) * HA.ψ(r,0,0, n=n), 0, Inf, maxevals=10^3)[1] + analytical = E(HA, n=n) * 2 + numerical = quadgk(r -> 4*π*r^2 * conj(ψ(HA,r,0,0, n=n)) * V(HA,r) * ψ(HA,r,0,0, n=n), 0, Inf, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%2d | %17.12f | %17.12f %s\n", n, analytical, numerical, acceptance ? "✔" : "✗") @@ -391,7 +391,7 @@ println(raw""" quadgk(phi -> quadgk(theta -> quadgk(r -> - r^2 * sin(theta) * conj(HA.ψ(r,theta,phi,n=n1,l=l1,m=m1)) * HA.ψ(r,theta,phi,n=n2,l=l2,m=m2) + r^2 * sin(theta) * conj(ψ(HA,r,theta,phi,n=n1,l=l1,m=m1)) * ψ(HA,r,theta,phi,n=n2,l=l2,m=m2) , 0, Inf, maxevals=50)[1] , 0, π, maxevals=4)[1] , 0, 2π, maxevals=8)[1] diff --git a/test/InfinitePotentialWell.jl b/test/InfinitePotentialWell.jl index daccbb3..6bd31a7 100644 --- a/test/InfinitePotentialWell.jl +++ b/test/InfinitePotentialWell.jl @@ -3,7 +3,7 @@ using Test using Printf using Markdown using QuadGK -IPW = antique(:InfinitePotentialWell, L=1.0, m=1.0, ℏ=1.0) +IPW = InfinitePotentialWell(L=1.0, m=1.0, ℏ=1.0) # <ψᵢ|ψⱼ> = ∫ψₙ*ψₙdx = δᵢⱼ @@ -27,7 +27,7 @@ println(raw""" for i in 1:10 for j in 1:10 analytical = (i == j ? 1 : 0) - numerical = quadgk(x -> conj(IPW.ψ(x, n=i, L=IPW.L)) * IPW.ψ(x, n=j, L=IPW.L), 0.0, IPW.L, maxevals=10^3)[1] + numerical = quadgk(x -> conj(ψ(IPW, x, n=i)) * ψ(IPW, x, n=j), 0.0, IPW.L, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%2d | %2d | %17.12f | %17.12f %s\n", i, j, analytical, numerical, acceptance ? "✔" : "✗") @@ -110,15 +110,16 @@ are given by the sum of 2 Taylor series: ```""") @testset "<ψₙ|H|ψₙ> = ∫ψₙ*Tψₙdx = Eₙ" begin - ψTψ(x; n=0, L=IPW.L, m=IPW.m, ℏ=IPW.ℏ, Δx=0.01) = -ℏ^2/(2*m)*conj(IPW.ψ(x,n=n,L=L))*(IPW.ψ(x+Δx,n=n,L=L)-2*IPW.ψ(x,n=n,L=L)+IPW.ψ(x-Δx,n=n,L=L))/Δx^2 + ψTψ(IPW, x; n=0, Δx=0.01) = -IPW.ℏ^2/(2*IPW.m)*conj(ψ(IPW,x,n=n))*(ψ(IPW,x+Δx,n=n)-2*ψ(IPW,x,n=n)+ψ(IPW,x-Δx,n=n))/Δx^2 println(" L | m | ℏ | n | analytical | numerical ") println("--- | --- | --- | -- | ----------------- | ----------------- ") for L in [0.1, 1.0] for m in [0.1, 1.0] for ℏ in [0.1, 1.0] for n in 1:10 - analytical = IPW.E(n=n, L=L, m=m, ℏ=ℏ) - numerical = quadgk(x -> ψTψ(x, n=n, L=L, m=m, ℏ=ℏ, Δx=L*0.0001), 0, L, maxevals=10^3)[1] + IPW = InfinitePotentialWell(L=L, m=m, ℏ=ℏ) + analytical = E(IPW, n=n) + numerical = quadgk(x -> ψTψ(IPW, x, n=n, Δx=L*0.0001), 0, L, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%.1f | %.1f | %.1f | %2d | %17.12f | %17.12f %s\n", L, m, ℏ, n, numerical, analytical, acceptance ? "✔" : "✗") @@ -153,8 +154,9 @@ Reference: println("--- | -- | ----------------- | ----------------- ") for L in [0.1, 0.5, 1.0, 7.0] for n in 1:1 + IPW = InfinitePotentialWell(L=L, m=1.0, ℏ=1.0) analytical = L/2 - numerical = quadgk(x -> conj(IPW.ψ(x, n=n, L=L)) * x * IPW.ψ(x, n=n, L=L), 0, L, maxevals=10^3)[1] + numerical = quadgk(x -> conj(ψ(IPW, x, n=n)) * x * ψ(IPW, x, n=n), 0, L, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%.1f | %2d | %17.12f | %17.12f %s\n", L, n, numerical, analytical, acceptance ? "✔" : "✗") @@ -186,8 +188,9 @@ Reference: println("--- | -- | ----------------- | ----------------- ") for L in [0.1, 0.5, 1.0, 7.0] for n in 1:1 + IPW = InfinitePotentialWell(L=L, m=1.0, ℏ=1.0) analytical = 2*L^2/π^3 * (π^3/6 - π/4) - numerical = quadgk(x -> conj(IPW.ψ(x, n=n, L=L)) * x^2 * IPW.ψ(x, n=n, L=L), 0, L, maxevals=10^3)[1] + numerical = quadgk(x -> conj(ψ(IPW, x, n=n)) * x^2 * ψ(IPW, x, n=n), 0, L, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%.1f | %2d | %17.12f | %17.12f %s\n", L, n, numerical, analytical, acceptance ? "✔" : "✗") @@ -274,13 +277,14 @@ are given by the sum of 2 Taylor series: ```""") @testset "<ψₙ|p|ψₙ> = ∫ψₙ*(-iℏd/dx)ψₙdx = 0" begin - ψpψ(x; n=0, L=IPW.L, m=IPW.m, ℏ=IPW.ℏ, Δx=0.01) = -im*ℏ*conj(IPW.ψ(x,n=n,L=L))*(IPW.ψ(x+Δx,n=n,L=L)-IPW.ψ(x-Δx,n=n,L=L))/2/Δx + ψpψ(IPW, x; n=0, Δx=0.01) = -im*IPW.ℏ*conj(ψ(IPW,x,n=n))*(ψ(IPW,x+Δx,n=n)-ψ(IPW,x-Δx,n=n))/2/Δx println(" L | n | analytical | numerical ") println("--- | -- | ----------------- | ----------------- ") for L in [0.1, 0.5, 1.0, 7.0] for n in 1:1 + IPW = InfinitePotentialWell(L=L, m=1.0, ℏ=1.0) analytical = 0 - numerical = abs(quadgk(x -> ψpψ(x, n=n, L=L, m=IPW.m, ℏ=IPW.ℏ, Δx=L*0.0001), 0, L, maxevals=10^3)[1]) + numerical = abs(quadgk(x -> ψpψ(IPW, x, n=n, Δx=L*0.0001), 0, L, maxevals=10^3)[1]) acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%.1f | %2d | %17.12f | %17.12f %s\n", L, n, numerical, analytical, acceptance ? "✔" : "✗") @@ -369,13 +373,14 @@ are given by the sum of 2 Taylor series: ```""") @testset "<ψₙ|p²|ψₙ> = ∫ψₙ*(-ℏ²d²/dx²)ψₙdx = π²ℏ²/L²" begin - ψp²ψ(x; n=0, L=IPW.L, m=IPW.m, ℏ=IPW.ℏ, Δx=0.01) = -ℏ^2*conj(IPW.ψ(x,n=n,L=L))*(IPW.ψ(x+Δx,n=n,L=L)-2*IPW.ψ(x,n=n,L=L)+IPW.ψ(x-Δx,n=n,L=L))/Δx^2 + ψp²ψ(IPW, x; n=0, Δx=0.01) = -IPW.ℏ^2*conj(ψ(IPW,x,n=n))*(ψ(IPW,x+Δx,n=n)-2*ψ(IPW,x,n=n)+ψ(IPW,x-Δx,n=n))/Δx^2 println(" L | n | analytical | numerical ") println("--- | -- | ----------------- | ----------------- ") for L in [0.1, 0.5, 1.0, 7.0] for n in 1:1 + IPW = InfinitePotentialWell(L=L, m=1.0, ℏ=1.0) analytical = π^2*IPW.ℏ^2/L^2 - numerical = quadgk(x -> ψp²ψ(x, n=n, L=L, m=IPW.m, ℏ=IPW.ℏ, Δx=L*0.0001), 0, L, maxevals=10^3)[1] + numerical = quadgk(x -> ψp²ψ(IPW, x, n=n, Δx=L*0.0001), 0, L, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%.1f | %2d | %17.12f | %17.12f %s\n", L, n, numerical, analytical, acceptance ? "✔" : "✗") diff --git a/test/MorsePotential.jl b/test/MorsePotential.jl index 80ece16..35226ef 100644 --- a/test/MorsePotential.jl +++ b/test/MorsePotential.jl @@ -7,7 +7,7 @@ using Symbolics using Latexify using LaTeXStrings using SpecialFunctions -MP = antique(:MorsePotential) +MP = MorsePotential() # Lₙ⁽ᵅ⁾(x) = x⁻ᵅeˣ/n! dⁿ/dxⁿ xⁿ⁺ᵅe⁻ˣ @@ -30,13 +30,13 @@ println(raw""" for α in 0:n # Rodriguesの公式の展開 @variables x - D = n==0 ? x->x : Differential(x)^n # dⁿ/dxⁿ - a = exp(x) * x^(-α) / factorial(n) # left - b = exp(-x) * x^(n+α) # right - c = a * D(b) # Rodrigues' formula - d = expand_derivatives(c) # expand dⁿ/dxⁿ - e = simplify(d, expand=true) # simplify - f = simplify(MP.Lαint(x, n=n, α=α), expand=true) # closed-form + D = n==0 ? x->x : Differential(x)^n # dⁿ/dxⁿ + a = exp(x) * x^(-α) / factorial(n) # left + b = exp(-x) * x^(n+α) # right + c = a * D(b) # Rodrigues' formula + d = expand_derivatives(c) # expand dⁿ/dxⁿ + e = simplify(d, expand=true) # simplify + f = simplify(L(MP, x, n=n, α=α), expand=true) # closed-form # latexify eq1 = latexify(e, env=:raw) eq2 = latexify(f, env=:raw) @@ -83,7 +83,7 @@ println(raw""" for i in 0:9 for j in 0:9 analytical = gamma(i+α+1)/factorial(i)*(i == j ? 1 : 0) - numerical = quadgk(x -> real(MP.L(x, n=i, α=α)) * real(MP.L(x, n=j, α=α)) * x^α * exp(-x), 0, Inf, maxevals=10^3)[1] + numerical = quadgk(x -> real(L(MP, x, n=i, α=α)) * real(L(MP, x, n=j, α=α)) * x^α * exp(-x), 0, Inf, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%4.2f | %2d | %2d | %17.12f | %17.12f %s\n", α, i, j, analytical, numerical, acceptance ? "✔" : "✗") @@ -114,7 +114,7 @@ println(raw""" for i in 0:9 for j in 0:9 analytical = (i == j ? 1 : 0) - numerical = quadgk(x -> conj(MP.ψ(x, n=i)) * MP.ψ(x, n=j), 0, Inf, maxevals=10^3)[1] + numerical = quadgk(x -> conj(ψ(MP, x, n=i)) * ψ(MP, x, n=j), 0, Inf, maxevals=10^3)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%2d | %2d | %17.12f | %17.12f %s\n", i, j, analytical, numerical, acceptance ? "✔" : "✗") @@ -195,13 +195,14 @@ are given by the sum of 2 Taylor series: ```""") @testset "<ψₙ|H|ψₙ> = ∫ψₙ*Hψₙdx = Eₙ" begin - ψHψ(r; n=0, rₑ=MP.rₑ, Dₑ=MP.Dₑ, k=MP.k, µ=MP.µ, ℏ=MP.ℏ, Δr=0.005) = MP.V(r,rₑ=rₑ,Dₑ=Dₑ,k=k)*MP.ψ(r,n=n,rₑ=rₑ,Dₑ=Dₑ,k=k,µ=µ,ℏ=ℏ)^2 - ℏ^2/(2*μ)*conj(MP.ψ(r,n=n,rₑ=rₑ,Dₑ=Dₑ,k=k,µ=µ,ℏ=ℏ))*(MP.ψ(r+Δr,n=n,rₑ=rₑ,Dₑ=Dₑ,k=k,µ=µ,ℏ=ℏ)-2*MP.ψ(r,n=n,rₑ=rₑ,Dₑ=Dₑ,k=k,µ=µ,ℏ=ℏ)+MP.ψ(r-Δr,n=n,rₑ=rₑ,Dₑ=Dₑ,k=k,µ=µ,ℏ=ℏ))/Δr^2 + ψHψ(MP, r; n=0, Δr=0.005) = V(MP,r)*ψ(MP,r,n=n)^2 - MP.ℏ^2/(2*MP.μ)*conj(ψ(MP,r,n=n))*(ψ(MP,r+Δr,n=n)-2*ψ(MP,r,n=n)+ψ(MP,r-Δr,n=n))/Δr^2 println(" k | n | analytical | numerical ") println("--- | -- | ----------------- | ----------------- ") for k in [0.1,0.2,0.3,MP.k] for n in 0:9 - analytical = MP.E(n=n,rₑ=MP.rₑ,Dₑ=MP.Dₑ,k=k,µ=MP.µ,ℏ=MP.ℏ) - numerical = quadgk(r -> ψHψ(r,n=n, rₑ=MP.rₑ, Dₑ=MP.Dₑ, k=k, µ=MP.µ, ℏ=MP.ℏ, Δr=0.0001), 0.0001, Inf, maxevals=10^4)[1] + MP = MorsePotential(k=k) + analytical = E(MP, n=n) + numerical = quadgk(r -> ψHψ(MP, r, n=n, Δr=0.0001), 0.0001, Inf, maxevals=10^4)[1] acceptance = iszero(analytical) ? isapprox(analytical, numerical, atol=1e-5) : isapprox(analytical, numerical, rtol=1e-5) @test acceptance @printf("%.1f | %2d | %17.12f | %17.12f %s\n", k, n, analytical, numerical, acceptance ? "✔" : "✗") diff --git a/test/runtests.jl b/test/runtests.jl index 5f20a69..75e28ef 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,9 +4,8 @@ using Suppressor @testset "Antique.jl" begin @suppress_out begin - include("./InfinitePotentialWell.jl") - include("./HarmonicOscillator.jl") - include("./MorsePotential.jl") - include("./HydrogenAtom.jl") + for model in Antique.models + include("./$(model).jl") + end end end