-
Notifications
You must be signed in to change notification settings - Fork 0
/
pusherFP.py
157 lines (121 loc) · 4.36 KB
/
pusherFP.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
# numpy
import numpy as np
np.random.seed(19680801)
from numpy.random import default_rng
rng = default_rng()
# import functions
from scipy.special import kv, iv, erf
from scipy.integrate import quad
from numpy import log, log10, sin, cos, exp, sqrt, pi
# interpolate
from scipy import interpolate
# physical constants
from scipy.constants import speed_of_light, fine_structure, hbar, elementary_charge, electron_mass
# root finding
from scipy.optimize import fsolve
from scipy import optimize
# plotting
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
# progress bar
from tqdm.notebook import tqdm
from tqdm import trange
from time import sleep
from tqdm import tqdm
# warnings
import warnings
warnings.simplefilter('ignore')
import pandas as pd
def arraycenter(x):
"""
returns centered array for histograms
"""
return [(x[i]+x[i+1])/2 for i in range(len(x)-1)]
"""
%-------------------------------------------------------------------------------
% FP pusher
%-------------------------------------------------------------------------------
"""
def g25(chi):
"""Niel2018 equation 25"""
res = 9*sqrt(3)/(8*pi) * quad(lambda v: 2*v**2 * kv(5/3,v)/(2+3*v*chi)**2 + 4*v* (3*v*chi)**2 *kv(2/3,v)/(2+3*v*chi)**4, 0, np.inf)[0]
return res
def h41(chi):
"""Niel2018 equation 41"""
res = 9*sqrt(3)/(4*pi) * quad(lambda v: 2*chi**3 * v**3 * kv(5/3,v)/(2+3*v*chi)**3 + 54* chi**5 * v**4 *kv(2/3,v)/(2+3*v*chi)**5, 0, np.inf)[0]
return res
def S39(chi,Kalpha):
"""Niel2018 equation 39"""
res = 2/3 * Kalpha * chi**2 * g25(chi); #fine_structure**2/taue
return res
def R40(chi,gm,Kalpha):
"""Niel2018 equation 40"""
res = 2/3 * Kalpha * gm * h41(chi);
return res
def pushFP(chi0, Nsmpl, tdim, nbins, g0, s0):
gmdist = g0 + s0*rng.standard_normal(Nsmpl)
gmdist[gmdist<1]=1 # make sure no particle has a physically inconsistent gamma
gmdist_dump1 = np.copy(gmdist) #[] dump1
gmdist_dump2 = np.copy(gmdist) * 0 #[] dump2
gmdist_dump3 = np.copy(gmdist) * 0 #[] dump3
# reference chi0
print("chi0 =",chi0)
#chi0 = 1e-1; #[] {1e-3,1e-2,1e-1,1e0}
if chi0==1e-3:
tmax = 20
elif chi0==1e-2:
tmax = 20
elif chi0==1e-1:
tmax = 5
elif chi0==1e0:
tmax = 3
else:
tmax = 5
#tmax = 5 #[1/omega_p] simulation duration {20,20,5,3}
# reference omega_p
omega_p = elementary_charge/electron_mass/1800 * 2.5e3 * (chi0/(1e-3))
#print("omega_p =",omega_p)
# dimensionless but physically relevant constant
Kalpha = fine_structure * electron_mass * speed_of_light**2 / (hbar * omega_p)
#print("Kalpha = alpha^2/tau_e/omega_p =",Kalpha)
# simulation parameters
dt = tmax/tdim #[1/omega_p]
print("dt =", dt)
# interpolation
gmlst = np.linspace(1,2*g0+3*s0,300); #[]
S_lst = np.array([S39(gm/1800*chi0,Kalpha) for gm in gmlst])
S_intrp = interpolate.interp1d(gmlst, S_lst) #interp1d / CubicSpline
R_lst = np.array([R40(gm/1800*chi0,gm,Kalpha) for gm in gmlst])
R_intrp = interpolate.interp1d(gmlst, R_lst)
mom1 = np.zeros(tdim); mom2 = np.zeros_like(mom1)
# for each time step
for t in tqdm(range(tdim)):
# interpolate S and R
S = S_intrp(gmdist);
R = R_intrp(gmdist);
# Gaussian random number
dW = sqrt(dt) * np.random.randn(Nsmpl)
# Niel2018 eq 42
dgamma = -S * dt + sqrt(R) * dW
gmdist = gmdist + dgamma
gmdist[gmdist < 1] = 1
gmdist[gmdist > 2*g0+3*s0] = g0+3*s0
# save average and spread
mom1[t] = np.mean(gmdist)
mom2[t] = np.std(gmdist)
# save distribution
if t == int(tdim/2):
gmdist_dump2 = np.copy(gmdist)
# save distribution
gmdist_dump3 = np.copy(gmdist)
# get histograms
# dump 1
gmdist_y,gmdist_x = np.histogram(gmdist_dump1,np.linspace(1,g0+4*s0,nbins))
gmdist1_y, gmdist1_x = gmdist_y, np.array(arraycenter(gmdist_x))
# dump 2
gmdist_y,gmdist_x = np.histogram(gmdist_dump2,np.linspace(1,g0+4*s0,nbins))
gmdist2_y, gmdist2_x = gmdist_y, np.array(arraycenter(gmdist_x))
# dump 3
gmdist_y,gmdist_x = np.histogram(gmdist_dump3,np.linspace(1,g0+4*s0,nbins))
gmdist3_y, gmdist3_x = gmdist_y, np.array(arraycenter(gmdist_x))
return gmdist1_y, gmdist1_x, gmdist2_y, gmdist2_x, gmdist3_y, gmdist3_x, mom1, mom2