# Usage

Here is a sample REPL session to draw a dispersion relationship by using this package.

First, construct FiniteSquareWell potential object.

julia> using ExtendedKronigPennyMatrixjulia> v0=10;julia> rho=0.5 # b/a;julia> pot=FiniteSquareWell(v0, rho)FiniteSquareWell(10, 0.5)

Use get_function method to acquire potential function.

julia> begin
pf = get_potential(pot)
a = 1
xs=-a:a/100:2a
plot(xs, pf.(xs), "k")
xlim(0,1)
xlabel(L"$x / a$")
ylabel(L"Energy / $E_0$")
title( L"$\rho =$"*string(rho))
end

Define Model by calling its constructor.

julia> Ka=0.0 # wavenumber multiplied by a;julia> model=Model(pot, Ka)Model{FiniteSquareWell}(FiniteSquareWell(10, 0.5), 0.0, 60, 121, Alternates(60), [5.0 3.183098861837907 … -0.0 -0.0; 3.183098861837907 9.0 … -0.05395082816674419 0.05218194855471979; … ; -0.0 -0.05395082816674419 … 14405.0 -0.0; -0.0 0.05218194855471979 … -0.0 14405.0])

The field hnm of model contains Hamiltonian matrix.

julia> typeof(model.hnm)Matrix{Float64} (alias for Array{Float64, 2})julia> size(model.hnm)(121, 121)julia> model.hnm[1:5,1:5]5×5 Matrix{Float64}:
5.0      3.1831    3.1831   -0.0      -0.0
3.1831   9.0      -0.0       3.1831   -1.06103
3.1831  -0.0       9.0      -1.06103   3.1831
-0.0      3.1831   -1.06103  21.0      -0.0
-0.0     -1.06103   3.1831   -0.0      21.0

Use LinearAlgebra.eigvals method to compute its energy eigenvalues. Refer to the LinearAlgebra standard library section in Julia documentation.

julia> using LinearAlgebrajulia> evs=eigvals(model.hnm);julia> evs[1:3]3-element Vector{Float64}:
1.9680655058115561
7.582164159224384
11.500508757008639

Draw dispersion curve by scanning Ka values between $[-\pi, \pi]$.

julia> using PyPlotjulia> clf()julia> begin
a = 1
xs=-a:a/100:2a
plot(xs .- 1/2, pf.(xs), "k")  # Holizontally shift to centerize the potential well
cm=get_cmap("tab10")
for Ka in (-18:18)/18*π
model=Model(pot, Ka)
ev = eigvals(model.hnm)
for i in 1:5
plot(Ka/ π, ev[i], ".", color=cm(i-1))
end
end
xlim(-1,1)
ylim(-2,32)
xlabel(L"$Ka / \pi$")
ylabel(L"Energy / $E_0$")
title( L"$\rho =$"*string(rho))
end