0
$\begingroup$

I am trying to calculate the spectrum of Bremmstrahlung, which involves calculating the Fourier transformed acceleration. I am solving a non-linear ODE to numerically calculate the acceleration in the time domain. After taking the Fourier transformation using Numpy's fft, the resultant spectrum looks highly non-smooth and "non-physical" . I cannot paste the entire code so I am posting what I think is relevant snippet. Can someone point out what I am doing wrong?

Note: My acceleration is a function of two variables (beta, and b, the impact parameter), and I want to plot it for the different b, and in order to factor out the beta term I am just summing over all the values of acceleration for the different beta, for a given impact parameter b.

Also my spectrum is Fourier transformed acceleration square times a constant factor (Larmor's formula) enter image description here enter image description here

'''
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)]
    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.fft(fft_input,norm=None))

        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()
$\endgroup$
  • $\begingroup$ I think there's some confusion about the order of the coefficients in the spectrum. Check np.fft.fftshift() to get your spectrum in the correct order. $\endgroup$ – AlexE Oct 15 at 14:49
  • $\begingroup$ @AlexE Thanks for the comment. fftshift shift the zero-frequency component to the center of the spectrum. I am not sure how that will help me in getting a smooth spectrum. I expect my Fourier transformed acceleration to rise up and then fall. $\endgroup$ – Prav001 Oct 15 at 18:00
  • $\begingroup$ I have a hard time to follow the outline in your question but just one thing: You said the spectra are highly non-smooth, but from the images that you attached I don't see that, cause they look pretty okayish to me. At least there are no wild noisy behavior, so I'm not sure what's your purpose here. $\endgroup$ – Alone Programmer Oct 17 at 13:45
  • $\begingroup$ I missed out on your use of fftfreq(). This should get you frequencies in a matching order, so you're probably fine without fftshift. Why are you abs()ing your frequencies for the plot? $\endgroup$ – AlexE Oct 17 at 14:57
  • $\begingroup$ @AlexE Well I have to plot the final spectrum on log-log plot so I need to take the absolute value. Anyway, since my data is real, I think all the information is contained in positive frequencies. I should have used rfft and rfftfreq to avoid confusion. $\endgroup$ – Prav001 Oct 17 at 18:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.