A = matrix([[1,1],[-1,1]])
D = [vector([0,0]), vector([1,0])]
@interact
def f(A = matrix([[1,1],[-1,1]]), D = '[[0,0],[1,0]]', k=(3..17)):
print "Det = ", A.det()
D = matrix(eval(D)).rows()
def Dn(k):
ans = []
for d in Tuples(D, k):
s = sum(A^n*d[n] for n in range(k))
ans.append(s)
return ans
G = points([v.list() for v in Dn(k)])
show(G, frame=True, axes=False)
Click to the left again to hide and once more to show the dynamic interactive window |
A = matrix([[0,0,2],[1,0,1],[0,1,-1]])
D = '[[0,0,0],[1,0,0]]'
def Dn(D,A,k):
ans = []
for d in Tuples(D, k):
s = sum(A^n*d[n] for n in range(k))
ans.append(s)
return ans
@interact
def f(A = matrix([[0,0,2],[1,0,1],[0,1,-1]]), D = '[[0,0,0],[1,0,0]]', k=(3..15), labels=True):
print "Det = ", A.det()
D = matrix(eval(D)).rows()
print "D:"
print D
G = point3d([v.list() for v in Dn(D,A,k)], size=8)#, opacity=.85)
if labels:
G += sum([text3d(str(v),v) for v in Dn(D,A,k)])
show(G, axes=False, frame=False)
Click to the left again to hide and once more to show the dynamic interactive window |
%cython
import numpy as np
cimport numpy as np
def mandelbrot_cython(float x0,float x1,float y0,float y1,
int N=200, int L=50, float R=3):
'''returns an array NxN to be plotted with matrix_plot
'''
cdef double complex c, z, I
cdef float deltax, deltay, R2 = R*R
cdef int h, j, k
cdef np.ndarray[np.uint16_t, ndim=2] m
m = np.zeros((N,N), dtype=np.uint16)
I = complex(0,1)
deltax = (x1-x0)/N
deltay = (y1-y0)/N
for j in range(N):
for k in range(N):
c = (x0+j*deltax)+ I*(y0+k*deltay)
z=0
h=0
while (h<L and
z.real**2 + z.imag**2 < R2):
z=z*z+c
h+=1
m[j,k]=h
return m
import pylab
x0_default = -2
y0_default = -1.5
side_default = 3.0
side = side_default
x0 = x0_default
y0 = y0_default
options = ['Reset','Upper Left', 'Upper Right', 'Stay', 'Lower Left', 'Lower Right']
@interact
def show_mandelbrot(option = selector(options, nrows = 2, width=8),
N = slider(100, 1000,100, 300),
L = slider(20, 300, 20, 60),
plot_size = slider(2,10,1,6),
auto_update = False):
global x0, y0, side
if option == 'Lower Right':
x0 += side/2
y0 += side/2
elif option == 'Upper Right':
y0 += side/2
elif option == 'Lower Left':
x0 += side/2
if option=='Reset':
side = side_default
x0 = x0_default
y0 = y0_default
elif option != 'Stay':
side = side/2
time m=mandelbrot_cython(x0 ,x0 + side ,y0 ,y0 + side , N, L )
# p = (matrix_plot(m) +
# line2d([(N/2,0),(N/2,N)], color='red', zorder=2) +
# line2d([(0,N/2),(N,N/2)], color='red', zorder=2))
# time show(p, figsize = (plot_size, plot_size))
pylab.clf()
pylab.imshow(m, cmap = pylab.cm.gray)
time pylab.savefig('mandelbrot.png')
Click to the left again to hide and once more to show the dynamic interactive window |
@interact
def mandel_plot(expo = slider(-10,10,0.1,2), \
formula = ['mandel','ff'],\
iterations=slider(1,100,1,30), \
zoom_x = range_slider(-2,2,0.01,(-2,1)), \
zoom_y = range_slider(-2,2,0.01,(-1.5,1.5))):
var('z c')
f(z,c) = z^expo + c
ff_m = fast_callable(f, vars=[z,c], domain=CDF)
# messing around with fast_callable
for i in range(int(iterations)/3):
f(z,c) = f(z,c)^expo+c
ff = fast_callable(f, vars=[z,c], domain=CDF)
def mandel(z):
c = z
for i in range(iterations):
z = ff_m(z,c)
if abs(z) > 2:
return z
return z
print 'z <- z^%s + c' % expo
# calling ff three times, otherwise it fast_callable exceeds a recursion limit
if formula is 'ff':
func = lambda z: ff(ff(ff(z,z),z),z)
elif formula is 'mandel':
func = mandel
complex_plot(func, zoom_x,zoom_y, plot_points=200, dpi=150).show(frame=True, aspect_ratio=1)
Click to the left again to hide and once more to show the dynamic interactive window |
%python
from numpy import zeros
def sierpinski(N):
'''Generates the Sierpiński Triangle fractal to N iterations using the Rule 90 elementary cellular automaton. N is in powers of 2 because these produce "whole" triangles.'''
M=zeros( (N,2*N+1), dtype=int)
M[0,N]=1
rule=[0, 1, 0, 1, 1, 0, 1, 0]
for j in range(1,N):
for k in range(N-j,N+j+1):
l = 4*M[j-1,k-1] + 2*M[j-1,k] + M[j-1,k+1]
M[j,k]=rule[ l ]
return M
@interact
def _(N=slider([2**a for a in range(0, 12)], label='Number of iterations',default=128), size = slider(1, 20, label= 'Size', step_size=1, default=9 )):
M = sierpinski(N)
plot_M = matrix_plot(M, cmap='binary')
plot_M.show( figsize=[size,size])
Click to the left again to hide and once more to show the dynamic interactive window |
"""
Barnsley fractal fern
Copyright David Joyner wdjoyner@gmail.com, 2008.
Licensed under the creative commons Attribution 3.0 license,
http://creativecommons.org/licenses/by/3.0/us/
"""
def f(i,P):
'''
Returns image under one of four coordinate mapping functions.
'''
x = P[0]
y = P[1]
if i==0:
return [ 0.65*x + 0.04*y - 0.10, -0.01*x + 0.65*y + 1.60]
if i==1:
return [-0.15*x + 0.28*y + 0.10, 0.26*x + 0.24*y + 0.44]
if i==2:
return [0.20*x - 0.26*y - 0.10, 0.23*x + 0.22*y + 0.60]
if i==3:
return [-0.10*x + 0.20*y + 0.10, 0.00*x + 0.26*y + 1.50]
def fern_leaf(num_pts, clr, pt_size=2, initial_pt = (0,0), center_pt = (0,0)):
"""
Returns a graphics objects with num_pts randomly plotted points
in the rgb color, clr, of size pt_size, with initial point initial_pt,
"centered" at center_pt.
EXAMPLES:
sage: P1 = fern_leaf(20000, (0,0.95,0.3))
sage: P2a = fern_leaf(1, (0.9, 0.1, 0), pt_size=50, center_pt=(-0.7,2))
sage: P2b = fern_leaf(2, (0.9, 0.1, 0.5), pt_size=50, center_pt=(1.2,4))
sage: P2c = fern_leaf(2, (0.9, 0.1, 0), pt_size=50, center_pt=(0.3,3.4))
sage: P2d = fern_leaf(2, (0.9, 0.1, 0), pt_size=50, center_pt=(0.25,4.7))
sage: P3 = fern_leaf(2, (0.7, 0.1, 0.2), pt_size=80, center_pt=(1.4,5))
sage: P4a = fern_leaf(2, (0.9, 0.1, 0.2), pt_size=80, center_pt=(1.5,3))
sage: P4b = fern_leaf(2, (0.9, 0.1, 0.6), pt_size=80, center_pt=(1.2,2.8))
sage: P4c = fern_leaf(1, (0.9, 0.1, 0.7), pt_size=80, center_pt=(-0.9,3))
sage: P5a = fern_leaf(2, (0.9, 0.1, 0.2), pt_size=10, center_pt=(1.0,5.8))
sage: P5b = fern_leaf(2, (0.9, 0.1, 0.2), pt_size=10, center_pt=(1.2,5.5))
sage: P5c = fern_leaf(2, (0.9, 0.1, 0.2), pt_size=10, center_pt=(1.5,4.8))
sage: (P2a+P2b+P2c+P2d+P3+P4a+P4b+P4c+P1+P5a+P5b+P5c).show(axes=False)
Here P1 is the branch or tree and P2-P5 are berries or ornaments:-).
"""
p0 = 0.01 # prob pick f0
p1 = 0.07 # prob pick f1
p2 = 0.07 # prob pick f2
p3 = 0.85 # prob pick f3
#p = [p0,1,p2,p3]
P = initial_pt
pts = [P]
for i in range(num_pts):
r = random()
#print i, r
if 0<r<p0:
P = tuple(f(3,P))
pts.append(P)
elif p0<r<p0+p1:
P = tuple(f(1,P))
pts.append(P)
elif p0+p1<r<p0+p1+p2:
P = tuple(f(2,P))
pts.append(P)
else:
P = tuple(f(0,P))
pts.append(P)
if center_pt != (0,0):
pts_centered = [center_pt]
for p in pts:
pts_centered.append((p[0]+center_pt[0],p[1]+center_pt[1]))
else:
pts_centered = pts
return point(pts_centered,rgbcolor=clr, axes=False, pointsize=pt_size)
def f_2008_12_5(i,P):
'''
One of four coordinate mapping functions.
'''
x = P[0]
y = P[1]
if i==0:
return [ 0.85*x + 0.04*y - 0.01, -0.04*x + 0.85*y + 1.60]
if i==1:
return [-0.15*x + 0.28*y + 0.01, 0.26*x + 0.24*y + 0.44]
if i==2:
return [0.20*x - 0.26*y - 0.01, 0.23*x + 0.22*y + 0.60]
if i==3:
return [-0.10*x + 0.20*y + 0.01, 0.00*x + 0.26*y + 1.50]
def f_2008_12_6(i,P):
'''
One of four coordinate mapping functions.
'''
x = P[0]
y = P[1]
if i==0:
return [ 0.85*x + 0.07*y - 0.10, -0.1*x + 0.85*y + 1.60]
if i==1:
return [-0.15*x + 0.28*y + 0.10, 0.26*x + 0.24*y + 0.44]
if i==2:
return [0.20*x - 0.26*y - 0.10, 0.23*x + 0.22*y + 0.60]
if i==3:
return [-0.10*x + 0.20*y + 0.10, 0.00*x + 0.26*y + 1.50]
def f_spikey(i,P):
'''
One of four coordinate mapping functions.
'''
x = P[0]
y = P[1]
if i==0:
return [ 0.45*x + 0.07*y - 0.10, -0.1*x + 0.85*y + 1.60]
if i==1:
return [-0.15*x + 0.28*y + 0.10, 0.26*x + 0.24*y + 0.44]
if i==2:
return [0.20*x - 0.26*y - 0.10, 0.23*x + 0.22*y + 0.60]
if i==3:
return [-0.10*x + 0.20*y + 0.10, 0.00*x + 0.26*y + 1.50]
|
|
P1 = fern_leaf(20000, (0,0.95,0.3))
P1.show()
|
|
def f(i,P): #f_2008_12_6(i,P):
'''
One of four coordinate mapping functions.
'''
x = P[0]
y = P[1]
if i==0:
return [ 0.85*x + 0.07*y - 0.10, -0.1*x + 0.85*y + 1.60]
if i==1:
return [-0.15*x + 0.28*y + 0.10, 0.26*x + 0.24*y + 0.44]
if i==2:
return [0.20*x - 0.26*y - 0.10, 0.23*x + 0.22*y + 0.60]
if i==3:
return [-0.10*x + 0.20*y + 0.10, 0.00*x + 0.26*y + 1.50]
|
|
P1 = fern_leaf(20000, (0,0.95,0.3))
P1.show()
|
|
def f(i,P): #f_spikey(i,P):
'''
One of four coordinate mapping functions.
'''
x = P[0]
y = P[1]
if i==0:
return [ 0.45*x + 0.07*y - 0.10, -0.1*x + 0.85*y + 1.60]
if i==1:
return [-0.15*x + 0.28*y + 0.10, 0.26*x + 0.24*y + 0.44]
if i==2:
return [0.20*x - 0.26*y - 0.10, 0.23*x + 0.22*y + 0.60]
if i==3:
return [-0.10*x + 0.20*y + 0.10, 0.00*x + 0.26*y + 1.50]
|
|
P1 = fern_leaf(20000, (0,0.95,0.3))
P1.show()
|
|
def f(i,P):
'''
One of four coordinate mapping functions.
'''
x = P[0]
y = P[1]
if i==0:
return [ 0.85*x + 0.07*y - 0.10, -0.1*x + 0.85*y + 1.60]
if i==1:
return [-0.15*x + 0.28*y + 0.10, 0.26*x + 0.24*y + 0.44]
if i==2:
return [0.20*x - 0.26*y - 0.10, 0.23*x + 0.22*y + 0.60]
if i==3:
return [-0.10*x + 0.20*y + 0.10, 0.00*x + 0.26*y + 1.50]
P1 = fern_leaf(20000, (0,0.95,0.3))
P2a = fern_leaf(1, (0.9, 0.1, 0), pt_size=50, center_pt=(-0.7,2))
P2b = fern_leaf(2, (0.9, 0.1, 0.5), pt_size=50, center_pt=(1.2,4))
P2c = fern_leaf(2, (0.9, 0.1, 0), pt_size=50, center_pt=(0.3,3.4))
P2d = fern_leaf(2, (0.9, 0.1, 0), pt_size=50, center_pt=(0.25,4.7))
P3 = fern_leaf(2, (0.7, 0.1, 0.2), pt_size=80, center_pt=(1.4,5))
P4a = fern_leaf(2, (0.9, 0.1, 0.2), pt_size=80, center_pt=(1.5,3))
P4b = fern_leaf(2, (0.9, 0.1, 0.6), pt_size=80, center_pt=(1.2,2.8))
P4c = fern_leaf(1, (0.9, 0.1, 0.7), pt_size=80, center_pt=(-0.9,3))
P5a = fern_leaf(2, (0.9, 0.1, 0.2), pt_size=10, center_pt=(1.0,5.8))
P5b = fern_leaf(2, (0.9, 0.1, 0.2), pt_size=10, center_pt=(1.2,5.5))
P5c = fern_leaf(2, (0.9, 0.1, 0.2), pt_size=10, center_pt=(1.5,4.8))
(P2a+P2b+P2c+P2d+P3+P4a+P4b+P4c+P1+P5a+P5b+P5c).show(axes=False)
|
|
# http://christopherolah.wordpress.com/tag/fractals/
|
|