Get the dispersion relation for slab geometry

We will here calculate the real and imaginary part of the dispersion relation given in Pécseli, H -Low Frequency Waves and Turbulence in Magnetized Laboratory Plasmas and in the Ionosphere, 2016



In [ ]:

    
from sympy import init_printing
from sympy import Eq, I
from sympy import re, im
from sympy import symbols
from sympy.solvers import solve
from IPython.display import display
from sympy import latex

om     = symbols('omega')
omI    = symbols('omega_i', real=True)
omStar = symbols('omega_S', real=True)

sigmaPar = symbols('sigma', positive=True)

b = symbols('b', real=True)

init_printing()

The pure drift wave

We are here checking if the dispersion relation in 6.4.1 - Pure drift-wave is easy to work with



In [ ]:

    
LHS = om*(om-omI)+I*sigmaPar*(om-omStar + b*(om-omI))
RHS = 0
eq = Eq(LHS, RHS)



In [ ]:

    
display(eq)



In [ ]:

    
sol1, sol2 = solve(eq, om)

display(sol1)
display(sol2)

The difference in the two solutions are just the sign of the square-root



In [ ]:

    
sol1Re = re(sol1)
sol1Im = im(sol1)

display(sol1Re)
display(sol1Im)

This is cumbersome to work with. Using something easier.

Resistive drift waves with $T_i=0$

We are here checking if the dispersion relation in 5.5 - Dispersion relation, is easy to work with.



In [ ]:

    
LHS = om**2 + I*sigmaPar*(om*(1+b)-omStar)
RHS = 0
eq = Eq(LHS, RHS)



In [ ]:

    
display(eq)



In [ ]:

    
sol1, sol2 = solve(eq, om)

display(sol1)
display(sol2)

The difference in the two solutions are just the sign of the square-root. The first solution gives the largest growth rate



In [ ]:

    
sol1Re = re(sol1.expand())
sol1Im = im(sol1.expand())

real = Eq(symbols("I"),sol1Im.simplify())
imag = Eq(symbols("R"),sol1Re.simplify())
display(real)
display(imag)

print(latex(real))
print(latex(imag))



In [ ]:

    
sol2Re = re(sol2.expand())
sol2Im = im(sol2.expand())

display(Eq(symbols("I"),sol2Im.simplify()))
display(Eq(symbols("R"),sol2Re.simplify()))

This also gives quite the mess...

Splitting to $\omega_R$ and $\omega_I$



In [ ]:

    
# NOTE: Do not confuse om_I with om_i
om_I = symbols('omega_I', real=True)
om_R = symbols('omega_R', real=True)



In [ ]:

    
LHSSplit = LHS.subs(om, om_R + I*om_I)
display(Eq(LHS, re(LHSSplit)+I*im(LHSSplit)))