DataVisTutorial/gen_fractal.py at main · AstroBrandt/DataVisTutorial · GitHub

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
import matplotlib.pyplot as plt
import scipy.fft as sfft
import numpy as np
import numpy.random as rnd
rnd.seed()

def gen_fractal(N, D):
    D = 2.3
    beta = 2*(4. - D)
    ndim= int(N)
    N2  = int(ndim/2)

    data = np.zeros((ndim, ndim, ndim), dtype=complex)
    fndim = float(ndim)

    ix = 0
    while ix < N2:
        iy = 0
        while iy < ndim:
            iz = 0
            while iz < ndim:
                phi = rnd.uniform(0, 2.0*np.pi)
                fix = min(float(ix), float(ndim-ix))
                fiy = min(float(iy), float(ndim-iy))
                fiz = min(float(iz), float(ndim-iz))
                freq = np.sqrt(fix*fix + fiy*fiy + fiz*fiz)
                amp = np.sqrt(float(1.0/(freq/ndim)**beta))

                data[ix,iy,iz] = complex(amp*np.cos(phi), amp*np.sin(phi))
                iz += 1
            iy += 1
        ix += 1
    ix = 0
    while ix < N2:
        iy = 0
        while iy < ndim:
            iz = 0
            while iz < ndim:
                data[(ndim-ix)%ndim,(ndim-iy)%ndim,(ndim-iz)%ndim] = np.conj(data[ix,iy,iz])
                iz += 1
            iy += 1
        ix += 1
    data[0,0,0]=complex(0.0,0.0)
    data[0,0,N2]=complex(abs(data[0,0,N2]),0.0)
    data[0,N2,0]=complex(abs(data[0,N2,0]),0.0)
    data[N2,0,0]=complex(abs(data[N2,0,0]),0.0)
    data[0,N2,N2]=complex(abs(data[0,N2,N2]),0.0)
    data[N2,0,N2]=complex(abs(data[N2,0,N2]),0.0)
    data[N2,N2,0]=complex(abs(data[N2,N2,0]),0.0)
    data[N2,N2,N2]=complex(abs(data[N2,N2,N2]),0.0)

    ffdata= sfft.fftn(data)
    ffdata= sfft.fftshift(ffdata)

    var=np.var(ffdata)
    ffdata=ffdata/np.sqrt(var)

    fbmcube=np.sqrt(ffdata**2).astype('float64')

    h = np.std(fbmcube.flatten())
    dens= np.exp(fbmcube/h)

    (xind, yind, zind) = np.unravel_index(dens.argmax(), dens.shape)
    xroll = int(xind - N2)
    yroll = int(yind - N2)
    zroll = int(zind - N2)
    dens = np.roll(dens, -xroll, axis=0)
    dens = np.roll(dens, -yroll, axis=1)
    dens = np.roll(dens, -zroll, axis=2)

    return dens