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.431570101700084e12 - 1.8551906133728024e7im, -4.4315701017000625e12 - 1.8551906134632695e7im, -4.431570101700047e12 - 1.8551906134349592e7im, -4.431557011709488e12 - 2.0627940049980577e7im, -4.431557011709475e12 - 2.062794004801128e7im, -4.431557011709474e12 - 2.0627940046806168e7im, -4.431545449453711e12 - 2.062794004937744e7im, -4.431545449453699e12 - 2.062794004786873e7im, -4.431545449453693e12 - 2.0627940047851562e7im, -4.431532359463206e12 - 1.855190613403436e7im  …  4.431532359463177e12 - 1.8551906135454748e7im, 4.431545449453714e12 - 2.062794004731323e7im, 4.431545449453722e12 - 2.06279400480085e7im, 4.431545449453743e12 - 2.062794005144938e7im, 4.431557011709484e12 - 2.0627940048331134e7im, 4.431557011709514e12 - 2.0627940047593117e7im, 4.431557011709526e12 - 2.0627940047838457e7im, 4.431570101700048e12 - 1.8551906135876447e7im, 4.431570101700061e12 - 1.8551906132385306e7im, 4.431570101700078e12 - 1.8551906135009766e7im]; [-4.431570151360942e12 - 1.860072693884575e7im, -4.43157015136092e12 - 1.8600726942454338e7im, -4.431570151360891e12 - 1.860072693988037e7im, -4.431557026923013e12 - 2.0682224099746305e7im, -4.431557026922984e12 - 2.0682224100723267e7im, -4.431557026922978e12 - 2.0682224101334333e7im, -4.431545434240246e12 - 2.068222410522461e7im, -4.431545434240231e12 - 2.0682224099172235e7im, -4.431545434240198e12 - 2.0682224100013733e7im, -4.431532309802333e12 - 1.8600726938232422e7im  …  4.431532309802347e12 - 1.8600726940387975e7im, 4.431545434240202e12 - 2.0682224101266354e7im, 4.431545434240221e12 - 2.0682224100937128e7im, 4.431545434240268e12 - 2.0682224098876953e7im, 4.431557026922972e12 - 2.068222410039997e7im, 4.431557026922979e12 - 2.068222410003985e7im, 4.431557026922996e12 - 2.068222410129693e7im, 4.431570151360881e12 - 1.860072694066275e7im, 4.431570151360894e12 - 1.8600726937062528e7im, 4.431570151360899e12 - 1.8600726939208984e7im]; … ; [-4.431574024906279e12 - 2.2408749778845437e7im, -4.43157402490627e12 - 2.24087497779541e7im, -4.431574024906255e12 - 2.240874977748241e7im, -4.431558213575569e12 - 2.4916380216080494e7im, -4.431558213575565e12 - 2.4916380216213316e7im, -4.431558213575546e12 - 2.491638021591405e7im, -4.431544247587649e12 - 2.4916380215430498e7im, -4.431544247587641e12 - 2.4916380215280555e7im, -4.431544247587635e12 - 2.4916380215271004e7im, -4.431528436256959e12 - 2.2408749778832797e7im  …  4.431528436256947e12 - 2.2408749778839e7im, 4.431544247587639e12 - 2.4916380215863112e7im, 4.431544247587664e12 - 2.4916380216311708e7im, 4.431544247587672e12 - 2.4916380215579025e7im, 4.431558213575551e12 - 2.4916380215834618e7im, 4.431558213575554e12 - 2.491638021595826e7im, 4.431558213575606e12 - 2.491638021504797e7im, 4.431574024906241e12 - 2.2408749778198242e7im, 4.43157402490625e12 - 2.2408749777162496e7im, 4.431574024906267e12 - 2.240874977734375e7im]; [-4.431574074567125e12 - 2.2457570584960938e7im, -4.431574074567098e12 - 2.2457570585989952e7im, -4.431574074567092e12 - 2.2457570582763672e7im, -4.431558228789078e12 - 2.497066426870724e7im, -4.431558228789065e12 - 2.4970664269293785e7im, -4.431558228789053e12 - 2.4970664266500782e7im, -4.431544232374218e12 - 2.4970664270507812e7im, -4.431544232374169e12 - 2.497066426835312e7im, -4.431544232374166e12 - 2.4970664266418457e7im, -4.431528386596139e12 - 2.2457570583017662e7im  …  4.431528386596122e12 - 2.2457570583129883e7im, 4.431544232374141e12 - 2.4970664268513232e7im, 4.431544232374152e12 - 2.497066426852417e7im, 4.431544232374186e12 - 2.4970664267468702e7im, 4.431558228789058e12 - 2.4970664269531276e7im, 4.431558228789061e12 - 2.497066426789093e7im, 4.43155822878912e12 - 2.497066426641944e7im, 4.4315740745671e12 - 2.2457570584716797e7im, 4.431574074567111e12 - 2.245757058544922e7im, 4.431574074567111e12 - 2.24575705826416e7im];;])

using CairoMakie

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