Case: Ion cyclotron emission

Ion cyclotron emission (ICE) driven by a ring ion beam distribution in a magnetized fusion device

using PlasmaBO
using PlasmaBO: gyrofrequency, Alfven_speed

using DelimitedFiles

# Magnetic field (Tesla)
B0 = 2.1

# Wave vector angle between k and B0
θ = deg2rad(89.5)

# Species table taken from the original BO MATLAB input
fpath = pkgdir(PlasmaBO, "test/ice_Irvine18.in")
readdlm(fpath)

4×8 Matrix{Any}:
   "qs(e)"   "ms(m_unit)"   "ns(m^-3)"  …   "vdsz/c"   "vdsr/c"   "ifv(1/0)"
  2         4.0            1.0e16          0.0        0.045      1
  1         2.0            9.98e18         0.0        0.0        1
 -1         0.0005447      1.0e19          0.0        0.0        1

Dispersion Curve Scan

tbl = readdlm(fpath, Float64; skipstart = 1)
species = map(eachrow(tbl)) do row
    Z, A, n_s, Tz_s, Tp_s, vdz_s, vdr_s = row[1:7]
    Maxwellian(n_s, Tz_s, Tp_s; vdz = vdz_s, vdr = vdr_s, Z , A)
end

ωn = gyrofrequency(B0, species[1])
vA = Alfven_speed(B0, species)
kn = ωn / vA

ks = (9.5:0.025:11.5) .* kn
results = solve(species, B0, ks, θ; N = 12, J = 4)

PlasmaBO.DispersionSolution{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, Float64, Matrix{Vector{ComplexF64}}}(93.30054330562844:0.24552774554112747:112.94276294891863, 1.562069680534925, Vector{ComplexF64}[[-4.431570101700072e12 - 1.855190613505257e7im, -4.43157010170007e12 - 1.855190613483269e7im, -4.431570101700049e12 - 1.8551906134210803e7im, -4.431557011709517e12 - 2.0627940043406587e7im, -4.431557011709508e12 - 2.062794004823205e7im, -4.43155701170948e12 - 2.062794004800417e7im, -4.431545449453785e12 - 2.0627940049560543e7im, -4.431545449453759e12 - 2.0627940047954217e7im, -4.431545449453696e12 - 2.0627940047713965e7im, -4.431532359463228e12 - 1.8551906135404333e7im  …  4.431532359463167e12 - 1.8551906134819284e7im, 4.431545449453702e12 - 2.0627940048085608e7im, 4.431545449453737e12 - 2.0627940048861165e7im, 4.431545449453751e12 - 2.0627940047313385e7im, 4.431557011709488e12 - 2.0627940048462257e7im, 4.431557011709488e12 - 2.0627940047668282e7im, 4.431557011709514e12 - 2.0627940044725567e7im, 4.431570101700056e12 - 1.855190613461667e7im, 4.431570101700069e12 - 1.8551906135385606e7im, 4.431570101700076e12 - 1.8551906133371823e7im]; [-4.431570151360947e12 - 1.860072693951416e7im, -4.431570151360918e12 - 1.860072694021328e7im, -4.431570151360914e12 - 1.860072693996179e7im, -4.431557026923019e12 - 2.0682224098560512e7im, -4.431557026923e12 - 2.0682224101250887e7im, -4.43155702692299e12 - 2.068222410185884e7im, -4.431545434240289e12 - 2.068222410107302e7im, -4.431545434240237e12 - 2.0682224100585938e7im, -4.431545434240227e12 - 2.0682224101379395e7im, -4.431532309802328e12 - 1.8600726939331055e7im  …  4.431532309802332e12 - 1.8600726939029694e7im, 4.431545434240213e12 - 2.0682224101455733e7im, 4.431545434240218e12 - 2.068222410039459e7im, 4.431545434240226e12 - 2.0682224099885337e7im, 4.431557026922984e12 - 2.068222410035515e7im, 4.431557026922995e12 - 2.0682224100404568e7im, 4.431557026923003e12 - 2.0682224101623535e7im, 4.431570151360892e12 - 1.8600726939086914e7im, 4.431570151360902e12 - 1.8600726938975785e7im, 4.43157015136092e12 - 1.860072694118377e7im]; … ; [-4.43157402490629e12 - 2.2408749780855328e7im, -4.431574024906268e12 - 2.2408749777956687e7im, -4.431574024906243e12 - 2.2408749777813967e7im, -4.431558213575577e12 - 2.4916380217041016e7im, -4.431558213575571e12 - 2.4916380215209354e7im, -4.431558213575533e12 - 2.49163802158575e7im, -4.431544247587658e12 - 2.491638021546048e7im, -4.431544247587645e12 - 2.49163802159729e7im, -4.43154424758764e12 - 2.4916380214701682e7im, -4.431528436256948e12 - 2.240874977709961e7im  …  4.43152843625696e12 - 2.240874977716065e7im, 4.431544247587606e12 - 2.491638021528017e7im, 4.431544247587629e12 - 2.4916380215296175e7im, 4.431544247587639e12 - 2.4916380215747047e7im, 4.431558213575538e12 - 2.4916380216951743e7im, 4.431558213575568e12 - 2.4916380215396117e7im, 4.431558213575569e12 - 2.491638021446719e7im, 4.431574024906242e12 - 2.24087497779541e7im, 4.431574024906246e12 - 2.2408749777442932e7im, 4.431574024906265e12 - 2.240874977716155e7im]; [-4.43157407456712e12 - 2.2457570583129883e7im, -4.43157407456709e12 - 2.2457570585205078e7im, -4.43157407456708e12 - 2.2457570583776765e7im, -4.431558228789059e12 - 2.4970664268961743e7im, -4.431558228789037e12 - 2.497066427054882e7im, -4.431558228789034e12 - 2.4970664268464573e7im, -4.431544232374177e12 - 2.4970664270355225e7im, -4.431544232374165e12 - 2.497066427085352e7im, -4.431544232374159e12 - 2.4970664268833887e7im, -4.43152838659614e12 - 2.2457570582315147e7im  …  4.431528386596135e12 - 2.245757058544922e7im, 4.431544232374152e12 - 2.4970664269206695e7im, 4.431544232374165e12 - 2.497066426679902e7im, 4.431544232374198e12 - 2.4970664270345785e7im, 4.431558228789025e12 - 2.4970664269531265e7im, 4.43155822878905e12 - 2.4970664269476555e7im, 4.431558228789074e12 - 2.4970664270992957e7im, 4.431574074567082e12 - 2.2457570584942393e7im, 4.431574074567094e12 - 2.2457570583251953e7im, 4.431574074567154e12 - 2.2457570582265157e7im];;])

using CairoMakie

f, (ax, ax2) = plot(results, kn, ωn)
ylims!(ax, [0, 15])
ylims!(ax2, [-0.5, 0.5])
f