Reaction Difusion Surface

Just had an amazing stroke of luck!!! found a working Reaction Diffusion code from an internet site
translating it from java to python was almost too easy

the site is:
http://cgjennings.ca/toybox/turingmorph/

and here are the results (on a plane for now)

and a closeup

and here is the python code, its pretty straightforward

enjoy

import math
import rhinoscriptsyntax as rs
import random as rnd
  
def InitMatrix (rndSeed):
    global A0
    global B0
    rnd.seed(rndSeed)
    for i in range(height):
        thisA =[]
        thisB =[]
        for j in range(width):
            thisA.append(rnd.random()*12+rnd.random()*2)
            thisB.append(rnd.random()*12+rnd.random()*2)
        A0.append(thisA)
        B0.append(thisB)

def CopyMatrix(src):
    res=[]
    for i in range(len(src)):
        list =[]
        for j in range(len(src[0])):
            list.append(src[i][j])
        res.append(list)
    return res

def NullMatrix():
    res=[]
    for i in range(height):
        list =[]
        for j in range(width):
            list.append(0)
        res.append(list)
    return res


def Solve(CA,CB,iterations) :
    An = NullMatrix()
    Bn = NullMatrix()
    global A0
    global B0

   
    ###uses Euler's method to solve the diff eqns
    for n in range(iterations):
        for i in range(height):
            ### treat the surface as a torus by wrapping at the edges
            iplus1 = i+1
            iminus1 = i-1
            if( i == 0 ) : iminus1 = height - 1
            if( i == height - 1 ) : iplus1 = 0
            for j in range (width):
                jplus1 = j+1
                jminus1 = j-1
                if( j == 0 ) : jminus1 = width - 1
                if( j == width - 1 ) : jplus1 = 0
                ###Component A
                DiA = CA * ( A0[iplus1][j] - 2.0 * A0[i][j] + A0[iminus1][j] + A0[i][jplus1] - 2.0 * A0[i][j] + A0[i][jminus1] )
                ReA = A0[i][j] * B0[i][j] - A0[i][j] - 12.0
                An[i][j] = A0[i][j] + 0.01 * (ReA + DiA)
                if( An[i][j] < 0.0 ) : An[i][j] = 0.0
                ###Component B
                DiB = CB * ( B0[iplus1][j] - 2.0 * B0[i][j] + B0[iminus1][j]   + B0[i][jplus1] - 2.0 * B0[i][j] + B0[i][jminus1] )
                ReB = 16.0 - A0[i][j] * B0[i][j]
                Bn[i][j] = B0[i][j] + 0.01 * (ReB + DiB)
                if( Bn[i][j] < 0.0 ) : Bn[i][j]=0.0
               
               
        A0=CopyMatrix(An)
        B0=CopyMatrix(Bn)

def ShowLines(which,spacing,Z):
    mat = []
    if which=='A' : mat=A0
    if which=='B' : mat=B0
   
    rs.EnableRedraw(False)
    for i in range(height):
        curve =[]
        for j in range(width):
            curve.append((i*spacing,j*spacing,mat[i][j]*Z))
        rs.AddCurve(curve)
    rs.EnableRedraw()

###main call

width = 200
height = 200
seed = 1
A0 =[]
B0 = []



InitMatrix(seed)
Solve(1,16,500)
ShowLines('A',3,0.5)