#### Simulation of a reaction-diffusion equation.
#### Here is the 1D version (2D too slow for VPython)
## This is the Gray-Scott equation, where 2 chemicals
## diffuse and interact. The interaction is nonlinear.
##
## Eqn:
##    d/dt U = D_u d^2/dx^2 U - U V^2 + F(1-U)
##    d/dt V = D_u d^2/dx^2 V + U V^2 - (F + k) V
##
## See web refs., e.g.
##
## The goal here is to try different values of the parameter
## F and k.

from visual import *
from random import *


print "Simulation of Gray-Scott reaction-diffusion equation in one"
print "dimension. The chemical concentrations of U and V are updated"
print "by spreading their concentrations (diffusion) and by reactions"
print "between the chemicals. The chemical U is constantly fed in, while"
print "both U and V decay, at rates set by the parameters F and k."
print "The concentration of V is indicated by the displacement of the"
print "blocks, while color indicates the concentration of U."
F = input("Choose a value for F (suggested around 0.03): ")
k = input("Choose a value for k (suggested around 0.05 to 0.06): ")

#### parameters: size of array of boxes, strength of
#### interactions, time step, and damping
N = 100
dt = 1
dx = 0.02
F = 0.025
k = 0.055
U = ones((N),Float)
V = zeros((N),Float)
delU = zeros((N),Float)
delV = zeros((N),Float)
DU = .00002/dx/dx
DV = .00001/dx/dx
for i in range(N/4,3*N/4):
        U[i] = random()*0.2+0.4
        V[i] = random()*0.2+0.2

##### function to take a number and turn into a color
##### zero gives dark magenta; large positive is red, large
##### negative is blue
def httocolor(z):   
    return (sqrt(sqrt(z)),0,(1-z))

# updateUV arguments: U, V arrays, DU,DV diffusion constants,
#                     parameters F and k.
def updateUV(U, V, DU, DV, F, k):
    for i in range(0,N):
        lf = (i+N-1)%N
        rt = (i+1)%N
        delU[i] = U[rt] + U[lf]
        delV[i] = V[rt] + V[lf]
    U += (DU * (delU - 2 * U) - U * V * V + F * (1. - U)) * dt
    V += (DV * (delV - 2 * V) + U * V * V - (F + k) * V) * dt

##### set up scene
scene.autoscale = 0
scene.center = (N/2,0,0)
scene.range = N*.75
scene.forward = (0,0,-1)
#scene.up = (0,0,1)


#### set up a list of lists of boxes, with initial velocity (zero)
bl = []
for i in range(N):
    bl.append(box(pos=(i,0,0)))


#### Go on simulation
timelabel = label(pos=(N,N/3,0))
count = 0
while 1:
    count += 1
    updateUV(U,V,DU,DV,F,k)
    timelabel.text = "Time = " + str(count)
    for i in range(N):
        bl[i].y = U[i]*2
        bl[i].color = httocolor(V[i])