PyPropogator/testcode.py at main · ELLI-E/PyPropogator · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#first run these dependencies to allow cells below to run
import numpy as np
from matplotlib import pyplot as plt
import resource.pulseshape as pulse
import resource.operators as op
import resource.ramanfunction as rf
from copy import deepcopy

#define numerical parameters
b2 = 0.001e6*1 #in fs^2m^-1
stepsize = 0.025*2 #in meters
gamma = 0.4*4
steps = 100 #z steps - s.t. total length = steps * stepsize in meters
timeNum = 6000 #number of time samples
t1 = -1000 #start time in fs
t2 = 1000 #end time in fs
samplingRate = timeNum/(t2-t1) #samples per femtosecond
t = np.linspace(t1,t2,timeNum)
gaussPulseInitial = pulse.GaussianPulse(t,50) #50fs pulse

#testing generalised GVD version of RK4IP
attenuation = 0
dispersionList = [-b2,0] # should give the same result as above when we only consider beta_2

#to visualize in 2d, save pulse at each z
gaussPulse = op.GeneralGVDRK4IP(dispersionList,attenuation,gamma,stepsize,samplingRate,gaussPulseInitial)
for i in range(steps-1):
    gaussPulse = op.GeneralGVDRK4IP(dispersionList,attenuation,gamma,stepsize,samplingRate,gaussPulse)

initialFT = np.fft.fft(gaussPulseInitial)
propFT = np.fft.fft(gaussPulse)
freqs = np.fft.fftfreq(timeNum,1/samplingRate)
#plot results
plt.plot(t,np.square(gaussPulseInitial),label = "input")
plt.plot(t,np.square(np.abs(gaussPulse)),label = "after propogation")
plt.legend(loc="best")
plt.ylabel(r"$I(T)$")
plt.xlabel(r"$T/T_0$")
plt.figure()
plt.plot(freqs,np.abs(initialFT),label = "input")
plt.plot(freqs,np.abs(propFT),label = "after propogation")
plt.legend(loc="best")
plt.ylabel(r"$I(\nu - \nu_0)$")
plt.xlabel(r"$\nu - \nu_0$")
plt.xlim(-0.2,0.2)
plt.show()
#testing generalised Full GNLSE version of RK4IP
attenuation = 0
dispersionList = [-b2,0] # should give the same result as above when we only consider beta_2
ramanFraction = 0 #0.18
centFrequency = 100 #pulse is normalised so instead define centFrequency as central angular frequency multiplied by duration of pulse - if pulse is NOT normalised use the real angular frequency
ramanCurve = rf.BlowWoodResponse2(t)
#to visualize in 2d, save pulse at each z
gaussPulse = op.GNLSERK4IP(dispersionList,attenuation,gamma,ramanCurve,ramanFraction,centFrequency,stepsize,samplingRate,gaussPulseInitial)
for i in range(steps-1):
    gaussPulse = op.GNLSERK4IP(dispersionList,attenuation,gamma,ramanCurve,ramanFraction,centFrequency,stepsize,samplingRate,gaussPulse)

initialFT = np.fft.fft(gaussPulseInitial)
propFT = np.fft.fft(gaussPulse)
freqs = np.fft.fftfreq(timeNum,1/samplingRate)
#plot results
plt.figure()
plt.plot(t,np.square(gaussPulseInitial),label = "input")
plt.plot(t,np.square(np.abs(gaussPulse)),label = "after propogation")
plt.legend(loc="best")
plt.ylabel(r"$I(T)$")
plt.xlabel(r"$T/T_0$")
plt.figure()
plt.plot(freqs,np.abs(initialFT),label = "input")
plt.plot(freqs,np.abs(propFT),label = "after propogation")
plt.legend(loc="best")
plt.ylabel(r"$I(\nu - \nu_0)$")
plt.xlabel(r"$\nu - \nu_0$")
plt.xlim(-0.2,0.2)
plt.show()