Я пытаюсь вычислить спектр тормозного излучения, который включает в себя вычисление ускорения, преобразованного Фурье. Теперь я решаю нелинейный ОДУ, чтобы численно рассчитать ускорение во временной области. Теперь, после преобразования Фурье с использованием fft от Numpy, результирующий спектр выглядит крайне негладким и «нефизическим». Я не могу вставить весь код, поэтому я публикую то, что я считаю релевантным фрагментом. Может кто-то указать, что я делаю не так?
Примечание. Мое ускорение является функцией двух переменных (бета и b, параметр воздействия), и я хочу построить его для другого b, и чтобы вычленить бета-член, я простосуммирование по всем значениям ускорения для другой бета-версии для заданного параметра удара b.
Также мой спектр представляет собой квадрат ускорения с преобразованием Фурье, умноженный на постоянный коэффициент (формула Лармора)
'''
Fourier Transformation of the acceleration
'''
N = 2**8
sampling_frequency = 10
a_w_normalized = [[] for a in range(len(impact_parameter))]
a_w =[[] for a in range(len(impact_parameter))]
intensity_normalized = [[] for a in range(len(impact_parameter))]
acc_summed_over_velocity = []
intensity_summed_over_velocity =[]
for index,b in enumerate(impact_parameter):
acc_sum =[0 for x in range(N//2 +1)]
intensity_sum = [0 for x in range(N)]
for j in range(len(velocity_z_component)):
#window_kaiser = signal.kaiser(N, 15)
#window_hann = signal.hann(N,sym=True)
window =1
fft_input = acc_normalized[index][j]*window
ft_acc_normalized = np.abs(np.fft.rfft(fft_input,norm=None))
minima_norm = signal.argrelextrema(ft_acc_normalized, np.less)[0]
acc_sum =np.add(ft_acc_normalized,acc_sum )
intensity_list = [power_spectrum_factor * (a ** 2) for a in ft_acc_normalized]
intensity_sum = np.add(intensity_list,intensity_sum)
intensity_normalized[index].append(intensity_list)
a_w_normalized[index].append(ft_acc_normalized)
acc_summed_over_velocity.append(acc_sum)
intensity_summed_over_velocity.append(intensity_sum *acceleration_factor**2)
#intensity_summed_over_velocity= ma.masked_less_equal(intensity_summed_over_velocity,1e-5)
'''
Plotting acceleration for selected values of impact paramters in time and frequency domain
'''
if plot is True:
plt.figure(figsize=(12, 8))
for i in range(len(impact_parameter)):
plt.plot(t,acc_timedomain_summed_over_velocity[i], label='b={:.3f}'.format(impact_parameter[i]), )
plt.legend()
plt.ylabel(r'$ a(\tilde t) $', fontsize=14)
plt.xlabel(r'$ \tilde t $', fontsize=14)
ax.spines['left'].set_position('center')
ax.spines['bottom'].set_position('center')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_smart_bounds(True)
ax.spines['bottom'].set_smart_bounds(True)
ax.xaxis.set_ticks_position('bottom')
ax.yaxis.set_ticks_position('left')
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
if screening is False:
plt.title('a(t) vs time without screening',fontsize=15)
plt.savefig('./Plots/acc_without_screening.png')
else:
plt.title('a(t) vs time with screening',fontsize=15)
plt.savefig('./Plots/acc_with_screening.png')
'''
Fourier transformed intensity for different impact paramater
'''
plt.figure(figsize=(12, 8))
freq_normalized = np.fft.fftfreq(N)*(2*np.pi*sampling_frequency)
for index,b in enumerate(plasma.impact_parameter):
plt.plot(np.abs(freq_normalized), intensity_summed_over_velocity[index], label=r'$ \tilde b={:.2f}$'.format(b), )
plt.xlabel(r'$ \tilde \omega $', fontsize=14)
plt.ylabel(r'$ I_{\omega} $', fontsize=16)
plt.legend(loc="lower left")
plt.xscale('log')
plt.yscale('log')
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
if screening is False:
plt.title('Single particle spectrum without screening')
plt.savefig('./Plots/spectrum_no_screening.png')
else:
plt.title('Single particle spectrum with screening')
plt.savefig('./Plots/spectrum_debye_screening.png')
plt.show()