In [1]:
%config InlineBackend.figure_format = 'retina'
import numpy as np
import matplotlib.pyplot as plt
import math as mt
In [2]:
def dispcold(rm, wce, th, wmin, wmax, nomega, kmax):
rm=1/rm
if th==90: th=89.9999
rd=th*mt.pi/180
sin2=mt.sin(rd)**2
cos2=mt.cos(rd)**2
omsave=np.empty((4,nomega*2))
isol=0
for iw in range(nomega):
om=wmin+(wmax-wmin)*iw/nomega
rr=1-rm/om/(om+rm*wce)-1/om/(om-wce)
ll=1-rm/om/(om-rm*wce)-1/om/(om+wce)
pp=1-(rm+1)/om/om
ss=(rr+ll)/2; dd=(rr-ll)/2
aa=ss*sin2+pp*cos2
bb=(ss**2-dd**2)*sin2+pp*ss*(1+cos2)
cc=pp*rr*ll
ff=bb**2-4*aa*cc
nnp=(bb+mt.sqrt(ff))/2/aa if ff>0 else -1
nnm=(bb-mt.sqrt(ff))/2/aa if ff>0 else -1
if nnp>0:
isol=isol+1
ak=mt.sqrt(nnp)*om; pol=(nnp-ss)/dd
ez=-nnp*mt.sqrt(cos2*sin2)/(pp-nnp*sin2); ey2=1/pol**2
esem=(mt.sqrt(sin2)+mt.sqrt(cos2)*ez)/mt.sqrt(1+ey2+ez**2)
omsave[0,isol]=ak
omsave[1,isol]=om
omsave[2,isol]=pol
omsave[3,isol]=(esem/mt.sqrt(1-esem**2))
#omsave[3,isol]=(esem/mt.sqrt(1-esem**2))
if nnm>0:
isol=isol+1
ak=mt.sqrt(nnm)*om; pol=(nnm-ss)/dd
ez=-nnm*mt.sqrt(cos2*sin2)/(pp-nnm*sin2); ey2=1/pol**2
esem=(mt.sqrt(sin2)+mt.sqrt(cos2)*ez)/mt.sqrt(1+ey2+ez**2)
omsave[0,isol]=ak
omsave[1,isol]=om
omsave[2,isol]=pol
omsave[3,isol]=(esem/mt.sqrt(1-esem**2))
plt.figure(figsize=(12,4))
plt.subplots_adjust(wspace=0.4, hspace=0.6)
plt.subplot(1,3,1)
plt.plot(omsave[0,:],omsave[1,:],'.')
if th==0:
plt.plot(omsave[0,:], np.zeros(nomega*2)+np.sqrt(1+rm), 'b-')
plt.axis([0, kmax, wmin, wmax])
plt.xlabel(r'$kc/\omega_{pe}$')
plt.ylabel(r'$\omega/\omega_{pe}$')
plt.subplot(1,3,2)
plt.plot(omsave[2,:],omsave[1,:],'.')
plt.title(r'$\longleftarrow\quad L\; |\; R \quad\longrightarrow$')
plt.axis([-10,10,wmin,wmax])
plt.xlabel(r'$iE_x/E_y$')
plt.subplot(1,3,3)
plt.plot(np.log10(omsave[3,:]),omsave[1,:],'.')
plt.title(r'$EM\quad\longleftrightarrow\quad ES$')
plt.axis([-5,5,wmin,wmax])
plt.xlabel(r'$Log\left|ES/EM\right|$')
plt.show()
In [3]:
rm=16; wce=2; th=0
wmin=0.01; wmax=10*mt.pi/8; nomega=2000
kmax=10*mt.pi/8
dispcold(rm, wce, th, wmin, wmax, nomega, kmax)
In [4]:
rm=16; wce=2; th=45
wmin=0.01; wmax=10*mt.pi/8; nomega=2000
kmax=10*mt.pi/8
dispcold(rm, wce, th, wmin, wmax, nomega, kmax)
In [5]:
rm=16; wce=2; th=90
wmin=0.01; wmax=10*mt.pi/8; nomega=2000
kmax=10*mt.pi/8
dispcold(rm, wce, th, wmin, wmax, nomega, kmax)
