In [27]:

    

import gnuplot as gp
import numpy as np
from numpy import random, fft
from functools import reduce, partial
from math import pi

def _wave_number(s):
    """compute wave numbers for fourier transforms"""
	N = s[0]
	i = np.indices(s)
	return np.where(i > N/2, i - N, i)

def _scale_filter(B, t):  
	"""returns discrete scale space filter, take care with units: 
		[res] = Mpc / pixel, [k] = 1 / Mpc, [t] = Mpc**2"""
	def f(K):
		return reduce(
			lambda x, y: x*y, 
			(np.exp(t / B.res**2 * (np.cos(k * B.res) - 1)) for k in K))
	return f

class Box:
    """store parameters for the box"""
	def __init__(self, dim, N, L):
		self.N = N
		self.L = L
		self.res = L/N
		self.dim = dim
		
		self.shape = (self.N,) * dim
		self.size = reduce(lambda x, y: x*y, self.shape)
	
		self.K = _wave_number(self.shape) * 2*pi / self.L
		self.k = np.sqrt((self.K**2).sum(axis=0))

		self.k_max = N*pi/L
		self.k_min = 2*pi/L

def _K_pow(k, n):
	"""raise |k| to the <n>-th power safely"""
	save = np.seterr(divide = 'ignore')
	a = np.where(k == 0, 0, k**n)
	np.seterr(**save)
	return a
        
B = Box(2, 512, 100.0)
# create some white noise
wn = random.normal(0, 1, B.shape)
# apply a power-spectrum
f = fft.ifftn(fft.fftn(wn) * np.sqrt(_K_pow(B.k, -2)) * _scale_filter(B, 0.1)(B.K)).real

gp.plot("unset key", "set size square", 
        "set xrange [-{0}:{0}] ; set yrange [-{0}:{0}]".format(B.L/2),
        gp.plot_data(gp.array(f, 
            "dx={0} dy={0} origin=(-{1}+{0}/2,-{1}+{0}/2) with image".format(
            B.L/B.N, B.L/2))))


    

In [39]:

    

from IPython.display import display, Image
import os


    

In [43]:

    

x = np.linspace(-B.L/2, B.L/2, B.N)
y = f[0]

g = gp.Gnuplot()
g("set term pdf font 'Bitstream Charter'",
  "set output 'tmp.pdf'",
  "set key bottom right box",
  gp.plot_data(gp.record(np.c_[x, y], "t'hoi' with lines")))
g.terminate()

os.system("convert tmp.pdf tmp.png")
display(Image(filename="tmp.png"))


    

In [ ]: