用 python 拟合等角螺线

飞蛾为什么要扑火？

暗梁闻语燕，夜烛见飞蛾。 飞蛾绕残烛，半夜人醉起。 

人类很早就注意到飞蛾扑火这一奇怪的现象，并且自作主张地赋予了飞蛾扑火很多含义，引申出为了理想和追求义无反顾、不畏牺牲的精神。但是，这种引申和比喻，征求过飞蛾的意见吗？

后来，生物学家又提出来昆虫趋光性这一假说来解释飞蛾扑火。不过，这个假说似乎也不成立。如果昆虫真的追逐光明，估计地球上早就没有昆虫了——它们应该齐刷刷整体移民到太阳或月亮上去了。

仔细观察飞蛾扑火，就会发现，昆虫们并不是笔直地飞向光源，而是绕着光源飞行，同时越来越接近光源，最终酿成了“惨案”。这一行为被解释成“失误”似乎更合理一点。既然火烛危险，那么飞蛾为什么要绕着火烛飞行呢？

最新的解释是，飞蛾在夜晚飞行时是依据月光和星光作为参照物进行导航的。星星和月亮离我们非常远，光到了地面上可以看成平行光，当飞蛾的飞行路径保持与光线方向成恒定夹角时，飞蛾就变成了直线飞行，如下图所示。

然而，当飞蛾遇到了火烛等危险光源时，还是按照以前的飞行方式，路径保持与光线方向成恒定夹角，以为依旧能飞成一条直线，结果悲剧了。此时它的飞行轨迹并不是一条直线，而是一条等角螺旋线，如下图所示。

可怜的飞蛾！亿万年进化出来的精准导航，在人工光源的干扰下竟如此不堪。

螺线及等角螺线

螺线家族很庞大，比如，阿基米德螺线、费马螺线、等角螺线、双曲螺线、连锁螺线、斐波那契螺线、欧拉螺线等等。等角螺线，又叫对数螺线，螺线家族的一员。

r = fixed*np.exp(theta/np.tan(fixed)) 

等角螺线在生活中也经常见到，比如，鹦鹉螺的花纹、玫瑰花瓣的排列，星系的悬臂，低气压云图等。

绘制等角螺线

给定中心点和固定角，一个等角螺线就被唯一地确定了。这个螺线可以绕很多圈，可以填满整个宇宙。但很多时候，我们往往只需要观察螺线上的一小部分，这时候就需要两个参数来约定：一个叫作 circle，表示你希望看到多少圈螺线，一个叫作 phase，表示螺线的可见部分向内（顺时针）或向外（逆时针螺线）旋转多少圈。

这是使用 matplotlib 绘制等角螺线的函数，其中固定角参数 fixed 做了一点处理：以度（°）为单位，以零为中心，大于零则为顺时针螺线，小于零则为逆时针螺线

import numpy as np import matplotlib.pyplot as plt from pylab import mpl mpl.rcParams['font.sans-serif'] = ['FangSong'] mpl.rcParams['axes.unicode_minus'] = False def plotSpiral(core, fixed, phase=0, circle=4): """绘制等角螺线 core - 等角螺线的中心坐标，tuple类型 fixed - 等角螺线的固定角度，单位：度（°）。fixed大于零则为顺时针螺线，小于零则为逆时针螺线 phase - 初始相位，单位：圈（360°）。对顺时针螺线，该数值越大，螺线越大，对逆时针螺线则相反 circle - 螺线可见部分的圈数，单位：圈（360°） """ plt.axis("equal") plt.plot([core[0]], [core[1]], c='red', marker='+', markersize=10) fixed_rad = np.radians(90 + fixed) theta = np.linspace(0, circle*2*np.pi, 361) + phase*2*np.pi r = fixed_rad*np.exp(theta/np.tan(fixed_rad)) x = r*np.cos(theta) + core[0] y = r*np.sin(theta) - core[1] plt.plot(x, y, c='blue') plt.show() 

拟合等角螺线

import numpy as np from scipy import optimize def fit_spiral(core, dots): """拟合等角螺线，返回定角fixed，初始相位phase""" fixed_ccw = 0.445*np.pi fixed_cw = 0.555*np.pi # 将dots拆分成x_list和y_list x_list, y_list = list(), list() for x, y in dots: x_list.append(x-core[0]) y_list.append(y-core[1]) # 计算距离 x = np.array(x_list) y = np.array(y_list) r = np.hypot(x,y) # 按照距离排序 sort_mask = np.argsort(r) x = x[sort_mask] y = y[sort_mask] r = r[sort_mask] # 计算角度 theta = np.arctan(y/x) theta[x<0] += np.pi # 确定顺序（CW-顺时针，CCW-逆时针） d = np.diff(theta) print(d) ccw = d[d>0].size > d[d<0].size print('ccw=',ccw) # 调整角度为升序（CCW）或降序（CW） if ccw: for i in range(1, theta.size): while theta[i] < theta[i-1]: theta[i] += 2*np.pi dtheta = theta[i] - theta[i-1] while r[i]/r[i-1] > 1.8*np.exp(dtheta/np.tan(fixed_ccw)): theta[i] += 2*np.pi dtheta = theta[i] - theta[i-1] else: for i in range(theta.size-1)[::-1]: while theta[i] < theta[i+1]: theta[i] += 2*np.pi dtheta = theta[i+1] - theta[i] while r[i+1]/r[i] > 1.8*np.exp(dtheta/np.tan(fixed_cw)): theta[i] += 2*np.pi dtheta = theta[i+1] - theta[i] # 定义拟合函数 def fmax(theta, fixed, phase): fixed = np.radians(90 + fixed) return fixed*np.exp((theta+phase*2*np.pi)/np.tan(fixed)) try: fita, fitb = optimize.curve_fit(fmax, theta, r, [2-int(ccw), 0], maxfev=10000) return fita except: return None core = (530, 496) dots = [(467,538), (448,675), (522,484), (513,451), (811,519)] result = fit_spiral(core, dots) if isinstance(result, np.ndarray): plotSpiral(core, result[0], phase=result[1], circle=4) else: print(u'拟合失败')