#* **************************************************************************
#*
#* CDTK, Chemical Dynamics Toolkit
#* A modular system for chemical dynamics applications and more
#*
#* Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016
#* Oriol Vendrell, DESY, <oriol.vendrell@desy.de>
#*
#* Copyright (C) 2017, 2018, 2019
#* Ralph Welsch, DESY, <ralph.welsch@desy.de>
#*
#* Copyright (C) 2020, 2021, 2022, 2023
#* Ludger Inhester, DESY, ludger.inhester@cfel.de
#*
#* This file is part of CDTK.
#*
#* CDTK is free software: you can redistribute it and/or modify
#* it under the terms of the GNU General Public License as published by
#* the Free Software Foundation, either version 3 of the License, or
#* (at your option) any later version.
#*
#* This program is distributed in the hope that it will be useful,
#* but WITHOUT ANY WARRANTY; without even the implied warranty of
#* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#* GNU General Public License for more details.
#*
#* You should have received a copy of the GNU General Public License
#* along with this program. If not, see <http://www.gnu.org/licenses/>.
#*
#* **************************************************************************
import sys
import math
import numpy as np
import CDTK.Tools.Conversion as cv
[docs]
class EField(object):
"""
Represents an oscillating field characterized by
E(t) = P A(t)
P: unit 3-vector, polarizarion axis of the field
A(t): function, field amplitude as a function of time
P and A must have been defined at execution time.
Units for A(t) must be au both for time and E-field.
"""
def __init__(self,**opts):
self.P = np.array([0.0,0.0,1.0],float)
[docs]
def A(self,t):
return 1.0
[docs]
def Envelope(self, t):
return 1.0
[docs]
def E(self,t):
return self.A(t) * np.asarray(self.P)
[docs]
def set_gaussian_pulse(self,omega,fwhm,amp,t0=0,phase=0,fluence=False):
"""
Gaussian shaped pulse with single photon energy
A(t) = amp * G(fwhm,t0) * cos( omega t + phase)
All input in au
omega -- photon energy
fwhm -- full width at half maximum
amp -- maximum field amplitude
t0 -- center of the pulse
phase -- phase of the cosine carrier, units of Pi
fluence -- if amp is given as fluence, then normalized gaussians have to be used
No output, modifies self.A and self.Envelope
"""
g = gaussian_envelope(fwhm,t0,fluence)
p = phase * math.pi
def a(t):
return g(t) * amp * np.cos( omega*t + p )
def env(t):
return g(t) * amp
self.A = a
self.Envelope = env
[docs]
def set_double_gaussian_pulse(self,omega,fwhm, fwhm2, ratio2, amp, t1=0., t2=0., phase=0,fluence=False):
"""
Gaussian shaped pulse with single photon energy
A(t) = amp * ((1-ratio2) * G(fwhm1,t0) + ratio2 * F(fwhm2,t0)) * cos( omega t + phase)
All input in au
omega -- photon energy
fwhm -- full width at half maximum
fwhm2 -- full width at half maximum
ratio2 -- fluence ratio of second pulse part
amp -- maximum field amplitude
t1 -- center of pulse1
t2 -- center of pulse2
phase -- phase of the cosine carrier, units of Pi
fluence -- if amp is given as fluence, then normalized gaussians have to be used
No output, modifies self.A and self.Envelope
"""
g1 = gaussian_envelope(fwhm, t1, fluence)
g2 = gaussian_envelope(fwhm2, t2, fluence)
amp1 = (1. - ratio2) * amp
amp2 = ratio2 * amp
p = phase * math.pi
def a(t):
return (amp1 * g1(t) + amp2 * g2(t)) * np.cos( omega*t + p )
def env(t):
return (amp1 * g1(t) + amp2 * g2(t))
self.A = a
self.Envelope = env
[docs]
def writeLog(self, times, filename='field.log'):
amplitudes = self.A(times)
envelope = self.Envelope(times)
f = open(filename, 'w')
for i in range(times.shape[0]):
f.write('%9.3f %12.6f %12.6f \n' % (times[i],
amplitudes[i],
envelope[i])
)
f.flush()
f.close()
[docs]
def plot(self,ti=0.0,tf=10.0,**opts):
"""
Use matplotlib to plot the time-profile of the pulse
ti -- float, intial time
tf -- float, final time
"""
try:
import matplotlib.pyplot as plt
except ImportError:
message = 'matplotlib module missing. plotting is not possible.\n'
sys.stdout.write(message)
tunit = opts.get('tunit','au') # time units
funit = opts.get('funit','au') # electric field units
npoints = opts.get('npoints',500) # number of points used for plotting
color = opts.get('color','r')
t = np.linspace(ti,tf,npoints)
au_ti = ti
au_tf = tf
if tunit == 'fs':
au_ti = ti*cv.fs2au
au_tf = tf*cv.fs2au
au_t = t*cv.fs2au
else:
au_ti = ti
au_tf = tf
au_t = t
amp = list( map(self.A,au_t) )
plt.plot(t,amp,color)
plt.show()
__k = 4.0*math.log(2.0)
__N = __k/math.pi
[docs]
def gaussian_envelope(fwhm,t0=0.0,norm=False):
"""
Return Gaussian function f(t)
fwhm -- float; full width at half maximum
t0 -- float; center of envelope; Def t0 = 0.0
norm -- bool; True: integral = 1, False: maximum = 1. Def. norm=False
"""
# N = __N/fwhm
N = 1.0/math.sqrt(2.0*math.pi*fwhm**2/(2.0*__k))
if norm:
def g(t):
return N*np.exp(-__k*((t-t0)/fwhm)**2)
else:
def g(t):
return np.exp(-__k*((t-t0)/fwhm)**2)
return g
if __name__ == '__main__':
ef = EField()
ef.P = np.array([0,0,1],float)
ef.set_gaussian_pulse(1.5*cv.ev2au, 0.05*cv.fs2au, 0.5, 2.8*cv.fs2au, 1.5)
ef.plot(ti=0,tf=10,tunit='fs')