diff --git a/examples/test_wl_del_freq.py b/examples/test_wl_del_freq.py new file mode 100644 index 0000000..af2534e --- /dev/null +++ b/examples/test_wl_del_freq.py @@ -0,0 +1,103 @@ +""" +This example is here to test new plot functions. This test is based on file +test_dispersion: +Example: Example of supercontinuum generation in anomalous dispersion regime +at a central wavelength of 835 nm in a 15 centimeter long fiber. +Comparision of results obtained with two dispersion input: +1. dispersion calculated from Taylor expansion +2. dispersion calculated from effective refractive indicies + +The python code based on MATLAB code published in +'Supercontinuum Generation in Optical Fibers' +by J. M. Dudley and J. R. Taylor, available at http://scgbook.info/. +""" + +import numpy as np +import matplotlib.pyplot as plt + +import gnlse + +import os + +if __name__ == '__main__': + setup = gnlse.GNLSESetup() + + # Numerical parameters + setup.resolution = 2**14 + setup.time_window = 12.5 # ps + setup.z_saves = 200 + + # Physical parameters + setup.wavelength = 835 # nm + setup.fiber_length = 0.15 # m + setup.nonlinearity = 0.11 # 1/W/m + setup.raman_model = gnlse.raman_blowwood + setup.self_steepening = True + + # The dispersion model is built from a Taylor expansion with coefficients + # given below. + loss = 0 + betas = np.array([ + -11.830e-3, 8.1038e-5, -9.5205e-8, 2.0737e-10, -5.3943e-13, 1.3486e-15, + -2.5495e-18, 3.0524e-21, -1.7140e-24 + ]) + setup.dispersion_model = gnlse.DispersionFiberFromTaylor(loss, betas) + + # Input pulse parameters + power = 10000 + # pulse duration [ps] + tfwhm = 0.05 + # hyperbolic secant + setup.pulse_model = gnlse.SechEnvelope(power, tfwhm) + solver = gnlse.GNLSE(setup) + solution = solver.run() + + plt.figure(figsize=(14, 8), facecolor='w', edgecolor='k') + + plt.subplot(4, 3, 1) + gnlse.plot_delay_vs_distance(solution, time_range=[-.5, 5], cmap="jet") + + plt.subplot(4, 3, 2) + gnlse.plot_frequency_vs_distance(solution, frequency_range=[-300, 200], + cmap="plasma") + + plt.subplot(4, 3, 3) + gnlse.plot_wavelength_vs_distance(solution, WL_range=[400, 1400]) + + plt.subplot(4, 3, 4) + gnlse.plot_delay_vs_distance_logarithmic(solution, time_range=[-.5, 5], + cmap="jet") + + plt.subplot(4, 3, 5) + gnlse.plot_frequency_vs_distance_logarithmic(solution, + frequency_range=[-300, 200], + cmap="plasma") + + plt.subplot(4, 3, 6) + gnlse.plot_wavelength_vs_distance_logarithmic(solution, + WL_range=[400, 1400]) + + plt.subplot(4, 3, 7) + gnlse.plot_delay_for_distance_slice(solution, time_range=[-.5, 5]) + + plt.subplot(4, 3, 8) + gnlse.plot_frequency_for_distance_slice(solution, + frequency_range=[-300, 200]) + + plt.subplot(4, 3, 9) + gnlse.plot_wavelength_for_distance_slice(solution, WL_range=[400, 1400]) + + plt.subplot(4, 3, 10) + ax = gnlse.plot_delay_for_distance_slice_logarithmic(solution, + time_range=[-.5, 5]) + + plt.subplot(4, 3, 11) + gnlse.plot_frequency_for_distance_slice_logarithmic(solution, + frequency_range=[-300, 200]) + + plt.subplot(4, 3, 12) + gnlse.plot_wavelength_for_distance_slice_logarithmic(solution, + WL_range=[400, 1400]) + + plt.tight_layout() + plt.show() diff --git a/gnlse/visualization.py b/gnlse/visualization.py index 34c6e7e..f368698 100644 --- a/gnlse/visualization.py +++ b/gnlse/visualization.py @@ -10,7 +10,8 @@ from gnlse.common import c -def plot_frequency_vs_distance_logarithmic(solver, ax=None, norm=None): +def plot_frequency_vs_distance_logarithmic(solver, ax=None, norm=None, + frequency_range=None, cmap="magma"): """Plotting results in logarithmic scale in frequency domain. Parameters @@ -42,7 +43,15 @@ def plot_frequency_vs_distance_logarithmic(solver, ax=None, norm=None): where=(np.abs(solver.AW)**2 > 0))) frequency = (solver.W - solver.w_0) / 2 / np.pi # frequency grid - ax.imshow(lIW, origin='lower', aspect='auto', cmap="magma", + if frequency_range is not None: + iis = np.logical_and(frequency >= frequency_range[0], + frequency <= frequency_range[1]) + # indices of interest + + frequency= frequency[iis] + lIW = lIW[:, iis] + + ax.imshow(lIW, origin='lower', aspect='auto', cmap=cmap, extent=[np.min(frequency), np.max(frequency), 0, np.max(solver.Z)], vmin=-40) ax.set_xlabel("Frequency [THz]") @@ -50,7 +59,8 @@ def plot_frequency_vs_distance_logarithmic(solver, ax=None, norm=None): return ax -def plot_frequency_vs_distance(solver, ax=None, norm=None): +def plot_frequency_vs_distance(solver, frequency_range=None, + ax=None, norm=None, cmap="magma"): """Plotting results in frequency domain. Linear scale. Parameters @@ -81,7 +91,15 @@ def plot_frequency_vs_distance(solver, ax=None, norm=None): np.abs(solver.AW)**2 / norm) frequency = (solver.W - solver.w_0) / 2 / np.pi # frequency grid - ax.imshow(IW, origin='lower', aspect='auto', cmap="magma", + if frequency_range is not None: + iis = np.logical_and(frequency >= frequency_range[0], + frequency <= frequency_range[1]) + # indices of interest + + frequency= frequency[iis] + IW = IW[:, iis] + + ax.imshow(IW, origin='lower', aspect='auto', cmap=cmap, extent=[np.min(frequency), np.max(frequency), 0, np.max(solver.Z)], vmin=0) ax.set_xlabel("Frequency [THz]") @@ -118,10 +136,9 @@ def plot_delay_for_distance_slice(solver, time_range=None, ax=None, time_range = [np.min(solver.t), np.max(solver.t)] if norm is None: - norm = np.max(np.abs(solver.At)**2) + norm = np.max(np.abs(solver.At[0][:])**2) - It = np.fliplr( - np.abs(solver.At)**2 / norm) + It = np.abs(solver.At)**2 / norm # indices of interest if no z_slice positions were given if z_slice is None: @@ -143,9 +160,138 @@ def plot_delay_for_distance_slice(solver, time_range=None, ax=None, return ax +def plot_wavelength_for_distance_slice(solver, WL_range=None, ax=None, + z_slice=None, norm=None): + """Plotting chosen slices of intensity + in linear scale in wavelength domain. + + Parameters + ---------- + solver : Solution + Model outputs in the form of a ``Solution`` object. + WL_range : list, (2, ) + Wavelength range. Set [400, 1350] as default. + ax : :class:`~matplotlib.axes.Axes` + :class:`~matplotlib.axes.Axes` instance for plotting + norm : float + Normalization factor for output spectrum. As default maximum of + square absolute of ``solver.AW`` variable is taken. + + Returns + ------- + ax : :class:`~matplotlib.axes.Axes` + Used :class:`~matplotlib.axes.Axes` instance. + """ + + if ax is None: + ax = plt.gca() + + if WL_range is None: + WL_range = [np.min(c / (solver.W / 2 / np.pi)), + np.max(c / (solver.W / 2 / np.pi))] + + if norm is None: + norm = np.max(np.abs(solver.AW)**2) + + IW = np.fliplr( + np.abs(solver.AW)**2 / norm) + WL = 2 * np.pi * c / solver.W # wavelength grid + WL_asc = np.flip(WL, ) # ascending order for interpolation + iio = np.logical_and(WL_asc > WL_range[0], + WL_asc < WL_range[1]) # indices in order + + WL_asc = WL_asc[iio] + IW = IW[:, iio] + + # indices of interest if no z_slice positions were given + if z_slice is None: + iis = [0, -1] + # indices of interest nearest to given z_slice positions + else: + iis = [np.nonzero( + np.min(np.abs(solver.Z - z)) == np.abs(solver.Z - z) + )[0][0] for z in z_slice] + + for i in iis: + label_i = "z = " + str(solver.Z[i]) + "m" + ax.plot(WL_asc, IW[i][:], label=label_i) + + ax.set_xlim([np.min(WL_asc), np.max(WL_asc)]) + ax.set_xlabel("Wavelength [nm]") + ax.set_ylabel("Normalized Spectral Density") + ax.legend() + return ax + + +def plot_wavelength_for_distance_slice_logarithmic(solver, WL_range=None, + ax=None, + z_slice=None, norm=None): + """Plotting chosen slices of intensity + in linear scale in wavelength domain. + + Parameters + ---------- + solver : Solution + Model outputs in the form of a ``Solution`` object. + WL_range : list, (2, ) + Wavelength range. Set [400, 1350] as default. + ax : :class:`~matplotlib.axes.Axes` + :class:`~matplotlib.axes.Axes` instance for plotting + norm : float + Normalization factor for output spectrum. As default maximum of + square absolute of ``solver.AW`` variable is taken. + + Returns + ------- + ax : :class:`~matplotlib.axes.Axes` + Used :class:`~matplotlib.axes.Axes` instance. + """ + + if ax is None: + ax = plt.gca() + + if WL_range is None: + WL_range = [np.min(c / (solver.W / 2 / np.pi)), + np.max(c / (solver.W / 2 / np.pi))] + + if norm is None: + norm = np.max(np.abs(solver.AW)**2) + + lIW = np.fliplr( + 10 * np.log10(np.abs(solver.AW)**2 / norm, + where=(np.abs(solver.AW)**2 > 0))) + WL = 2 * np.pi * c / solver.W # wavelength grid + WL_asc = np.flip(WL, ) # ascending order for interpolation + iio = np.logical_and(WL_asc > WL_range[0], + WL_asc < WL_range[1]) # indices in order + + WL_asc = WL_asc[iio] + lIW = lIW[:, iio] + + # indices of interest if no z_slice positions were given + if z_slice is None: + iis = [0, -1] + # indices of interest nearest to given z_slice positions + else: + iis = [np.nonzero( + np.min(np.abs(solver.Z - z)) == np.abs(solver.Z - z) + )[0][0] for z in z_slice] + + for i in iis: + label_i = "z = " + str(solver.Z[i]) + "m" + ax.plot(WL_asc, lIW[i][:], label=label_i) + + ax.set_xlim([np.min(WL_asc), np.max(WL_asc)]) + ax.set_ylim(-40) + ax.set_xlabel("Wavelength [nm]") + ax.set_ylabel("Normalized Spectral Density") + ax.legend() + return ax + + def plot_delay_for_distance_slice_logarithmic(solver, time_range=None, ax=None, z_slice=None, norm=None): - """Plotting intensity in logarithmic scale in time domain. + """Plotting chosen slices of intensity in linear scale in time domain. Parameters ---------- @@ -248,7 +394,7 @@ def plot_frequency_for_distance_slice(solver, frequency_range=None, ax=None, ax.set_xlim(frequency_range) ax.set_xlabel("Frequency [Thz]") - ax.set_ylabel("Normalized Power") + ax.set_ylabel("Normalized Spectral Density") ax.legend() return ax @@ -306,14 +452,15 @@ def plot_frequency_for_distance_slice_logarithmic(solver, frequency_range=None, ax.plot((solver.W - solver.w_0) / 2 / np.pi, lIW[i][:], label=label_i) ax.set_xlim(frequency_range) + ax.set_ylim(-40) ax.set_xlabel("Frequency [Thz]") - ax.set_ylabel("Normalized Power") + ax.set_ylabel("Normalized Spectral Density") ax.legend() return ax def plot_delay_vs_distance_logarithmic(solver, time_range=None, ax=None, - norm=None): + norm=None, cmap="magma"): """Plotting intensity in logarithmic scale in time domain. Parameters @@ -346,14 +493,15 @@ def plot_delay_vs_distance_logarithmic(solver, time_range=None, ax=None, where=(np.abs(solver.At)**2 > 0)) ax.pcolormesh(solver.t, solver.Z, lIT, shading="auto", vmin=-40, - cmap="magma") + cmap=cmap) ax.set_xlim(time_range) ax.set_xlabel("Delay [ps]") ax.set_ylabel("Distance [m]") return ax -def plot_delay_vs_distance(solver, time_range=None, ax=None, norm=None): +def plot_delay_vs_distance(solver, time_range=None, ax=None, norm=None, + cmap="magma"): """Plotting normalized intensity in linear scale in time domain. Parameters @@ -385,7 +533,7 @@ def plot_delay_vs_distance(solver, time_range=None, ax=None, norm=None): lIT = np.abs(solver.At)**2 / norm ax.pcolormesh(solver.t, solver.Z, lIT, shading="auto", vmin=0, - cmap="jet") + cmap=cmap) ax.set_xlim(time_range) ax.set_xlabel("Delay [ps]") ax.set_ylabel("Distance [m]") @@ -393,7 +541,7 @@ def plot_delay_vs_distance(solver, time_range=None, ax=None, norm=None): def plot_wavelength_vs_distance(solver, WL_range=None, ax=None, - norm=None): + norm=None, cmap="magma"): """Plotting results in linear scale in wavelength domain. Parameters @@ -438,7 +586,7 @@ def plot_wavelength_vs_distance(solver, WL_range=None, ax=None, newWL = np.linspace(np.min(WL_asc), np.max(WL_asc), IW.shape[1]) toshow = interpolator(newWL, solver.Z) - ax.imshow(toshow, origin='lower', aspect='auto', cmap="jet", + ax.imshow(toshow, origin='lower', aspect='auto', cmap=cmap, extent=[np.min(WL_asc), np.max(WL_asc), 0, np.max(solver.Z)], vmin=0) ax.set_xlabel("Wavelength [nm]") @@ -447,7 +595,7 @@ def plot_wavelength_vs_distance(solver, WL_range=None, ax=None, def plot_wavelength_vs_distance_logarithmic(solver, WL_range=None, - ax=None, norm=None): + ax=None, norm=None, cmap="magma"): """Plotting results in logarithmic scale in wavelength domain. Parameters @@ -493,7 +641,7 @@ def plot_wavelength_vs_distance_logarithmic(solver, WL_range=None, newWL = np.linspace(np.min(WL_asc), np.max(WL_asc), lIW.shape[1]) toshow = interpolator(newWL, solver.Z) - ax.imshow(toshow, origin='lower', aspect='auto', cmap="magma", + ax.imshow(toshow, origin='lower', aspect='auto', cmap=cmap, extent=[np.min(WL_asc), np.max(WL_asc), 0, np.max(solver.Z)], vmin=-40) ax.set_xlabel("Wavelength [nm]")