From f9bd6310f182b62f9d373f4cbbba11b0af312b9d Mon Sep 17 00:00:00 2001
From: Alessio Ferrari <alessio.ferrari@polito.it>
Date: Tue, 11 Jun 2019 13:38:39 +0200
Subject: [PATCH] introcude NLI solver

---
 gnpy/core/science_utils.py | 105 +++++++++++++++++++++++++++++++++++++
 1 file changed, 105 insertions(+)

diff --git a/gnpy/core/science_utils.py b/gnpy/core/science_utils.py
index e015444a..8fdcece8 100644
--- a/gnpy/core/science_utils.py
+++ b/gnpy/core/science_utils.py
@@ -340,3 +340,108 @@ class RamanSolver:
                 dpdz[f_ind][z_ind] = dpdz_element
 
         return np.vstack(dpdz)
+
+class NliSolver:
+    """ This class implements the NLI models.
+        Model and method can be specified in `self.nli_params.method`.
+        List of implemented methods:
+        'gn_model_analytic': brute force triple integral solution
+        'GGN_spectrally_separated_xpm_spm': XPM plus SPM
+    """
+
+    def __init__(self, nli_params=None, fiber_params=None):
+        """ Initialize the fiber object with its physical parameters
+        """
+        self.fiber_params = fiber_params
+        self.nli_params = nli_params
+        self.srs_profile = None
+
+    @property
+    def fiber_params(self):
+        return self.___fiber_params
+
+    @fiber_params.setter
+    def fiber_params(self, fiber_params):
+        self.___fiber_params = fiber_params
+
+    @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 nli_params(self):
+        return self.__nli_params
+
+    @nli_params.setter
+    def nli_params(self, nli_params):
+        """
+        :param model_params: namedtuple containing the parameters used to compute the NLI.
+        """
+        self.__nli_params = nli_params
+
+    def alpha0(self, f_eval=193.5e12):
+        if not hasattr(self.fiber_params.loss_coef, 'alpha_power'):
+            alpha0 = self.fiber_params.loss_coef
+        else:
+            alpha_interp = interp1d(self.fiber_params.loss_coef['frequency'],
+                                    self.fiber_params.loss_coef['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.nli_params.nli_method_name.lower():
+            carrier_nli = self._gn_analytic(carrier, *carriers)
+        else:
+            raise ValueError(f'Method {self.nli_params.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
+        beta2 = self.fiber_params.beta2
+        gamma = self.fiber_params.gamma
+        length = self.fiber_params.length
+        effective_length = (1 - np.exp(-2 * alpha * length)) / (2 * alpha)
+        asymptotic_length = 1 / (2 * alpha)
+
+        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_params.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