Отказ от ответственности: Слово "pythonic" не имеет большого смысла для меня;так что это просто два решения, если им нужно имя, я бы предложил «erinaceidaeic» или «axolotlable».
В этой проблеме есть физический компонент, который будет определять, что "ширина "значит.Я мог бы представить, что было бы интересно определить параметр w моды Эрмита-Гаусса.Однако, для цели этого вопроса как вопроса программирования, скажем так, вы хотите найти полную ширину на половине максимума (FWHM).
Вот две функции FWHM.
import numpy as np
#########
# FWHM function
#########
def cFWHM(x,y):
""" returns coarse full width at half maximum, and the two
xcoordinates of the first and last values above the half maximum """
where, = np.where(y >= y.max()/2.)
maxi = x[where[-1]]
mini = x[where[0]]
return maxi-mini, mini, maxi
def fFWHM(x,y):
""" returns interpolated full width at half maximum, and the two
xcoordinates at the (interpolated) half maximum """
def find_roots(x,y):
s = np.abs(np.diff(np.sign(y))).astype(bool)
return x[:-1][s] + np.diff(x)[s]/(np.abs(y[1:][s]/y[:-1][s])+1)
z = find_roots(x,y-y.max()/2)
return z.max()-z.min(), z.min(), z.max()
Первая, как упоминалось в вопросе, найдет минимальную и максимальную координаты вдоль осей, где значение y больше или равно половине максимального значения y в данных.Разница между ними заключается в ширине.Это дает разумные результаты для достаточно плотных точек.
Если требуется большая точность, можно вместо этого найти нули y-ymax/2 путем интерполяции между точками данных.(Решение взято из моего ответа на этот вопрос ).
Полный пример:
import matplotlib.pyplot as plt
#########
# Generate some data
#########
def H(n, x):
# Get the nth hermite polynomial, evaluated at x
coeff = np.zeros(n)
coeff[-1] = 1
return np.polynomial.hermite.hermval(x, coeff)
def E(x,y,w,l,m, E0=1):
# get the hermite-gaussian TEM_l,m mode in the focus (z=0)
return E0*H(l,np.sqrt(2)*x/w)*H(m,np.sqrt(2)*y/w)*np.exp(-(x**2+y**2)/w**2)
g = np.linspace(-4.5e-3,4.5e-3,901)
X,Y = np.meshgrid(g,g)
f = E(X,Y,.9e-3, 5,7)**2
# Intensity profiles along x and y direction
Int_x = np.sum(f, axis=0)
Int_y = np.sum(f, axis=1)
#########
# Plotting
#########
fig = plt.figure(figsize=(8,4.5))
ax = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(2,2,2)
ax3 = fig.add_subplot(2,2,4)
dx = np.diff(X[0])[0]; dy = np.diff(Y[:,0])[0]
extent = [X[0,0]-dx/2., X[0,-1]+dx/2., Y[0,0]-dy/2., Y[-1,0]+dy/2.]
ax.imshow(f, extent=extent)
ax2.plot(g,Int_x)
ax3.plot(g,Int_y)
width, x1, x2 = cFWHM(g,Int_x) # compare to fFWHM(g,Int_x)
height, y1, y2 = cFWHM(g,Int_y)
ax2.plot([x1, x2],[Int_x.max()/2.]*2, color="crimson", marker="o")
ax3.plot([y1, y2],[Int_y.max()/2.]*2, color="crimson", marker="o")
annkw = dict(xytext=(0,10),
textcoords="offset pixels", color="crimson", ha="center")
ax2.annotate(width, xy=(x1+width/2, Int_x.max()/2.), **annkw)
ax3.annotate(height, xy=(y1+height/2, Int_y.max()/2.), **annkw)
plt.tight_layout()
plt.show()
А вот визуальное сравнение двух функций.Максимальная ошибка использования первой вместо второй функции - это разница между последующими значениями данных.В этом случае это, вероятно, будет разрешение камеры.(Однако обратите внимание, что истинная ошибка, конечно, должна учитывать ошибку при определении максимального значения и, следовательно, вероятно, будет намного больше.)
