Z = IMPEDANCE = TYPE_Z = 'Z' Y = ADMITTANCE = TYPE_Y = 'Y' S = SCATTERING = TYPE_S = 'S' T = SCATTERING_TRANSFER = TYPE_T = 'T' H = HYBRID = TYPE_H = 'H' G = INVERSE_HYBRID = TYPE_G = 'G' ABCD = TRANSMISSION = TYPE_ABCD = 'ABCD' Z = IMPEDANCE = TYPE_Z = 'Z'
Y = ADMITTANCE = TYPE_Y = 'Y'
S = SCATTERING = TYPE_S = 'S'
T = SCATTERING_TRANSFER = TYPE_T = 'T'
H = HYBRID = TYPE_H = 'H'
G = INVERSE_HYBRID = TYPE_G = 'G'
ABCD = TRANSMISSION = TYPE_ABCD = 'ABCD'
        P = type(cname, cbases, cdict)
class NPortMatrixBase(np.ndarray):
class NPortBase(NPortMatrixBase):
    def __new__(cls, matrix, type, z0=None):
    def _z0test(self, type, z0):
    def __repr__(self):
    def __array_finalize__(self, obj):
    def __reduce__(self):
    def __setstate__(self,state):
    def ports(self):
    def parameters(self):
    def get_parameter(self, port1, port2):
    def convert_z0test(self, type, z0):
    def __metaclass__(cname, cbases, cdict):
        def make_op(name, func):
            def op(this, other):
    def __new__(cls, freqs, matrices, type, z0=None):
    def __array_finalize__(self, obj):
    def __getitem__(self, index):
    def __getslice__(self, start, end):
    def __reduce__(self):
    def __setstate__(self,state):
    def parameters(self):
    def get_parameter(self, port1, port2):
    def get_element(self, port1, port2):
    def at(self, freqs):
    def average(self, n): class NPortMatrixBase(np.ndarray):
    """Base class representing an n-port matrix (Z, Y, S, T, G, H or ABCD)
    
    :type matrix: *n* by *n* array
    :param type: matrix type
    :type type: :data:`Z`, :data:`Y`, :data:`S`, :data:`T`, :data:`G`,
                :data:`H` or :data:`ABCD`
    :param z0: normalizing impedance (only :class:`S` and :class:`T`)
    :type z0: :class:`float`

    """
    # TODO: implement type checking for operations (as in NPortBase)
    def __new__(cls, matrix, type, z0=None):
        # http://docs.scipy.org/doc/numpy/user/basics.subclassing.html
        obj = matrix.view(cls)
        if type not in (Z, Y, S, T, H, G, ABCD):
            raise ValueError("illegal matrix type specified")
        obj.type = type
        obj.z0 = obj._z0test(type, z0)
        return obj

    def _z0test(self, type, z0):
        """Verify whether a normalizing impedance may be specified for the given
        n-port parameter type
        
        :param type: n-port parameter type
        :type type: :data:`Z`, :data:`Y`, :data:`S`, :data:`T`, :data:`G`,
                    :data:`H` or :data:`ABCD`
        :param z0: normalizing impedance
        :type z0: :class:`float`
        
        """
        if type in (S, T):
            if z0 is None:
                z0 = 50.0
        elif z0 is not None:
            raise ValueError("the specified n-port representation (%s) does "
                             "not require a reference impedance" % type)
        return z0

    def __repr__(self):
        array_repr = repr(self.view(np.ndarray))[6:-1].replace('\n     ', '\n')
        attrs = 'type=%s' % self.type
        if self.type in (S, T):
            attrs += ', z0=%g' % self.z0
        return '%s(%s)\n %s' % (self.__class__.__name__,
                                attrs,
                                array_repr)
    
    def __array_finalize__(self, obj):
        if obj is None: return
        self.type = getattr(obj, 'type', None)
        self.z0 = getattr(obj, 'z0', None)

    def __reduce__(self):
        object_state = list(np.ndarray.__reduce__(self))
        subclass_state = (self.type, self.z0, )
        object_state[2] = (object_state[2], subclass_state)
        return tuple(object_state)
    
    def __setstate__(self,state):
        nd_state, own_state = state
        np.ndarray.__setstate__(self, nd_state)
        self.type, self.z0 = own_state

    @property
    def ports(self):
        """The number of ports of this :class:`NPortMatrixBase`
        
        :rtype: :class:`int`
        
        """
        raise NotImplementedError

    @property
    def parameters(self):
        """Return an iterator over the parameters in row-major order
        
        :returns: iterator over parameter
        :rtype: iterator of :class:`complex`
        
        """
        return self.flat

    def get_parameter(self, port1, port2):
        """Return the parameter as specified by the indices `port1` and `port2`
        as an ndarray
        
        :param port1: first index
        :type port1: :class:`int`
        :param port2: second index
        :type port2: :class:`int`
        :returns: parameter at indices `port1` and `port2`
        :rtype: :class:`complex`
        
        """
        return self[port1 - 1, port2 - 1]

    def convert_z0test(self, type, z0):
        """Check supplied `z0` for conversion and set default value depending on
        the type
        
        :param type: target type for conversion
        :type type: :data:`Z`, :data:`Y`, :data:`S`, :data:`T`, :data:`G`,
                    :data:`H` or :data:`ABCD`
        :param z0: target normalization impedance (only :data:`S` and :data:`T`)
        :type z0: :class:`float`
        :returns: `z0` or default value when `z0 == None`
        :rtype: :class:`float`
        
        """
        if type in (S, T):
            if z0 is None:
                if self.type in (S, T):
                    z0 = self.z0
                else:
                    z0 = 50.0
        elif z0 is not None:
            raise ValueError("the specified n-port representation (%s) does"
                             " not require a reference impedance" % type)
        return z0 class NPortBase(NPortMatrixBase):
    """Base class representing an n-port across a list of frequencies

    :param freqs: list of frequency samples
    :type freqs: :class:`list`
    :param matrices: list of matrix elements
    :type matrix: list of 2 by 2 arrays in which each element is an 
                  *n* by *n* array
    :param type: matrix type
    :type type: :data:`Z`, :data:`Y`, :data:`S`, :data:`T`, :data:`G`, :data:`H`
                or :data:`ABCD`
    :param z0: normalizing impedance (only :data:`S` and :data:`T`)
    :type z0: :class:`float`

    """
    matrix_cls = NPortMatrixBase
    
    def __metaclass__(cname, cbases, cdict):
        # http://thread.gmane.org/gmane.comp.python.general/373788/focus=373795
        P = type(cname, cbases, cdict)
        def make_op(name, func):
            def op(this, other):
                if issubclass(type(other), NPortBase):
                    if other.type == this.type:
                        if other.z0 != this.z0:
                            raise ValueError("operands have different Z0")
                        result_type = this.type
                        result_z0 = this.z0
                    else:
                        raise ValueError("operands have different types")
                    min_freq = max(this.freqs[0], other.freqs[0])
                    max_freq = min(this.freqs[-1], other.freqs[-1])
                    result_freqs = list(set(this.freqs).union(set(other.freqs)))
                    result_freqs.sort()
                    result_freqs = \
                        np.array(result_freqs[result_freqs.index(min_freq):
                                              result_freqs.index(max_freq)+1])
                    this_matrices = this.at(result_freqs)
                    other_matrices = other.at(result_freqs)
                    result_matrices = func(this_matrices, other_matrices)
                else:
                    result_freqs = this.freqs
                    result_matrices = func(this, other)
                    result_type = this.type
                    result_z0 = this.z0
                subclass = type(this)
                return subclass(result_freqs, result_matrices, result_type,
                                result_z0)
            op.__name__ = name
            return op
        for name, func in np.ndarray.__dict__.items():
            if callable(func) and name in ['__add__', '__sub__', '__mul__',
                '__div__', '__radd__', '__rsub__', '__rmul__', '__rdiv__']:
                setattr(P, name, make_op(name, func))
        return P

    def __new__(cls, freqs, matrices, type, z0=None):
        if len(freqs) != len(matrices):
            raise ValueError("the list of frequencies and the list of "
                             "matrices should have equal lenghts")
        if matrices[0].shape[0] != matrices[0].shape[1]:
            raise ValueError("the matrices should be square")
        obj = super(NPortBase, cls).__new__(cls, matrices, type, z0)
        obj.freqs = np.asarray(freqs)
        return obj

    def __array_finalize__(self, obj):
        super(NPortBase, self).__array_finalize__(obj)
        self.freqs = getattr(obj, 'freqs', None)

    def __getitem__(self, index):
        if type(index) == int:
            return self.matrix_cls(np.asarray(self)[index], self.type, self.z0)
        else:
            return np.asarray(self).__getitem__(index)

    def __getslice__(self, start, end):
        return self.__class__(self.freqs[start:end],
                              np.asarray(self)[start:end], self.type, self.z0)

    def __reduce__(self):
        object_state = list(super(NPortBase, self).__reduce__())
        subclass_state = (self.freqs, )
        object_state[2] = (object_state[2], subclass_state)
        return tuple(object_state)
    
    def __setstate__(self,state):
        super_state, own_state = state
        super(NPortBase, self).__setstate__(super_state)
        self.freqs, = own_state

    @property
    def parameters(self):
        """Return an iterator over the parameters in row-major order
        
        :returns: iterator over parameters
        :rtype: iterator over :class:`ndarray`
        
        """
        for row in range(self.ports):
            for column in range(self.ports):
                yield np.asarray(self[:, row, column])

    def get_parameter(self, port1, port2):
        """Return the parameter as specified by the indices `port1` and `port2`
        as an ndarray
        
        :param port1: first index
        :type port1: :class:`int`
        :param port2: second index
        :type port2: :class:`int`
        :returns: parameter at indices `port1` and `port2`
        :rtype: :class:`ndarray`
        
        """
        return np.asarray(self[:, port1 - 1, port2 - 1])

    def get_element(self, port1, port2):
        """Return the submatrices made up of the element as specified by the
        indices `port1` and `port2`
        
        :param port1: first index
        :type port1: :class:`int`
        :param port2: second index
        :type port2: :class:`int`
        :returns: submatrices at indices `port1` and `port2`
        :rtype: :class:`ndarray`
        
        """
        subclass = type(self)
        return subclass(self.freqs,
            np.asarray(self[:, (port1 - 1, ), :][:, :, (port2 - 1, )]),
            self.type, self.z0)

    def at(self, freqs):
        """Return the linearly interpolated n-port data at the given
        * list of frequencies (`freqs` is iterable), or
        * at a single frequency (`freqs` is a value)
        
        :param freqs: frequency point or list of frequencies at which to return
                      the n-port matrices
        :type freqs: :class:`float` or iterable of :class:`float`s
        :rtype: `type(self)`
        
        """
        # check whether freqs is iterable
        try:
            it = iter(freqs)
            single = False
        except TypeError:
            single = True
            freqs = [freqs]
        subclass = type(self)
        func = interpolate.interp1d(self.freqs, self, axis=0)
        interpolated = func(freqs)
        interpolated_nport = subclass(freqs, interpolated, self.type, self.z0)
        if single:
            return interpolated_nport[0]
        else:
            return interpolated_nport
            
    def average(self, n):
        """Take a moving average over `n` frequency samples
        
        :param n: number of samples to average over
        :type n: :class:`int`
        :rtype: `type(self)`
        
        """
        averaged = np.zeros_like(np.asarray(self))
        for i in range(len(self)):
            for j in range(- int(np.floor(n/2)), int(np.ceil(n/2))):
                index = i + j
                if index < 0:
                    index = 0
                elif index >= len(self):
                    index = len(self) - 1
                averaged[i] += self[index] / n        
        # numpy convolve
        #~ ones = np.ones(n) / n
        #~ averaged = np.zeros_like(np.asarray(self))
        #~ for i in range(self.ports):
            #~ for j in range(self.ports):
                #~ averaged[:, i, j] = np.convolve(ones, self[:, i, j])
        # or scipy convolve
        #~ averaged = scipy.signal.convolve(np.ones((n, ) + self.shape[1:])/n, np.asarray(self), 'same')
        subclass = type(self)
        return subclass(self.freqs, averaged, self.type, self.z0) def rad(arg):
    """Return the phase in radians of the elements of the array"""
    return np.angle(arg) class TwoPortMatrix(NPortMatrix):
    """Class representing a 2-port matrix (Z, Y, S, T, G, H or ABCD)
    
    :param matrix: matrix elements
    :type matrix: 2 by 2 array
    :param type: matrix type
    :type type: :data:`Z`, :data:`Y`, :data:`S`, :data:`T`, :data:`G`, :data:`H`
                or :data:`ABCD`
    :param z0: normalizing impedance (only :data:`S` and :data:`T`)
    :type z0: :class:`float`

    """
    def convert(self, type, z0=None):
        """Convert to another 2-port matrix representation
        
        :param type: new matrix type
        :type type: :data:`Z`, :data:`Y`, :data:`S`, :data:`T`, :data:`G`,
                    :data:`H` or :data:`ABCD`
        :param z0: normalizing impedance (only :data:`S` and :data:`T`)
        :type z0: :class:`float`
        :rtype: :class:`NPortMatrix`

        """
        # references:
        #  * http://qucs.sourceforge.net/tech/node98.html
        #           "Transformations of n-Port matrices"
        z0 = self.convert_z0test(type, z0)

        if type not in (Z, Y, S, T, H, G, ABCD):
            raise TypeError("Unknown 2-port parameter type")

        # TODO: check for singularities (condition number)
        result = None
        if self.type == SCATTERING:
            s11 = self[0, 0]
            s12 = self[0, 1]
            s21 = self[1, 0]
            s22 = self[1, 1]
            if type == HYBRID:
                h11 = ((1 + s11) * (1 + s22) - s12 * s21) * self.z0
                h12 = 2 * s12
                h21 = -2 * s21
                h22 = ((1 - s11) * (1 - s22) - s12 * s21) / self.z0
                d = (1 - s11) * (1 + s22) + s12 * s21
                result = np.asarray([[h11, h12], [h21, h22]]) / d
            elif type == INVERSE_HYBRID:
                g11 = ((1 - s11) * (1 - s22) - s12 * s21) / self.z0
                g12 = -2 * s12
                g21 = 2 * s21
                g22 = ((1 + s11) * (1 + s22) - s12 * s21) * self.z0
                d = (1 + s11) * (1 - s22) + s12 * s21
                result = np.asarray([[g11, g12], [g21, g22]]) / d
        elif self.type == SCATTERING_TRANSFER:
            t11 = self[0, 0]
            t12 = self[0, 1]
            t21 = self[1, 0]
            t22 = self[1, 1]
            if type == HYBRID:
                h11 = (- t11 + t12 - t21 + t22) * self.z0
                h12 = - 2 * (t11 * t22 - t12 * t21)
                h21 = - 2
                h22 = (t11 + t12 - t21 - t22) / self.z0
                d = - t11 + t12 + t21 - t22
                result = np.asarray([[h11, h12], [h21, h22]]) / d
            elif type == INVERSE_HYBRID:
                g11 = (t11 + t12 - t21 - t22) / self.z0
                g12 = 2 * (t11 * t22 - t12 * t21)
                g21 = 2
                g22 = (- t11 + t12 - t21 + t22) * self.z0
                d = t11 + t12 + t21 + t22
                result = np.asarray([[g11, g12], [g21, g22]]) / d
        elif self.type == IMPEDANCE:
            z11 = self[0, 0]
            z12 = self[0, 1]
            z21 = self[1, 0]
            z22 = self[1, 1]
            dz = z11 * z22 - z12 * z21
            if type == HYBRID:
                result = np.asarray([[dz, z12], [-z21, 1]]) / z22
            elif type == INVERSE_HYBRID:
                result = np.asarray([[1, -z12], [z21, dz]]) / z11
        elif self.type == ADMITTANCE:
            y11 = self[0, 0]
            y12 = self[0, 1]
            y21 = self[1, 0]
            y22 = self[1, 1]
            dy = y11 * y22 - y12 * y21
            if type == HYBRID:
                result = np.asarray([[1, -y12], [y21, dy]]) / y11
            elif type == INVERSE_HYBRID:
                result = np.asarray([[dy, y12], [-y21, 1]]) / y22
        if self.type == TRANSMISSION:
            a11 = self[0, 0]
            a12 = self[0, 1]
            a21 = self[1, 0]
            a22 = self[1, 1]
            da = a11 * a22 - a12 * a21
            if type == HYBRID:
                result = np.asarray([[a12, da], [-1, a21]]) / a22
            elif type == INVERSE_HYBRID:
                result = np.asarray([[a21, -da], [1, a12]]) / a11
        elif self.type == HYBRID:
            h11 = self[0, 0]
            h12 = self[0, 1]
            h21 = self[1, 0]
            h22 = self[1, 1]
            dh = h11 * h22 - h12 * h21
            if type == SCATTERING:
                h11_ = h11 / z0
                h22_ = h22 * z0
                s11 = (h11_ - 1) * (1 + h22_) - h12 * h21
                s12 = 2 * h12
                s21 = -2 * h21
                s22 = (h11_ + 1) * (1 - h22_) + h12 * h21
                d = (h11_ + 1) * (1 + h22_) - h12 * h21
                result = np.asarray([[s11, s12], [s21, s22]]) / d
            elif type == SCATTERING_TRANSFER:
                h11_ = h11 / z0
                h22_ = h22 * z0
                t11 = + (h11_ + 1) * (1 - h22_) + h12 * h21
                t12 = - (h11_ + 1) * (1 + h22_) + h12 * h21
                t21 = + (h11_ - 1) * (1 - h22_) + h12 * h21
                t22 = - (h11_ - 1) * (1 + h22_) + h12 * h21
                result = np.asarray([[t11, t12], [t21, t22]]) / (2 * h21)
            elif type == IMPEDANCE:
                result = np.asarray([[dh, h12], [-h21, 1]]) / h22
            elif type == ADMITTANCE:
                result = np.asarray([[1, -h12], [h21, dh]]) / h11
            elif type == TRANSMISSION:
                result = np.asarray([[-dh, -h11], [-h22, -1]]) / h21
            elif type == HYBRID:
                result = self
            elif type == INVERSE_HYBRID:
                result = np.asarray([[h22, -h12], [-h21, h11]]) / dh
        elif self.type == INVERSE_HYBRID:
            g11 = self[0, 0]
            g12 = self[0, 1]
            g21 = self[1, 0]
            g22 = self[1, 1]
            dg = g11 * g22 - g12 * g21
            if type == SCATTERING:
                g11_ = g11 * z0
                g22_ = g22 / z0
                s11 = (1 - g11_) * (g22_ + 1) + g12 * g21
                s12 = -2 * g12
                s21 = 2 * g21
                s22 = (1 + g11_) * (g22_ - 1) - g12 * g21
                d = (1 + g11_) * (g22_ + 1) - g12 * g21
                result = np.asarray([[s11, s12], [s21, s22]]) / d
            elif type == SCATTERING_TRANSFER:
                g11_ = g11 * z0
                g22_ = g22 / z0
                t11 = + (g11_ + 1) * (1 - g22_) + g12 * g21
                t12 = + (g11_ + 1) * (1 + g22_) - g12 * g21
                t21 = - (g11_ - 1) * (1 - g22_) - g12 * g21
                t22 = - (g11_ - 1) * (1 + g22_) + g12 * g21
                result = np.asarray([[t11, t12], [t21, t22]]) / (2 * g21)
            elif type == IMPEDANCE:
                result = np.asarray([[1, -g12], [g21, dg]]) / g11
            elif type == ADMITTANCE:
                result = np.asarray([[dg, g12], [-g21, 1]]) / g22
            elif type == TRANSMISSION:
                result = np.asarray([[1, g22], [g11, dg]]) / g21
            elif type == HYBRID:
                result = np.asarray([[g22, -g12], [-g21, g11]]) / dg
            elif type == INVERSE_HYBRID:
                result = self
                
        if result is not None:
            return NPortMatrix(result, type, z0)
        else:
            nportmatrix = super(TwoPortMatrix, self)
            return nportmatrix.twonportmatrix().convert(type, z0).nportmatrix()

    def stability_k(self):
        """Returns the Rollet stability factor K and Delta
        
        :rtype: tuple
        
        """
        return stability_k(self)

    def stability_mu(self):
        """Returns the mu stability factor
        
        :rtype: :class:`float`
        
        """
        return stability_mu(self)

    def conditional_stability_mu(self, r1, r2):
        """Returns the conditional stability factor mu as defined in
        
        "Discussion and new proofs of the conditional stability criteria for 
        multidevice microwave amplifiers" by M. Balsi, G. Scotti, P. Tommasino
        and A. Trifiletti in "IEE Proceedings - Microwaves, Antennas and 
        Propagation", vol. 153, no. 2, pp. 177-181, 2006

        :param r1: maximum allowed value for the absolute value of the input
                   reflection coefficient (0 < r1 <= 1)
        :type r1: :class:`float`
        :param r2: maximum allowed value for the absolute value of the input
                   reflection coefficient (0 < r2 <= 1)
        :type r2: :class:`float`
        :rtype: :class:`float`
        
        """
        return conditional_stability_mu(self, r1, r2)

    def is_stable_k(self):
        k, delta = self.stability_k()
        return k > 1 and np.abs(delta) < 1            
    
    def is_stable_mu(self):
        return self.stability_mu() > 1

    def is_conditionally_stable_mu(self, r1, r2):
        return self.conditional_stability_mu(r1, r2) > 1

    def stability_circle_source(self):
        return stability_circle_source(self)

    def stability_circle_load(self):
        return stability_circle_load(self) class TwoPort(NPort):
    matrix_cls = TwoPortMatrix

    def convert(self, type, z0=None):
        # for two-ports, we can operate on arrays, speeding things up a lot
        # TODO: add other conversions (now handled by NPort.convert())
        z0 = self.convert_z0test(type, z0)
        
        if self.type == ABCD and type == S:
            # via scattering transfer parameters
            t = self.convert(T, z0)
            result = t.convert(S)
        elif self.type == T and type == S:
            t11, t12, t21, t22 = self.parameters
            t11i = 1 / t11
            s11 = t21 * t11i
            s12 = t22 - t21 * t11i * t12
            s21 = t11i
            s22 = - t11i * t12
            matrices = np.column_stack(((s11, s12, s21, s22))).reshape((-1,2,2))
            result = self.__class__(self.freqs, matrices, S, self.z0)
            if z0 != self.z0:
                result = result.renormalize(z0)
        elif self.type == S and type == T:
            # first renormalize
            if z0 != self.z0:
                result = self.renormalize(z0).convert(T, z0)
            else:
                s11, s12, s21, s22 = self.parameters
                s21i = 1 / s21
                t11 = s21i
                t12 = - s21i * s22
                t21 = s11 * s21i
                t22 = s12 - s11 * s21i * s22
                matrices = np.column_stack(((t11, t12,
                                             t21, t22))).reshape((-1,2,2))
                t = self.__class__(self.freqs, matrices, T, z0)
                result = t.renormalize(z0)
        elif self.type == ABCD and type == T:
            a, b, c, d = self.parameters
            b0 = b / z0
            c0 = c * z0
            t11 = a + b0 + c0 + d
            t12 = a - b0 + c0 - d
            t21 = a + b0 - c0 - d
            t22 = a - b0 - c0 + d
            matrices = np.column_stack(((t11, t12, t21, t22))).reshape((-1,2,2))
            result = 0.5 * self.__class__(self.freqs, matrices, T, z0)
        else:
            result = super(TwoPort, self).convert(type, z0)
            
        return result

    def stability_k(self):
        """Returns the Rollet stability factor K and Delta
        
        :rtype: tuple of :class:`ndarray`s
        
        """
        return stability_k(self)

    def stability_mu(self):
        """Returns the mu stability factor
        
        :rtype: :class:`ndarray`
        
        """
        return stability_mu(self)

    def conditional_stability_mu(self, r1, r2):
        """Returns the conditional stability factor mu as defined in
        
        "Discussion and new proofs of the conditional stability criteria for 
        multidevice microwave amplifiers" by M. Balsi, G. Scotti, P. Tommasino
        and A. Trifiletti in "IEE Proceedings - Microwaves, Antennas and 
        Propagation", vol. 153, no. 2, pp. 177-181, 2006

        :param r1: maximum allowed value for the absolute value of the input
                   reflection coefficient (0 < r1 <= 1)
        :type r1: :class:`float`
        :param r2: maximum allowed value for the absolute value of the input
                   reflection coefficient (0 < r2 <= 1)
        :type r2: :class:`float`
        :rtype: :class:`ndarray`
        
        """
        return conditional_stability_mu(self, r1, r2)

    def is_stable_k(self):
        for twoportmatrix in self:
            if not twoportmatrix.is_stable_k():
                return False
        return True
    
    def is_stable_mu(self):
        for twoportmatrix in self:
            if not twoportmatrix.is_stable_mu():
                return False
        return True

    def is_conditionally_stable_mu(self, r1, r2):
        for twoportmatrix in self:
            if not twoportmatrix.is_conditionally_stable_mu(r1, r2):
                return False
        return True

    def stability_circle_source(self):
        return stability_circle_source(self)

    def stability_circle_load(self):
        return stability_circle_load(self)

    def transmission_line(self, length, reciprocal=True):
        """Returns a :class:`TransmissionLine` represented by this
        :class:`TwoPort`
        
        :param length: physical length of the transmission line in meters
        :type length: :class:`float`
        :param reciprocal: specifies whether this TwoPort represents a reciprocal
                           transmission line (speeds up calculations)
        :type twonport: :class:`bool`
        
        """
        a, b, c, d = self.convert(TRANSMISSION).parameters
        sum = a + d
        diff = a - d
        ad_bc = a * d - b * c
        sq = unwrap_sqrt(sum**2 - 4 * ad_bc)

        exp_mgl_forward = 2 / (sum + sq)
        gamma_forward = - unwrap_log(exp_mgl_forward) / length
        z0_forward = (diff + sq) / (2 * c)

        if reciprocal:
            return GammaZ0TransmissionLine(self.freqs, gamma_forward,
                                           z0_forward)
        else:
            exp_mgl_backward = (sum - sq) / 2
            gamma_backward = - unwrap_log(exp_mgl_backward) / length
            z0_backward = (sq - diff) / (2 * c)
            return  GammaZ0TransmissionLine(self.freqs, gamma_forward,
                                           z0_forward, gamma_backward,
                                           z0_backward) class TwoNPortMatrix(NPortMatrixBase):
    """Class representing a 2n-port matrix (Z, Y, S, T or ABCD)

    :param freqs: list of frequency samples
    :type freqs: :class:`list`
    :param matrix: matrix elements
    :type matrix: 2 by 2 array where each element is an *n* by *n* array
    :param type: matrix type
    :type type: :data:`Z`, :data:`Y`, :data:`S`, :data:`T` or :data:`ABCD`
    :param z0: normalizing impedance (only :data:`S` and :data:`T`)
    :type z0: :class:`float`

    """
    # See also "block matrices" in the Matrix Cookbook
    def __new__(cls, matrix, type, z0=None):
        matrix = np.asarray(matrix, dtype=complex)
        obj = NPortMatrixBase.__new__(cls, matrix, type, z0)
        if len(obj.shape) != 4:
            raise ValueError("the matrix should be four-dimensional")
        if obj.shape[0] != 2 or obj.shape[0] != obj.shape[1]:
            raise ValueError("the matrices should be 2x2 matrices")
        if obj.shape[2] != obj.shape[3]:
            raise ValueError("the submatrices should be square")
        return obj

    @property
    def ports(self):
        """The number of ports of this :class:`TwoNPortMatrix`
        
        :rtype: :class:`int`
        
        """
        return 2 * self.shape[2]

    def nportmatrix(self):
        """Convert this :class:`TwoNPortMatrix` to an :class:`NPortMatrix`
        
        :rtype: :class:`NPortMatrix`
        
        """
        matrix = np.vstack((np.hstack((self[0,0], self[0,1])),
            np.hstack((self[1,0], self[1,1]))))
        return NPortMatrix(matrix, self.type, self.z0)

    def renormalize(self, z0):
        """Renormalize the 2n-port parameters to `z0`
        
        :param z0: new normalizing impedance
        :type z0: :class:`float`
        
        """
        assert self.type in (S, T)
        if z0 == self.z0:
            result = self
        elif self.type == S:
            result = self.nportmatrix().renormalize(z0).twonportmatrix()
        elif self.type == T:
            result = self.convert(S).renormalize(z0).convert(T)
        return result

    def convert(self, type, z0=None):
        """Convert from one 2n-port matrix representation to another
        
        :param type: new matrix type
        :type type: :data:`Z`, :data:`Y`, :data:`S`, :data:`T` or :data:`ABCD`
        :param z0: normalizing impedance (only :data:`S` and :data:`T`)
        :type z0: :class:`float`

        """
        z0 = self.convert_z0test(type, z0)
        idty = np.identity(self.shape[2], dtype=complex)
        
        # TODO: check for singularities
        if type == self.type:
            if type in (S, T):
                return self.renormalize(z0)
            else:
                return self                
        elif self.type == S and type == ABCD:
            s11 = np.asmatrix(self[0,0])
            s12 = np.asmatrix(self[0,1])
            s21 = np.asmatrix(self[1,0])
            s22 = np.asmatrix(self[1,1])
            s21i = np.linalg.inv(s21)
            i_plus_s11__times_s21i = (idty + s11) * s21i
            i_minus_s11__times_s21i = (idty - s11) * s21i
            i_plus_s22 = (idty + s22)
            i_minus_s22 = (idty - s22)
            a = (i_plus_s11__times_s21i * i_minus_s22 + s12)
            b = (i_plus_s11__times_s21i * i_plus_s22 - s12) * self.z0
            c = (i_minus_s11__times_s21i * i_minus_s22 - s12) / self.z0
            d = (i_minus_s11__times_s21i * i_plus_s22 + s12)
            result = 0.5 * self.__class__([[a, b], [c, d]], ABCD, z0)
        elif self.type == ABCD and type == S:
            # via scattering transfer parameters
            # TODO
            t = self.convert(T, z0)
            result = t.convert(S)
        elif self.type == S and type == T:
            # first renormalize
            if z0 != self.z0:
                result = self.renormalize(z0).convert(T, z0)
            else:
                s11 = np.asmatrix(self[0,0])
                s12 = np.asmatrix(self[0,1])
                s21 = np.asmatrix(self[1,0])
                s22 = np.asmatrix(self[1,1])
                s21i = np.linalg.inv(s21)
                t11 = s21i
                t12 = - s21i * s22
                t21 = s11 * s21i
                t22 = s12 - s11 * s21i * s22
                t = self.__class__([[t11, t12], [t21, t22]], T, z0)
                result = t.renormalize(z0)
        elif self.type == T and type == S:
            t11 = np.asmatrix(self[0,0])
            t12 = np.asmatrix(self[0,1])
            t21 = np.asmatrix(self[1,0])
            t22 = np.asmatrix(self[1,1])
            t11i = np.linalg.inv(t11)
            s11 = t21 * t11i
            s12 = t22 - t21 * t11i * t12
            s21 = t11i
            s22 = - t11i * t12
            result = self.__class__([[s11, s12], [s21, s22]], S, self.z0)
            if z0 != self.z0:
                result = result.renormalize(z0)
        elif self.type == T and type == ABCD:
            t11 = np.asmatrix(self[0,0])
            t12 = np.asmatrix(self[0,1])
            t21 = np.asmatrix(self[1,0])
            t22 = np.asmatrix(self[1,1])
            a = t11 + t12 + t21 + t22
            b = (t11 - t12 + t21 - t22) * self.z0
            c = (t11 + t12 - t21 -t22) / self.z0
            d = t11 - t12 -t21 + t22
            result = 0.5 * self.__class__([[a, b], [c, d]], ABCD, z0)
        elif self.type == ABCD and type == T:
            a = np.asmatrix(self[0,0])
            b0 = np.asmatrix(self[0,1]) / z0
            c0 = np.asmatrix(self[1,0]) * z0
            d = np.asmatrix(self[1,1])
            t11 = a + b0 + c0 + d
            t12 = a - b0 + c0 - d
            t21 = a + b0 - c0 - d
            t22 = a - b0 - c0 + d
            result = 0.5 * self.__class__([[t11, t12], [t21, t22]], T, z0)
        elif self.type == ABCD and type == Z:
            a = np.asmatrix(self[0,0])
            b = np.asmatrix(self[0,1])
            c = np.asmatrix(self[1,0])
            d = np.asmatrix(self[1,1])
            ci = np.linalg.inv(c)
            z11 = a * ci
            z12 = a * ci * d - b
            z21 = ci
            z22 = ci * d
            result = self.__class__([[z11, z12], [z21, z22]], Z, z0)
        elif self.type == Z and type == ABCD:
            z11 = np.asmatrix(self[0,0])
            z12 = np.asmatrix(self[0,1])
            z21 = np.asmatrix(self[1,0])
            z22 = np.asmatrix(self[1,1])
            z21i = np.linalg.inv(z21)
            a = z11 * z21i
            b = z11 * z21i * z22 - z12
            c = z21i
            d = z21i * z22
            result = self.__class__([[a, b], [c, d]], ABCD, z0)
        elif self.type == ABCD and type == Y:
            a = np.asmatrix(self[0,0])
            b = np.asmatrix(self[0,1])
            c = np.asmatrix(self[1,0])
            d = np.asmatrix(self[1,1])
            bi = np.linalg.inv(b)
            y11 = d * bi
            y12 = c - d * bi * a
            y21 = - bi
            y22 = bi * a
            result = self.__class__([[y11, y12], [y21, y22]], Y, z0)
        elif self.type == Y and type == ABCD:
            y11 = np.asmatrix(self[0,0])
            y12 = np.asmatrix(self[0,1])
            y21 = np.asmatrix(self[1,0])
            y22 = np.asmatrix(self[1,1])
            y21i = np.linalg.inv(y21)
            a = - y21i * y22
            b = - y21i
            c = y12 - y11 * y21i * y22
            d = - y11 * y21i
            result = self.__class__([[a, b], [c, d]], ABCD, z0)
        elif self.type == Z and type == T:
            z11 = np.asmatrix(self[0,0]) / z0
            z12 = np.asmatrix(self[0,1]) / z0
            z21 = np.asmatrix(self[1,0]) / z0
            z22 = np.asmatrix(self[1,1]) / z0
            z21i = np.linalg.inv(z21)
            i_plus_z11__times_z21i = (idty + z11) * z21i
            z11_minus_i__times_z21i = (z11 - idty) * z21i
            i_plus_z22 = (idty + z22)
            i_minus_z22 = (idty - z22)
            t11 = i_plus_z11__times_z21i * i_plus_z22 - z12
            t12 = i_plus_z11__times_z21i * i_minus_z22 + z12
            t21 = z11_minus_i__times_z21i * i_plus_z22 - z12
            t22 = z11_minus_i__times_z21i * i_minus_z22 + z12
            result = 0.5 * self.__class__([[t11, t12], [t21, t22]], T, z0)
        elif self.type == T and type == Z:
            # via scattering parameters
            s = self.convert(S).nportmatrix()
            z = s.convert(Z)
            return z.twonportmatrix()            
        elif self.type == Y and type == T:
            y11 = np.asmatrix(self[0,0]) * z0
            y12 = np.asmatrix(self[0,1]) * z0
            y21 = np.asmatrix(self[1,0]) * z0
            y22 = np.asmatrix(self[1,1]) * z0
            y21i = np.linalg.inv(y21)
            i_plus_y11__times_y21i = (idty + y11) * y21i
            i_minus_y11__times_y21i = (idty - y11) * y21i
            i_plus_y22 = (idty + y22)
            i_minus_y22 = (idty - y22)
            t11 = - i_plus_y11__times_y21i * i_plus_y22 + y12
            t12 = i_plus_y11__times_y21i * i_minus_y22 + y12
            t21 = - i_minus_y11__times_y21i * i_plus_y22 - y12
            t22 = i_minus_y11__times_y21i * i_minus_y22 - y12
            result = 0.5 * self.__class__([[t11, t12], [t21, t22]], T, z0)
        elif self.type == T and type == Y:
            # via scattering parameters
            s = self.convert(S).nportmatrix()
            y = s.convert(Y)
            return y.twonportmatrix()
        else:
            nportmatrix = self.nportmatrix()
            # catch infinite recursion for TwoPortMatrix
            nportmatrix.__class__ = NPortMatrix
            result = nportmatrix.convert(type, z0)
            return result.twonportmatrix()

        return result class TwoNPort(NPortBase):
    """Class representing a 2n-port across a list of frequencies

    :param matrices: list of matrix elements
    :type matrix: list of 2 by 2 arrays in which each element is an
                  *n* by *n* array
    :param type: matrix type
    :type type: :data:`Z`, :data:`Y`, :data:`S`, :data:`T` or :data:`ABCD`
    :param z0: normalizing impedance (only :data:`S` and :data:`T`)
    :type z0: :class:`float`

    """
    matrix_cls = TwoNPortMatrix

    def __new__(cls, freqs, matrices, type, z0=None):
        matrices = np.asarray(matrices, dtype=complex)
        obj = NPortBase.__new__(cls, freqs, matrices, type, z0)
        if len(obj[0].shape) != 4:
            raise ValueError("the matrices should be four-dimensional")
        if obj[0].shape[0] != 2 or obj[0].shape[0] != obj[0].shape[1]:
            raise ValueError("the matrices should be a 2x2 matrix")
        if obj[0].shape[2] != obj[0].shape[3]:
            raise ValueError("the submatrices should be square")
        return obj

    def __getitem__(self, arg):
        if type(arg) == int:
            return TwoNPortMatrix(np.asarray(self).__getitem__(arg), self.type,
                                  self.z0)
        else:
            return np.asarray(self).__getitem__(arg)

    @property
    def ports(self):
        """The number of ports of this :class:`TwoNPort`
        
        :rtype: :class:`int`
        
        """
        return 2 * self[0].shape[2]

    def add(self, freq, matrix):
        """Return a :class:`TwoNPort` with the specified frequency sample added.

        If `matrix` is a :class:`TwoNPortMatrix`, its elements will be converted
        to this :class:`TwoNPort`'s type and characteristic impedance. If
        `matrix` is a complex array, it is assumed the elements are in the
        format of this :class:`TwoNPort`.
        
        :param freq: frequency at which to insert `matrix`
        :type freq: :class:`float`
        :param matrix: matrix to insert at `freq`
        :type matrix: :class:`TwoNPortMatrix` or a complex array
        :rtype: :class:`TwoNPort`

        """
        if type(matrix) == TwoNPortMatrix:
            if matrix.type != self.type or matrix.z0 != self.z0:
                matrix = matrix.convert(self.type, self.z0)
        index = self.freqs.searchsorted(freq)
        freqs = np.insert(self.freqs, index, freq)
        matrices = np.insert(self.matrices, index, matrix, 0)
        return TwoNPort(freqs, matrices, self.type, self.z0)

    def nport(self):
        """Convert this :class:`TwoNPort` to an :class:`NPort`
        
        :rtype: :class:`NPort`

        """
        nportmatrices = []
        for matrix in self:
            twonportmatrix = TwoNPortMatrix(matrix, self.type, self.z0)
            nportmatrices.append(twonportmatrix.nportmatrix())
        return NPort(self.freqs, nportmatrices, self.type, self.z0)

    def renormalize(self, z0):
        """Renormalize the 2n-port parameters to `z0`
        
        :param z0: new normalizing impedance
        :type z0: :class:`float`
        :rtype: :class:`TwoNPort`
        
        """
        assert self.type in (S, T)
        if z0 == self.z0:
            result = self
        else:
            renormalized = []
            for matrix in self:
                twonportmatrix = TwoNPortMatrix(matrix, self.type, self.z0)
                renormalized.append(twonportmatrix.renormalize(z0))
            result = self.__class__(self.freqs, renormalized, self.type, z0)
        return result

    def convert(self, type, z0=None):
        """Convert to another 2n-port matrix representation
        
        :param type: new matrix type
        :type type: :data:`Z`, :data:`Y`, :data:`S`, :data:`T` or :data:`ABCD`
        :param z0: normalizing impedance (only :data:`S` and :data:`T`)
        :type z0: :class:`float`

        """
        converted = []
        for matrix in self:
            twonportmatrix = TwoNPortMatrix(matrix, self.type, self.z0)
            converted.append(twonportmatrix.convert(type, z0))
        return TwoNPort(self.freqs, converted, type, z0)

    #~ def invert(self):
        #~ """Return a TwoNPort described by the inverse of this TwoNPort's
        #~ matrices
        
        #~ """
        #~ inverted = []
        #~ for matrix in self:
            #~ nportmatrix = NPortMatrix(matrix, self.type, self.z0)
            #~ inverted.append(np.linalg.inv(nportmatrix))
        #~ return NPort(self.freqs, inverted, self.type, self.z0)
import numpy as np
from .base import Z, Y, S, T, H, G, ABCD
from .base import IMPEDANCE, ADMITTANCE, SCATTERING
from .base import NPortMatrixBase, NPortBase
from .parameter import rad
from .twoport import TwoPortMatrix, TwoPort
from .twonport import TwoNPortMatrix, TwoNPort

from __future__ import division




class NPortMatrix(NPortMatrixBase):
    """Class representing an n-port matrix (Z, Y or S; for 2-ports also T, G,
    H or ABCD)

    :param matrix: matrix elements
    :type matrix: *n* by *n* array
    :param type: matrix type
    :type type: :data:`Z`, :data:`Y`, :data:`S`, :data:`T`, :data:`G`, :data:`H`
                or :data:`ABCD`
    :param z0: normalizing impedance (only :class:`S` and :class:`T`)
    :type z0: :class:`float`
    
    """
    def __new__(cls, matrix, type, z0=None):
        matrix = np.asarray(matrix, dtype=complex)
        obj = NPortMatrixBase.__new__(cls, matrix, type, z0)
        if len(obj.shape) != 2:
            raise ValueError("the matrix should be two-dimensional")
        if obj.shape[0] != obj.shape[1]:
            raise ValueError("the matrix should be square")
        if obj.shape[0] == 2:
            obj.__class__ = TwoPortMatrix
        return obj

    @property
    def ports(self):
        """The number of ports of this :class:`NPortMatrix`
        
        :rtype: :class:`int`
        
        """
        return self.shape[0]

    def power(self, n):
        """Return this :class:`NPortMatrix` raised to the `n`\ th power
        
        :rtype: :class:`NPortMatrix`
        
        """
        return np.linalg.matrix_power(self, n)

    def twonportmatrix(self, inports=None, outports=None):
        """Return the 2n-port matrix represented by this n-port matrix
        
        :param inports: the list of ports that make up the inputs of the 2n-port
        :type inports: class:`tuple` or :class:`list`
        :param outports: the list of ports that make up the outputs
        :type outports: :class:`tuple` or :class:`list`
        :rtype: :class:`TwoNPortMatrix`
        
        """
        if self.ports % 2 != 0:
            raise TypeError("this NPortMatrix' number of ports is not a "
                            "multiple of 2")
        n = int(self.ports / 2)
        if inports is None and outports is None:
            inports = range(n)
            outports = range(n, 2*n)
            matrix = self
        else:
            # check whether the given sets of ports are valid
            assert inports is not None and outports is not None
            assert len(inports) == n
            assert len(outports) == n
            allports = set(inports).union(set(outports))
            assert len(allports) == 2*n
            assert min(allports) == 1 and max(allports) == 2*n
            inports = [port - 1 for port in inports]
            outports = [port - 1 for port in outports]

            ports = inports + outports

            # shuffle rows and columns to obtain 2n-port format
            matrix = []
            for row in ports:
                matrix.append(self[row])
            matrix = np.asarray(matrix).T
            matrix2 = []
            for column in ports:
                matrix2.append(matrix[column])
            matrix = np.asarray(matrix2).T

        # TODO: rewrite using reshape()?
        matrix = np.asmatrix(matrix)
        x11 = matrix[0:n  , 0:n  ]
        x12 = matrix[0:n  , n:2*n]
        x21 = matrix[n:2*n, 0:n  ]
        x22 = matrix[n:2*n, n:2*n]
        matrix = np.array([[x11, x12], [x21, x22]])
        return TwoNPortMatrix(matrix, self.type, self.z0)

    def reverse(self):
        """Return an :class:`NPortMatrix` with the ports reversed. This simply
        flips the matrix horizontally and vertically.
        
        :rtype: :class:`NPortMatrix`
        
        """
        return self[::-1, ::-1]

    def renormalize(self, z0):
        """Renormalize the n-port parameters to `z0`
        
        :param z0: new normalizing impedance
        :type z0: :class:`float`
        :rtype: :class:`NPortMatrix`
        
        """
        # http://qucs.sourceforge.net/tech/node98.html
        # "Renormalization of S-parameters to different port impedances"

        if self.type not in (S, T):