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.30054329941338:0.24552774552477205:112.94276294139515, 1.562069680534925, Vector{ComplexF64}[[-4.431570101700075e12 - 1.855190613244629e7im, -4.431570101700067e12 - 1.8551906133163072e7im, -4.431570101700028e12 - 1.855190613408616e7im, -4.431557011709507e12 - 2.062794004681134e7im, -4.431557011709477e12 - 2.0627940046658434e7im, -4.431557011709469e12 - 2.062794004653199e7im, -4.431545449453727e12 - 2.06279400481596e7im, -4.4315454494537e12 - 2.0627940046884015e7im, -4.431545449453698e12 - 2.0627940046569824e7im, -4.431532359463172e12 - 1.8551906132834185e7im  …  4.431532359463154e12 - 1.8551906132080078e7im, 4.431545449453695e12 - 2.0627940047058407e7im, 4.4315454494537e12 - 2.0627940046984132e7im, 4.431545449453726e12 - 2.0627940046526838e7im, 4.431557011709482e12 - 2.06279400462265e7im, 4.431557011709485e12 - 2.0627940046517104e7im, 4.431557011709488e12 - 2.0627940047902927e7im, 4.431570101700022e12 - 1.855190613338525e7im, 4.431570101700027e12 - 1.855190613305664e7im, 4.431570101700048e12 - 1.855190613251046e7im]; [-4.431570151360895e12 - 1.8600726939056396e7im, -4.431570151360891e12 - 1.860072693896772e7im, -4.431570151360883e12 - 1.8600726938079223e7im, -4.43155702692301e12 - 2.068222409949025e7im, -4.431557026923006e12 - 2.0682224098266196e7im, -4.431557026922971e12 - 2.068222409939257e7im, -4.431545434240248e12 - 2.0682224099306107e7im, -4.431545434240228e12 - 2.068222410021921e7im, -4.431545434240207e12 - 2.068222409980899e7im, -4.431532309802341e12 - 1.8600726938720703e7im  …  4.431532309802325e12 - 1.860072693853542e7im, 4.431545434240227e12 - 2.0682224099835217e7im, 4.431545434240234e12 - 2.068222409930747e7im, 4.431545434240284e12 - 2.068222409858606e7im, 4.431557026922997e12 - 2.0682224100201726e7im, 4.431557026923009e12 - 2.068222409930248e7im, 4.431557026923038e12 - 2.0682224100760806e7im, 4.431570151360884e12 - 1.8600726938537598e7im, 4.431570151360896e12 - 1.8600726938361667e7im, 4.431570151360927e12 - 1.860072693866293e7im]; … ; [-4.431574024906275e12 - 2.240874977635449e7im, -4.431574024906273e12 - 2.240874977734375e7im, -4.431574024906251e12 - 2.2408749776245117e7im, -4.431558213575583e12 - 2.491638021390364e7im, -4.431558213575572e12 - 2.4916380215523038e7im, -4.431558213575566e12 - 2.4916380215050504e7im, -4.431544247587669e12 - 2.4916380213623047e7im, -4.431544247587646e12 - 2.49163802131958e7im, -4.431544247587616e12 - 2.4916380217285156e7im, -4.431528436257006e12 - 2.2408749776678182e7im  …  4.431528436256956e12 - 2.240874977646785e7im, 4.431544247587649e12 - 2.4916380214265347e7im, 4.431544247587661e12 - 2.491638021435379e7im, 4.431544247587685e12 - 2.4916380214425642e7im, 4.431558213575538e12 - 2.491638021465744e7im, 4.431558213575568e12 - 2.4916380214038454e7im, 4.431558213575593e12 - 2.4916380213420868e7im, 4.431574024906261e12 - 2.2408749776611328e7im, 4.431574024906271e12 - 2.240874977650365e7im, 4.431574024906272e12 - 2.2408749775563605e7im]; [-4.431574074567117e12 - 2.245757058266449e7im, -4.431574074567111e12 - 2.245757058198595e7im, -4.431574074567091e12 - 2.245757058215332e7im, -4.431558228789055e12 - 2.497066426670551e7im, -4.431558228789051e12 - 2.497066426733589e7im, -4.431558228789051e12 - 2.497066426690674e7im, -4.431544232374193e12 - 2.4970664265716553e7im, -4.431544232374187e12 - 2.4970664266750723e7im, -4.431544232374158e12 - 2.497066426679775e7im, -4.431528386596122e12 - 2.2457570583498973e7im  …  4.431528386596134e12 - 2.245757058236636e7im, 4.431544232374158e12 - 2.4970664267274607e7im, 4.431544232374161e12 - 2.4970664267339855e7im, 4.431544232374186e12 - 2.4970664267047137e7im, 4.431558228789071e12 - 2.497066426588159e7im, 4.431558228789078e12 - 2.4970664267065268e7im, 4.431558228789086e12 - 2.4970664267511394e7im, 4.431574074567089e12 - 2.245757058249668e7im, 4.431574074567089e12 - 2.245757058236754e7im, 4.4315740745671e12 - 2.2457570580828488e7im];;])

using CairoMakie

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