Usage

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

First, construct FiniteSquareWell potential object.

julia> using ExtendedKronigPennyMatrix
julia> 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 LinearAlgebra
julia> 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 PyPlot
julia> 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