gn model introduced in the science utils module

gnpy/core/science_utils.py

@@ -1,3 +1,8 @@
+import progressbar
+import numpy as np
+from operator import attrgetter
+from scipy.interpolate import interp1d
+
 from gnpy.core.utils import load_json
 import scipy.constants as ph
 from scipy.integrate import solve_bvp
@@ -339,4 +344,113 @@ class RamanSolver:
                 dpdz_element = prop_direct[f_ind] * (-alphap_fiber[f_ind] + raman_gain - raman_loss) * power_sample
                 dpdz[f_ind][z_ind] = dpdz_element
 
-        return np.vstack(dpdz)
+        return np.vstack(dpdz)
+
+class NLI:
+    """ This class implements the NLI models.
+        Model and method can be specified in `self.model_parameters.method`.
+        List of implemented methods:
+        'gn_analytic': brute force triple integral solution
+        'GGN_spectrally_separated_xpm_spm': XPM plus SPM
+    """
+
+    def __init__(self, fiber_information=None):
+        """ Initialize the fiber object with its physical parameters
+        """
+        self.fiber_information = fiber_information
+        self.srs_profile = None
+        self.model_parameters = None
+
+    @property
+    def fiber_information(self):
+        return self.__fiber_information
+
+    @fiber_information.setter
+    def fiber_information(self, fiber_information):
+        self.__fiber_information = fiber_information
+
+    @property
+    def srs_profile(self):
+        return self.__srs_profile
+
+    @srs_profile.setter
+    def srs_profile(self, srs_profile):
+        self.__srs_profile = srs_profile
+
+    @property
+    def model_parameters(self):
+        return self.__model_parameters
+
+    @model_parameters.setter
+    def model_parameters(self, model_params):
+        """
+        :param model_params: namedtuple containing the parameters used to compute the NLI.
+        """
+        self.__model_parameters = model_params
+
+    def alpha0(self, f_eval=193.5e12):
+        if len(self.fiber_information.attenuation_coefficient.alpha_power) == 1:
+            alpha0 = self.fiber_information.attenuation_coefficient.alpha_power[0]
+        else:
+            alpha_interp = interp1d(self.fiber_information.attenuation_coefficient.frequency,
+                                    self.fiber_information.attenuation_coefficient.alpha_power)
+            alpha0 = alpha_interp(f_eval)
+        return alpha0
+
+    def compute_nli(self, carrier, *carriers):
+        """ Compute NLI power generated by the WDM comb `*carriers` on the channel under test `carrier`
+        at the end of the fiber span.
+        """
+        if 'gn_model_analytic' == self.model_parameters.method.lower():
+            carrier_nli = self._gn_analytic(carrier, *carriers)
+        else:
+            raise ValueError(f'Method {self.model_parameters.method_nli} not implemented.')
+
+        return carrier_nli
+
+    # Methods for computing spectrally separated GN
+    def _gn_analytic(self, carrier, *carriers):
+        """ Computes the nonlinear interference power on a single carrier.
+        The method uses eq. 120 from arXiv:1209.0394.
+        :param carrier: the signal under analysis
+        :param carriers: the full WDM comb
+        :return: carrier_nli: the amount of nonlinear interference in W on the under analysis
+        """
+
+        alpha = self.alpha0() / 2
+        length = self.fiber_information.length
+        effective_length = (1 - np.exp(-2 * alpha * length)) / (2 * alpha)
+        asymptotic_length = 1 / (2 * alpha)
+
+        beta2 = self.fiber_information.beta2
+        gamma = self.fiber_information.gamma
+
+        g_nli = 0
+        for interfering_carrier in carriers:
+            g_interfearing = interfering_carrier.power.signal / interfering_carrier.baud_rate
+            g_signal = carrier.power.signal / carrier.baud_rate
+            g_nli += g_interfearing**2 * g_signal * self._psi(carrier, interfering_carrier)
+        g_nli *= (16.0 / 27.0) * (gamma * effective_length)**2 /\
+                 (2 * np.pi * abs(beta2) * asymptotic_length)
+
+        carrier_nli = carrier.baud_rate * g_nli
+        return carrier_nli
+
+    def _psi(self, carrier, interfering_carrier):
+        """ Calculates eq. 123 from arXiv:1209.0394.
+        """
+        alpha = self.alpha0() / 2
+        beta2 = self.fiber_information.beta2
+
+        asymptotic_length = 1 / (2 * alpha)
+        if carrier.channel_number == interfering_carrier.channel_number:  # SPM
+            psi = np.arcsinh(0.5 * np.pi**2 * asymptotic_length
+                              * abs(beta2) * carrier.baud_rate**2)
+        else:  # XPM
+            delta_f = carrier.frequency - interfering_carrier.frequency
+            psi = np.arcsinh(np.pi**2 * asymptotic_length * abs(beta2) *
+                             carrier.baud_rate * (delta_f + 0.5 * interfering_carrier.baud_rate))
+            psi -= np.arcsinh(np.pi**2 * asymptotic_length * abs(beta2) *
+                              carrier.baud_rate * (delta_f - 0.5 * interfering_carrier.baud_rate))
+
+        return psi