# -*- coding: utf-8 -*-
from benchpress.benchmarks import util
import numpy as np
bench = util.Benchmark("Quasi Crystal simulation", "k*stripes*image_size*iterations")
def main():
k = bench.args.size[0] # number of plane waves
stripes = bench.args.size[1] # number of stripes per wave
N = bench.args.size[2] # image size in pixels
ite = bench.args.size[3] # iterations
phases = np.arange(0, 2 * np.pi, 2 * np.pi / ite)
image = np.empty((N, N), dtype=bench.dtype)
d = np.arange(-N / 2, N / 2, dtype=bench.dtype)
xv, yv = np.meshgrid(d, d)
theta = np.arctan2(yv, xv)
r = np.log(np.sqrt(xv * xv + yv * yv))
r[np.isinf(r) == True] = 0
tcos = theta * np.cos(np.arange(0, np.pi, np.pi / k))[:, np.newaxis, np.newaxis]
rsin = r * np.sin(np.arange(0, np.pi, np.pi / k))[:, np.newaxis, np.newaxis]
inner = (tcos - rsin) * stripes
cinner = np.cos(inner)
sinner = np.sin(inner)
bench.start()
for phase in phases:
image[:] = np.sum(cinner * np.cos(phase) - sinner * np.sin(phase), axis=0) + k
bench.flush()
bench.stop()
bench.pprint()
if __name__ == "__main__":
main()