Subphase rotation of 2D sinewave

# (c) Computational Alchemy
# Works on Python 2.7
power_x = 10
power_y = 10
steps = 2000

filename = "trig-1"

import math, cmath, numpy
import colorsys
from PIL import Image
import os.path

def BaseSelect(base, n):
  accumulator = 1
  bit = 0
  while n > 0:
    if (n & 1) == 1:
      accumulator = accumulator * base[bit]
    n = n >> 1
    bit = bit + 1
  return accumulator

#make palette
palette = [0] * 256
for i in range(256):
  c = colorsys.hsv_to_rgb(math.pi*i/256, 1.0, 1.0)
  palette[i] = (int(c[0]*256), int(c[1]*256), int(c[2]*256))

print "Calculate number sequencies..."

W = 2**power_x
H = 2**power_y
base_x = [cmath.exp(2j*math.pi*2**(i+1)/2**power_x) for i in range(power_x)]
base_y = [cmath.exp(2j*math.pi*2**(i+1)/2**power_y) for i in range(power_y)]
print "base_x=", base_x
print "base_y=", base_y
n_x = [BaseSelect(base_x, i) for i in range(2**power_x)]
n_y = [BaseSelect(base_y, i) for i in range(2**power_y)]

print "Find boundaries..."

max_int = 1.e+10
min_int = -max_int
real_min = max_int
real_max = min_int
imag_min = max_int
imag_max = min_int
mag_min = max_int
mag_max = min_int
arg_min = max_int
arg_max = min_int
for y in range(H):
  for x in range(W):
    n = (n_x[x] * n_y[y])
    real_min = min(real_min, n.real)
    real_max = max(real_max, n.real)
    imag_min = min(imag_min, n.imag)
    imag_max = max(imag_max, n.imag)
    r, phi = cmath.polar(n)
    mag_min = min(mag_min, r)
    mag_max = max(mag_max, r)
    arg_min = min(arg_min, phi)
    arg_max = max(arg_max, phi)


print "Generating images..."
templates = [
  (real_min, real_max, lambda n: n.real, "real"),
  #(imag_min, imag_max, lambda n: n.imag, "imag"),
  #(mag_min, mag_max, lambda n: cmath.polar(n)[0], "mag"),
  #(arg_min, arg_max, lambda n: cmath.polar(n)[1], "arg"),
]

angle = cmath.exp(2j*math.pi/steps)

for step in range(steps):
  print "*** Step", step
  n_x = [BaseSelect(base_x, i) for i in range(2**power_x)]
  n_y = [BaseSelect(base_y, i) for i in range(2**power_y)]
  for template in templates:
    min_n, max_n, num_part, name = template
    magnitude = max_n - min_n
    print "Generating image", name, "..."
    print "min_n=", min_n
    print "max_n=", max_n
    print "magnitude=", magnitude

    im = Image.new("RGB", (W, H), (0,0,0))
    pix = im.load()

    for y in range(H):
      for x in range(W):
        pix[x, y] = palette[int(255*(num_part(n_x[x] * n_y[y])-min_n)/(magnitude))]

    dstfilename, _ = os.path.splitext(filename)
    dstfilename = "%s-%s-%04d.png" % (dstfilename, name, step)
    print "Saving image to", dstfilename
    im.save(dstfilename, "PNG")

  base_x = [p*angle for p in base_x]
  base_y = [p*angle for p in base_y]