Drawing random fields, problem on monopol and highest frequency
Jakob and I encountered a problem with .power_synthesize.
from nifty import *
import numpy as np
x1 = RGSpace(200, zerocenter=False, distances=.1)
k1 = RGRGTransformation.get_codomain(x1)
# setting spaces
npix = np.array([200]) # number of pixels
total_volume = 1. # total length
# setting signal parameters
lambda_s = .2 # signal correlation length
sigma_s = 1.5 # signal variance
# calculating parameters
k_0 = 4./(2*np.pi*lambda_s)
a_s = sigma_s**2.*lambda_s*total_volume
# creation of spaces
x1 = RGSpace(npix, distances=total_volume/npix,
zerocenter=False)
k1 = RGRGTransformation.get_codomain(x1)
p1 = PowerSpace(harmonic_partner=k1, logarithmic=False)
# # creating Power Operator with given spectrum
spec = (lambda k: a_s/(1+(k/k_0)**2)**2)
p_field = Field(p1, val=spec, dtype=np.float64)
# creating FFT-Operator and Response-Operator with Gaussian convolution
FFTOp = FFTOperator(domain=x1, target=k1,
domain_dtype=np.float64,
target_dtype=np.complex128)
# drawing a random field
sp = Field(p1, val=lambda z: spec(z) ** (1. / 2))
no_spec= 10000
sps = Field(p1, val=0.)
for ii in xrange(no_spec):
sk = sp.power_synthesize(real_signal=True, mean=0.)
sps += Field(p1, val=sk.power_analyze()**2, dtype=np.float64)
mean_power = sps/10000.
First and last mode of mean_power are wrong by a factor of 2.