numpy.imag

Here are the examples of the python api numpy.imag taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.

98 Examples 7

3 Source : _ode.py
with GNU General Public License v3.0
from adityaprakash-bobby

    def _wrap(self, t, y, *f_args):
        f = self.cf(*((t, y[::2] + 1j * y[1::2]) + f_args))
        # self.tmp is a real-valued array containing the interleaved
        # real and imaginary parts of f.
        self.tmp[::2] = real(f)
        self.tmp[1::2] = imag(f)
        return self.tmp

    def _wrap_jac(self, t, y, *jac_args):

3 Source : _ode.py
with GNU General Public License v3.0
from adityaprakash-bobby

    def set_initial_value(self, y, t=0.0):
        """Set initial conditions y(t) = y."""
        y = asarray(y)
        self.tmp = zeros(y.size * 2, 'float')
        self.tmp[::2] = real(y)
        self.tmp[1::2] = imag(y)
        return ode.set_initial_value(self, self.tmp, t)

    def integrate(self, t, step=False, relax=False):

3 Source : pyhht.py
with GNU General Public License v3.0
from computer-animation-perception-group

def plothilbert(imfs):
    for i in range(imfs.shape[0]):
        h=hilbert(imfs[i,:])
        plot(real(h),imag(h))
    show()

def symmetrydemo():

3 Source : csdm.py
with BSD 3-Clause "New" or "Revised" License
from deepanshs

    def imag(self):
        """Return a csdm object with only the imaginary part of the dependent variable
        components."""
        return np.imag(self)

    @property

3 Source : gradients.py
with MIT License
from DTUWindEnergy

def cs(f, vector_interdependence=False, argnum=0, step=1e-20):
    def cs_gradient(res, _, step):
        return np.imag(res) / np.imag(step)
    return _step_grad(f, argnum, cs_gradient, step * 1j, vector_interdependence)


def autograd(f, vector_interdependence=False, argnum=0):

3 Source : test_preprocessing.py
with MIT License
from ECSHackWeek

def test_ignoreBelowX():
    filtered_freq, filtered_Z = ignoreBelowX(frequencies, Z_correct)

    assert (np.imag(filtered_Z) == Zi_np).all()


def test_cropFreq_maxonly():

3 Source : ltv.py
with BSD 3-Clause "New" or "Revised" License
from edwardcwang

    def _print_debug_msg(result):
        res_imag = np.imag(result).flatten()
        res_real = np.real(result).flatten()
        res_ratio = np.abs(res_imag / (res_real + 1e-18))
        idx = np.argmax(res_ratio)
        print('max imag/real ratio: %.4g, imag = %.4g, real = %.4g' %
              (res_ratio[idx], res_imag[idx], res_real[idx]))

    def __call__(self, t, tau, debug=False):

3 Source : wavefront.py
with MIT License
from ehpor

	def imag(self):
		'''The imaginary part of the wavefront as function of 2D position on the plane.
		'''
		return np.imag(self.electric_field)

	@property

3 Source : jax_backend_test.py
with Apache License 2.0
from google

def test_randn_non_zero_imag(dtype):
  backend = jax_backend.JaxBackend()
  a = backend.randn((4, 4), dtype=dtype)
  assert np.linalg.norm(np.imag(a)) != 0.0


@pytest.mark.parametrize("dtype", [np.complex64, np.complex128])

3 Source : jax_backend_test.py
with Apache License 2.0
from google

def test_random_uniform_non_zero_imag(dtype):
  backend = jax_backend.JaxBackend()
  a = backend.random_uniform((4, 4), dtype=dtype)
  assert np.linalg.norm(np.imag(a)) != 0.0


@pytest.mark.parametrize("dtype", np_dtypes)

3 Source : numpy_backend_test.py
with Apache License 2.0
from google

def test_randn_non_zero_imag(dtype):
  backend = numpy_backend.NumPyBackend()
  a = backend.randn((4, 4), dtype=dtype, seed=10)
  assert np.linalg.norm(np.imag(a)) != 0.0


@pytest.mark.parametrize("dtype", [np.complex64, np.complex128])

3 Source : numpy_backend_test.py
with Apache License 2.0
from google

def test_random_uniform_non_zero_imag(dtype):
  backend = numpy_backend.NumPyBackend()
  a = backend.random_uniform((4, 4), dtype=dtype, seed=10)
  assert np.linalg.norm(np.imag(a)) != 0.0


@pytest.mark.parametrize("dtype", np_dtypes)

3 Source : plotting.py
with The Unlicense
from johanna-rock

def phase_by_rd(rd):
    rd_imag = np.imag(rd)
    rd_imag[np.logical_or(np.isnan(rd_imag), np.isinf(rd_imag))] = 0
    rd_real = np.real(rd)
    rd_real[np.logical_or(np.isnan(rd_real), np.isinf(rd_real))] = 0
    rd_phase = np.arctan(rd_imag.astype('float') / rd_real.astype('float'))
    return rd_phase


def plot_row_column_cuts(rd_log_mag_signals, signal_labels, obj_col, obj_row, obj_id, packet_id, phase, is_rd, is_mag):

3 Source : integration.py
with GNU General Public License v3.0
from joselado

def recursive_asr_imag(f,a,b,eps,whole):
    "Recursive implementation of adaptive Simpson's rule."
    c = (a+b) / 2.0
    left = simpsons_rule(f,a,c)
    right = simpsons_rule(f,c,b)
    if np.max(np.abs(np.imag(left + right - whole)))   <  = 15*eps:
        return left + right + (left + right - whole)/15.0
    return recursive_asr(f,a,c,eps/2.0,left) + recursive_asr(f,c,b,eps/2.0,right)



def simpson(f,xlim=[0.,1.],eps=0.1):

3 Source : _ode.py
with MIT License
from ktraunmueller

    def _wrap(self, t, y, *f_args):
        f = self.cf(*((t, y[::2] + 1j * y[1::2]) + f_args))
        self.tmp[::2] = real(f)
        self.tmp[1::2] = imag(f)
        return self.tmp

    def _wrap_jac(self, t, y, *jac_args):

3 Source : _ode.py
with MIT License
from ktraunmueller

    def _wrap_jac(self, t, y, *jac_args):
        jac = self.cjac(*((t, y[::2] + 1j * y[1::2]) + jac_args))
        self.jac_tmp[1::2, 1::2] = self.jac_tmp[::2, ::2] = real(jac)
        self.jac_tmp[1::2, ::2] = imag(jac)
        self.jac_tmp[::2, 1::2] = -self.jac_tmp[1::2, ::2]
        return self.jac_tmp

    @property

3 Source : _ode.py
with MIT License
from ktraunmueller

    def set_initial_value(self, y, t=0.0):
        """Set initial conditions y(t) = y."""
        y = asarray(y)
        self.tmp = zeros(y.size * 2, 'float')
        self.tmp[::2] = real(y)
        self.tmp[1::2] = imag(y)
        if self.cjac is not None:
            self.jac_tmp = zeros((y.size * 2, y.size * 2), 'float')
        return ode.set_initial_value(self, self.tmp, t)

    def integrate(self, t, step=0, relax=0):

3 Source : exchangeCL2.py
with BSD 2-Clause "Simplified" License
from mailhexu

    def A_to_Jtensor(self):
        for key, val in self.JJ.items():
            # key:(R, iatom, jatom)
            R, iatom, jatom = key
            ispin = self.ispin(iatom)
            jspin = self.ispin(jatom)
            keyspin = (R, ispin, jspin)
            is_nonself = not (R == (0, 0, 0) and iatom == jatom)
            Jij = np.imag(val) / np.sign(
                np.dot(self.spinat[iatom], self.spinat[jatom]))
            Jorbij = np.imag(self.Jorb[key]) / np.sign(
                np.dot(self.spinat[iatom], self.spinat[jatom]))
            if is_nonself:
                self.exchange_Jdict[keyspin] = Jij
                self.Jiso_orb[
                    keyspin] = self.simplify_orbital_contributions(Jorbij, iatom, jatom)

    def get_rho_e(self, rho_up, rho_dn):

3 Source : exchangeCL2.py
with BSD 2-Clause "Simplified" License
from mailhexu

    def integrate(self, method="simpson"):
        if method == "trapezoidal":
            integrate = trapezoidal_nonuniform
        elif method == 'simpson':
            integrate = simpson_nonuniform
        self.rho_up = np.imag(integrate(self.contour.path, self.rho_up_list))
        self.rho_dn = np.imag(integrate(self.contour.path, self.rho_dn_list))
        for R, ijpairs in self.R_ijatom_dict.items():
            for iatom, jatom in ijpairs:
                self.Jorb[(R, iatom, jatom)] = integrate(
                    self.contour.path, self.Jorb_list[(R, iatom, jatom)])
                self.JJ[(R, iatom,
                         jatom)] = integrate(self.contour.path,
                                             self.JJ_list[(R, iatom, jatom)])

    def get_AijR_rhoR(self, e):

3 Source : test_adap_basis.py
with MIT License
from MesoscienceLab

def const_gw_sysbath(nsite, e_lambda, gamma, temp, gamma_mark):
    # define exponential parameters
    (g_0, w_0) = bcf_convert_sdl_to_exp(e_lambda, gamma, 0.0, temp)
    # define parameter lists
    gw_sysbath = []
    lop_list = []
    for i in range(nsite):
        loperator = sparse.coo_matrix(([1], ([i], [i])), shape=(nsite, nsite))
        gw_sysbath.append([g_0, w_0])
        lop_list.append(loperator)
        gw_sysbath.append([-1j * np.imag(g_0), gamma_mark])
        lop_list.append(loperator)
    return gw_sysbath, lop_list


def linear_chain(nsite, sb_params, V, maxhier=3, seed=None):

3 Source : _array.py
with MIT License
from michaelnowotny

def imag(a: ndarray):
    """
    Return the imaginary part of the complex argument.
    """

    return _unary_function(a, af_func=af.imag, np_func=np.imag)


def isempty(num: ndarray) -> bool:

3 Source : kinematics.py
with GNU General Public License v2.0
from mikecokina

def calculate_mode_derivatives(displacement, angular_frequency):
    """
    Calculates derivatives of angular displacement.

    :param displacement: numpy.array; complex
    :param angular_frequency: np.float;
    :return: numpy.array;
    """
    return angular_frequency * np.imag(displacement)


# _______________________acceleration coordinates_______________________
def calculate_mode_second_derivatives(displacement, angular_frequency):

3 Source : disp-solver.py
with MIT License
from natj

def fastest_growing_mode(k, params):
    omega = ES1d(k, params) 

    omega1 = np.imag(omega[0])
    if isnan(omega1):
        omega1 = 0.0 + 0.0j
    if isinf(omega1):
        omega1 = 0.0 + 0.0j

    return -omega1


def find_max_growth(params):

3 Source : _ode.py
with GNU Affero General Public License v3.0
from nccgroup

    def set_initial_value(self, y, t=0.0):
        """Set initial conditions y(t) = y."""
        y = asarray(y)
        self.tmp = zeros(y.size * 2, 'float')
        self.tmp[::2] = real(y)
        self.tmp[1::2] = imag(y)
        return ode.set_initial_value(self, self.tmp, t)

    def integrate(self, t, step=0, relax=0):

3 Source : functions.py
with GNU General Public License v3.0
from nrc-cnrc

def imag(x):
    """
    returns x.imag
    """
    return _bcallg(np.imag,x)

def conj(x):

3 Source : test_dsp.py
with MIT License
from pyfar

def test_zero_phase():
    """Test zero phase generation."""
    # generate test signal and zero phase version
    signal = pf.Signal([0, 0, 0, 2], 44100)
    signal_zero = dsp.zero_phase(signal)
    # assert type and id
    assert isinstance(signal_zero, pf.Signal)
    assert id(signal) != id(signal_zero)
    # assert freq data
    assert np.any(np.abs(np.imag(signal.freq)) > 1e-15)
    assert np.all(np.abs(np.imag(signal_zero.freq)) == 0)
    # assert time data
    npt.assert_allclose(signal_zero.time, np.atleast_2d([2, 0, 0, 0]))


def test_zero_phase_assertion():

3 Source : NTMpy.py
with MIT License
from udcm-su

    def ref2delta(self,refindex,lambdavac): 
        """
        Use the refractive index and compute the optical penetration depth
        This is used for Lambert Beer´s absorption law. 
        """
        lambdavac_m = lambdavac*1e-9 
        #corp away the two layers of air and only consider target film layers
        refindex = refindex[1:-1]
        #If there is no imaginary part we avoid dividing over 0
        for i in range(0,len(refindex)): 
            if np.imag(refindex[i]) == 0: 
                refindex[i] = refindex[i] + 1e-9j        
        deltap = (4*np.pi/lambdavac_m*np.imag(refindex))**(-1)
        return(deltap)
           
    def Gaussian(self,xmg,tmg,lam,A,sigma2,x0,customtime = None):

3 Source : strategy.py
with GNU General Public License v3.0
from ufvceiec

    def calculate_conn(self, data_intervals, i, j, sample_rate, channels, bands):
        _, Pxx = signal.welch(data_intervals[i], fs=sample_rate)
        _, Pyy = signal.welch(data_intervals[j], fs=sample_rate)
        f, Pxy = signal.csd(data_intervals[i],data_intervals[j],fs=sample_rate)
        icoh = np.imag(Pxy)/(np.sqrt(Pxx*Pyy))
        
        delta, theta, alpha, beta, gamma = frequency_bands(f, icoh)
        
        return delta.mean(), theta.mean(), alpha.mean(), beta.mean(), gamma.mean()
    
class Corr_cross_correlation_Estimator(Cross_correlation_rescaled):

0 Source : _ode.py
with GNU General Public License v3.0
from adityaprakash-bobby

    def _wrap_jac(self, t, y, *jac_args):
        # jac is the complex Jacobian computed by the user-defined function.
        jac = self.cjac(*((t, y[::2] + 1j * y[1::2]) + jac_args))

        # jac_tmp is the real version of the complex Jacobian.  Each complex
        # entry in jac, say 2+3j, becomes a 2x2 block of the form
        #     [2 -3]
        #     [3  2]
        jac_tmp = zeros((2 * jac.shape[0], 2 * jac.shape[1]))
        jac_tmp[1::2, 1::2] = jac_tmp[::2, ::2] = real(jac)
        jac_tmp[1::2, ::2] = imag(jac)
        jac_tmp[::2, 1::2] = -jac_tmp[1::2, ::2]

        ml = getattr(self._integrator, 'ml', None)
        mu = getattr(self._integrator, 'mu', None)
        if ml is not None or mu is not None:
            # Jacobian is banded.  The user's Jacobian function has computed
            # the complex Jacobian in packed format.  The corresponding
            # real-valued version has every other column shifted up.
            jac_tmp = _transform_banded_jac(jac_tmp)

        return jac_tmp

    @property

0 Source : mmio.py
with GNU General Public License v3.0
from adityaprakash-bobby

    def _write(self, stream, a, comment='', field=None, precision=None,
               symmetry=None):
        if isinstance(a, list) or isinstance(a, ndarray) or \
           isinstance(a, tuple) or hasattr(a, '__array__'):
            rep = self.FORMAT_ARRAY
            a = asarray(a)
            if len(a.shape) != 2:
                raise ValueError('Expected 2 dimensional array')
            rows, cols = a.shape

            if field is not None:

                if field == self.FIELD_INTEGER:
                    if not can_cast(a.dtype, 'intp'):
                        raise OverflowError("mmwrite does not support integer "
                                            "dtypes larger than native 'intp'.")
                    a = a.astype('intp')
                elif field == self.FIELD_REAL:
                    if a.dtype.char not in 'fd':
                        a = a.astype('d')
                elif field == self.FIELD_COMPLEX:
                    if a.dtype.char not in 'FD':
                        a = a.astype('D')

        else:
            if not isspmatrix(a):
                raise ValueError('unknown matrix type: %s' % type(a))
            
            rep = 'coordinate'
            rows, cols = a.shape

        typecode = a.dtype.char

        if precision is None:
            if typecode in 'fF':
                precision = 8
            else:
                precision = 16
        if field is None:
            kind = a.dtype.kind
            if kind == 'i':
                if not can_cast(a.dtype, 'intp'):
                    raise OverflowError("mmwrite does not support integer "
                                        "dtypes larger than native 'intp'.")
                field = 'integer'
            elif kind == 'f':
                field = 'real'
            elif kind == 'c':
                field = 'complex'
            elif kind == 'u':
                field = 'unsigned-integer'
            else:
                raise TypeError('unexpected dtype kind ' + kind)

        if symmetry is None:
            symmetry = self._get_symmetry(a)

        # validate rep, field, and symmetry
        self.__class__._validate_format(rep)
        self.__class__._validate_field(field)
        self.__class__._validate_symmetry(symmetry)

        # write initial header line
        stream.write(asbytes('%%MatrixMarket matrix {0} {1} {2}\n'.format(rep,
            field, symmetry)))

        # write comments
        for line in comment.split('\n'):
            stream.write(asbytes('%%%s\n' % (line)))

        template = self._field_template(field, precision)
        # write dense format
        if rep == self.FORMAT_ARRAY:
            # write shape spec
            stream.write(asbytes('%i %i\n' % (rows, cols)))

            if field in (self.FIELD_INTEGER, self.FIELD_REAL, self.FIELD_UNSIGNED):
                if symmetry == self.SYMMETRY_GENERAL:
                    for j in range(cols):
                        for i in range(rows):
                            stream.write(asbytes(template % a[i, j]))
                            
                elif symmetry == self.SYMMETRY_SKEW_SYMMETRIC:
                    for j in range(cols):
                        for i in range(j + 1, rows):
                            stream.write(asbytes(template % a[i, j]))
                            
                else:
                    for j in range(cols):
                        for i in range(j, rows):
                            stream.write(asbytes(template % a[i, j]))

            elif field == self.FIELD_COMPLEX:

                if symmetry == self.SYMMETRY_GENERAL:
                    for j in range(cols):
                        for i in range(rows):
                            aij = a[i, j]
                            stream.write(asbytes(template % (real(aij),
                                                             imag(aij))))
                else:
                    for j in range(cols):
                        for i in range(j, rows):
                            aij = a[i, j]
                            stream.write(asbytes(template % (real(aij),
                                                             imag(aij))))

            elif field == self.FIELD_PATTERN:
                raise ValueError('pattern type inconsisted with dense format')

            else:
                raise TypeError('Unknown field type %s' % field)

        # write sparse format
        else:
            coo = a.tocoo()  # convert to COOrdinate format

            # if symmetry format used, remove values above main diagonal
            if symmetry != self.SYMMETRY_GENERAL:
                lower_triangle_mask = coo.row >= coo.col
                coo = coo_matrix((coo.data[lower_triangle_mask],
                                 (coo.row[lower_triangle_mask],
                                  coo.col[lower_triangle_mask])),
                                 shape=coo.shape)

            # write shape spec
            stream.write(asbytes('%i %i %i\n' % (rows, cols, coo.nnz)))

            template = self._field_template(field, precision-1)

            if field == self.FIELD_PATTERN:
                for r, c in zip(coo.row+1, coo.col+1):
                    stream.write(asbytes("%i %i\n" % (r, c)))
            elif field in (self.FIELD_INTEGER, self.FIELD_REAL, self.FIELD_UNSIGNED):
                for r, c, d in zip(coo.row+1, coo.col+1, coo.data):
                    stream.write(asbytes(("%i %i " % (r, c)) +
                                         (template % d)))
            elif field == self.FIELD_COMPLEX:
                for r, c, d in zip(coo.row+1, coo.col+1, coo.data):
                    stream.write(asbytes(("%i %i " % (r, c)) +
                                         (template % (d.real, d.imag))))
            else:
                raise TypeError('Unknown field type %s' % field)


def _is_fromfile_compatible(stream):

0 Source : _decomp_ldl.py
with GNU General Public License v3.0
from adityaprakash-bobby

def ldl(A, lower=True, hermitian=True, overwrite_a=False, check_finite=True):
    """ Computes the LDLt or Bunch-Kaufman factorization of a symmetric/
    hermitian matrix.

    This function returns a block diagonal matrix D consisting blocks of size
    at most 2x2 and also a possibly permuted unit lower triangular matrix
    ``L`` such that the factorization ``A = L D L^H`` or ``A = L D L^T``
    holds. If ``lower`` is False then (again possibly permuted) upper
    triangular matrices are returned as outer factors.

    The permutation array can be used to triangularize the outer factors
    simply by a row shuffle, i.e., ``lu[perm, :]`` is an upper/lower
    triangular matrix. This is also equivalent to multiplication with a
    permutation matrix ``P.dot(lu)`` where ``P`` is a column-permuted
    identity matrix ``I[:, perm]``.

    Depending on the value of the boolean ``lower``, only upper or lower
    triangular part of the input array is referenced. Hence a triangular
    matrix on entry would give the same result as if the full matrix is
    supplied.

    Parameters
    ----------
    a : array_like
        Square input array
    lower : bool, optional
        This switches between the lower and upper triangular outer factors of
        the factorization. Lower triangular (``lower=True``) is the default.
    hermitian : bool, optional
        For complex-valued arrays, this defines whether ``a = a.conj().T`` or
        ``a = a.T`` is assumed. For real-valued arrays, this switch has no
        effect.
    overwrite_a : bool, optional
        Allow overwriting data in ``a`` (may enhance performance). The default
        is False.
    check_finite : bool, optional
        Whether to check that the input matrices contain only finite numbers.
        Disabling may give a performance gain, but may result in problems
        (crashes, non-termination) if the inputs do contain infinities or NaNs.

    Returns
    -------
    lu : ndarray
        The (possibly) permuted upper/lower triangular outer factor of the
        factorization.
    d : ndarray
        The block diagonal multiplier of the factorization.
    perm : ndarray
        The row-permutation index array that brings lu into triangular form.

    Raises
    ------
    ValueError
        If input array is not square.
    ComplexWarning
        If a complex-valued array with nonzero imaginary parts on the
        diagonal is given and hermitian is set to True.

    Examples
    --------
    Given an upper triangular array `a` that represents the full symmetric
    array with its entries, obtain `l`, 'd' and the permutation vector `perm`:

    >>> import numpy as np
    >>> from scipy.linalg import ldl
    >>> a = np.array([[2, -1, 3], [0, 2, 0], [0, 0, 1]])
    >>> lu, d, perm = ldl(a, lower=0) # Use the upper part
    >>> lu
    array([[ 0. ,  0. ,  1. ],
           [ 0. ,  1. , -0.5],
           [ 1. ,  1. ,  1.5]])
    >>> d
    array([[-5. ,  0. ,  0. ],
           [ 0. ,  1.5,  0. ],
           [ 0. ,  0. ,  2. ]])
    >>> perm
    array([2, 1, 0])
    >>> lu[perm, :]
    array([[ 1. ,  1. ,  1.5],
           [ 0. ,  1. , -0.5],
           [ 0. ,  0. ,  1. ]])
    >>> lu.dot(d).dot(lu.T)
    array([[ 2., -1.,  3.],
           [-1.,  2.,  0.],
           [ 3.,  0.,  1.]])

    Notes
    -----
    This function uses ``?SYTRF`` routines for symmetric matrices and
    ``?HETRF`` routines for Hermitian matrices from LAPACK. See [1]_ for
    the algorithm details.

    Depending on the ``lower`` keyword value, only lower or upper triangular
    part of the input array is referenced. Moreover, this keyword also defines
    the structure of the outer factors of the factorization.

    .. versionadded:: 1.1.0

    See also
    --------
    cholesky, lu

    References
    ----------
    .. [1] J.R. Bunch, L. Kaufman, Some stable methods for calculating
       inertia and solving symmetric linear systems, Math. Comput. Vol.31,
       1977. DOI: 10.2307/2005787

    """
    a = atleast_2d(_asarray_validated(A, check_finite=check_finite))
    if a.shape[0] != a.shape[1]:
        raise ValueError('The input array "a" should be square.')
    # Return empty arrays for empty square input
    if a.size == 0:
        return empty_like(a), empty_like(a), np.array([], dtype=int)

    n = a.shape[0]
    r_or_c = complex if iscomplexobj(a) else float

    # Get the LAPACK routine
    if r_or_c is complex and hermitian:
        s, sl = 'hetrf', 'hetrf_lwork'
        if np.any(imag(diag(a))):
            warn('scipy.linalg.ldl():\nThe imaginary parts of the diagonal'
                 'are ignored. Use "hermitian=False" for factorization of'
                 'complex symmetric arrays.', ComplexWarning, stacklevel=2)
    else:
        s, sl = 'sytrf', 'sytrf_lwork'

    solver, solver_lwork = get_lapack_funcs((s, sl), (a,))
    lwork = _compute_lwork(solver_lwork, n, lower=lower)
    ldu, piv, info = solver(a, lwork=lwork, lower=lower,
                            overwrite_a=overwrite_a)
    if info   <   0:
        raise ValueError('{} exited with the internal error "illegal value '
                         'in argument number {}". See LAPACK documentation '
                         'for the error codes.'.format(s.upper(), -info))

    swap_arr, pivot_arr = _ldl_sanitize_ipiv(piv, lower=lower)
    d, lu = _ldl_get_d_and_l(ldu, pivot_arr, lower=lower, hermitian=hermitian)
    lu, perm = _ldl_construct_tri_factor(lu, swap_arr, pivot_arr, lower=lower)

    return lu, d, perm


def _ldl_sanitize_ipiv(a, lower=True):

0 Source : test_filter_design.py
with GNU General Public License v3.0
from adityaprakash-bobby

    def test_high_order(self):
        # high even order, 'phase'
        z, p, k = bessel(24, 100, analog=True, output='zpk')
        z2 = []
        p2 = [
             -9.055312334014323e+01 + 4.844005815403969e+00j,
             -8.983105162681878e+01 + 1.454056170018573e+01j,
             -8.837357994162065e+01 + 2.426335240122282e+01j,
             -8.615278316179575e+01 + 3.403202098404543e+01j,
             -8.312326467067703e+01 + 4.386985940217900e+01j,
             -7.921695461084202e+01 + 5.380628489700191e+01j,
             -7.433392285433246e+01 + 6.388084216250878e+01j,
             -6.832565803501586e+01 + 7.415032695116071e+01j,
             -6.096221567378025e+01 + 8.470292433074425e+01j,
             -5.185914574820616e+01 + 9.569048385258847e+01j,
             -4.027853855197555e+01 + 1.074195196518679e+02j,
             -2.433481337524861e+01 + 1.207298683731973e+02j,
             ]
        k2 = 9.999999999999989e+47
        assert_array_equal(z, z2)
        assert_allclose(sorted(p, key=np.imag),
                        sorted(np.union1d(p2, np.conj(p2)), key=np.imag))
        assert_allclose(k, k2, rtol=1e-14)

        # high odd order, 'phase'
        z, p, k = bessel(23, 1000, analog=True, output='zpk')
        z2 = []
        p2 = [
             -2.497697202208956e+02 + 1.202813187870698e+03j,
             -4.126986617510172e+02 + 1.065328794475509e+03j,
             -5.304922463809596e+02 + 9.439760364018479e+02j,
             -9.027564978975828e+02 + 1.010534334242318e+02j,
             -8.909283244406079e+02 + 2.023024699647598e+02j,
             -8.709469394347836e+02 + 3.039581994804637e+02j,
             -8.423805948131370e+02 + 4.062657947488952e+02j,
             -8.045561642249877e+02 + 5.095305912401127e+02j,
             -7.564660146766259e+02 + 6.141594859516342e+02j,
             -6.965966033906477e+02 + 7.207341374730186e+02j,
             -6.225903228776276e+02 + 8.301558302815096e+02j,
             -9.066732476324988e+02]
        k2 = 9.999999999999983e+68
        assert_array_equal(z, z2)
        assert_allclose(sorted(p, key=np.imag),
                        sorted(np.union1d(p2, np.conj(p2)), key=np.imag))
        assert_allclose(k, k2, rtol=1e-14)

        # high even order, 'delay' (Orchard 1965 "The Roots of the
        # Maximally Flat-Delay Polynomials" Table 1)
        z, p, k = bessel(31, 1, analog=True, output='zpk', norm='delay')
        p2 = [-20.876706,
              -20.826543 + 1.735732j,
              -20.675502 + 3.473320j,
              -20.421895 + 5.214702j,
              -20.062802 + 6.961982j,
              -19.593895 + 8.717546j,
              -19.009148 + 10.484195j,
              -18.300400 + 12.265351j,
              -17.456663 + 14.065350j,
              -16.463032 + 15.889910j,
              -15.298849 + 17.746914j,
              -13.934466 + 19.647827j,
              -12.324914 + 21.610519j,
              -10.395893 + 23.665701j,
              - 8.005600 + 25.875019j,
              - 4.792045 + 28.406037j,
              ]
        assert_allclose(sorted(p, key=np.imag),
                        sorted(np.union1d(p2, np.conj(p2)), key=np.imag))

        # high odd order, 'delay'
        z, p, k = bessel(30, 1, analog=True, output='zpk', norm='delay')
        p2 = [-20.201029 + 0.867750j,
              -20.097257 + 2.604235j,
              -19.888485 + 4.343721j,
              -19.572188 + 6.088363j,
              -19.144380 + 7.840570j,
              -18.599342 + 9.603147j,
              -17.929195 + 11.379494j,
              -17.123228 + 13.173901j,
              -16.166808 + 14.992008j,
              -15.039580 + 16.841580j,
              -13.712245 + 18.733902j,
              -12.140295 + 20.686563j,
              -10.250119 + 22.729808j,
              - 7.901170 + 24.924391j,
              - 4.734679 + 27.435615j,
              ]
        assert_allclose(sorted(p, key=np.imag),
                        sorted(np.union1d(p2, np.conj(p2)), key=np.imag))

    def test_refs(self):

0 Source : test_filter_design.py
with GNU General Public License v3.0
from adityaprakash-bobby

    def test_bandpass(self):
        z, p, k = butter(8, [0.25, 0.33], 'band', output='zpk')
        z2 = [1, 1, 1, 1, 1, 1, 1, 1,
              -1, -1, -1, -1, -1, -1, -1, -1]
        p2 = [
            4.979909925436156e-01 + 8.367609424799387e-01j,
            4.979909925436156e-01 - 8.367609424799387e-01j,
            4.913338722555539e-01 + 7.866774509868817e-01j,
            4.913338722555539e-01 - 7.866774509868817e-01j,
            5.035229361778706e-01 + 7.401147376726750e-01j,
            5.035229361778706e-01 - 7.401147376726750e-01j,
            5.307617160406101e-01 + 7.029184459442954e-01j,
            5.307617160406101e-01 - 7.029184459442954e-01j,
            5.680556159453138e-01 + 6.788228792952775e-01j,
            5.680556159453138e-01 - 6.788228792952775e-01j,
            6.100962560818854e-01 + 6.693849403338664e-01j,
            6.100962560818854e-01 - 6.693849403338664e-01j,
            6.904694312740631e-01 + 6.930501690145245e-01j,
            6.904694312740631e-01 - 6.930501690145245e-01j,
            6.521767004237027e-01 + 6.744414640183752e-01j,
            6.521767004237027e-01 - 6.744414640183752e-01j,
            ]
        k2 = 3.398854055800844e-08
        assert_array_equal(z, z2)
        assert_allclose(sorted(p, key=np.imag),
                        sorted(p2, key=np.imag), rtol=1e-13)
        assert_allclose(k, k2, rtol=1e-13)

        # bandpass analog
        z, p, k = butter(4, [90.5, 110.5], 'bp', analog=True, output='zpk')
        z2 = np.zeros(4)
        p2 = [
            -4.179137760733086e+00 + 1.095935899082837e+02j,
            -4.179137760733086e+00 - 1.095935899082837e+02j,
            -9.593598668443835e+00 + 1.034745398029734e+02j,
            -9.593598668443835e+00 - 1.034745398029734e+02j,
            -8.883991981781929e+00 + 9.582087115567160e+01j,
            -8.883991981781929e+00 - 9.582087115567160e+01j,
            -3.474530886568715e+00 + 9.111599925805801e+01j,
            -3.474530886568715e+00 - 9.111599925805801e+01j,
            ]
        k2 = 1.600000000000001e+05
        assert_array_equal(z, z2)
        assert_allclose(sorted(p, key=np.imag), sorted(p2, key=np.imag))
        assert_allclose(k, k2, rtol=1e-15)

    def test_bandstop(self):

0 Source : test_filter_design.py
with GNU General Public License v3.0
from adityaprakash-bobby

    def test_bandstop(self):
        z, p, k = butter(7, [0.45, 0.56], 'stop', output='zpk')
        z2 = [-1.594474531383421e-02 + 9.998728744679880e-01j,
              -1.594474531383421e-02 - 9.998728744679880e-01j,
              -1.594474531383421e-02 + 9.998728744679880e-01j,
              -1.594474531383421e-02 - 9.998728744679880e-01j,
              -1.594474531383421e-02 + 9.998728744679880e-01j,
              -1.594474531383421e-02 - 9.998728744679880e-01j,
              -1.594474531383421e-02 + 9.998728744679880e-01j,
              -1.594474531383421e-02 - 9.998728744679880e-01j,
              -1.594474531383421e-02 + 9.998728744679880e-01j,
              -1.594474531383421e-02 - 9.998728744679880e-01j,
              -1.594474531383421e-02 + 9.998728744679880e-01j,
              -1.594474531383421e-02 - 9.998728744679880e-01j,
              -1.594474531383421e-02 + 9.998728744679880e-01j,
              -1.594474531383421e-02 - 9.998728744679880e-01j]
        p2 = [-1.766850742887729e-01 + 9.466951258673900e-01j,
              -1.766850742887729e-01 - 9.466951258673900e-01j,
               1.467897662432886e-01 + 9.515917126462422e-01j,
               1.467897662432886e-01 - 9.515917126462422e-01j,
              -1.370083529426906e-01 + 8.880376681273993e-01j,
              -1.370083529426906e-01 - 8.880376681273993e-01j,
               1.086774544701390e-01 + 8.915240810704319e-01j,
               1.086774544701390e-01 - 8.915240810704319e-01j,
              -7.982704457700891e-02 + 8.506056315273435e-01j,
              -7.982704457700891e-02 - 8.506056315273435e-01j,
               5.238812787110331e-02 + 8.524011102699969e-01j,
               5.238812787110331e-02 - 8.524011102699969e-01j,
              -1.357545000491310e-02 + 8.382287744986582e-01j,
              -1.357545000491310e-02 - 8.382287744986582e-01j]
        k2 = 4.577122512960063e-01
        assert_allclose(sorted(z, key=np.imag), sorted(z2, key=np.imag))
        assert_allclose(sorted(p, key=np.imag), sorted(p2, key=np.imag))
        assert_allclose(k, k2, rtol=1e-14)

    def test_ba_output(self):

0 Source : test_filter_design.py
with GNU General Public License v3.0
from adityaprakash-bobby

    def test_bandpass(self):
        z, p, k = cheby1(8, 1, [0.3, 0.4], 'bp', output='zpk')
        z2 = [1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1]
        p2 = [3.077784854851463e-01 + 9.453307017592942e-01j,
              3.077784854851463e-01 - 9.453307017592942e-01j,
              3.280567400654425e-01 + 9.272377218689016e-01j,
              3.280567400654425e-01 - 9.272377218689016e-01j,
              3.677912763284301e-01 + 9.038008865279966e-01j,
              3.677912763284301e-01 - 9.038008865279966e-01j,
              4.194425632520948e-01 + 8.769407159656157e-01j,
              4.194425632520948e-01 - 8.769407159656157e-01j,
              4.740921994669189e-01 + 8.496508528630974e-01j,
              4.740921994669189e-01 - 8.496508528630974e-01j,
              5.234866481897429e-01 + 8.259608422808477e-01j,
              5.234866481897429e-01 - 8.259608422808477e-01j,
              5.844717632289875e-01 + 8.052901363500210e-01j,
              5.844717632289875e-01 - 8.052901363500210e-01j,
              5.615189063336070e-01 + 8.100667803850766e-01j,
              5.615189063336070e-01 - 8.100667803850766e-01j]
        k2 = 5.007028718074307e-09
        assert_array_equal(z, z2)
        assert_allclose(sorted(p, key=np.imag),
                        sorted(p2, key=np.imag), rtol=1e-13)
        assert_allclose(k, k2, rtol=1e-13)

    def test_bandstop(self):

0 Source : test_filter_design.py
with GNU General Public License v3.0
from adityaprakash-bobby

    def test_bandstop(self):
        z, p, k = cheby1(7, 1, [0.5, 0.6], 'stop', output='zpk')
        z2 = [-1.583844403245361e-01 + 9.873775210440450e-01j,
              -1.583844403245361e-01 - 9.873775210440450e-01j,
              -1.583844403245361e-01 + 9.873775210440450e-01j,
              -1.583844403245361e-01 - 9.873775210440450e-01j,
              -1.583844403245361e-01 + 9.873775210440450e-01j,
              -1.583844403245361e-01 - 9.873775210440450e-01j,
              -1.583844403245361e-01 + 9.873775210440450e-01j,
              -1.583844403245361e-01 - 9.873775210440450e-01j,
              -1.583844403245361e-01 + 9.873775210440450e-01j,
              -1.583844403245361e-01 - 9.873775210440450e-01j,
              -1.583844403245361e-01 + 9.873775210440450e-01j,
              -1.583844403245361e-01 - 9.873775210440450e-01j,
              -1.583844403245361e-01 + 9.873775210440450e-01j,
              -1.583844403245361e-01 - 9.873775210440450e-01j]
        p2 = [-8.942974551472813e-02 + 3.482480481185926e-01j,
              -8.942974551472813e-02 - 3.482480481185926e-01j,
               1.293775154041798e-01 + 8.753499858081858e-01j,
               1.293775154041798e-01 - 8.753499858081858e-01j,
               3.399741945062013e-02 + 9.690316022705607e-01j,
               3.399741945062013e-02 - 9.690316022705607e-01j,
               4.167225522796539e-04 + 9.927338161087488e-01j,
               4.167225522796539e-04 - 9.927338161087488e-01j,
              -3.912966549550960e-01 + 8.046122859255742e-01j,
              -3.912966549550960e-01 - 8.046122859255742e-01j,
              -3.307805547127368e-01 + 9.133455018206508e-01j,
              -3.307805547127368e-01 - 9.133455018206508e-01j,
              -3.072658345097743e-01 + 9.443589759799366e-01j,
              -3.072658345097743e-01 - 9.443589759799366e-01j]
        k2 = 3.619438310405028e-01
        assert_allclose(sorted(z, key=np.imag),
                        sorted(z2, key=np.imag), rtol=1e-13)
        assert_allclose(sorted(p, key=np.imag),
                        sorted(p2, key=np.imag), rtol=1e-13)
        assert_allclose(k, k2, rtol=1e-15)

    def test_ba_output(self):

0 Source : basic.py
with GNU General Public License v3.0
from adityaprakash-bobby

def clpmn(m, n, z, type=3):
    """Associated Legendre function of the first kind for complex arguments.

    Computes the associated Legendre function of the first kind of order m and
    degree n, ``Pmn(z)`` = :math:`P_n^m(z)`, and its derivative, ``Pmn'(z)``.
    Returns two arrays of size ``(m+1, n+1)`` containing ``Pmn(z)`` and
    ``Pmn'(z)`` for all orders from ``0..m`` and degrees from ``0..n``.

    Parameters
    ----------
    m : int
       ``|m|   <  = n``; the order of the Legendre function.
    n : int
       where ``n >= 0``; the degree of the Legendre function.  Often
       called ``l`` (lower case L) in descriptions of the associated
       Legendre function
    z : float or complex
        Input value.
    type : int, optional
       takes values 2 or 3
       2: cut on the real axis ``|x| > 1``
       3: cut on the real axis ``-1  <  x  <  1`` (default)

    Returns
    -------
    Pmn_z : (m+1, n+1) array
       Values for all orders ``0..m`` and degrees ``0..n``
    Pmn_d_z : (m+1, n+1) array
       Derivatives for all orders ``0..m`` and degrees ``0..n``

    See Also
    --------
    lpmn: associated Legendre functions of the first kind for real z

    Notes
    -----
    By default, i.e. for ``type=3``, phase conventions are chosen according
    to [1]_ such that the function is analytic. The cut lies on the interval
    (-1, 1). Approaching the cut from above or below in general yields a phase
    factor with respect to Ferrer's function of the first kind
    (cf. `lpmn`).

    For ``type=2`` a cut at ``|x| > 1`` is chosen. Approaching the real values
    on the interval (-1, 1) in the complex plane yields Ferrer's function
    of the first kind.

    References
    ----------
    .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
           Functions", John Wiley and Sons, 1996.
           https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
    .. [2] NIST Digital Library of Mathematical Functions
           http://dlmf.nist.gov/14.21

    """
    if not isscalar(m) or (abs(m) > n):
        raise ValueError("m must be  < = n.")
    if not isscalar(n) or (n  <  0):
        raise ValueError("n must be a non-negative integer.")
    if not isscalar(z):
        raise ValueError("z must be scalar.")
    if not(type == 2 or type == 3):
        raise ValueError("type must be either 2 or 3.")
    if (m  <  0):
        mp = -m
        mf, nf = mgrid[0:mp+1, 0:n+1]
        with ufuncs.errstate(all='ignore'):
            if type == 2:
                fixarr = where(mf > nf, 0.0,
                               (-1)**mf * gamma(nf-mf+1) / gamma(nf+mf+1))
            else:
                fixarr = where(mf > nf, 0.0, gamma(nf-mf+1) / gamma(nf+mf+1))
    else:
        mp = m
    p, pd = specfun.clpmn(mp, n, real(z), imag(z), type)
    if (m  <  0):
        p = p * fixarr
        pd = pd * fixarr
    return p, pd


def lqmn(m, n, z):

0 Source : column.py
with MIT License
from AlexShkarin

    def imag(self):
        return ListDataColumn([np.imag(e) for e in self._column])
    def conjugate(self):

0 Source : column.py
with MIT License
from AlexShkarin

    def imag(self):
        return LinearDataColumn(self.length,np.imag(self.start),np.imag(self.step))
    def conjugate(self):

0 Source : coatingthermal.py
with MIT License
from ark0015

def getCoatDopt(ifo, T, dL, dCap=0.5):
    """Coating layer optical thicknesses to match desired transmission

    ifo = gwinc IFO model, or vector of 3 refractive indexes [nS, nL, nH]
    nS, nL, nH = refractive index of substrate, low and high-n materials
                 nL is used for the even layers, nH for odd layers
                 this algorithm may fail if nS > nL
    opticName = name of the Optic struct to use for T, dL and dCap
    T = power transmission of coating
    dL = optical thickness of low-n layers
         high-n layers have dH = 0.5 - dL
    dCap = first layer (low-n) thickness (default 0.5)
    dOpt = optical thickness vector Nlayer x 1

    """
    ##############################################
    def getTrans(ifo, Ndblt, dL, dH, dCap, dTweak):

        # the optical thickness vector
        dOpt = zeros(2 * Ndblt)
        dOpt[0] = dCap
        dOpt[1::2] = dH
        dOpt[2::2] = dL

        N = dTweak.size
        T = zeros(N)
        for n in range(N):
            dOpt[-1] = dTweak[n]
            r = getCoatRefl(ifo, dOpt)[0]
            T[n] = 1 - abs(r ** 2)

        return T

    ##############################################
    def getTweak(ifo, T, Ndblt, dL, dH, dCap, dScan, Nfit):

        # tweak bottom layer
        Tn = getTrans(ifo, Ndblt, dL, dH, dCap, dScan)
        pf = polyfit(dScan, Tn - T, Nfit)
        rts = roots(pf)
        if not any(imag(rts) == 0 & (rts > 0)):
            dTweak = None
            Td = 0
            return dTweak, Td
        dTweak = real(min(rts[(imag(rts) == 0) & (rts > 0)]))

        # compute T for this dTweak
        Td = getTrans(ifo, Ndblt, dL, dH, dCap, np.array([dTweak]))

        return dTweak, Td

        # plot for debugging
        #   plot(dScan, [Tn - T, polyval(pf, dScan)], dTweak, Td - T, 'ro')
        #   grid on
        #   legend('T exact', 'T fit', 'T at dTweak')
        #   title(sprintf('%d doublets', Ndblt))
        #   pause(1)

    # get IFO model stuff (or create it for other functions)
    pS = ifo.Materials.Substrate
    pC = ifo.Materials.Coating

    nS = pS.RefractiveIndex
    nL = pC.Indexlown
    nH = pC.Indexhighn

    ########################
    # find number of quarter-wave layers required, as first guess
    nR = nH / nL
    a1 = (2 - T + 2 * sqrt(1 - T)) / (nR * nH * T)
    Ndblt = int(ceil(log(a1) / (2 * log(nR))))

    # search through number of doublets to find Ndblt
    # which gives T lower than required
    dH = 0.5 - dL
    Tn = getTrans(ifo, Ndblt, dL, dH, dCap, np.array([dH]))
    while Tn   <   T and Ndblt > 1:
        # strange, but T is too low... remove doublets
        Ndblt = Ndblt - 1
        Tn = getTrans(ifo, Ndblt, dL, dH, dCap, np.array([dH]))
    while Tn > T and Ndblt  <  1e3:
        # add doublets until T > tN
        Ndblt = Ndblt + 1
        Tn = getTrans(ifo, Ndblt, dL, dH, dCap, np.array([dH]))

    ########################
    # tweak bottom layer
    delta = 0.01
    dScan = np.arange(0, 0.25 + delta, delta)
    dTweak = getTweak(ifo, T, Ndblt, dL, dH, dCap, dScan, 5)[0]

    if not dTweak:
        if nS > nL:
            raise Exception("Coating tweak layer not sufficient since nS > nL.")
        else:
            raise Exception("Coating tweak layer not found... very strange.")

    # now that we are close, get a better result with a linear fit
    delta = 0.001
    dScan = np.linspace(dTweak - 3 * delta, dTweak + 3 * delta, 7)
    dTweak, Td = getTweak(ifo, T, Ndblt, dL, dH, dCap, dScan, 3)

    # negative values are bad
    if dTweak  <  0.01:
        dTweak = 0.01

    # check the result
    if abs(log(Td / T)) > 1e-3:
        print("Exact coating tweak layer not found... %g%% error." % abs(log(Td / T)))

    ########################
    # return dOpt vector
    dOpt = zeros(2 * Ndblt)
    dOpt[0] = dCap
    dOpt[1::2] = dH
    dOpt[2::2] = dL
    dOpt[-1] = dTweak

    return dOpt


def getCoatRefl(ifo, dOpt):

0 Source : coatingthermal.py
with MIT License
from ark0015

def getCoatRefl2(nIn, nOut, nLayer, dOpt):
    """Coating reflection and phase derivatives

    nIn = refractive index of input medium (e.g., vacuum = 1)
    nOut = refractive index of output medium (e.g., SiO2 = 1.45231 @ 1064nm)
    nLayer = refractive index of each layer, ordered input to output (N x 1)
    dOpt   = optical thickness / lambda of each layer
           = geometrical thickness * refractive index / lambda

    rCoat = amplitude reflectivity of coating (complex) = rbar(0)
    dcdp = d reflection phase / d round-trip layer phase
    rbar = amplitude reflectivity of coating from this layer down
    r = amplitude reflectivity of this interface (r(1) is nIn to nLayer(1))

    see GWINC getCoatRefl and getCoatTOPhase for more information
    (see also T080101)

    """

    # Z-dir (1 => away from the substrate, -1 => into the substrate)
    zdir = 1

    # vector of all refractive indexs
    nAll = np.concatenate(([nIn], nLayer, [nOut]))

    # reflectivity of each interface
    r = (nAll[:-1] - nAll[1:]) / (nAll[:-1] + nAll[1:])

    # combine reflectivities
    rbar = zeros(r.size, dtype=complex)
    ephi = zeros(r.size, dtype=complex)

    # round-trip phase in each layer
    ephi[0] = 1
    ephi[1:] = exp(4j * zdir * pi * np.asarray(dOpt))

    rbar[-1] = ephi[-1] * r[-1]
    for n in range(len(dOpt), 0, -1):
        # accumulate reflectivity
        rbar[n - 1] = ephi[n - 1] * (r[n - 1] + rbar[n]) / (1 + r[n - 1] * rbar[n])

    # reflectivity derivatives
    dr_dphi = ephi[:-1] * (1 - r[:-1] ** 2) / (1 + r[:-1] * rbar[1:]) ** 2
    dr_dphi = (1j * zdir * rbar[1:]) * np.multiply.accumulate(dr_dphi)

    # shift rbar index
    rCoat = rbar[0]
    rbar = rbar[1:]

    # phase derivatives
    dcdp = -imag(dr_dphi / rCoat)  ### Where did this minus come from???

    return rCoat, dcdp, rbar, r

0 Source : suspensionthermal.py
with MIT License
from ark0015

def susptherm(f, ifo):
    """Suspention thermal noise.

    'Temp' must either be set for each stage individually, or globally
    in ifo.Suspension.Temp.  The latter will be preferred if
    specified, so if you wish to use per-stage tempurature you must
    remove Suspension.Temp.

    Assumes suspension transfer functions and V-H coupling have been
    pre-calculated and populated into the relevant `ifo` struct
    fields.

    """
    # Assign Physical Constants
    kB = const.kB

    # and vertical to beamline coupling angle
    theta = ifo.Suspension.VHCoupling.theta

    noise = zeros((1, f.size))

    # if the temperature is uniform along the suspension
    if "Temp" in ifo.Suspension:
        ##########################################################
        # Suspension TFs
        ##########################################################

        hForce = ifo.Suspension.hForce
        vForce = ifo.Suspension.vForce

        ##########################################################
        # Thermal Noise Calculation
        ##########################################################

        # convert to beam line motion
        #  theta is squared because we rotate by theta into the suspension
        #  basis, and by theta to rotate back to the beam line basis
        dxdF = hForce + theta ** 2 * vForce

        # thermal noise (m^2/Hz) for one suspension
        w = 2 * pi * f
        noise = 4 * kB * ifo.Suspension.Temp * abs(imag(dxdF)) / w

    # if the temperature is set for each suspension stage
    else:
        ##########################################################
        # Suspension TFs
        ##########################################################

        hForce = ifo.Suspension.hForce_singlylossy[:, :]
        vForce = ifo.Suspension.vForce_singlylossy[:, :]

        ##########################################################
        # Thermal Noise Calculation
        ##########################################################

        dxdF = zeros(hForce.shape, dtype=complex)
        for n, stage in enumerate(reversed(ifo.Suspension.Stage)):
            # add up the contribution from each stage

            # convert to beam line motion.  theta is squared because
            # we rotate by theta into the suspension basis, and by
            # theta to rotate back to the beam line basis
            dxdF[n, :] = hForce[n, :] + theta ** 2 * vForce[n, :]

            # thermal noise (m^2/Hz) for one suspension
            w = 2 * pi * f
            noise += 4 * kB * stage.Temp * abs(imag(dxdF[n, :])) / w

    # 4 masses, turn into gravitational wave strain
    noise *= 4 * ifo.gwinc.dhdl_sqr

    return np.squeeze(noise)

0 Source : suspension.py
with MIT License
from ark0015

def suspQuad(f, ifo, material="Silica"):
    """Suspension for quadruple pendulum

    `f` is frequency vector, `ifo` is IFO model.  `material` specifies
    material used for test mass suspension stage.  steel used for all
    other stages.  Violin modes are included.

    ifo.Suspension.FiberType should be: 0=round, 1=ribbons.

    hForce, vForce = transfer functions from the force on the TM to TM
    motion these should have the correct losses for the mechanical
    system such that the thermal noise is:

    dxdF = force on TM along beam line to position of TM along beam line
         = hForce + theta^2 * vForce
         = admittance / (i * w)

    where theta = ifo.Suspension.VHCoupling.theta.

    Since this is just suspension thermal noise, the TM internal modes
    and coating properties should not be included.

    hTable, vTable = TFs from support motion to TM motion

    Ah = horizontal equations of motion
    Av = vertical equations of motion

    Adapted from code by Morag Casey (Matlab) and Geppo Cagnoli
    (Maple).  Modification for the different temperatures between the
    stages by K.Arai.

    """
    g = const.g

    sus = ifo.Suspension

    # bottom stage fiber Type
    FiberType = sus.FiberType
    assert FiberType in FIBER_TYPES

    ##############################
    # Parameter Assignment

    def scarray():
        return np.zeros(len(sus.Stage))

    ds_w = scarray()
    dil0 = scarray()
    mass = scarray()
    length = scarray()
    kv0 = scarray()
    r_w = scarray()
    t_b = scarray()
    N_w = scarray()
    Temp = scarray()
    # wire material properties
    alpha_w = scarray()
    beta_w = scarray()
    rho_w = scarray()
    C_w = scarray()
    K_w = scarray()
    Y_w = scarray()
    phi_w = scarray()
    # blade material properties
    alpha_b = scarray()
    beta_b = scarray()
    rho_b = scarray()
    C_b = scarray()
    K_b = scarray()
    Y_b = scarray()
    phi_b = scarray()

    # NOTE: reverse suspension list so that the last stage in the list
    # is the test mass
    last_stage = len(sus.Stage) - 1
    for n, stage in enumerate(reversed(sus.Stage)):
        ##############################
        # main suspension parameters

        mass[n] = stage.Mass
        length[n] = stage.Length
        dil0[n] = stage.Dilution
        kv0[n] = stage.K

        if np.isnan(stage.WireRadius):
            if "FiberRadius" in stage:
                r_w[n] = stage.FiberRadius
            elif FiberType == "Ribbon":
                r_w[n] = sus.Ribbon.Width
            else:
                r_w[n] = sus.Fiber.Radius
        else:
            r_w[n] = stage.WireRadius

        # blade thickness
        t_b[n] = stage.Blade
        # number of support wires
        N_w[n] = stage.NWires

        if "Temp" in stage:
            Temp[n] = stage.Temp
        else:
            Temp[n] = sus.Temp

        ##############################
        # support wire material parameters

        if "WireMaterial" in stage:
            WireMaterial = stage.WireMaterial
        elif n == last_stage:
            WireMaterial = material
        else:
            WireMaterial = "C70Steel"

        logging.debug("stage {} wires: {}".format(n, WireMaterial))
        wireMat = sus[WireMaterial]

        alpha_w[n] = wireMat.Alpha  # coeff. thermal expansion
        beta_w[n] = wireMat.dlnEdT  # temp. dependence Youngs modulus
        rho_w[n] = wireMat.Rho  # mass density
        C_w[n] = wireMat.C  # heat capacity
        K_w[n] = wireMat.K  # W/(m kg)
        Y_w[n] = wireMat.Y  # Young's modulus
        phi_w[n] = wireMat.Phi  # loss angle

        # surface loss dissipation depth
        if "Dissdepth" in wireMat:
            ds_w[n] = wireMat.Dissdepth
        else:
            # otherwise ignore surface effects
            ds_w[n] = 0

        ##############################
        # support blade material parameters

        if "BladeMaterial" in stage:
            BladeMaterial = stage.BladeMaterial
        else:
            BladeMaterial = "MaragingSteel"

        logging.debug("stage {} blades: {}".format(n, BladeMaterial))
        bladeMat = sus[BladeMaterial]

        alpha_b[n] = bladeMat.Alpha  # coeff. thermal expansion
        beta_b[n] = bladeMat.dlnEdT  # temp. dependence Youngs modulus
        rho_b[n] = bladeMat.Rho  # mass density
        C_b[n] = bladeMat.C  # heat capacity
        K_b[n] = bladeMat.K  # W/(m kg)
        Y_b[n] = bladeMat.Y  # Young's modulus
        phi_b[n] = bladeMat.Phi  # loss angle

    # weight support by lower stages
    Mg = g * np.flipud(np.cumsum(np.flipud(mass)))

    # Correction for the pendulum restoring force
    kh0 = Mg / length  # N/m, horiz. spring constant, stage n

    ##############################
    # Thermoelastic Calculations for wires and blades

    # wire geometry
    tension = Mg / N_w  # Tension
    xsect = pi * r_w ** 2  # cross-sectional area
    xII = r_w ** 4 * pi / 4  # x-sectional moment of inertia
    mu_h = 4 / r_w  # surface to volume ratio, horizontal
    mu_v = 2 / r_w  # surface to volume ratio, vertical (wire)

    # horizontal TE time constant, wires
    # The constant 7.37e-2 is 1/(4*q0^2) from eq 12, C. Zener 10.1103/PhysRev.53.90 (1938)
    tau_h = 7.37e-2 * 4 * (rho_w * C_w * xsect) / (pi * K_w)

    # vertical TE time constant, blades
    tau_v = (rho_b * C_b * t_b ** 2) / (K_b * pi ** 2)

    # vertical delta, blades
    delta_v = Y_b * alpha_b ** 2 * Temp / (rho_b * C_b)

    # deal with ribbon geometry for last stage
    if FiberType == "Ribbon":
        W = sus.Ribbon.Width
        t = sus.Ribbon.Thickness
        xsect[-1] = W * t  # cross-sectional area
        xII[-1] = (W * t ** 3) / 12  # x-sectional moment of inertia
        mu_v[-1] = 2 * (W + t) / (W * t)
        mu_h[-1] = mu_v[-1] * (3 * N_w[-1] * W + t) / (N_w[-1] * W + t)
        tau_h[-1] = (rho_w[-1] * C_w[-1] * t ** 2) / (K_w[-1] * pi ** 2)

    # horizontal delta, wires
    delta_h = (
        (alpha_w - tension * beta_w / (xsect * Y_w)) ** 2 * Y_w * Temp / (rho_w * C_w)
    )

    # deal with tapered geometry for last stage
    if FiberType == "Tapered":
        r_end = sus.Fiber.EndRadius

        # recompute these for
        xsectEnd = pi * r_end ** 2  # cross-sectional area (for delta_h)
        xII[-1] = pi * r_end ** 4 / 4  # x-sectional moment of inertia
        mu_h[-1] = 4 / r_end  # surface to volume ratio, horizontal

        # use this xsect for thermo-elastic noise
        delta_h[-1] = (
            (alpha_w[-1] - tension[-1] * beta_w[-1] / (xsectEnd * Y_w[-1])) ** 2
            * Y_w[-1]
            * Temp[-1]
            / (rho_w[-1] * C_w[-1])
        )

    # bending length, and dilution factors
    d_bend = sqrt(Y_w * xII / tension)
    dil = length / d_bend
    # dil(~isnan(dil0)) = dil0(~isnan(dil0))
    dil = np.where(~np.isnan(dil0), dil0, dil)

    ##############################
    # Loss Calculations for wires and blades

    # these calculations use the frequency vector
    w = 2 * pi * f

    phih = np.zeros([len(sus.Stage), len(w)], dtype=complex)
    kh = np.zeros([len(sus.Stage), len(w)], dtype=complex)
    phiv = np.zeros([len(sus.Stage), len(w)], dtype=complex)
    kv = np.zeros([len(sus.Stage), len(w)], dtype=complex)

    for n, stage in enumerate(sus.Stage):
        # horizontal loss factor, wires
        phih[n, :] = phi_w[n] * (1 + mu_h[n] * ds_w[n]) + delta_h[n] * tau_h[n] * w / (
            1 + w ** 2 * tau_h[n] ** 2
        )

        # complex spring constant, horizontal
        kh[n, :] = kh0[n] * (1 + 1j * phih[n, :] / dil[n])

        # vertical loss factor, blades
        phiv[n, :] = phi_b[n] + delta_v[n] * tau_v[n] * w / (1 + w ** 2 * tau_v[n] ** 2)

        # complex spring constant, vertical
        kv[n, :] = kv0[n] * (1 + 1j * phiv[n, :])

    ##############################
    # last suspension stage
    # Equations from "GG" (maybe?)
    #   Suspensions thermal noise in the LIGO gravitational wave detector
    #   Gabriela Gonzalez, Class. Quantum Grav. 17 (2000) 4409?4435
    #
    # Note:
    #  I in GG = xII
    #  rho in GG = rho_w * xsect
    #  delta in GG = d_bend
    ##############################

    ### Vertical (bounce) ###
    # loss factor, last stage suspension, vertical (no blades)
    phiv[-1, :] = phi_w[-1] * (1 + mu_v[-1] * ds_w[-1])

    # vertical Young's modulus
    Y_v = Y_w[-1] * (1 + 1j * phiv[-1, :])

    # vertical spring constant, last stage
    k_z = sqrt(rho_w[-1] / Y_v) * w
    kv4 = N_w[-1] * xsect[-1] * Y_v * k_z / (tan(k_z * length[-1]))

    # deal with tapered geometry for last stage
    if FiberType == "Tapered" and "EndLength" in sus.Fiber:
        l_end = 2 * sus.Fiber.EndLength
        l_mid = length[-1] - l_end

        kv_mid = N_w[-1] * xsect[-1] * Y_v * k_z / (tan(k_z * l_mid))
        kv_end = N_w[-1] * xsectEnd * Y_v * k_z / (tan(k_z * l_end))
        kv4 = kv_mid * kv_end / (kv_mid + kv_end)

    if np.isnan(kv0[-1]):
        kv[-1, :] = kv4  # no blades
    else:
        kv[-1, :] = kv[-1, :] * kv4 / (kv[-1, :] + kv4)  # with blades

    ### Horizontal (pendulum and violins) ###
    # horizontal Young's modulus
    Y_h = Y_w[-1] * (1 + 1j * phih[-1, :])

    # simplification factors for later calculations
    ten4 = tension[-1]  # T in GG
    k4 = sqrt(rho_w[-1] * xsect[-1] / ten4) * w  # k in GG
    d_bend4 = sqrt(Y_h * xII[-1] / ten4)  # complex d_bend(4)
    dk4 = k4 * d_bend4

    # simp3a is inherited from the previous version of this suspension
    # thermal noise calculation (part of simp3).
    #
    # I'm not sure where this comes from, but it differs from
    # 1 by less than 1e-6 for frequencies below 100Hz (and less than
    # 1e-3 up to 10kHz).  What's more, the units appear to be wrong.
    # (Missing factor of rho * length?)
    # [mevans June 2015]
    #
    # simp3a = sqrt(1 + d_bend4 .* xsect(4) .* w.^2 / ten4)
    simp3a = 1

    coskl = simp3a * cos(k4 * length[-1])
    sinkl = sin(k4 * length[-1])

    # numerator, horiz spring constant, last stage
    #   numerator of K_xx in eq 9 of GG
    #     = T k (cos(k L) + k delta sin(k L))
    #   for w -> 0, this reduces to N_w * T * k
    kh4num = (
        N_w[-1] * ten4 * k4 * simp3a * (simp3a ** 2 + dk4 ** 2) * (coskl + dk4 * sinkl)
    )

    # denominator, horiz spring constant, last stage
    #   D after equation 8 in GG
    #   D = sin(k L) - 2 k delta cos(k L)
    #   for w -> 0, this reduces to k (L - 2 delta)
    kh4den = (simp3a ** 2 - dk4 ** 2) * sinkl - 2 * dk4 * coskl

    # horizontal spring constant, last stage
    #   K_xx in eq 9 of GG
    kh[-1, :] = kh4num / kh4den

    ###############################################################
    # Equations of motion for the system
    ###############################################################

    # Calculate TFs turning on the loss of each stage one by one
    hForce = Struct()
    vForce = Struct()
    hForce.singlylossy = np.zeros([len(sus.Stage), len(w)], dtype=complex)
    vForce.singlylossy = np.zeros([len(sus.Stage), len(w)], dtype=complex)

    for n in range(len(mass)):
        # horizontal
        # only the imaginary part of the specified stage is used.
        k = real(kh) + 1j * imag(
            [
                kh[0, :] * (n == 0),
                kh[1, :] * (n == 1),
                kh[2, :] * (n == 2),
                kh[3, :] * (n == 3),
            ]
        )
        # calculate TFs
        hForce.singlylossy[n, :] = tst_force_to_tst_displ(k, mass, f)

        # vertical
        # only the imaginary part of the specified stage is used
        k = real(kv) + 1j * imag(
            [
                kv[0, :] * (n == 0),
                kv[1, :] * (n == 1),
                kv[2, :] * (n == 2),
                kv[3, :] * (n == 3),
            ]
        )
        # calculate TFs
        vForce.singlylossy[n, :] = tst_force_to_tst_displ(k, mass, f)

    # calculate horizontal TFs with all losses on
    hForce.fullylossy = tst_force_to_tst_displ(kh, mass, f)
    hTable = top_displ_to_tst_displ(kh, mass, f)

    # calculate vertical TFs with all losses on
    vForce.fullylossy = tst_force_to_tst_displ(kv, mass, f)
    vTable = top_displ_to_tst_displ(kv, mass, f)

    return hForce, vForce, hTable, vTable


def suspBQuad(f, ifo):

0 Source : core.py
with GNU General Public License v3.0
from block-hczhai

    def imag(self):
        return np.imag(self)

    @staticmethod

0 Source : core.py
with GNU General Public License v3.0
from block-hczhai

    def imag(self):
        return np.imag(self)

    def quick_deflate(self):

0 Source : core.py
with GNU General Public License v3.0
from block-hczhai

    def imag(self):
        return np.imag(self)

    def deflate(self, cutoff=0):

0 Source : flat.py
with GNU General Public License v3.0
from block-hczhai

    def imag(self):
        return np.imag(self)

    def cast_assign(self, other):

0 Source : basic.py
with MIT License
from buds-lab

def clpmn(m, n, z, type=3):
    """Associated Legendre function of the first kind for complex arguments.

    Computes the associated Legendre function of the first kind of order m and
    degree n, ``Pmn(z)`` = :math:`P_n^m(z)`, and its derivative, ``Pmn'(z)``.
    Returns two arrays of size ``(m+1, n+1)`` containing ``Pmn(z)`` and
    ``Pmn'(z)`` for all orders from ``0..m`` and degrees from ``0..n``.

    Parameters
    ----------
    m : int
       ``|m|   <  = n``; the order of the Legendre function.
    n : int
       where ``n >= 0``; the degree of the Legendre function.  Often
       called ``l`` (lower case L) in descriptions of the associated
       Legendre function
    z : float or complex
        Input value.
    type : int, optional
       takes values 2 or 3
       2: cut on the real axis ``|x| > 1``
       3: cut on the real axis ``-1  <  x  <  1`` (default)

    Returns
    -------
    Pmn_z : (m+1, n+1) array
       Values for all orders ``0..m`` and degrees ``0..n``
    Pmn_d_z : (m+1, n+1) array
       Derivatives for all orders ``0..m`` and degrees ``0..n``

    See Also
    --------
    lpmn: associated Legendre functions of the first kind for real z

    Notes
    -----
    By default, i.e. for ``type=3``, phase conventions are chosen according
    to [1]_ such that the function is analytic. The cut lies on the interval
    (-1, 1). Approaching the cut from above or below in general yields a phase
    factor with respect to Ferrer's function of the first kind
    (cf. `lpmn`).

    For ``type=2`` a cut at ``|x| > 1`` is chosen. Approaching the real values
    on the interval (-1, 1) in the complex plane yields Ferrer's function
    of the first kind.

    References
    ----------
    .. [1] Zhang, Shanjie and Jin, Jianming. "Computation of Special
           Functions", John Wiley and Sons, 1996.
           https://people.sc.fsu.edu/~jburkardt/f_src/special_functions/special_functions.html
    .. [2] NIST Digital Library of Mathematical Functions
           https://dlmf.nist.gov/14.21

    """
    if not isscalar(m) or (abs(m) > n):
        raise ValueError("m must be  < = n.")
    if not isscalar(n) or (n  <  0):
        raise ValueError("n must be a non-negative integer.")
    if not isscalar(z):
        raise ValueError("z must be scalar.")
    if not(type == 2 or type == 3):
        raise ValueError("type must be either 2 or 3.")
    if (m  <  0):
        mp = -m
        mf, nf = mgrid[0:mp+1, 0:n+1]
        with ufuncs.errstate(all='ignore'):
            if type == 2:
                fixarr = where(mf > nf, 0.0,
                               (-1)**mf * gamma(nf-mf+1) / gamma(nf+mf+1))
            else:
                fixarr = where(mf > nf, 0.0, gamma(nf-mf+1) / gamma(nf+mf+1))
    else:
        mp = m
    p, pd = specfun.clpmn(mp, n, real(z), imag(z), type)
    if (m  <  0):
        p = p * fixarr
        pd = pd * fixarr
    return p, pd


def lqmn(m, n, z):

0 Source : func_dynamic.py
with Apache License 2.0
from cfsengineering

def longi_root_identification(A_longi):
    """Identifies the root of the longitudinal charcateristic equation

    Args:
        A_longi (2D matrix) : State space matrix in CONCISE form

    Returns:
        If roots are not complex conjugate the finction returns:
            None
        If roots are well complex conjugate the function returns:
            sp1 : first complex root of the short period mode
            sp2 : complex conjugate root of the short period mode
            ph1 : first complex root of the phugoid mode
            ph2 : complex conjugate root of the phugoid mode
    """

    # Compute eigenvalues and eigenvectors of longi CONCISE equations at trim
    eg_value_longi, eg_vector_longi = linalg.eig(A_longi)

    # TODO: maybe use as test
    # delta = a S^4 + bS^3 + cS^2 + dS  + e  , delta_longi  = [a,b,c,d,e]     from page 423 Cook
    # delta_longi = [1, \
    #                        -(m_q+x_u+z_w),\
    #                        (m_q*z_w - m_w*z_q) + (m_q*x_u - m_u*x_q) + (x_u*z_w - x_w*z_u) - m_theta,\
    #                        (m_theta*x_u-m_u*x_theta)+(m_theta*z_w-m_w*z_theta)+m_q*(x_w*z_u-x_u*z_w)+x_q*(m_u*z_w-m_w*z_u)+z_q*(m_w*x_u-m_u*x_w),\
    #                        m_theta*(x_w*z_u-x_u*z_w)+x_theta*(m_u*z_w-m_w*z_u)+z_theta*(m_w*x_u-m_u*x_w)]
    # roots_longi = np.roots(delta_longi)
    # print('\n Eigenvalues Longitudinal :')
    # print(eg_value_longi)
    # print('\n DEN TF roots Longitudinal :')
    # print(roots_longi)
    # print('\n Eigenvectors Longitudinal :')
    # print(eg_vector_longi)

    count = 0
    real_root = []
    complex_root = []
    for root in eg_value_longi:
        if isinstance(root, complex):
            count += 1
            complex_root.append(root)
        else:
            real_root.append(root)

    # count number of complex numbers
    if count == 2:  #  if phugoid not stable (real number):
        if real_root[0]   <  = real_root[0]:
            ph1 = real_root[0]
            ph2 = real_root[1]
        else:
            ph1 = real_root[1]
            ph2 = real_root[0]
        if np.imag(complex_root[0])  <  np.imag(complex_root[1]):
            sp1 = complex_root[0]
            sp2 = complex_root[1]
        else:
            sp1 = complex_root[1]
            sp2 = complex_root[0]
    else:  #  if phugoid and short period both complex:
        root1 = eg_value_longi[0]
        root2 = eg_value_longi[1]
        root3 = eg_value_longi[2]
        root4 = eg_value_longi[3]
        if (np.imag(root1)) ** 2 == (np.imag(root2)) ** 2:
            couple1 = [root1, root2]
            couple2 = [root3, root4]
        elif (np.imag(root1)) ** 2 == (np.imag(root3)) ** 2:
            couple1 = [root1, root3]
            couple2 = [root2, root4]
        else:
            couple1 = [root1, root4]
            couple2 = [root2, root3]

        if abs(np.imag(root1)) > abs(np.imag(couple2[0])) ** 2:
            sp1 = root1
            sp2 = couple1[1]
            ph1 = couple2[0]
            ph2 = couple2[1]
        else:
            ph1 = root1
            ph2 = couple1[1]
            sp1 = couple2[0]
            sp2 = couple2[1]

    # If short period roots are not correctly identified all longi roots are not identified and roots are not physically possible, and programm has to stop
    if sp1 != np.conj(sp2) or ph1 != np.conj(ph2):
        return (None, eg_value_longi)
    else:
        # Calculate eigen vectors magnitude :
        eg_vector_longi_magnitude = []
        for row in eg_vector_longi:
            new_row = []
            for vector in row:
                vector_magnitude = abs(vector)  # Module of complex number
                new_row.append(format(vector_magnitude, ".3g"))
            eg_vector_longi_magnitude.append(new_row)
        return (sp1, sp2, ph1, ph2, eg_value_longi, eg_vector_longi, eg_vector_longi_magnitude)


def longi_mode_characteristic(sp1, sp2, ph1, ph2, load_factor):

0 Source : func_dynamic.py
with Apache License 2.0
from cfsengineering

def direc_root_identification(A_direc):
    """identifies the root of the lateral charcateristic equation

    Args:
        roots_direc_list : a list of four roots

    Returns:
        If dutch roll roots are not complex conjugate, and spiral and roll not identified, the function returns:
            None
        If roots are well identified:
            roll : the root coresponding to the roll equation
            spiral : the root coresponding to the spiral equation
            dr1 : first complex root of the dutch roll mode
            dr1 : complex conjugate root of the dutch roll mode
    """

    # Compute eigenvalues and eigenvectors of lateral CONCISE equations at trim
    eg_value_direc, eg_vector_direc = linalg.eig(A_direc)

    # TODO: from book, maybe use as test
    # # In cincise form,   delta = aS^5 + bS^4 + cS^3+ dS^2 + eS + f ,  delta_direc =[a,b,c,d,e,f]
    # delta_direc = [1,\
    #                         -(l_p+n_r+y_v),
    #                         (l_p*n_r-l_r*n_p)+(n_r*y_v-n_v*y_r)+(l_p*y_v-l_v*y_p)-(l_phi+n_psi),\
    #                         (l_p*n_psi-l_psi*n_p)+(l_phi*n_r-l_r*n_phi)+l_v*(n_r*y_p-n_p*y_r-y_phi)+n_v*(l_p*y_r-l_r*y_p-y_psi)+y_v*(l_r*n_p-l_p*n_r+l_phi+n_psi),\
    #                         (l_phi*n_psi-l_psi*n_psi)+l_v*((n_r*y_phi-n_phi*y_r)+(n_psi*y_p-n_p*y_psi))+n_v*((l_phi*y_r-l_r*y_phi)+(l_p*y_psi-l_psi*y_p))+y_v*((l_r*n_phi-l_phi*n_r)+(l_psi*n_p-l_p*n_psi)),\
    #                         l_v*(n_psi*y_phi-n_phi*y_psi)+n_v*(l_phi*y_psi-l_psi*y_phi)+y_v*(l_psi*n_phi-l_phi*n_psi)]

    roll = None
    spiral = None
    dr1 = None
    dr2 = None
    reals = []

    for root in eg_value_direc:
        if root.imag != 0:  # if complex root
            if np.imag(root) > 0:  # If imaginary part >0
                dr1 = root  # Dutch Roll Root1
            else:
                dr2 = root  # Dutch Roll Root2
        else:  # If real root
            reals.append(root.real)

    absolute = [abs(value) for value in reals]
    absolute.sort(reverse=True)

    for elem in reals:
        if abs(elem) == absolute[0]:
            roll = elem
        elif abs(elem) == absolute[1]:
            spiral = elem
        else:
            pass

    if (
        dr1 == None or dr2 == None or roll == None or spiral == None
    ):  # Check that all variable have been assigned a value
        return (None, eg_value_direc)
    elif dr1 != np.conj(dr2):  # Check that dr roots are complex conjugate
        return (None, eg_value_direc)
    else:
        # Calculate eigen vectors magnitude :
        eg_vector_direc_magnitude = []
        for row in eg_vector_direc:
            new_row = []
            for vector in row:
                vector_magnitude = abs(vector)  # Module of complex number
                new_row.append(format(vector_magnitude, ".4g"))
            eg_vector_direc_magnitude.append(new_row)
        return (roll, spiral, dr1, dr2, eg_value_direc, eg_vector_direc, eg_vector_direc_magnitude)


def direc_mode_characteristic(roll, spiral, dr1, dr2):

See More Examples