I am trying to obtain harmonics using fft and fft shift of ay_list of the following program. I dont understand why I am not able to obtain first harmonics at omega/omega0. Can some one plea ehelp me to visualize the first harmonics in frequency axis. I have tried lot but could land into right answe.
Here, ay_list is the acceleration of the electron. In our case omega0 is w1. and I guess omega is xfa that we have obtained after applying fftshift. I request to please help me with the problem.
Thanking you, Kind regards
import numpy as np # for math
import matplotlib.pyplot as plt # for graphs
#import mpld3;mpld3.enable_notebook() # interactive graphs
import seaborn; seaborn.set() # prettier graphs
from scipy.fftpack import fft, fftfreq, fftshift
import matplotlib.pyplot as plt
from numpy import zeros
##### Constants:
q = 1.602e-19 # Coulombs Charge of electron
c = 2.99792458e8 # m/s Speed of light
eo = 8.8541e-12 # C^2/(Nm^2) Permittivity of vacuum
me = 9.109e-31 # kg Mass of electron
ke = 8.985551e9 #N m^2 C-2 Coulomb's constant
pulseFWHM=40.e-15
"""############## Arrays to store value ###################################"""
Zero_array = 40000
Ey_list = np.zeros(Zero_array) #store the y-values for later
y_list = np.zeros(Zero_array) #store the y-values for later
ay_list = np.zeros(Zero_array)
vy_list= np.zeros(Zero_array)
constant = 1
Velocity_mu = 0
Velocity_sigma = 2.18*1e6/20# mean and standard deviation
Velocity_s_y = np.random.normal(Velocity_mu, Velocity_sigma, constant)
Velocity_s_x = np.random.normal(Velocity_mu, Velocity_sigma, constant)
Position_mu = 0
Position_sigma = 0.05e-9/20# mean and standard deviation
Position_s_y = np.random.normal(Position_mu, Position_sigma, constant)
Position_s_x = np.random.normal(Position_mu, Position_sigma, constant)
def E_and_P(t=np.linspace(-20*800e-9/c,20*800e-9/c,100*40),
I1=1e14,I2=1e14,tb=0,
wavelength1=800e-9,wavelength2=400e-9,
ellip1=0,ellip2=0,
delay=0,
time_step=.0001e-15, k = 0, vy = 0):
w1 = c/wavelength1 * 2. * np.pi # Angular frequency of field 1
Eo1 = np.sqrt(2*I1*10**4/(c*eo)) # Electric field in V/m
# Scale Eo to get same intensity for all ellipticity values
Eo1 = -Eo1/np.sqrt(ellip1**2+1)
# Calculate components of E-fields
Factor= 15
Ey1 = Factor*ellip1*Eo1*np.sin(w1*t) * np.exp(-t**2/(2*(pulseFWHM / 2.35482)**2))
# Add fields to make circular/elliptical
Ey = Ey1
x,y,z = Position_s_x[tb],Position_s_y[tb],0
vx,vy,vz = (Velocity_s_x[tb],Velocity_s_y[tb],0)
for index in range(0,len(t)):
r = np.sqrt(x**2+y**2+z**2)
EyG = Ey[index]
Fcy = -q**2 * ke * y/np.sqrt(r**2+ 5.40e-11**2)**3
ay = EyG*(-q)/me + Fcy/me
vy = vy + ay*time_step
y = y + vy * time_step
#####################################################################
Ey_list[index] = EyG ####################################################################
y_list[index] = y
ay_list[index] = ay
vy_list[index] = vy ####################################################################
matrix_ay[k][index]= ay
return ay_list
def main_a():
k = -1
vy = 0
for tb in range(1):
k = k+1
ay_list=
E_and_P(t=np.linspace(-20*800e9/c,20*800e9/c,1000*40),I1=1e14,I2=1e14,tb=tb,wavelength1=800e-9,wavelength2=400e-9, ellip1=1,ellip2=-1,delay=0,time_step=.001e-15, k = k,vy = vy)
plt.subplot(2,1,1)
plt.plot(ay_list)
N = 40000
f_ay = np.fft.fft(ay_list)
ayplot = np.fft.fftshift(f_ay)
xfa = np.fft.fftfreq(N, 15*2.66e-15/40000)
xfa = np.fft.fftshift(xfa)
plt.subplot(2,1,2)
plt.semilogy(xfa, 1.0/N * np.abs(ayplot))
#plt.xlim(0,90e15)
if __name__=="__main__":
main_a()