当前位置：首页 > news >正文

先睹为快_Mandelbrot集

news 来源：原创 2024/5/18 22:41:10

本文引用自作者编写的下述图书; 本文允许以个人学习、教学等目的引用、讲授或转载，但需要注明原作者"海洋饼干叔
叔"；本文不允许以纸质及电子出版为目的进行抄摘或改编。
1.《Python编程基础及应用》，陈波，刘慧君，高等教育出版社。免费授课视频 Python编程基础及应用
2.《Python编程基础及应用实验教程》, 陈波，熊心志，张全和，刘慧君，赵恒军，高等教育出版社Python编程基础及应用实验教程
3. 《简明C及C++语言教程》，陈波，待出版书稿。免费授课视频

先睹为快_Mandelbrot集

简易版本的Mandelbrot集

"""
Draw mandelbrot set with matplotlib.
Matrix generated by Pure Python code.
chenbo@cqu.edu.cn,   Dec 2018.

Usage:
    key-up, key-down:  switch color map
"""



from matplotlib import pyplot as plt
from matplotlib import cm
import numpy as np

def getEscapeTime(c):
    "计算参数c的逃逸时间，该逃逸速度将用作点的颜色"
    z = 0
    for i in range(100):
        if abs(z) > 2:
            return i
        z = z*z + c
    return i

def computeMandelbrot(xCenter,yCenter,semiWidth,N):
    xFrom,xTo,yFrom,yTo = xCenter-semiWidth,xCenter+semiWidth,\
                          yCenter-semiWidth,yCenter+semiWidth
    y,x = np.ogrid[yFrom:yTo:N*1j,xFrom:xTo:N*1j]
    c = x + y*1j
    print("c.shape:",c.shape,"x.shape:",x.shape,"y.shape:",y.shape)
    return np.frompyfunc(getEscapeTime,1,1)(c).astype(np.float)

def drawMandelbrot(ax,xCenter,yCenter,semiWidth,N,cmap):
    "(xCenter,yCenter)-中心点，semiWidth-矩形半宽，N*N像素."
    ax.set_axis_off()
    ds = computeMandelbrot(xCenter,yCenter,semiWidth,N)
    ax.imshow(ds,cmap=cmap)

def refresh():
    c = cm.get_cmap(para.cmaps[para.idxColorMap%len(para.cmaps)])
    drawMandelbrot(para.ax0, para.x, para.y,
                   semiWidth=0.2, N=600,cmap=c)
    drawMandelbrot(para.ax1, para.x, para.y,
                   semiWidth=0.2 ** 3, N=600,cmap=c)
    para.fig.canvas.draw()

def on_key_release(event):
    if event.key == 'up':
        para.idxColorMap+=1
    elif event.key == 'down':
        para.idxColorMap-=1
    else:
        return
    refresh()

class Para:
    pass

para = Para()
para.x, para.y = 0.27322626, 0.595153338
para.idxColorMap = 0
para.cmaps = ['rainbow', 'jet', 'nipy_spectral', 'gist_ncar','flag',
              'prism', 'ocean', 'gist_earth', 'terrain', 'gist_stern',
              'gnuplot', 'gnuplot2', 'CMRmap', 'cubehelix', 'brg',
              'gist_rainbow']
para.fig = plt.figure(figsize=(12,6),dpi=100)
para.fig.canvas.mpl_connect('key_release_event',on_key_release)
plt.subplots_adjust(0,0,1,1,0.0,0)
para.ax0 = plt.subplot(121)
para.ax1 = plt.subplot(122)
refresh()
plt.show()

在这里插入图片描述

交互版本的Mandelbrot集

"""
Draw mandelbrot set with matplotlib.
Escape time computed by Cython.
chenbo@cqu.edu.cn,   Dec 2018.

Usage:
    left-click mouse: Zoom in
    right-click mouse: Zoom out
    key escape:  Reset
    key up,key down: Switch color map.
"""

from matplotlib import pyplot as plt
from matplotlib import cm
import numpy as np
from MandelbrotComp import getEscapeTime

def computeMandelbrot(xCenter,yCenter,semiWidth,N):
    xFrom,xTo,yFrom,yTo = xCenter-semiWidth,xCenter+semiWidth,\
                          yCenter-semiWidth,yCenter+semiWidth
    y,x = np.ogrid[yFrom:yTo:N*1j,xFrom:xTo:N*1j]
    c = x + y*1j
    print("c.shape:",c.shape,"x.shape:",x.shape,"y.shape:",y.shape)
    return np.frompyfunc(getEscapeTime,1,1)(c).astype(np.float)

def drawMandelbrot(ax,xCenter,yCenter,semiWidth,N,cmap):
    "(xCenter,yCenter)-中心点，semiWidth-矩形半宽，N*N像素."
    ax.set_axis_off()
    ds = computeMandelbrot(xCenter,yCenter,semiWidth,N)
    ax.imshow(ds,cmap=cmap)

def refresh():
    print("x,y =",para.x,para.y,"semiWidth =", para.semiWidth)
    c = cm.get_cmap(para.cmaps[para.idxColorMap%len(para.cmaps)])
    drawMandelbrot(para.ax, para.x, para.y,
                   para.semiWidth, N=700,cmap=c)
    para.fig.canvas.draw()

def on_key_release(event):
    if event.key == 'up':
        para.idxColorMap+=1
    elif event.key == 'down':
        para.idxColorMap-=1
    elif event.key == 'escape':
        para.x, para.y = -0.5, 0.0
        para.semiWidth = 1.5
    else:
        return
    refresh()

def on_button_release(event):
    para.x = (para.x - para.semiWidth) + \
             2*para.semiWidth*event.xdata/para.figWidth
    para.y = (para.y - para.semiWidth) + \
             2*para.semiWidth*event.ydata/para.figHeight
    if event.button == 1:
        para.semiWidth /= 3.0
    elif event.button == 3:
        para.semiWidth *= 3.0
    refresh()

class Para:
    pass

para = Para()
para.x, para.y = -0.5, 0.0
para.semiWidth = 1.5
para.idxColorMap = 0
para.cmaps = ['rainbow', 'jet', 'nipy_spectral', 'gist_ncar','flag',
              'prism', 'ocean', 'gist_earth', 'terrain', 'gist_stern',
              'gnuplot', 'gnuplot2', 'CMRmap', 'cubehelix', 'brg',
              'gist_rainbow']
para.fig = plt.figure(figsize=(7,7),dpi=100)
para.figWidth,para.figHeight = 700,700
para.fig.canvas.mpl_connect('key_release_event',on_key_release)
para.fig.canvas.mpl_connect('button_release_event',on_button_release)
plt.subplots_adjust(0,0,1,1,0.0,0)
para.ax = plt.subplot(111)
refresh()
plt.show()

这个需要Cython的支持，使用Cython语言编写一个核心函数，然后将其编译成Python模块（MandelbrotComp ）来调用。

def getEscapeTime(complex c):
    "计算参数c的逃逸时间，该逃逸速度将用作点的颜色"
    cdef complex z = 0
    cdef int i
    for i in range(100):
        if z.real * z.real + z.imag*z.imag > 4:
            return i
        z = z*z + c
    return i