In [1]:
from math import exp, pi
import numpy as np
from numpy import fft, random
from scipy.integrate import quad
import cft
import gnuplot as gp
#pylab.rcParams['figure.figsize'] = (10.0, 8.0)
def EH(k):
Omega0 = 0.32 ; h = 0.69 ; Theta = 2.7255/2.7
ns = 0.96 ; A = 1122670 ; e = exp(1)
q = k * Theta**2/(Omega0*h)
L0 = np.log(2*e + 1.8*q)
C0 = 14.2 + 731.0/(1 + 62.5*q)
T0 = L0 / (L0 + C0 * q**2)
return A * k**ns * T0**2
class PowerSpectrum(cft.Filter):
def __init__(self, f):
k = cft.Power_law(1)
cft.Filter.__init__(self, lambda K: f(k(K)))
class FromAbs(cft.Filter):
def __init__(self, f):
k = cft.Power_law(1)
cft.Filter.__init__(self, lambda K: f(k(K)))
P = PowerSpectrum(EH)
def Gaussian_window(R):
def f(k):
return np.exp(-k**2 * R**2 / 2)
return f
def Tophat_window(R):
def f(k):
y = k*R
return 3*cft._K_pow(y, -3) * (np.sin(y) - y*np.cos(y))
return f
class PS_norm:
def __init__(self, P, W, dim = 3):
self.P = P
self.W = W
self.d = dim
def integrant(self):
def f(k):
return self.P(k) * self.W(k)**2 * k**(self.d-1) * 1./(2*pi**(self.d-1))
return f
def __call__(self, k_min, k_max):
f = self.integrant()
n, err = quad(f, k_min, k_max)
return np.sqrt(n)
In [115]:
N = 64 ; dim = 3
B = cft.Box(dim, N, 20.0)
data = cft.garfield(B, P)
data_f = fft.fftn(data)
power = (data_f * data_f.conj()).real
idx = np.argsort(B.k.flat)
ps = (power/B.size).flat[idx].reshape([N, N**(dim-1)]).mean(axis=1)
km = B.k.flat[idx].reshape([N, N**(dim-1)]).mean(axis=1)
S = cft.Scale(B, 8.0 / sqrt(5))
data_s = fft.ifftn(data_f * S(B.K))
def W(K):
y = (K**2).sum(axis=0)**(1./2)*8
return 3*cft._K_pow(y,-3) * (np.sin(y) - y*np.cos(y))
th0 = np.sqrt(P(B.K).mean())
th1 = np.sqrt((P * S**2)(B.K).mean())
th2 = np.sqrt(PS_norm(EH)(8.0, 0, B.k_max))
th3 = np.sqrt((P * cft.Filter(W)**2)(B.K).mean())
print(th1/data_s.std(), th0/data.std(), th1/th3, th2, th3)
#import gnuplot as gp
#np.savetxt("r.dat", np.c_[km, ps, EH(km)])
#gp.plot("set log xy", "plot 'r.dat' u 1:2 w l, '' u 1:3 w l")
In [5]:
#P = PowerSpectrum(EH)
#P = cft.Power_law(EH)
result = []
for L in np.linspace(20, 500.0, 40):
dim = 3
#W = Gaussian_window(8.0/np.sqrt(5))
#EH = lambda k: k**2
W = Tophat_window(8.0)
m = 2
B1 = cft.Box(dim, 32*m, L)
#B2 = cft.Box(dim, 48*m, L) ; B3 = cft.Box(dim, 64*m, L)
n1 = PS_norm(EH, W, dim)(B1.k_min, B1.k_max)
#n2 = PS_norm(EH, W)(B2.k_min, B2.k_max)
#n3 = PS_norm(EH, W)(B3.k_min, B3.k_max)
#n1 = n2 = n3 = 1
result.append([L,
np.sqrt((EH(B1.k) * W(B1.k)**2).mean()) / B1.res**(dim/2),
np.sqrt((EH(B1.k) * W(B1.k)**2).mean()) / B1.res**(dim/2) / n1,
n1])
#np.sqrt((EH(B2.k) * W(B2.k)).mean()) / n2 / B2.res**(dim/2),
#np.sqrt((EH(B3.k) * W(B3.k)).mean()) / n3 / B3.res**(dim/2)])
np.savetxt('r0.dat', np.array(result))
gp.plot("plot 'r0.dat' u 1:2 w l, '' u 1:3 w l, '' u 1:4 w l")
In [6]:
B = cft.Box(3, 128, 200.0)
QN = np.indices(B.shape)
def testsigma(N, L):
B = cft.Box(3, N, L)
data = cft.garfield(B, P)
Q = QN*B.res
q = np.sqrt(np.sum(((Q+B.L/2)%B.L-B.L/2)**2, axis=0))
disc = np.where(q < 8.0, 1, 0)
disc_f = fft.fftn(disc) * cft.Scale(B, B.res)(B.K)
disc_f /= disc_f[0,0]
f = fft.ifftn(fft.fftn(data) * disc_f).real
print('.', end='', file=sys.stderr, flush=True)
return f.std()
In [7]:
result2 = np.array([[L, testsigma(128, L)] for L in np.linspace(100.0, 1000.0, 20)])
In [11]:
plot(result2[:,0], result2[:,1]*result2[:,0]**2) #/(result2[4:,0]/128)**(1./2))
Out[11]:
In [126]:
result2
Out[126]:
In [115]:
def testnorm(N, L):
k_min = 2*pi/L ; k_max = N*pi/L
return sqrt(PS_norm(EH)(8, k_min, k_max))
result3 = np.array([[L, testnorm(128, L)] for L in np.linspace(100.0, 1000.0, 20)])
In [116]:
plot(result3[:,0], result3[:,1])
Out[116]:
In [42]:
import thulsa
f = open('../c.density.init.conan', 'rb')
h,_ = thulsa.read_array(f)
_,d = thulsa.read_array(f)
L = eval(h['size'])
N = eval(h['N'])
dim = eval(h['dim'])
shape = (N,)*dim
B = cft.Box(dim, N, L)
d = d.reshape(shape)
smd = fft.ifftn(fft.fftn(d) * Gaussian_window(8.0/np.sqrt(5))(B.k)).real
print(smd.std())
In [ ]: