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
3
Source : _ode.py
with GNU General Public License v3.0
from adityaprakash-bobby
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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