Here is the algorithm shown in the paper,
reference: Zenitani2015PoP
Here I try to reproduce the algorithm suggested in Zenitani2015PoP to load relativistic particles in PIC simulations.
%config InlineBackend.figure_format='retina'
from spacepy import pycdf
import math as mt
import matplotlib.pyplot as plt
import numpy as np
import scipy.fft as sf
import datetime
import bisect
import matplotlib.patches as patches
import matplotlib.dates as mdates
from IPython.display import display,Math
import matplotlib as mpl
from ai import cdas
import datetime
import pandas as pd
from matplotlib.gridspec import GridSpec
import matplotlib.animation as animation
from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, mark_inset
from matplotlib.ticker import ScalarFormatter
color='#2c2c2c'
#color='#fffff0'
mpl.rcParams['font.family' ]='Arial'
mpl.rcParams['mathtext.default' ]='rm'
mpl.rcParams['font.size' ]=14
mpl.rcParams['figure.dpi' ]=120
mpl.rcParams['axes.linewidth' ]=2
mpl.rcParams['xtick.major.size' ]=5
mpl.rcParams['xtick.minor.size' ]=4
mpl.rcParams['ytick.major.size' ]=5
mpl.rcParams['ytick.minor.size' ]=3
mpl.rcParams['xtick.major.width']=1
mpl.rcParams['ytick.major.width']=1
mpl.rcParams['xtick.direction' ]='out'
mpl.rcParams['ytick.direction' ]='out'
mpl.rcParams['axes.facecolor' ]='None'
mpl.rcParams['figure.facecolor' ]='None'
mpl.rcParams['axes.edgecolor' ]=color
mpl.rcParams['xtick.color' ]=color
mpl.rcParams['ytick.color' ]=color
mpl.rcParams['axes.labelcolor' ]=color
mpl.rcParams['text.color' ]=color
mpl.rcParams['patch.edgecolor' ]=color
mpl.rcParams['savefig.bbox' ]='tight'
mpl.rcParams['animation.html' ]='jshtml'
###color scheme for lines: https://www.nature.com/articles/nmeth.1618 ###
from cycler import cycler
line_cycler = (cycler(color =[color, '#e69f00', '#56b4e9', '#009e73', '#f0e442', '#0072b2', '#d55e00', '#cc79a7']))
#line_cycler = (cycler(color =[color, color, color, color])+
# cycler(linestyle=["-" , "--" , "-." , ":" ]))
mpl.rcParams['axes.prop_cycle']=line_cycler
#########################################################################
Here is the algorithm shown in the paper,
def RM(Gamma0):
nop=200000
T=1
up=np.zeros((nop,3))
#Gamma0=10
beta0 =np.sqrt(1-1/Gamma0**2)
for ip in range(nop):
uu =0
eta=0
while eta**2-uu**2 < 1:
x1,x2,x3,x4=np.random.rand(4)
uu =-T*np.log(x1*x2*x3)
eta=-T*np.log(x1*x2*x3*x4)
x5,x6,x7=np.random.rand(3)
up[ip,0]= uu*(2*x5-1)
up[ip,1]=2*uu*np.sqrt(x5*(1-x5))*np.cos(2*np.pi*x6)
up[ip,2]=2*uu*np.sqrt(x5*(1-x5))*np.sin(2*np.pi*x6)
gam=np.sqrt(1+uu**2)
vx=up[ip,0]/gam
if(-beta0*vx>x7): up[ip,0]=-up[ip,0]
up[ip,0]=Gamma0*(up[ip,0]+beta0*np.sqrt(1+uu**2))
return up
An analitical form of a boosted relativistic Maxwell distributrion can be written, $$ f\left( {u_x^\prime } \right) = \int_0^\infty {\int_0^{2\pi } {{f^\prime }\left( {{{\mathbf{u}}^\prime }} \right)u_ \bot ^\prime d\phi du_ \bot ^\prime = \frac{{N\left(\Gamma \sqrt {1 + u_x^{\prime 2}} + T\right)}}{{2{\Gamma ^2}{K_2}\left( {1/T} \right)}}\exp \left( { - \frac{{\Gamma \left( {\sqrt {1 + u_x^{\prime 2}} - \beta u_x^\prime } \right)}}{T}} \right),} } $$ where $K_2(x)$ is the modified Bessel function of the second kind.
### three cases for loading relativistic particles ###
up1=RM(1 )
up2=RM(1.1)
up3=RM(10)
from scipy.special import kn
ux=np.linspace(-200,200,20000)
T=1
G1=1
G2=1.1
G3=10
beta1=np.sqrt(1-1/G1**2)
beta2=np.sqrt(1-1/G2**2)
beta3=np.sqrt(1-1/G3**2)
f1=(G1*np.sqrt(1+ux**2)+T)/(2*G1**3*kn(2,1/T))*np.exp(-G1*(np.sqrt(1+ux**2)-beta1*ux)/T)
f2=(G2*np.sqrt(1+ux**2)+T)/(2*G2**3*kn(2,1/T))*np.exp(-G2*(np.sqrt(1+ux**2)-beta2*ux)/T)
f3=(G3*np.sqrt(1+ux**2)+T)/(2*G3**3*kn(2,1/T))*np.exp(-G3*(np.sqrt(1+ux**2)-beta3*ux)/T)
#plt.hist(up1[:,0],bins=200,edgecolor='C0',alpha=0.1,density=True)
#plt.hist(up2[:,0],bins=200,edgecolor='black',alpha=0.25,density=True)
fp1,binx1=np.histogram(up1[:,0],bins= 20,density=True)
fp2,binx2=np.histogram(up2[:,0],bins= 20,density=True)
fp3,binx3=np.histogram(up3[:,0],bins=100,density=True)
xx1=0.5*(binx1[1:]+binx1[:-1])
xx2=0.5*(binx2[1:]+binx2[:-1])
xx3=0.5*(binx3[1:]+binx3[:-1])
plt.plot(xx1,fp1,'C0s',zorder=5,label=r'$\Gamma=1$',markersize=5)
plt.plot(xx2,fp2,'o',zorder=5,label=r'$\Gamma=1.1$', markersize=5, markeredgewidth=1, markeredgecolor='C0',color='none')
plt.plot(xx3,fp3,'C0x',zorder=5,label=r'$\Gamma=10$',markersize=5)
plt.plot(ux,f1,'C0-',lw=0.5)
plt.plot(ux,f2,'C0-',lw=0.5)
plt.plot(ux,f3,'C0-',lw=0.5)
plt.yscale('log')
plt.xlim(-20,50)
plt.ylim(1e-5,1e0)
plt.legend()
plt.xlabel(r'$u_x^\prime$')
plt.ylabel(r'$\frac{f^\prime(u_x^\prime)}{\Gamma N}$')
plt.show()
