import numpy as np
import matplotlib.pyplot as plt
def GetSASE(CentralEnergy, dE_FWHM, dt_FWHM, samples=0, Axis=True):
h=4.135667662 #in eV*fs
dE=dE_FWHM/2.355 #in eV, converts to sigma
dt=dt_FWHM/2.355 #in fs, converts to sigma
if samples == 0:
samples=int(400.*dt*CentralEnergy/h)
else:
if (samples < 400.*dt*CentralEnergy/h):
print("Number of samples is a little small, proceeding anyway. Got", samples, "prefer more than",400.*dt*CentralEnergy/h)
EnAxis=np.linspace(0.,20.*CentralEnergy,num=samples)
EnInput=np.zeros(samples, dtype=np.complex64)
#for i in range(samples):
EnInput=np.exp(-(EnAxis-CentralEnergy)**2/2./dE**2+2*np.pi*1j*np.random.random(size=samples))
En_FFT=np.fft.fft(EnInput)
TAxis=np.fft.fftfreq(samples,d=(20.*CentralEnergy)/samples)*h
TOutput=np.exp(-TAxis**2/2./dt**2)*En_FFT
EnOutput=np.fft.ifft(TOutput)
if (Axis):
return EnAxis, EnOutput, TAxis, TOutput
else:
return EnOutput, TOutput
# set the main parameters here:
CentralEnergy=80. # in eV
bandwidth=0.5 # bandwidth in %
dt_FWHM=30. # FWHM of the temporal duration on average
dE_FWHM=CentralEnergy/100 *bandwidth # calculate bandwidth of the spectrum in eV
# calculate 3 SASE pulses
EnAxis, EnOutput, TAxis, TOutput = GetSASE(CentralEnergy=CentralEnergy, dE_FWHM=dE_FWHM, dt_FWHM=dt_FWHM)
EnAxis2, EnOutput2, TAxis2, TOutput2 = GetSASE(CentralEnergy=CentralEnergy, dE_FWHM=dE_FWHM, dt_FWHM=dt_FWHM)
EnAxis3, EnOutput3, TAxis3, TOutput3 = GetSASE(CentralEnergy=CentralEnergy, dE_FWHM=dE_FWHM, dt_FWHM=dt_FWHM)
# plot spectrum
ax1 = plt.subplot(1, 2, 1)
plt.plot(EnAxis,np.absolute(EnOutput),EnAxis2,np.absolute(EnOutput2),EnAxis3,np.absolute(EnOutput3) )
plt.xlim(CentralEnergy-2.*dE_FWHM,CentralEnergy+2.*dE_FWHM)
plt.title('Average pulse duration: %.1f fs' % dt_FWHM )
ax1.set_xlabel('Photon energy in eV')
ax1.set_ylabel('spectral intensity')
# plot time structure
ax1 =plt.subplot(1, 2, 2)
plt.plot(TAxis,np.absolute(TOutput),TAxis2,np.absolute(TOutput2), TAxis3,np.absolute(TOutput3))
plt.xlim(-2.*dt_FWHM,+2.*dt_FWHM)
ax1.set_xlabel('time in fs')
ax1.set_ylabel('pulse amplitude')
plt.show()