chris/microring-dynamics
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Formatting

Browse Source
This commit is contained in:
Chris Wong 2018-12-17 09:55:58 +08:00
parent 8d65c8622f
commit 6f212468af
6 changed files with 139 additions and 98 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
microring_dynamics.py
View File

@ -7,8 +7,7 @@ K = 0.1
t = np.sqrt(T) + 0.0j t = np.sqrt(T) + 0.0j
k = 0.0 - 1j * np.sqrt(K) k = 0.0 - 1j * np.sqrt(K)
Coupler = np.array([[t, k], Coupler = np.array([[t, k], [k.conjugate(), -t.conjugate()]])
[k.conjugate(), -t.conjugate()]])
num_cycle = 1000 num_cycle = 1000
a = np.zeros((num_cycle, 2), dtype=complex) a = np.zeros((num_cycle, 2), dtype=complex)

9
ref/simple_anim.py
View File

@ -11,22 +11,23 @@ import matplotlib.animation as animation
fig, ax = plt.subplots() fig, ax = plt.subplots()
x = np.arange(0, 2*np.pi, 0.01) x = np.arange(0, 2 * np.pi, 0.01)
line, = ax.plot(x, np.sin(x)) line, = ax.plot(x, np.sin(x))
def init(): # only required for blitting to give a clean slate. def init(): # only required for blitting to give a clean slate.
line.set_ydata([np.nan] * len(x)) line.set_ydata([np.nan] * len(x))
return line, return (line,)
def animate(i): def animate(i):
line.set_ydata(np.sin(x + i / 100)) # update the data. line.set_ydata(np.sin(x + i / 100)) # update the data.
return line, return (line,)
ani = animation.FuncAnimation( ani = animation.FuncAnimation(
fig, animate, init_func=init, interval=2, blit=True, save_count=50) fig, animate, init_func=init, interval=2, blit=True, save_count=50
)
# To save the animation, use e.g. # To save the animation, use e.g.
# #

12
ref/strip_chart.py
View File

@ -21,7 +21,7 @@ class Scope(object):
self.ydata = [0] self.ydata = [0]
self.line = Line2D(self.tdata, self.ydata) self.line = Line2D(self.tdata, self.ydata)
self.ax.add_line(self.line) self.ax.add_line(self.line)
self.ax.set_ylim(-.1, 1.1) self.ax.set_ylim(-0.1, 1.1)
self.ax.set_xlim(0, self.maxt) self.ax.set_xlim(0, self.maxt)
def update(self, y): def update(self, y):
@ -36,18 +36,19 @@ class Scope(object):
self.tdata.append(t) self.tdata.append(t)
self.ydata.append(y) self.ydata.append(y)
self.line.set_data(self.tdata, self.ydata) self.line.set_data(self.tdata, self.ydata)
return self.line, return (self.line,)
def emitter(p=0.03): def emitter(p=0.03):
'return a random value with probability p, else 0' "return a random value with probability p, else 0"
while True: while True:
v = np.random.rand(1) v = np.random.rand(1)
if v > p: if v > p:
yield 0. yield 0.0
else: else:
yield np.random.rand(1) yield np.random.rand(1)
# Fixing random state for reproducibility # Fixing random state for reproducibility
np.random.seed(19680801) np.random.seed(19680801)
@ -56,7 +57,6 @@ fig, ax = plt.subplots()
scope = Scope(ax) scope = Scope(ax)
# pass a generator in "emitter" to produce data for the update func # pass a generator in "emitter" to produce data for the update func
ani = animation.FuncAnimation(fig, scope.update, emitter, interval=10, ani = animation.FuncAnimation(fig, scope.update, emitter, interval=10, blit=True)
blit=True)
plt.show() plt.show()

114
split_step_fourier_method.py
View File

@ -12,58 +12,57 @@ import matplotlib.pyplot as plt
import time import time
# ---Specify input parameters # ---Specify input parameters
distance = 150. # Enter fiber length (in units of L_c)= distance = 150.0 # Enter fiber length (in units of L_c)=
# Normalized 2nd-dispersion: kappa=beta2*f^2*L/2): # Normalized 2nd-dispersion: kappa=beta2*f^2*L/2):
# +ve for normal,-ve for anomalous*) # +ve for normal,-ve for anomalous*)
kappa = -0.001 kappa = -0.001
sigma = 0. # Normalized 3rd-dispersion: sigma=beta3*f^3*L/6 sigma = 0.0 # Normalized 3rd-dispersion: sigma=beta3*f^3*L/6
G = 1. # small signal gain of Ramam Amp: G=g*L G = 1.0 # small signal gain of Ramam Amp: G=g*L
Is = 1. # gain saturation parameter Is = 1.0 # gain saturation parameter
alpha = 0.4 # Normalized fiber amplitude absorption coeff: alpha=l*L alpha = 0.4 # Normalized fiber amplitude absorption coeff: alpha=l*L
# Nonlinear parameter n=') # Nonlinear parameter n=')
# sqrt(L_D/L_NL)=sqrt(gamma*P0*T0^2/|beta2|) or QT: n=kappa^0.5 # sqrt(L_D/L_NL)=sqrt(gamma*P0*T0^2/|beta2|) or QT: n=kappa^0.5
n = 2. ** 0.5 n = 2.0 ** 0.5
# ---Specify filter parameters # ---Specify filter parameters
bdwidth = 2. * np.pi * 6. bdwidth = 2.0 * np.pi * 6.0
delta = 2. * np.pi * 0.5 delta = 2.0 * np.pi * 0.5
a = np.log(np.sqrt(0.95)) a = np.log(np.sqrt(0.95))
perta = 0.3 perta = 0.3
pertfsr = 0.17 pertfsr = 0.17
T = 0.2 T = 0.2
t = np.sqrt(T) t = np.sqrt(T)
r = 1.j * np.sqrt(1. - T) r = 1.0j * np.sqrt(1.0 - T)
# ---Specify input parameters # ---Specify input parameters
mshape = -1. # m=0 for sech,m>0 for super-Gaussian= mshape = -1.0 # m=0 for sech,m>0 for super-Gaussian=
chirp0 = 0. # % input pulse chirp (default value) chirp0 = 0.0 # % input pulse chirp (default value)
# P = 1/(gamma*L); (P is the ref peak power); # P = 1/(gamma*L); (P is the ref peak power);
# uu = A/sqrt(P); # uu = A/sqrt(P);
# z = z0/L_c; (z0 is the real length, L_c is the cavity length); # z = z0/L_c; (z0 is the real length, L_c is the cavity length);
# tau = f*t; (t is the time of reference traveling frame); # tau = f*t; (t is the time of reference traveling frame);
# ---set simulation parameters # ---set simulation parameters
nt = 2 ** 13 # % FFT points (powers of 2) nt = 2 ** 13 # % FFT points (powers of 2)
Tmax = 100. # (half) window size Tmax = 100.0 # (half) window size
stepno = 1 * round(20 * distance * n ** 2) # No.of z steps to stepno = 1 * round(20 * distance * n ** 2) # No.of z steps to
dz = distance / stepno # step size in z dz = distance / stepno # step size in z
dtau = (2. * Tmax) / nt # step size in tau dtau = (2.0 * Tmax) / nt # step size in tau
Twin = 5. Twin = 5.0
fmax = (1. / (2. * Tmax)) * nt / 2. fmax = (1.0 / (2.0 * Tmax)) * nt / 2.0
fwin = 5. fwin = 5.0
filterz = 0.5 filterz = 0.5
plotz = 5 plotz = 5
# ---tau and omega arrays # ---tau and omega arrays
tau = np.arange(-nt / 2., nt / 2.) * dtau # temporal grid tau = np.arange(-nt / 2.0, nt / 2.0) * dtau # temporal grid
# [(0:nt/2-1) (-nt/2:-1)] # [(0:nt/2-1) (-nt/2:-1)]
omega = (np.pi / Tmax) * \ omega = (np.pi / Tmax) * np.append(np.arange(0.0, nt / 2.0), np.arange(-nt / 2.0, 0.0))
np.append(np.arange(0.0, nt / 2.), np.arange(-nt / 2., 0.0))
# frequency grid # frequency grid
delaytau = dtau * np.arange(-round(Twin / dtau), round(Twin / dtau) + 1) delaytau = dtau * np.arange(-round(Twin / dtau), round(Twin / dtau) + 1)
@ -71,21 +70,22 @@ delaytau = dtau * np.arange(-round(Twin / dtau), round(Twin / dtau) + 1)
# Input Field profile # Input Field profile
if mshape == 0: if mshape == 0:
# ;% soliton # ;% soliton
uu = np.exp(-0.5j * chirp0 * tau ** 2.) / np.cosh(tau) uu = np.exp(-0.5j * chirp0 * tau ** 2.0) / np.cosh(tau)
elif mshape > 0: elif mshape > 0:
# super-Gaussian # super-Gaussian
uu = np.exp(-0.5 * (1. + 1.j * chirp0) * tau ** (2. * mshape)) uu = np.exp(-0.5 * (1.0 + 1.0j * chirp0) * tau ** (2.0 * mshape))
else: else:
# White noise # White noise
uu = (np.random.randn(nt) + 1.j * np.random.randn(nt)) * np.sqrt(0.5) uu = (np.random.randn(nt) + 1.0j * np.random.randn(nt)) * np.sqrt(0.5)
# temp = np.fft.fftshift(np.fft.ifft(uu) * (nt * dtau) / np.sqrt(2. * np.pi)) # temp = np.fft.fftshift(np.fft.ifft(uu) * (nt * dtau) / np.sqrt(2. * np.pi))
tempomega = np.fft.fftshift(omega) tempomega = np.fft.fftshift(omega)
# ---store dispersive phase shifts to speedup code # ---store dispersive phase shifts to speedup code
# % nonlinear phase factor # % nonlinear phase factor
dispersion = np.exp((-alpha + 1.j * kappa * omega ** 2. + dispersion = np.exp(
1.j * sigma * omega ** 3.) * dz) (-alpha + 1.0j * kappa * omega ** 2.0 + 1.0j * sigma * omega ** 3.0) * dz
)
# comb filter type # comb filter type
# original comb filter + BPF # original comb filter + BPF
@ -93,23 +93,27 @@ dispersion = np.exp((-alpha + 1.j * kappa * omega ** 2. +
# (1 - r ** 2 * np.exp(-1.j * (omega + delta) + a)) # (1 - r ** 2 * np.exp(-1.j * (omega + delta) + a))
# perturbated comb filter + BPF # perturbated comb filter + BPF
filtert = np.exp(-omega ** 2. / bdwidth ** 2. - filtert = (
perta * np.sin(0.5 * omega / pertfsr) ** 2) * \ np.exp(-omega ** 2.0 / bdwidth ** 2.0 - perta * np.sin(0.5 * omega / pertfsr) ** 2)
(t ** 2) / (1.0 - r ** 2 * np.exp(-1.j * (omega + delta) + a)) * (t ** 2)
/ (1.0 - r ** 2 * np.exp(-1.0j * (omega + delta) + a))
)
fig1 = plt.figure() fig1 = plt.figure()
plt.plot(tempomega / (2. * np.pi), plt.plot(
np.fft.fftshift(10. * np.log10(np.abs(filtert) ** 2.))) tempomega / (2.0 * np.pi), np.fft.fftshift(10.0 * np.log10(np.abs(filtert) ** 2.0))
)
plt.title("Perturbated comb filter + BPF") plt.title("Perturbated comb filter + BPF")
plt.xlim(-fwin, fwin) plt.xlim(-fwin, fwin)
plt.xlim(-fwin, fwin) plt.xlim(-fwin, fwin)
plt.ylim(-30., 10.) plt.ylim(-30.0, 10.0)
plt.show() plt.show()
# %*********[Beginning of MAIN Loop]*********** # %*********[Beginning of MAIN Loop]***********
# % scheme:1/2N\[Rule]D\[Rule]1/2N;first half step nonlinear # % scheme:1/2N\[Rule]D\[Rule]1/2N;first half step nonlinear
temp = uu * np.exp((1.j * np.abs(uu) ** 2. + temp = uu * np.exp(
G / (1.0 + np.abs(uu) ** 2. / Is)) * dz / 2.) (1.0j * np.abs(uu) ** 2.0 + G / (1.0 + np.abs(uu) ** 2.0 / Is)) * dz / 2.0
)
# % note hhz/2 # % note hhz/2
start_time = time.time() start_time = time.time()
@ -136,7 +140,9 @@ for i in range(stepno):
ftemp = np.fft.ifft(temp) * dispersion ftemp = np.fft.ifft(temp) * dispersion
uu = np.fft.fft(ftemp) uu = np.fft.fft(ftemp)
temp = uu * np.exp((1.j * np.abs(uu) ** 2. + G / (1.0 + np.abs(uu) ** 2 / Is)) * dz) temp = uu * np.exp(
(1.0j * np.abs(uu) ** 2.0 + G / (1.0 + np.abs(uu) ** 2 / Is)) * dz
)
z = z + dz z = z + dz
if round((z % plotz) / dz) == 0 or round(((z % plotz) - plotz) / dz) == 0: if round((z % plotz) / dz) == 0 or round(((z % plotz) - plotz) / dz) == 0:
@ -144,7 +150,7 @@ for i in range(stepno):
print("Z: " + str(z)) print("Z: " + str(z))
# fig2.suptitle("i = " + str(i) + ", z = " + str(z) + " Time: " + str(time_used)) # fig2.suptitle("i = " + str(i) + ", z = " + str(z) + " Time: " + str(time_used))
line1.set_data(tau, np.abs(temp) ** 2.) line1.set_data(tau, np.abs(temp) ** 2.0)
ax1.relim() ax1.relim()
ax1.autoscale_view(True, True, True) ax1.autoscale_view(True, True, True)
ax1.set_xlim([-Twin, Twin]) ax1.set_xlim([-Twin, Twin])
@ -152,7 +158,9 @@ for i in range(stepno):
ftemp0 = np.fft.fftshift(ftemp * (nt * dtau) / np.sqrt(2 * np.pi)) ftemp0 = np.fft.fftshift(ftemp * (nt * dtau) / np.sqrt(2 * np.pi))
line2.set_data(tempomega / (2. * np.pi), 10. * np.log10(np.abs(ftemp0) ** 2.)) line2.set_data(
tempomega / (2.0 * np.pi), 10.0 * np.log10(np.abs(ftemp0) ** 2.0)
)
# #
ax2.relim() ax2.relim()
ax2.autoscale_view(True, True, True) ax2.autoscale_view(True, True, True)
@ -165,14 +173,20 @@ for i in range(stepno):
ax3.set_xlim([-Tmax, Tmax]) ax3.set_xlim([-Tmax, Tmax])
ax3.set_title("Time domain (Full Scale)") ax3.set_title("Time domain (Full Scale)")
line4.set_data(tempomega / (2. * np.pi), 10. * np.log10(np.abs(ftemp0) ** 2.)) line4.set_data(
tempomega / (2.0 * np.pi), 10.0 * np.log10(np.abs(ftemp0) ** 2.0)
)
ax4.relim() ax4.relim()
ax4.autoscale_view(True, True, True) ax4.autoscale_view(True, True, True)
ax4.set_xlim([-fmax, fmax]) ax4.set_xlim([-fmax, fmax])
ax4.set_title("Spectrum (Full Scale)") ax4.set_title("Spectrum (Full Scale)")
autocorr0 = np.fft.fftshift( autocorr0 = np.fft.fftshift(
np.fft.ifft(np.fft.fft(np.abs(temp) ** 2) * np.conjugate(np.fft.fft(np.abs(temp) ** 2)))) np.fft.ifft(
np.fft.fft(np.abs(temp) ** 2)
* np.conjugate(np.fft.fft(np.abs(temp) ** 2))
)
)
line5.set_data(tau, np.abs(autocorr0) / max(np.abs(autocorr0))) line5.set_data(tau, np.abs(autocorr0) / max(np.abs(autocorr0)))
ax5.relim() ax5.relim()
@ -187,7 +201,23 @@ plt.show()
# Exporting results # Exporting results
fname = "result" # file name fname = "result" # file name
np.savetxt(fname + "_time.csv", (tau, np.real(uu), np.imag(uu), np.abs(uu) ** 2), delimiter=",") np.savetxt(
np.savetxt(fname + "_freq.csv", (tempomega / (2. * np.pi), np.real(temp), np.imag(temp), fname + "_time.csv", (tau, np.real(uu), np.imag(uu), np.abs(uu) ** 2), delimiter=","
np.abs(temp) ** 2, np.fft.fftshift(np.abs(filtert) ** 2.)), delimiter=",") )
np.savetxt(fname + "_autocorr.csv", (tau, np.abs(autocorr0), np.abs(autocorr0) / max(np.abs(autocorr0))), delimiter=",") np.savetxt(
fname + "_freq.csv",
(
tempomega / (2.0 * np.pi),
np.real(temp),
np.imag(temp),
np.abs(temp) ** 2,
np.fft.fftshift(np.abs(filtert) ** 2.0),
),
delimiter=",",
)
np.savetxt(
fname + "_autocorr.csv",
(tau, np.abs(autocorr0), np.abs(autocorr0) / max(np.abs(autocorr0))),
delimiter=",",
)

96
split_step_fourier_method_v0.py
View File

@ -12,63 +12,62 @@ import matplotlib.pyplot as plt
import time import time
# ---Specify input parameters # ---Specify input parameters
distance = 150. # Enter fiber length (in units of L_c)= distance = 150.0 # Enter fiber length (in units of L_c)=
# Normalized 2nd-dispersion: kappa=beta2*f^2*L/2): # Normalized 2nd-dispersion: kappa=beta2*f^2*L/2):
# +ve for normal,-ve for anomalous*) # +ve for normal,-ve for anomalous*)
kappa = -0.001 kappa = -0.001
sigma = 0. # Normalized 3rd-dispersion: sigma=beta3*f^3*L/6 sigma = 0.0 # Normalized 3rd-dispersion: sigma=beta3*f^3*L/6
G = 1. # small signal gain of Ramam Amp: G=g*L G = 1.0 # small signal gain of Ramam Amp: G=g*L
Is = 1. # gain saturation parameter Is = 1.0 # gain saturation parameter
alpha = 0.4 # Normalized fiber amplitude absorption coeff: alpha=l*L alpha = 0.4 # Normalized fiber amplitude absorption coeff: alpha=l*L
# Nonlinear parameter n=') # Nonlinear parameter n=')
# sqrt(L_D/L_NL)=sqrt(gamma*P0*T0^2/|beta2|) or QT: n=kappa^0.5 # sqrt(L_D/L_NL)=sqrt(gamma*P0*T0^2/|beta2|) or QT: n=kappa^0.5
n = 2. ** 0.5 n = 2.0 ** 0.5
# ---Specify filter parameters # ---Specify filter parameters
bdwidth = 2. * np.pi * 6. bdwidth = 2.0 * np.pi * 6.0
delta = 2. * np.pi * 0.5 delta = 2.0 * np.pi * 0.5
a = np.log(np.sqrt(0.95)) a = np.log(np.sqrt(0.95))
perta = 0.3 perta = 0.3
pertfsr = 0.17 pertfsr = 0.17
T = 0.2 T = 0.2
t = np.sqrt(T) t = np.sqrt(T)
r = 1.j * np.sqrt(1. - T) r = 1.0j * np.sqrt(1.0 - T)
# ---Specify input parameters # ---Specify input parameters
mshape = -1. # m=0 for sech,m>0 for super-Gaussian= mshape = -1.0 # m=0 for sech,m>0 for super-Gaussian=
chirp0 = 0. # % input pulse chirp (default value) chirp0 = 0.0 # % input pulse chirp (default value)
# P = 1/(gamma*L); (P is the ref peak power); # P = 1/(gamma*L); (P is the ref peak power);
# uu = A/sqrt(P); # uu = A/sqrt(P);
# z = z0/L_c; (z0 is the real length, L_c is the cavity length); # z = z0/L_c; (z0 is the real length, L_c is the cavity length);
# tau = f*t; (t is the time of reference traveling frame); # tau = f*t; (t is the time of reference traveling frame);
# ---set simulation parameters # ---set simulation parameters
nt = 2 ** 13 # % FFT points (powers of 2) nt = 2 ** 13 # % FFT points (powers of 2)
Tmax = 100. # (half) window size Tmax = 100.0 # (half) window size
stepno = 1 * round(20 * distance * n ** 2) # No.of z steps to stepno = 1 * round(20 * distance * n ** 2) # No.of z steps to
dz = distance / stepno # step size in z dz = distance / stepno # step size in z
dtau = (2. * Tmax) / nt # step size in tau dtau = (2.0 * Tmax) / nt # step size in tau
Twin = 5. Twin = 5.0
fmax = (1. / (2. * Tmax)) * nt / 2. fmax = (1.0 / (2.0 * Tmax)) * nt / 2.0
fwin = 5. fwin = 5.0
filterz = 0.5 filterz = 0.5
plotz = 5 plotz = 5
# ---tau and omega arrays # ---tau and omega arrays
tau = np.arange(-nt / 2., nt / 2.) * dtau # temporal grid tau = np.arange(-nt / 2.0, nt / 2.0) * dtau # temporal grid
# [(0:nt/2-1) (-nt/2:-1)] # [(0:nt/2-1) (-nt/2:-1)]
omega = (np.pi / Tmax) * \ omega = (np.pi / Tmax) * np.append(np.arange(0.0, nt / 2.0), np.arange(-nt / 2.0, 0.0))
np.append(np.arange(0.0 , nt / 2.), np.arange(-nt / 2., 0.0))
# frequency grid # frequency grid
@ -78,23 +77,24 @@ delaytau = dtau * np.arange(-round(Twin / dtau), round(Twin / dtau) + 1)
# Input Field profile # Input Field profile
if mshape == 0: if mshape == 0:
# ;% soliton # ;% soliton
uu = np.exp(-0.5j * chirp0 * tau ** 2.) / np.cosh(tau) uu = np.exp(-0.5j * chirp0 * tau ** 2.0) / np.cosh(tau)
elif mshape > 0: elif mshape > 0:
# super-Gaussian # super-Gaussian
uu = np.exp(-0.5 * (1. + 1.j * chirp0) * tau ** (2. * mshape)) uu = np.exp(-0.5 * (1.0 + 1.0j * chirp0) * tau ** (2.0 * mshape))
else: else:
# White noise # White noise
uu = (np.random.randn(nt) + 1.j * np.random.randn(nt)) * np.sqrt(0.5) uu = (np.random.randn(nt) + 1.0j * np.random.randn(nt)) * np.sqrt(0.5)
temp = np.fft.fftshift(np.fft.ifft(uu) * (nt * dtau) / np.sqrt(2.* np.pi)) temp = np.fft.fftshift(np.fft.ifft(uu) * (nt * dtau) / np.sqrt(2.0 * np.pi))
tempomega = np.fft.fftshift(omega) tempomega = np.fft.fftshift(omega)
# ---store dispersive phase shifts to speedup code # ---store dispersive phase shifts to speedup code
# % nonlinear phase factor # % nonlinear phase factor
dispersion = np.exp((-alpha + 1.j * kappa * omega **2. + \ dispersion = np.exp(
1.j * sigma * omega ** 3.) * dz) (-alpha + 1.0j * kappa * omega ** 2.0 + 1.0j * sigma * omega ** 3.0) * dz
)
# comb filter type # comb filter type
# original comb filter + BPF # original comb filter + BPF
@ -102,22 +102,26 @@ dispersion = np.exp((-alpha + 1.j * kappa * omega **2. + \
# (1 - r ** 2 * np.exp(-1.j * (omega + delta) + a)) # (1 - r ** 2 * np.exp(-1.j * (omega + delta) + a))
# perturbated comb filter + BPF # perturbated comb filter + BPF
filtert = np.exp(-omega ** 2. / bdwidth ** 2. - \ filtert = (
perta * np.sin(0.5 * omega / pertfsr) ** 2) * \ np.exp(-omega ** 2.0 / bdwidth ** 2.0 - perta * np.sin(0.5 * omega / pertfsr) ** 2)
(t ** 2) / (1.0 - r ** 2 * np.exp(-1.j*(omega + delta) + a)) * (t ** 2)
/ (1.0 - r ** 2 * np.exp(-1.0j * (omega + delta) + a))
)
plt.figure() plt.figure()
plt.plot(tempomega / (2. * np.pi), plt.plot(
np.fft.fftshift(10. * np.log10(np.abs(filtert) ** 2.))) tempomega / (2.0 * np.pi), np.fft.fftshift(10.0 * np.log10(np.abs(filtert) ** 2.0))
)
plt.xlim(-fwin, fwin) plt.xlim(-fwin, fwin)
plt.xlim(-fwin, fwin) plt.xlim(-fwin, fwin)
plt.ylim(-30., 10.) plt.ylim(-30.0, 10.0)
plt.show() plt.show()
# %*********[Beginning of MAIN Loop]*********** # %*********[Beginning of MAIN Loop]***********
# % scheme:1/2N\[Rule]D\[Rule]1/2N;first half step nonlinear # % scheme:1/2N\[Rule]D\[Rule]1/2N;first half step nonlinear
temp = uu * np.exp((1.j * np.abs(uu) ** 2. + temp = uu * np.exp(
G / (1.0 + np.abs(uu) ** 2. / Is )) * dz / 2.) (1.0j * np.abs(uu) ** 2.0 + G / (1.0 + np.abs(uu) ** 2.0 / Is)) * dz / 2.0
)
# % note hhz/2 # % note hhz/2
start_time = time.time() start_time = time.time()
@ -131,13 +135,15 @@ z = 0
plt.ion() plt.ion()
fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2) fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2)
for i in range(stepno): for i in range(stepno):
if round((z % 1 - filterz)/dz) == 0: if round((z % 1 - filterz) / dz) == 0:
ftemp = np.fft.ifft(temp) * filtert * dispersion ftemp = np.fft.ifft(temp) * filtert * dispersion
else: else:
ftemp = np.fft.ifft(temp) * dispersion ftemp = np.fft.ifft(temp) * dispersion
uu = np.fft.fft(ftemp) uu = np.fft.fft(ftemp)
temp = uu*np.exp((1.j * np.abs(uu)**2. + G / (1.0 + np.abs(uu)**2 / Is)) * dz) temp = uu * np.exp(
(1.0j * np.abs(uu) ** 2.0 + G / (1.0 + np.abs(uu) ** 2 / Is)) * dz
)
z = z + dz z = z + dz
if round((z % plotz) / dz) == 0 or round(((z % plotz) - plotz) / dz) == 0: if round((z % plotz) / dz) == 0 or round(((z % plotz) - plotz) / dz) == 0:
@ -145,29 +151,33 @@ for i in range(stepno):
print(i) print(i)
ax1.clear() ax1.clear()
ax1.plot(tau, np.abs(temp)**2.) ax1.plot(tau, np.abs(temp) ** 2.0)
ax1.set_xlim([-Twin, Twin]) ax1.set_xlim([-Twin, Twin])
ax1.set_title("i = " + str(i) + ", z = " + str(z)) ax1.set_title("i = " + str(i) + ", z = " + str(z))
ftemp0 = np.fft.fftshift(ftemp * (nt * dtau) / np.sqrt(2 * np.pi)) ftemp0 = np.fft.fftshift(ftemp * (nt * dtau) / np.sqrt(2 * np.pi))
ax2.clear() ax2.clear()
ax2.plot(tempomega / (2. * np.pi), 10. * np.log10(np.abs(ftemp0)**2.)) ax2.plot(tempomega / (2.0 * np.pi), 10.0 * np.log10(np.abs(ftemp0) ** 2.0))
ax2.set_xlim([-fwin, fwin]) ax2.set_xlim([-fwin, fwin])
ax2.set_title(" Time = " + str(time_used)) ax2.set_title(" Time = " + str(time_used))
ax3.clear() ax3.clear()
ax3.plot(tau, np.abs(temp)**2) ax3.plot(tau, np.abs(temp) ** 2)
ax3.set_xlim([-Tmax, Tmax]) ax3.set_xlim([-Tmax, Tmax])
ax3.set_title("i = " + str(i) + ", z = " + str(z)) ax3.set_title("i = " + str(i) + ", z = " + str(z))
ax4.clear() ax4.clear()
ax4.plot(tempomega / (2. * np.pi), 10. * np.log10(np.abs(ftemp0)**2.)) ax4.plot(tempomega / (2.0 * np.pi), 10.0 * np.log10(np.abs(ftemp0) ** 2.0))
ax4.set_xlim([-fmax, fmax]) ax4.set_xlim([-fmax, fmax])
ax4.set_title(" Time = " + str(time_used)) ax4.set_title(" Time = " + str(time_used))
autocorr0 = np.fft.fftshift( autocorr0 = np.fft.fftshift(
np.fft.ifft(np.fft.fft(np.abs(temp)**2) * np.conjugate(np.fft.fft(np.abs(temp)**2)))) np.fft.ifft(
np.fft.fft(np.abs(temp) ** 2)
* np.conjugate(np.fft.fft(np.abs(temp) ** 2))
)
)
ax5.clear() ax5.clear()
ax5.plot(tau, autocorr0 / max(autocorr0)) ax5.plot(tau, autocorr0 / max(autocorr0))

3
tmp.py
View File

@ -1 +1,2 @@
print(1.3%1.0) print(1.3 % 1.0)