alpha0 computation moved from NLI solver to Fiber Params

gnpy/core/science_utils.py

@@ -140,6 +140,21 @@ class FiberParams():
     def temperature(self):
         return self._temperature
 
+    def alpha0(self, f_ref=193.5e12):
+        """ It returns the zero element of the series expansion of attenuation coefficient alpha(f) in the
+        reference frequency f_ref
+
+        :param f_ref: reference frequency of series expansion [Hz]
+        :return: alpha0: power attenuation coefficient in f_ref [Neper/m]
+        """
+        if not hasattr(self.loss_coef, 'alpha_power'):
+            alpha0 = self.loss_coef
+        else:
+            alpha_interp = interp1d(self.loss_coef['frequency'],
+                                    self.loss_coef['alpha_power'])
+            alpha0 = alpha_interp(f_ref)
+        return alpha0
+
 def load_sim_params(path_sim_params):
     sim_params = load_json(path_sim_params)
     return SimParams(params=sim_params)
@@ -208,27 +223,27 @@ def propagate_raman_fiber(fiber, *carriers):
         yield carrier._replace(power=pwr)
 
 def frequency_resolution(carrier, carriers, sim_params, fiber_params):
-    def _get_freq_res_k_phi(delta_count, grid_size, loss_coef, delta_z, beta2, k_tol, phi_tol):
+    def _get_freq_res_k_phi(delta_count, grid_size, alpha0, delta_z, beta2, k_tol, phi_tol):
         res_phi = _get_freq_res_phase_rotation(delta_count, grid_size, delta_z, beta2, phi_tol)
-        res_k = _get_freq_res_dispersion_attenuation(delta_count, grid_size, loss_coef, beta2, k_tol)
+        res_k = _get_freq_res_dispersion_attenuation(delta_count, grid_size, alpha0, beta2, k_tol)
         res_dict = {'res_phi': res_phi, 'res_k': res_k}
         method = min(res_dict, key=res_dict.get)
         return res_dict[method], method, res_dict
 
-    def _get_freq_res_dispersion_attenuation(delta_count, grid_size, loss_coef, beta2, k_tol):
-        return k_tol * abs(loss_coef) / abs(beta2) / (1 + delta_count) / (4 * np.pi ** 2 * grid_size)
+    def _get_freq_res_dispersion_attenuation(delta_count, grid_size, alpha0, beta2, k_tol):
+        return k_tol * abs(alpha0) / abs(beta2) / (1 + delta_count) / (4 * np.pi ** 2 * grid_size)
 
     def _get_freq_res_phase_rotation(delta_count, grid_size, delta_z, beta2, phi_tol):
         return phi_tol / abs(beta2) / (1 + delta_count) / delta_z / (4 * np.pi ** 2 * grid_size)
 
     grid_size = sim_params.nli_params.wdm_grid_size
-    loss_coef = fiber_params.loss_coef
     delta_z = sim_params.raman_params.space_resolution
+    alpha0 = fiber_params.alpha0()
     beta2 = fiber_params.beta2
     k_tol = sim_params.nli_params.dispersion_tolerance
     phi_tol = sim_params.nli_params.phase_shift_tollerance
     f_pump_resolution, method_f_pump, res_dict_pump = \
-        _get_freq_res_k_phi(0, grid_size, loss_coef, delta_z, beta2, k_tol, phi_tol)
+        _get_freq_res_k_phi(0, grid_size, alpha0, delta_z, beta2, k_tol, phi_tol)
     f_cut_resolution = {}
     method_f_cut = {}
     res_dict_cut = {}
@@ -236,7 +251,7 @@ def frequency_resolution(carrier, carriers, sim_params, fiber_params):
         delta_number = cut_carrier.channel_number - carrier.channel_number
         delta_count = abs(delta_number)
         f_res, method, res_dict = \
-            _get_freq_res_k_phi(delta_count, grid_size, loss_coef, delta_z, beta2, k_tol, phi_tol)
+            _get_freq_res_k_phi(delta_count, grid_size, alpha0, delta_z, beta2, k_tol, phi_tol)
         f_cut_resolution[f'delta_{delta_number}'] = f_res
         method_f_cut[delta_number] = method
         res_dict_cut[delta_number] = res_dict
@@ -605,15 +620,6 @@ class NliSolver:
         """
         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.
@@ -681,7 +687,7 @@ class NliSolver:
         :param carriers: the full WDM comb
         :return: carrier_nli: the amount of nonlinear interference in W on the carrier under analysis
         """
-        alpha = self.alpha0() / 2
+        alpha = self.fiber_params.alpha0() / 2
         beta2 = self.fiber_params.beta2
         gamma = self.fiber_params.gamma
         length = self.fiber_params.length
@@ -701,7 +707,7 @@ class NliSolver:
     def _psi(self, carrier, interfering_carrier):
         """ Calculates eq. 123 from arXiv:1209.0394.
         """
-        alpha = self.alpha0() / 2
+        alpha = self.fiber_params.alpha0() / 2
         beta2 = self.fiber_params.beta2
         asymptotic_length = 1 / (2 * alpha)
 
@@ -746,7 +752,7 @@ class NliSolver:
         :return: generalized_psi
         """
         # Fiber parameters
-        alpha0 = self.alpha0(f_eval)
+        alpha0 = self.fiber_params.alpha0(f_eval)
         beta2 = self.fiber_params.beta2
         beta3 = self.fiber_params.beta3
         f_ref_beta = self.fiber_params.f_ref_beta
@@ -779,7 +785,7 @@ class NliSolver:
         :return: generalized_psi
         """
         # Fiber parameters
-        alpha0 = self.alpha0(f_eval)
+        alpha0 = self.fiber_params.alpha0(f_eval)
         beta2 = self.fiber_params.beta2
         beta3 = self.fiber_params.beta3
         f_ref_beta = self.fiber_params.f_ref_beta