Introduce computation of the chromatic dispersion

Change-Id: I3ee039154568d4255444fa8db5e89945851010f4

gnpy/core/elements.py

@@ -22,7 +22,7 @@ instance as a result.
 
 from numpy import abs, arange, array, divide, errstate, ones
 from numpy import interp, mean, pi, polyfit, polyval, sum
-from scipy.constants import h
+from scipy.constants import h, c
 from collections import namedtuple
 
 from gnpy.core.utils import lin2db, db2lin, arrange_frequencies, snr_sum
@@ -82,6 +82,12 @@ class Transceiver(_Node):
         self.snr = None
         self.passive = False
         self.baud_rate = None
+        self.chromatic_dispersion = None
+
+    def _calc_cd(self, spectral_info):
+        """ Updates the Transceiver property with the CD of the received channels. CD in ps/nm.
+        """
+        self.chromatic_dispersion = [carrier.chromatic_dispersion * 1e3 for carrier in spectral_info.carriers]
 
     def _calc_snr(self, spectral_info):
         with errstate(divide='ignore'):
@@ -141,7 +147,8 @@ class Transceiver(_Node):
                 f'osnr_ase_01nm={self.osnr_ase_01nm!r}, '
                 f'osnr_ase={self.osnr_ase!r}, '
                 f'osnr_nli={self.osnr_nli!r}, '
-                f'snr={self.snr!r})')
+                f'snr={self.snr!r}, '
+                f'chromatic_dispersion={self.chromatic_dispersion!r})')
 
     def __str__(self):
         if self.snr is None or self.osnr_ase is None:
@@ -151,16 +158,19 @@ class Transceiver(_Node):
         osnr_ase = round(mean(self.osnr_ase), 2)
         osnr_ase_01nm = round(mean(self.osnr_ase_01nm), 2)
         snr_01nm = round(mean(self.snr_01nm), 2)
+        cd = mean(self.chromatic_dispersion)
 
         return '\n'.join([f'{type(self).__name__} {self.uid}',
 
                           f'  OSNR ASE (0.1nm, dB):      {osnr_ase_01nm:.2f}',
                           f'  OSNR ASE (signal bw, dB):  {osnr_ase:.2f}',
                           f'  SNR total (signal bw, dB): {snr:.2f}',
-                          f'  SNR total (0.1nm, dB): {snr_01nm:.2f}'])
+                          f'  SNR total (0.1nm, dB):     {snr_01nm:.2f}',
+                          f'  CD (ps/nm):                {cd:.2f}'])
 
     def __call__(self, spectral_info):
         self._calc_snr(spectral_info)
+        self._calc_cd(spectral_info)
         return spectral_info
 
 
@@ -377,6 +387,20 @@ class Fiber(_Node):
         """
         return self.alpha(f_ref * ones(1))[0]
 
+    def chromatic_dispersion(self, freq=193.5e12):
+        """ Returns accumulated chromatic dispersion (CD).
+
+        :param freq: the frequency at which the chromatic dispersion is computed
+        :return: chromatic dispersion: the accumulated dispersion [s/m]
+        """
+        beta2 = self.params.beta2
+        beta3 = self.params.beta3
+        ref_f = self.params.ref_frequency
+        length = self.params.length
+        beta = beta2 + 2 * pi * beta3 * (freq - ref_f)
+        dispersion = -beta * 2 * pi * ref_f**2 / c
+        return dispersion * length
+
     def _gn_analytic(self, carrier, *carriers):
         """Computes the nonlinear interference power on a single carrier.
         The method uses eq. 120 from `arXiv:1209.0394 <https://arxiv.org/abs/1209.0394>`__.
@@ -423,7 +447,8 @@ class Fiber(_Node):
             pwr = pwr._replace(signal=pwr.signal / self.params.lin_attenuation / attenuation,
                                nli=(pwr.nli + carrier_nli) / self.params.lin_attenuation / attenuation,
                                ase=pwr.ase / self.params.lin_attenuation / attenuation)
-            yield carrier._replace(power=pwr)
+            chromatic_dispersion = carrier.chromatic_dispersion + self.chromatic_dispersion(carrier.frequency)
+            yield carrier._replace(power=pwr, chromatic_dispersion=chromatic_dispersion)
 
     def update_pref(self, pref):
         self.pch_out_db = round(pref.p_spani - self.loss, 2)
@@ -461,6 +486,9 @@ class RamanFiber(Fiber):
 
     def propagate(self, *carriers):
         for propagated_carrier in propagate_raman_fiber(self, *carriers):
+            chromatic_dispersion = propagated_carrier.chromatic_dispersion + \
+                                   self.chromatic_dispersion(propagated_carrier.frequency)
+            propagated_carrier = propagated_carrier._replace(chromatic_dispersion=chromatic_dispersion)
             yield propagated_carrier
 
 

gnpy/core/info.py

@@ -17,8 +17,16 @@ class Power(namedtuple('Power', 'signal nli ase')):
     """carriers power in W"""
 
 
-class Channel(namedtuple('Channel', 'channel_number frequency baud_rate roll_off power')):
+class Channel(namedtuple('Channel', 'channel_number frequency baud_rate roll_off power chromatic_dispersion')):
-    pass
+    """ Class containing the parameters of a WDM signal.
+
+        :param channel_number: channel number in the WDM grid
+        :param frequency: central frequency of the signal (Hz)
+        :param baud_rate: the symbol rate of the signal (Baud)
+        :param roll_off: the roll off of the signal. It is a pure number between 0 and 1
+        :param power (gnpy.core.info.Power): power of signal, ASE noise and NLI (W)
+        :param chromatic_dispersion: chromatic dispersion (s/m)
+    """
 
 
 class Pref(namedtuple('Pref', 'p_span0, p_spani, neq_ch ')):
@@ -42,6 +50,7 @@ def create_input_spectral_information(f_min, f_max, roll_off, baud_rate, power,
         pref=Pref(pref, pref, lin2db(nb_channel)),
         carriers=[
             Channel(f, (f_min + spacing * f),
-                    baud_rate, roll_off, Power(power, 0, 0)) for f in range(1, nb_channel + 1)
+                    baud_rate, roll_off, Power(power, 0, 0), 0) for f in range(1, nb_channel + 1)
-        ])
+        ]
+    )
     return si

tests/invocation/transmission_main_example

@@ -11,7 +11,8 @@ Transceiver Site_A
   OSNR ASE (0.1nm, dB):      inf
   OSNR ASE (signal bw, dB):  inf
   SNR total (signal bw, dB): inf
-  SNR total (0.1nm, dB): inf
+  SNR total (0.1nm, dB):     inf
+  CD (ps/nm):                0.00
 Fiber          Span1
   type_variety:                SSMF
   length (km):                 80.00
@@ -37,7 +38,8 @@ Transceiver Site_B
   OSNR ASE (0.1nm, dB):      32.03
   OSNR ASE (signal bw, dB):  27.95
   SNR total (signal bw, dB): 26.27
-  SNR total (0.1nm, dB): 30.35
+  SNR total (0.1nm, dB):     30.35
+  CD (ps/nm):                1336.00
 
 Transmission result for input power = 0.00 dBm:
   Final SNR total (0.1 nm): 30.35 dB

tests/invocation/transmission_main_example__raman

@@ -11,7 +11,8 @@ Transceiver Site_A
   OSNR ASE (0.1nm, dB):      inf
   OSNR ASE (signal bw, dB):  inf
   SNR total (signal bw, dB): inf
-  SNR total (0.1nm, dB): inf
+  SNR total (0.1nm, dB):     inf
+  CD (ps/nm):                0.00
 RamanFiber          Span1
   type_variety:                SSMF
   length (km):                 80.00
@@ -37,7 +38,8 @@ Transceiver Site_B
   OSNR ASE (0.1nm, dB):      32.65
   OSNR ASE (signal bw, dB):  28.57
   SNR total (signal bw, dB): 26.48
-  SNR total (0.1nm, dB): 30.56
+  SNR total (0.1nm, dB):     30.56
+  CD (ps/nm):                1336.00
 
 Transmission result for input power = 0.00 dBm:
   Final SNR total (0.1 nm): 30.56 dB

tests/test_propagation.py

@@ -11,7 +11,7 @@ from gnpy.core.network import build_network
 from gnpy.tools.json_io import load_network, load_equipment
 from pathlib import Path
 from networkx import dijkstra_path
-from numpy import mean
+from numpy import mean, ones
 
 #network_file_name = 'tests/test_network.json'
 network_file_name = Path(__file__).parent.parent / 'tests/LinkforTest.json'
@@ -61,7 +61,8 @@ def propagation(input_power, con_in, con_out, dest):
     print(f'pw: {input_power} conn in: {con_in} con out: {con_out}',
           f'OSNR@0.1nm: {round(mean(sink.osnr_ase_01nm),2)}',
           f'SNR@bandwitdth: {round(mean(sink.snr),2)}')
-    return sink, nf
+    return sink, nf, path
+
 
 
 test = {'a': (-1, 1, 0), 'b': (-1, 1, 1), 'c': (0, 1, 0), 'd': (1, 1, 1)}
@@ -74,13 +75,13 @@ def test_snr(osnr_test, dest):
     pw = test[osnr_test][0]
     conn_in = test[osnr_test][1]
     conn_out = test[osnr_test][2]
-    sink, nf = propagation(pw, conn_in, conn_out, dest)
+    sink, nf, _ = propagation(pw, conn_in, conn_out, dest)
     osnr = round(mean(sink.osnr_ase), 3)
     nli = 1.0 / db2lin(round(mean(sink.snr), 3)) - 1.0 / db2lin(osnr)
     pw = expected[osnr_test][0]
     conn_in = expected[osnr_test][1]
     conn_out = expected[osnr_test][2]
-    sink, exp_nf = propagation(pw, conn_in, conn_out, dest)
+    sink, exp_nf, _ = propagation(pw, conn_in, conn_out, dest)
     expected_osnr = round(mean(sink.osnr_ase), 3)
     expected_nli = 1.0 / db2lin(round(mean(sink.snr), 3)) - 1.0 / db2lin(expected_osnr)
     # compare OSNR taking into account nf change of amps
@@ -89,6 +90,24 @@ def test_snr(osnr_test, dest):
     assert osnr_diff < 0.01 and nli_diff < 0.01
 
 
+@pytest.mark.parametrize("dest", ['trx B', 'trx F'])
+@pytest.mark.parametrize("cd_test", ['a', 'b', 'c', 'd'])
+def test_chromatic_dispersion(cd_test, dest):
+    pw = test[cd_test][0]
+    conn_in = test[cd_test][1]
+    conn_out = test[cd_test][2]
+    sink, _, path = propagation(pw, conn_in, conn_out, dest)
+
+    chromatic_dispersion = sink.chromatic_dispersion
+
+    num_ch = len(chromatic_dispersion)
+    expected_cd = 0
+    for el in path:
+        expected_cd += el.params.dispersion * el.params.length if isinstance(el, Fiber) else 0
+    expected_cd = expected_cd * ones(num_ch) * 1e3
+    assert chromatic_dispersion == pytest.approx(expected_cd)
+
+
 if __name__ == '__main__':
     from logging import getLogger, basicConfig, INFO
     logger = getLogger(__name__)