1. 程式人生 > 程式設計 >用python擬合等角螺線的實現示例

用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,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'擬合失敗')

擬合效果如下圖:

在這裡插入圖片描述

以上就是本文的全部內容,希望對大家的學習有所幫助,也希望大家多多支援我們。