pixelserver2/apps/swifthohenberg.py

64 lines
1.5 KiB
Python
Raw Normal View History

2024-10-28 17:10:03 +00:00
#!/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)