In [14]:

    

import numpy as np
import cft
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import collections
from scipy.spatial import ConvexHull

from adhesion import get_convex_hull, plot_regular_triangulation, plot_power_diagram


    

In [15]:

    

plt.rcParams['figure.figsize'] = (15, 10)
plt.rcParams['lines.solid_capstyle'] = 'round'

def plot_image(box, data, title):
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.set_title(title)
    ax.imshow(data, extent=[-box.L/2, box.L/2, -box.L/2, box.L/2], interpolation='nearest', cmap='RdYlBu')
    plt.show()


    

In [16]:

    

box = cft.Box(dim=2, N=32, L=25.)
P = cft.Power_law(n=-0.75) * cft.Scale(box, 0.5)
delta_0 = cft.garfield(B=box, P=P, seed=44)
delta_0 /= delta_0.std()
pot_0 = np.fft.ifftn(np.fft.fftn(delta_0) * cft.Potential()(box.K)).real

ch, selection, valid = get_convex_hull(box, pot_0, 2.0)

xlim = None
ylim = None

plot_regular_triangulation(ch, selection, xlim, ylim)
plt.show()


    

In [17]:

    

plot_power_diagram(box, ch, valid, xlim, ylim, point_scale=10)
plt.show()


    

In [22]:

    

xlim = None #[0, 15]
ylim = None #[-10, 5]

box = cft.Box(dim=2, N=32, L=25.)
P = cft.Power_law(n=-0.75) * cft.Scale(box, 0.5)
delta_0 = cft.garfield(B=box, P=P, seed=44)
delta_0 /= delta_0.std()
pot_0 = np.fft.ifftn(np.fft.fftn(delta_0) * cft.Potential()(box.K)).real

plt.rcParams['figure.figsize'] = (7, 10)
fig, axes = plt.subplots(3, 2) # , sharex=True, sharey=True)
for ax in axes.flat:
    ax.get_xaxis().set_ticks([])
    ax.get_yaxis().set_ticks([])
ch, selection, valid = get_convex_hull(box, pot_0, 0.0001)
plot_regular_triangulation(ch, selection, xlim, ylim, ax=axes[0,0])
plot_power_diagram(box, ch, selection, xlim, ylim, ax=axes[0,1], point_scale=20)
ch, selection, valid = get_convex_hull(box, pot_0, 0.75)
plot_regular_triangulation(ch, selection, xlim, ylim, ax=axes[1,0])
plot_power_diagram(box, ch, selection, xlim, ylim, ax=axes[1,1], point_scale=20)
ch, selection, valid = get_convex_hull(box, pot_0, 1.5)
plot_regular_triangulation(ch, selection, xlim, ylim, ax=axes[2,0])
plot_power_diagram(box, ch, selection, xlim, ylim, ax=axes[2,1], point_scale=20)
fig.tight_layout()
axes[0,0].set_ylabel('D=0')
axes[1,0].set_ylabel('D=0.5')
axes[2,0].set_ylabel('D=1.0')
plt.show()


    

In [23]:

    

fig.savefig('time-evolution.pdf', bbox_inches='tight')


    

In [19]:

    

from numpy import fft

xlim = [-80, 80]
ylim = [-80, 80]

box = cft.Box(dim=2, N=128, L=200.)
seed = 43

fig, axes = plt.subplots(3, 2, sharex=True, sharey=True, figsize=(7,10))
#for ax in axes.flat:
#    ax.get_xaxis().set_ticks([])
#    ax.get_yaxis().set_ticks([])
def sigma_8(box, delta_0):
    f = fft.ifftn(fft.fftn(delta_0) * cft.Scale(box, 8.0)(box.K)).real
    return f.std()

P = cft.Power_law(n=1) # * cft.Scale(box, 1.0)
delta_0 = cft.garfield(B=box, P=P, seed=seed)
delta_0 /= sigma_8(box, delta_0)
pot_0 = np.fft.ifftn(np.fft.fftn(delta_0) * cft.Potential()(box.K)).real
ch, selection, valid = get_convex_hull(box, pot_0, 0.7)
plot_regular_triangulation(ch, selection, xlim, ylim, ax=axes[0,0])
plot_power_diagram(box, ch, valid, xlim, ylim, ax=axes[0,1], point_scale=2, line_scale=0.3, plot_grid=False)

P = cft.Power_law(n=0) # * cft.Scale(box, 1.0)
delta_0 = cft.garfield(B=box, P=P, seed=seed)
delta_0 /= sigma_8(box, delta_0)
pot_0 = np.fft.ifftn(np.fft.fftn(delta_0) * cft.Potential()(box.K)).real
ch, selection, valid = get_convex_hull(box, pot_0, 0.7)
plot_regular_triangulation(ch, selection, xlim, ylim, ax=axes[1,0])
plot_power_diagram(box, ch, valid, xlim, ylim, ax=axes[1,1], point_scale=2, line_scale=0.3, plot_grid=False)

P = cft.Power_law(n=-1) # * cft.Scale(box, 1.0)
delta_0 = cft.garfield(B=box, P=P, seed=seed)
delta_0 /= sigma_8(box, delta_0)
pot_0 = np.fft.ifftn(np.fft.fftn(delta_0) * cft.Potential()(box.K)).real
ch, selection, valid = get_convex_hull(box, pot_0, 0.7)
plot_regular_triangulation(ch, selection, xlim, ylim, ax=axes[2,0])
plot_power_diagram(box, ch, valid, xlim, ylim, ax=axes[2,1], point_scale=2, line_scale=0.3, plot_grid=False)
fig.tight_layout()
axes[0,0].set_ylabel('n=1')
axes[1,0].set_ylabel('n=0')
axes[2,0].set_ylabel('n=-1')

# plt.savefig('ps-index-evolution.pdf', bbox_inches='tight')
plt.show()


    

In [13]:

    

fig.savefig('ps-index-evolution.pdf', bbox_inches='tight')


    

In [ ]: