当前位置:首页 » python教程 » 正文

用python拟合等角螺线的实现示例

456 次  2020-12-25  分类 : python教程

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

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

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

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

在这里插入图片描述

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

在这里插入图片描述

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

螺线及等角螺线

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

早在2000多年以前,古希腊数学家阿基米德就对螺旋线进行了研究。公元1638年,著名数学家笛卡尔首先描述了对数螺旋线(等角螺旋线),并列出了螺旋线的解析式。这种螺旋线有很多特点,其中最突出的一点就是它的形状,无论你把它放大或缩小它都不会有任何的改变。就像我们不能把角放大或缩小一样。

在这里插入图片描述

用极坐标分析法分析飞蛾扑火的飞行轨迹,可知,轨迹线上任意一点的切线与该点与原点的连线之间的夹角是固定的,这就是等角螺线得名的由来。因为分析过程使用了对数,所以等角螺线又叫对数螺线。我不太会用LaTeX写数学公式,所以就用 python 的方法写出螺线方程。其中,fixed 表示螺线固定角,大于 pi/2 则为顺时针螺线,小于 pi/2 则为逆时针螺线。theta 表示旋转弧度,r 表示距离中心点距离。

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()

下图展示了逆时针等角螺线各个参数的意义:

在这里插入图片描述

下图展示了顺时针等角螺线各个参数的意义:

在这里插入图片描述

拟合等角螺线

在台风定位时,需要手动确定台风中心位置,并标识出台风螺线轨迹上的部分点,然后逆合出螺线方程。如下图所示,蓝色十字为台风中心点,5个黄色圆点是手工标注的台风螺线轨迹上的点。

在这里插入图片描述

以下为拟合函数

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'拟合失败')

拟合效果如下图:

在这里插入图片描述

以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持python博客。

来源:python博客 欢迎分享!

本文链接:https://www.94e.cn/info/5159

标签:numpy  matplotlib  

<< 上一篇 下一篇 >>

  Powered By python教程网   鲁ICP备18013710号
python博客 - 小白学python最友好的网站!

扫一扫,关注微信公众号 扫一扫,关注微信公众号