文章目录飞蛾为什么要扑火？螺线及等角螺线绘制等角螺线拟合等角螺线
飞蛾为什么要扑火？暗梁闻语燕，夜烛见飞蛾。飞蛾绕残烛，半夜人醉起。
人类很早就注意到飞蛾扑火这一奇怪的现象，并且自作主张地赋予了飞蛾扑火很多含义，引申出为了理想和追求义无反顾、不畏牺牲的精神。但是，这种引申和比喻，征求过飞蛾的意见吗？
后来，生物学家又提出来昆虫趋光性这一假说来解释飞蛾扑火。不过，这个假说似乎也不成立。如果昆虫真的追逐光明，估计地球上早就没有昆虫了——它们应该齐刷刷整体移民到太阳或月亮上去了。
仔细观察飞蛾扑火，就会发现，昆虫们并不是笔直地飞向光源，而是绕着光源飞行，同时越来越接近光源，最终酿成了“惨案”。这一行为被解释成“失误”似乎更合理一点。既然火烛危险，那么飞蛾为什么要绕着火烛飞行呢？
最新的解释是，飞蛾在夜晚飞行时是依据月光和星光作为参照物进行导航的。星星和月亮离我们非常远，光到了地面上可以看成平行光，当飞蛾的飞行路径保持与光线方向成恒定夹角时，飞蛾就变成了直线飞行，如下图所示。
然而，当飞蛾遇到了火烛等危险光源时，还是按照以前的飞行方式，路径保持与光线方向成恒定夹角，以为依旧能飞成一条直线，结果悲剧了。此时它的飞行轨迹并不是一条直线，而是一条等角螺旋线，如下图所示。
可怜的飞蛾！亿万年进化出来的精准导航，在人工光源的干扰下竟如此不堪。
螺线及等角螺线
螺线家族很庞大，比如，阿基米德螺线、费马螺线、等角螺线、双曲螺线、连锁螺线、斐波那契螺线、欧拉螺线等等。等角螺线，又叫对数螺线，螺线家族的一员。
早在2000多年以前，古希腊数学家阿基米德就对螺旋线进行了研究。公元1638年，著名数学家笛卡尔首先描述了对数螺旋线（等角螺旋线），并列出了螺旋线的解析式。这种螺旋线有很多特点，其中最突出的一点就是它的形状，无论你把它放大或缩小它都不会有任何的改变。就像我们不能把角放大或缩小一样。用极坐标分析法分析飞蛾扑火的飞行轨迹，可知，轨迹线上任意一点的切线与该点与原点的连线之间的夹角是固定的，这就是等角螺线得名的由来。因为分析过程使用了对数，所以等角螺线又叫对数螺线。我不太会用LaTeX写数学公式，所以就用python的方法写出螺线方程。其中，fixed表示螺线固定角，大于pi/2则为顺时针螺线，小于pi/2则为逆时针螺线。theta表示旋转弧度，r表示距离中心点距离。
r=fixed*np.exp(theta/np.tan(fixed))
等角螺线在生活中也经常见到，比如，鹦鹉螺的花纹、玫瑰花瓣的排列，星系的悬臂，低气压云图等。
绘制等角螺线
给定中心点和固定角，一个等角螺线就被唯一地确定了。这个螺线可以绕很多圈，可以填满整个宇宙。但很多时候，我们往往只需要观察螺线上的一小部分，这时候就需要两个参数来约定：一个叫作circle，表示你希望看到多少圈螺线，一个叫作phase，表示螺线的可见部分向内（顺时针）或向外（逆时针螺线）旋转多少圈。
这是使用matplotlib绘制等角螺线的函数，其中固定角参数fixed做了一点处理：以度（°）为单位，以零为中心，大于零则为顺时针螺线，小于零则为逆时针螺线
importnumpyasnpimportmatplotlib.pyplotaspltfrompylabimportmplmpl.rcParams['font.sans-serif']=['FangSong']mpl.rcParams['axes.unicode_minus']=FalsedefplotSpiral(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.pir=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()
下图展示了逆时针等角螺线各个参数的意义：下图展示了顺时针等角螺线各个参数的意义：
拟合等角螺线
在台风定位时，需要手动确定台风中心位置，并标识出台风螺线轨迹上的部分点，然后逆合出螺线方程。如下图所示，蓝色十字为台风中心点，5个黄色圆点是手工标注的台风螺线轨迹上的点。以下为拟合函数
importnumpyasnpfromscipyimportoptimizedeffit_spiral(core,dots):"""拟合等角螺线，返回定角fixed，初始相位phase"""fixed_ccw=0.445*np.pifixed_cw=0.555*np.pi#将dots拆分成x_list和y_listx_list,y_list=list(),list()forx,yindots: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].sizeprint('ccw=',ccw)#调整角度为升序（CCW）或降序（CW）ifccw:foriinrange(1,theta.size):whiletheta[i]<theta[i-1]:theta[i]+=2*np.pidtheta=theta[i]-theta[i-1]whiler[i]/r[i-1]>1.8*np.exp(dtheta/np.tan(fixed_ccw)):theta[i]+=2*np.pidtheta=theta[i]-theta[i-1]else:foriinrange(theta.size-1)[::-1]:whiletheta[i]<theta[i+1]:theta[i]+=2*np.pidtheta=theta[i+1]-theta[i]whiler[i+1]/r[i]>1.8*np.exp(dtheta/np.tan(fixed_cw)):theta[i]+=2*np.pidtheta=theta[i+1]-theta[i]#定义拟合函数deffmax(theta,fixed,phase):fixed=np.radians(90+fixed)returnfixed*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)returnfitaexcept:returnNonecore=(530,496)dots=[(467,538),(448,675),(522,484),(513,451),(811,519)]result=fit_spiral(core,dots)ifisinstance(result,np.ndarray):plotSpiral(core,result[0],phase=result[1],circle=4)else:print(u'拟合失败')
拟合效果如下图：
