diff --git a/fgspectra/cross.py b/fgspectra/cross.py
index 58ae819..0dbb25f 100644
--- a/fgspectra/cross.py
+++ b/fgspectra/cross.py
@@ -135,6 +135,8 @@ def eval(self, sed_kwargs={}, cl_kwargs={}):
         """
         f_nu = self._sed(**sed_kwargs)
         cl = self._cl(**cl_kwargs)
+        if len(f_nu.shape) == len(cl.shape) == 1:
+            return np.outer(f_nu, f_nu)[:, :, np.newaxis] * cl
         if f_nu.shape[0] != cl.shape[-1] or (f_nu.shape[0] == 1 and cl.shape[-1] == 1):
             f_nu = f_nu[np.newaxis]
 
diff --git a/fgspectra/frequency.py b/fgspectra/frequency.py
index c0deab3..c4918fc 100644
--- a/fgspectra/frequency.py
+++ b/fgspectra/frequency.py
@@ -9,12 +9,15 @@
 and BeFoRe (David Alonso and Ben Thorne).
 """
 
-import inspect
-import types
 import numpy as np
 from scipy import constants
+from abc import abstractmethod
 from .model import Model
 
+try:
+    from numpy import trapezoid
+except ImportError:
+    from numpy import trapz as trapezoid
 
 T_CMB = 2.72548
 H_OVER_KT_CMB = constants.h * 1e9 / constants.k / T_CMB
@@ -39,71 +42,57 @@ def _rj2cmb(nu):
     return (np.expm1(x) / x) ** 2 / np.exp(x)
 
 
-def _bandpass_integration():
-    """Bandpass integrated version of the caller
+class FreqModel(Model):
+    def eval_bandpass(self, **kw):
+        """Bandpass integrated version of eval()
 
-    The caller should have
+        The inheriting class's eval() function should have
 
-        if isinstance(nu, list):
-            return _bandpass_integration()
+            if isinstance(nu, list):
+                return self.eval_bandpass(**kw)
 
-    at the very beginning.
-    This function
+        - to pass integrated list frequencies back to eval() one by one
 
-    * iterates over the ``nu`` argument of the caller
-      (while keeping all the other arguments fixed)
-    * splits each element of the iteration in ``nu_band, transmittance``
-    * integrates the caller function over the bandpass.
-      ``np.trapz(caller(nu_band) * transmittance, nu_band)``
-      Note that no normalization nor unit conversion is done to the
-      transmittance
-    * stacks the output of the iteration (the frequency dimension is the last)
-      and returns it
+        * iterates over the ``nu`` argument of the caller
+          (while keeping all the other arguments fixed)
+        * splits each element of the iteration in ``nu_band, transmittance``
+        * integrates the caller function over the bandpass.
+          ``np.trapezoid(self.eval(nu_band) * transmittance, nu_band)``
+          Note that no normalization nor unit conversion is done to the
+          transmittance
+        * stacks the output of the iteration (the frequency dimension is the last)
+          and returns it
 
-    """
-    # This piece of code is fairly complicated, we did because:
-    # 1) We want to call eval on each element of the nu list (i.e. we iterate
-    #    over the bandpasses) but we don't want to define a new eval_bandpass
-    #    function for every class
-    # 2) We don't want to use a decorator because it breaks the signature
-    #    handling of eval and the modification of its defaults.
-    #    _bandpass_integration does from the insides of eval the same thing that
-    #    a decorator would do from the outside. This is achieved through the
-    #    following pretty ugly kludge
-    # Simpler code that achieve the same result is welcome
-
-    # You are here because this function was called inside eval before any other
-    # variable was defined.
-    # We now retrieve the keyword arguments that were passed to eval because we
-    # have to use them for the evaluation of eval on each bandpass
-    # It assumes that _bandpass_integration was called inside
-    # f(self, **kw) -- f is typically the eval method.
-    frame = inspect.currentframe().f_back
-    kw = frame.f_locals
-    self = kw["self"]
-    del kw["self"]  # self was in the locals but is not a keyword argument
-
-    # We create a copy of eval itself, we'll call it for each bandpass
-    f = types.FunctionType(frame.f_code, frame.f_globals)
-    # Store the nu-transmittance list because the nu keyword argumnt has to be
-    # modified with the frequencies of each bandpass
-    nus_transmittances = kw["nu"]
-
-    # Get the shape of the output from the result of the first bandpass
-    kw["nu"] = nus_transmittances[0][0]
-    res = np.trapz(f(self, **kw) * nus_transmittances[0][1], kw["nu"])
-    # Append the frequency dimension and put res in its first entry
-    res = res[..., np.newaxis] * np.array([1.0] + [0.0] * (len(nus_transmittances) - 1))
-
-    # Fill the remaining entries by iterating over the rest of the bandpasses
-    for i_band, (nu, transmittance) in enumerate(nus_transmittances[1:], 1):
-        kw["nu"] = nu
-        res[..., i_band] = np.trapz(f(self, **kw) * transmittance, nu)
-
-    return res
-
-
-class PowerLaw(Model):
+        """
+        # This piece of code is fairly complicated, we did because:
+        # 1) We want to call eval on each element of the nu list (i.e. we iterate
+        #    over the bandpasses) but we don't want to define a new eval_bandpass
+        #    function for every class
+        # 2) We don't want to use a decorator because it breaks the signature
+        #    handling of eval and the modification of its defaults.
+
+        # You are here because this function was called inside eval before any other
+        # variable was defined.
+
+        # Store the nu-transmittance list because the nu keyword argument has to be
+        # modified with the frequencies of each bandpass
+        nus_transmittances = kw.pop("nu")
+        res = None
+        # Fill the entries by iterating over the bandpasses
+        for i_band, (nu, transmittance) in enumerate(nus_transmittances):
+            integral = trapezoid(self.eval(nu=nu, **kw) * transmittance, nu)
+            if res is None:
+                res = np.empty(integral.shape + (len(nus_transmittances),))
+            res[..., i_band] = integral
+
+        return res
+
+    @abstractmethod
+    def eval(self, **kwargs):
+        pass
+
+
+class PowerLaw(FreqModel):
     r"""Power Law
 
     .. math:: f(\nu) = (\nu / \nu_0)^{\beta}
@@ -147,7 +136,7 @@ def eval(self, nu=None, beta=None, nu_0=None):
 
         """
         if isinstance(nu, list):
-            return _bandpass_integration()
+            return self.eval_bandpass(nu=nu, beta=beta, nu_0=nu_0)
 
         beta = np.array(beta)[..., np.newaxis]
         nu_0 = np.array(nu_0)[..., np.newaxis]
@@ -160,7 +149,7 @@ class Synchrotron(PowerLaw):
     pass
 
 
-class ModifiedBlackBody(Model):
+class ModifiedBlackBody(FreqModel):
     r"""Modified black body in K_RJ
 
     .. math:: f(\nu) = (\nu / \nu_0)^{\beta + 1} / (e^x - 1)
@@ -190,7 +179,7 @@ def eval(self, nu=None, nu_0=None, temp=None, beta=None):
             dimensions of `beta` and `temp`.
         """
         if isinstance(nu, list):
-            return _bandpass_integration()
+            return self.eval_bandpass(nu=nu, nu_0=nu_0, temp=temp, beta=beta)
 
         beta = np.array(beta)[..., np.newaxis]
         temp = np.array(temp)[..., np.newaxis]
@@ -201,12 +190,12 @@ def eval(self, nu=None, nu_0=None, temp=None, beta=None):
 
 
 class CIB(ModifiedBlackBody):
-    """Alias of :class:`ModifiedBlackBOdy`"""
+    """Alias of :class:`ModifiedBlackBody`"""
 
     pass
 
 
-class ThermalSZ(Model):
+class ThermalSZ(FreqModel):
     r"""Thermal Sunyaev-Zel'dovich in K_CMB
 
     This class implements the
@@ -229,12 +218,12 @@ def eval(self, nu=None, nu_0=None):
         T_CMB (optional) : float
         """
         if isinstance(nu, list):
-            return _bandpass_integration()
+            return self.eval_bandpass(nu=nu, nu_0=nu_0)
 
         return ThermalSZ.f(nu) / ThermalSZ.f(nu_0)
 
 
-class FreeFree(Model):
+class FreeFree(FreqModel):
     r"""Free-free
 
     .. math:: f(\nu) = EM * ( 1 + log( 1 + (\nu_{ff} / \nu)^{3/\pi} ) )
@@ -273,18 +262,18 @@ def eval(self, nu=None, EM=None, Te=None):
 
         """
         if isinstance(nu, list):
-            return _bandpass_integration()
+            return self.eval_bandpass(nu=nu, EM=EM, Te=Te)
 
         EM = np.array(EM)[..., np.newaxis]
         Te = np.array(Te)[..., np.newaxis]
-        Teff = (Te / 1.0e3) ** (1.5)
+        Teff = (Te / 1.0e3) ** 1.5
         nuff = 255.33e9 * Teff
         gff = 1.0 + np.log(1.0 + (nuff / nu) ** (np.sqrt(3) / np.pi))
         print("warning: I need to check the units on this")
         return EM * gff
 
 
-class ConstantSED(Model):
+class ConstantSED(FreqModel):
     """Frequency-independent component."""
 
     def eval(self, nu=None, amp=1.0):
@@ -305,14 +294,14 @@ def eval(self, nu=None, amp=1.0):
             Note that the last dimension is guaranteed to be the frequency.
         """
         if isinstance(nu, list):
-            return _bandpass_integration()
+            return self.eval_bandpass(nu=nu, amp=amp)
 
         amp = np.array(amp)[..., np.newaxis]
         return amp * np.ones_like(np.array(nu))
 
 
-class FreeSED(Model):
-    """Frequency-dependent component for which every entries of the SED are specifified."""
+class FreeSED(FreqModel):
+    """Frequency-dependent component for which every entry of the SED is specified."""
 
     def eval(self, nu=None, sed=None):
         """Evaluation of the SED
@@ -338,11 +327,11 @@ def eval(self, nu=None, sed=None):
             except:
                 print("SED and nu must have the same shape.")
         if isinstance(nu, list):
-            return _bandpass_integration()
+            return self.eval_bandpass(nu=nu, sed=sed)
         return np.asarray(sed)
 
 
-class Join(Model):
+class Join(FreqModel):
     """Join several SED models together"""
 
     def __init__(self, *seds, **kwargs):
diff --git a/fgspectra/model.py b/fgspectra/model.py
index cf9fc28..f1e770f 100644
--- a/fgspectra/model.py
+++ b/fgspectra/model.py
@@ -13,7 +13,7 @@ class Model(ABC):
     evaluation means is defined by the `eval` method that any child class has to
     override.
 
-    If the `eval` method of a hypotetical ``Child`` class calls the `eval`
+    If the `eval` method of a hypothetical ``Child`` class calls the `eval`
     method of other `Model`s, it is likely that ``Child`` has also to override
     `set_defaults`, `defaults` and `_get_repr`
     """
diff --git a/fgspectra/power.py b/fgspectra/power.py
index ce0303e..691c6e3 100644
--- a/fgspectra/power.py
+++ b/fgspectra/power.py
@@ -251,9 +251,9 @@ def eval(self, kwseq=None):
 
         Parameters
         ----------
-        *argss
-            The length of `argss` has to be equal to the number of SEDs joined.
-            ``argss[i]`` is the argument list of the ``i``-th SED.
+        kwseq
+            The length of `kwseq` has to be equal to the number of ps joined.
+            ``kwseq[i]`` is the argument list of the ``i``-th ps.
 
         """
         spectra = np.array(
@@ -362,7 +362,7 @@ class SZxCIB_Reichardt2012(PowerSpectraAndCorrelation):
     """PowerSpectrum for SZxCIB (Dunkley et al. 2013)."""
 
     def __init__(self, **kwargs):
-        """Intialize object with parameters."""
+        """Initialize object with parameters."""
         power_spectra = [
             PowerSpectrumFromFile(_get_power_file("tsz_150_bat")),
             PowerLaw(),
diff --git a/tests/test_frequency.py b/tests/test_frequency.py
index d31171a..c25fe47 100644
--- a/tests/test_frequency.py
+++ b/tests/test_frequency.py
@@ -1,5 +1,4 @@
-import fgspectra.frequency as fgf
-from fgspectra.model import Model
+from fgspectra.frequency import FreqModel
 import numpy as np
 from numpy.testing import assert_allclose as aac
 
@@ -9,10 +8,10 @@ def test_bandpass_integration():
 
     transmittance = np.array([[1.0, 2.0, 4.0], [1.0, 3.0, 1.0]])
 
-    class Mock(Model):
+    class Mock(FreqModel):
         def eval(self, nu=None, par=None):
             if isinstance(nu, list):
-                return fgf._bandpass_integration()
+                return self.eval_bandpass(nu=nu, par=par)
             return np.ones_like(nu) * np.array(par)[..., np.newaxis]
 
     mock = Mock()