Page MenuHomec4science
Files F60152920

surface_1d.py
No OneTemporary
Actions

File Metadata

Created
Sat, Apr 27, 21:55

surface_1d.py
View Options

#!/usr/bin/env python
# @file
# @section LICENSE
#
# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
from __future__ import division
import tamaas as tm
from numpy.fft import fft, fftshift, fftfreq
import matplotlib.pyplot as plt
tm.initialize()
def draw_triangle(pos, w, ax, style, alpha=1.5, base=10):
h = alpha * w
ax.loglog([base**pos[0], base**(pos[0]+w)],
[base**pos[1], base**pos[1]], style)
ax.loglog([base**pos[0], base**pos[0]],
[base**pos[1], base**(pos[1]+h)], style)
ax.loglog([base**(pos[0]+w), base**pos[0]],
[base**pos[1], base**(pos[1]+h)], style)
# ax.text(base**(pos[0]+w/3), base**(pos[1]+h/3), "$-{}$".format(alpha),
# horizontalalignment='center',
# verticalalignment='center')
ax.text(base**(pos[0]+w/2), base**(pos[1]-0.2), '$1$',
horizontalalignment='center',
verticalalignment='top')
ax.text(base**(pos[0]-0.15), base**(pos[1]+h/2), '${}$'.format(alpha),
horizontalalignment='right',
verticalalignment='center')
iso = tm.Isopowerlaw1D()
helper = tm.ParamHelper(iso)
params = {
"Q0": 16,
"Q1": 16,
"Q2": 512,
"Hurst": 0.8
}
helper.set(params)
gen = tm.SurfaceGeneratorFilter1D()
gen.setFilter(iso)
gen.setSizes([2048])
# gen.setRandomSeed(0) # uncomment for fixed seed
surf = gen.buildSurface()
plt.plot(surf)
surf_fft = fft(surf)
psd = surf_fft*surf_fft.conj()
psd = fftshift(psd)
plt.figure()
freqs = fftshift(fftfreq(psd.size, d=1/psd.size))
plt.plot(freqs[psd.size/2+1:], psd.real[psd.size/2+1:])
draw_triangle([5, -4], 2, plt.gca(), '-k', 2*(params["Hurst"]+1), 2)
plt.yscale('log')
plt.xscale('log', basex=2)
plt.ylim([10**(-2), 10**8])
plt.show()
tm.finalize()

Event Timeline

c4science · Help