用python拟合等角螺线的实现示例
人类很早就注意到飞蛾扑火这一奇怪的现象,并且自作主张地赋予了飞蛾扑火很多含义,引申出为了理想和追求义无反顾、不畏牺牲的精神。但是,这种引申和比喻,征求过飞蛾的意见吗?
后来,生物学家又提出来昆虫趋光性这一假说来解释飞蛾扑火。不过,这个假说似乎也不成立。如果昆虫真的追逐光明,估计地球上早就没有昆虫了——它们应该齐刷刷整体移民到太阳或月亮上去了。
仔细观察飞蛾扑火,就会发现,昆虫们并不是笔直地飞向光源,而是绕着光源飞行,同时越来越接近光源,最终酿成了“惨案”。这一行为被解释成“失误”似乎更合理一点。既然火烛危险,那么飞蛾为什么要绕着火烛飞行呢?
最新的解释是,飞蛾在夜晚飞行时是依据月光和星光作为参照物进行导航的。星星和月亮离我们非常远,光到了地面上可以看成平行光,当飞蛾的飞行路径保持与光线方向成恒定夹角时,飞蛾就变成了直线飞行,如下图所示。
螺线及等角螺线
螺线家族很庞大,比如,阿基米德螺线、费马螺线、等角螺线、双曲螺线、连锁螺线、斐波那契螺线、欧拉螺线等等。等角螺线,又叫对数螺线,螺线家族的一员。
早在2000多年以前,古希腊数学家阿基米德就对螺旋线进行了研究。公元1638年,著名数学家笛卡尔首先描述了对数螺旋线(等角螺旋线),并列出了螺旋线的解析式。这种螺旋线有很多特点,其中最突出的一点就是它的形状,无论你把它放大或缩小它都不会有任何的改变。就像我们不能把角放大或缩小一样。
用极坐标分析法分析飞蛾扑火的飞行轨迹,可知,轨迹线上任意一点的切线与该点与原点的连线之间的夹角是固定的,这就是等角螺线得名的由来。因为分析过程使用了对数,所以等角螺线又叫对数螺线。我不太会用LaTeX写数学公式,所以就用python的方法写出螺线方程。其中,fixed表示螺线固定角,大于pi/2则为顺时针螺线,小于pi/2则为逆时针螺线。theta表示旋转弧度,r表示距离中心点距离。
r=fixed*np.exp(theta/np.tan(fixed))
等角螺线在生活中也经常见到,比如,鹦鹉螺的花纹、玫瑰花瓣的排列,星系的悬臂,低气压云图等。
绘制等角螺线
给定中心点和固定角,一个等角螺线就被唯一地确定了。这个螺线可以绕很多圈,可以填满整个宇宙。但很多时候,我们往往只需要观察螺线上的一小部分,这时候就需要两个参数来约定:一个叫作circle,表示你希望看到多少圈螺线,一个叫作phase,表示螺线的可见部分向内(顺时针)或向外(逆时针螺线)旋转多少圈。
这是使用matplotlib绘制等角螺线的函数,其中固定角参数fixed做了一点处理:以度(°)为单位,以零为中心,大于零则为顺时针螺线,小于零则为逆时针螺线
importnumpyasnp importmatplotlib.pyplotasplt frompylabimportmpl mpl.rcParams['font.sans-serif']=['FangSong'] mpl.rcParams['axes.unicode_minus']=False defplotSpiral(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()
下图展示了逆时针等角螺线各个参数的意义:
下图展示了顺时针等角螺线各个参数的意义:
拟合等角螺线
在台风定位时,需要手动确定台风中心位置,并标识出台风螺线轨迹上的部分点,然后逆合出螺线方程。如下图所示,蓝色十字为台风中心点,5个黄色圆点是手工标注的台风螺线轨迹上的点。
以下为拟合函数
importnumpyasnp fromscipyimportoptimize deffit_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() 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].size print('ccw=',ccw) #调整角度为升序(CCW)或降序(CW) ifccw: foriinrange(1,theta.size): whiletheta[i]1.8*np.exp(dtheta/np.tan(fixed_ccw)): theta[i]+=2*np.pi dtheta=theta[i]-theta[i-1] else: foriinrange(theta.size-1)[::-1]: whiletheta[i] 1.8*np.exp(dtheta/np.tan(fixed_cw)): theta[i]+=2*np.pi dtheta=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) returnfita except: returnNone core=(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'拟合失败')
拟合效果如下图:
以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持毛票票。
声明:本文内容来源于网络,版权归原作者所有,内容由互联网用户自发贡献自行上传,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任。如果您发现有涉嫌版权的内容,欢迎发送邮件至:czq8825#qq.com(发邮件时,请将#更换为@)进行举报,并提供相关证据,一经查实,本站将立刻删除涉嫌侵权内容。