forked from JoeMcEwen/FAST-PT
-
Notifications
You must be signed in to change notification settings - Fork 1
Expand file tree
/
Copy pathRG_RK4_example.py
More file actions
49 lines (37 loc) · 1.02 KB
/
RG_RK4_example.py
File metadata and controls
49 lines (37 loc) · 1.02 KB
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
'''
Example script to run the RK4 method for RG results.
We compare these results to the same resulst from the Cotper code.
'''
import numpy as np
from RG_RK4 import RG_RK4
# load the Copter data
d=np.loadtxt('RGcopter_1_500.dat')
k=d[:,0]; P=d[:,1]; copter=d[:,2]
# set the STS parameters here
# this combo seems to work well for k_max=10, 2000 grid points
# if you encounter instabilities you may want to fiddle with these values.
P_window=np.array([.2,.2])
C_window=.75
step=.1
max=1
n_pad=500
P_rg=RG_RK4('test_RK4',k,P,step,max,n_pad,P_window,C_window)
import matplotlib.pyplot as plt
ax=plt.subplot(121)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel(r'$k$', size=20)
ax.set_ylabel(r'$P(k)$', size=20)
ax.plot(k,P, label='linear')
ax.plot(k,copter, label='RG copter')
ax.plot(k, P_rg, label='RG FAST-PT')
plt.grid()
plt.legend(loc=3)
ax=plt.subplot(122)
ax.set_xscale('log')
ax.set_ylim(.8,1.2)
ax.set_xlabel(r'$k$', size=20)
ax.plot(k,P_rg/copter, label='ratio')
plt.grid()
plt.legend(loc=3)
plt.show()