Initial commit

convergence.py

@@ -0,0 +1,61 @@
+#!/usr/bin/env python3
+import numpy as np
+import time
+import os
+import sys
+import colorsys
+Nx = int(sys.argv[1])
+Ny = int(sys.argv[2])
+N = max(Nx, Ny)
+
+def CpxRandNormal(min, max, N):
+    mu = (min+max)/2
+    sigma = (max-min)/2
+    return np.random.normal(mu, sigma, (N, N))+1.0j*np.random.normal(mu, sigma, (N, N))
+def CpxRand(min, max, N):
+    return np.random.uniform(min, max, (N, N))+1.0j*np.random.uniform(min, max, (N, N))
+
+A = CpxRandNormal(-100, 100, N)
+B = CpxRand(-1, 1, N)*0.1
+C = CpxRand(-1, 1, N)
+
+A_ex = np.linalg.inv(np.eye(N)-B)@C
+
+dir = -1
+solve = True
+t = 0
+while True:
+    t += 0.003
+    if solve:
+        A = 0.9*A+0.1*(B@A+C)
+    else:
+        A += CpxRandNormal(0, 0.008, N)
+    diff = A-A_ex
+    #angles = (np.angle(diff)+np.pi)/(2*np.pi
+    #diff = np.abs(diff)
+    #diff = 1-diff/(1+diff)
+    a = np.abs(np.real(diff))
+    a = a/(a+1)
+    b = np.abs(np.imag(diff))
+    b = b/(b+1)
+    img = np.zeros((Ny, Nx, 3))
+    for x in range(Nx):
+        for y in range(Ny):
+            c = colorsys.hsv_to_rgb((t+0.2+a[x, y]*0.2)%1, 1, 0.5*b[x, y]+0.5)
+            img[y,x,:] = np.array(c)*255
+    #img[:,:,0] = np.minimum(np.abs(diff*255).astype(int), 255)[:Ny,:Nx]
+    out = img.reshape((Nx*Ny*3,)).astype(np.uint8)
+    #A += np.random.uniform(-0.01, 0.01, (N, N))
+
+    #print(len(out))
+    os.write(1, out.tobytes())
+
+    if solve and np.sum(np.abs(diff)) <= 0.001*N*N:
+        solve = False
+    elif not solve and np.sum(np.abs(diff)) >= 1*N*N:
+        solve = True
+
+    #os.write(2, b"frame")
+    #print(angles, diff)
+
+    time.sleep(0.01)

pendlum.py

@@ -0,0 +1,81 @@
+#!/usr/bin/env python3
+import numpy as np
+import sys
+import os
+import time
+import colorsys
+
+Nx = int(sys.argv[1])
+Ny = int(sys.argv[2])
+N = 100
+iterations = 2
+try:
+    N = int(sys.argv[3])
+except:
+    pass
+
+def rk4(f, x, dt):
+    k1 = f(x)
+    k2 = f(x+k1*dt/2)
+    k3 = f(x+k2*dt/2)
+    k4 = f(x+k3*dt)
+    return x+dt*(k1+2*k2+2*k3+k4)/6.0
+dt = 0.08 
+m = 1
+l = 1
+g = 1
+
+def sq(x): return x*x
+
+def rhs(x):
+    phi1 = x[0]
+    p1 = x[1]
+    phi2 = x[2]
+    p2 = x[3]
+
+    dphi1 = 6/(m*l*l)*(2*p1-3*np.cos(phi1-phi2)*p2)/(16-9*sq(np.cos(phi1-phi2))) 
+    dphi2 = 6/(m*l*l)*(8*p2-3*np.cos(phi1-phi2)*p1)/(16-9*sq(np.cos(phi1-phi2)))
+
+    dp1 = -0.5*m*l*l*(dphi1*dphi2*np.sin(phi1-phi2)+3*g/l*np.sin(phi1))
+    dp2 = -0.5*m*l*l*(-dphi1*dphi2*np.sin(phi1-phi2)+g/l*np.sin(phi2))
+    return np.array([dphi1, dp1, dphi2, dp2])
+
+s = np.array([np.pi+0.1, 0, np.pi-0.1, 0])
+
+angles = np.linspace(np.pi, np.pi/3, N)+np.random.uniform(-np.pi/6/N, np.pi/6/N, N)
+angles = np.random.uniform(-np.pi/6, np.pi/3, N)+np.pi
+states = [np.array([np.pi, 0, a, 0]) for a in angles]
+
+def clamp(x, min_, max_):
+    return  max(min_, min(max_, x))
+
+buffer = bytearray(b"\x00"*Nx*Ny*3)
+
+def setpixel(x, y, r, g, b):
+    xi = int(Nx*(x+1.2)/2.4)
+    yi = int(Ny*(y+1.2)/2.4)
+    if xi < 0 or xi >= Nx or yi < 0 or yi >= Ny:
+        return
+    idx = xi+Nx*yi
+    buffer[3*idx+0] = r
+    buffer[3*idx+1] = g
+    buffer[3*idx+2] = b
+timestep = 0
+timefac = 0.03
+while True:
+    timestep += 1
+    for s, i in zip(states, range(len(states))):
+        phi1 = s[0]
+        phi2 = s[2]
+        x1 = np.sin(phi1)*l
+        y1 = np.cos(phi1)*l
+        x2 = np.sin(phi2)*l+x1
+        y2 = np.cos(phi2)*l+y1
+        h = 0.2*i/N+0*timestep*timefac/iterations
+        r, g, b = colorsys.hsv_to_rgb(h%1, 1, 1)
+        setpixel(x2, y2, int(r*255), int(g*255), int(b*255))
+        states[i] = rk4(rhs, s, dt/iterations)
+    if timestep > 10*iterations and timestep % iterations == 0:
+        #time.sleep(0.01)
+        os.write(1, buffer)
+

quadratic.py

@@ -0,0 +1,80 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import os
+import sys
+import time
+import colorsys
+
+Nx = int(sys.argv[1])
+Ny = int(sys.argv[2])
+iterations = 10
+
+class Slider:
+    def __init__(self, lo, hi, step, pos=0):
+        self.lo = lo
+        self.hi = hi
+        self.stepSize = step
+        self.pos = 0
+        self.dir = 1
+    
+    def step(self):
+        if self.pos >= self.hi:
+            self.dir = -1
+        elif self.pos <= self.lo:
+            self.dir = 1
+        
+        self.pos += self.dir * self.stepSize
+        return self.pos
+
+Sa = Slider(-10, 10, 0.1/iterations, 0) #1.1
+Sb = Slider(-10, 10, 0.2/iterations, 1) #-2.1
+Sc = Slider(-10, 10, 0.2/iterations, 2) #-0.33
+Sd = Slider(-10, 10, 0.1/iterations, 3) #1.7
+Se = Slider(-10, 10, 0.2/iterations, 4)
+Sf = Slider(-10, 10, 0.01/iterations, 5)
+
+x, y = np.meshgrid(
+    np.linspace(-4, 4, Nx),
+    np.linspace(-4, 4, Ny))
+
+color = 1.0
+
+img = np.zeros( [Ny, Nx, 3] )
+
+cr = 1.0
+cg = 0.5
+cb = 0.25
+
+step = 0
+while True:
+    a = Sa.step()
+    b = Sb.step()
+    c = Sc.step()
+    d = Sd.step()
+    e = Se.step()
+    f = Sf.step()
+
+    curve = (np.abs(a*x**2 + b*x*y + c*y**2 + d*x + e*y + f) <= 1.5)*1
+
+    cr, cg, cb = colorsys.hsv_to_rgb(color, 1, 1)
+
+
+    s = np.shape(curve)
+
+    #img[:,:,0] = curve
+    img[:,:,0] = np.where(curve > 0, curve*cr*255, img[:,:,0])
+    img[:,:,1] = np.where(curve > 0, curve*cg*255, img[:,:,1])
+    img[:,:,2] = np.where(curve > 0, curve*cb*255, img[:,:,2])
+
+    step += 1
+    if step % iterations == 0:
+        out = img.reshape((Nx*Ny*3,)).astype(np.uint8)
+        os.write(1, out.tobytes())
+        time.sleep(0.01)
+
+    color = (color+0.01/iterations ) % 1
+
+    cb = (cb + 0.05) % 1
+
+

swifthohenberg.py

@@ -0,0 +1,63 @@
+#!/usr/bin/env python3
+import sys
+import os
+import numpy as np
+Nx = int(sys.argv[1])
+Ny = int(sys.argv[2])
+param = sys.argv[3]
+
+buffer = bytearray(b" "*(3*Nx*Ny))
+
+
+Lx = 32 #128 #Physikalische Laenge des Gebietes in x-Richtung
+Ly =32 #128#Physikalische Laenge des Gebietes in y-Richtung
+
+x, y = np.meshgrid(np.arange(Nx)* Lx/Nx, np.arange(Ny)* Ly/Ny) #x-Array
+kx, ky = np.meshgrid(np.fft.fftfreq(Nx,Lx/(Nx*2.0*np.pi)), np.fft.fftfreq(Ny,Ly/(Ny*2.0*np.pi)))
+ksq = kx*kx + ky*ky
+
+c = 0.0+(np.random.random((Ny,Nx))-0.5)*0.1
+eps=0.3
+delta=0.0
+#eps = 0.1
+#delta = 1.0
+
+t = 0.0
+dt = 0.01
+T_End = 300000 
+N_t = int(T_End / dt)
+
+plotEveryNth = 100
+ck = np.fft.fft2(c) #FFT(c)
+# Lineare Terme
+def rhs_lin(ksq):
+    result=eps-(1.0-ksq)**2
+    return result
+
+Eu=1.0/(1.0-dt*rhs_lin(ksq))
+i = 0
+def lerp_sat(a, t, b):
+    v = (1-t)*a+t*b
+    v = int(v)
+    if v < 0:
+        v = 0
+    if v > 255:
+        v = 255
+    return v
+
+while True:
+    i+= 1
+    ck=Eu*(ck+dt*(delta*np.fft.fft2(np.fft.ifft2(ck)**2)-np.fft.fft2(np.fft.ifft2(ck)**3)))
+    c=np.fft.ifft2(ck)
+    eps = 0.1+0.2*np.cos(i/10000)
+    delta = np.sin(i/10000)
+    if(i % plotEveryNth == 0):
+        myc = c.real
+        myc = (myc-np.min(myc))/(np.max(myc)-np.min(myc))
+        for px in range(Nx):
+            for py in range(Ny):
+                idx = 3*(px+Nx*py)
+                buffer[idx+0] = lerp_sat(0xff, myc[py,px], 0x00)
+                buffer[idx+1] = lerp_sat(0x00, myc[py,px], 0xff)
+                buffer[idx+2] = 0
+        os.write(1, buffer)