Python numpy.polyval() Examples

The following are 30 code examples for showing how to use numpy.polyval(). These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example.

You may check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module numpy , or try the search function .

Example 1
Project: pyberny   Author: jhrmnn   File: Math.py    License: Mozilla Public License 2.0 6 votes vote down vote up
def fit_cubic(y0, y1, g0, g1):
    """Fit cubic polynomial to function values and derivatives at x = 0, 1.

    Returns position and function value of minimum if fit succeeds. Fit does
    not succeeds if

    1. polynomial doesn't have extrema or
    2. maximum is from (0,1) or
    3. maximum is closer to 0.5 than minimum
    """
    a = 2 * (y0 - y1) + g0 + g1
    b = -3 * (y0 - y1) - 2 * g0 - g1
    p = np.array([a, b, g0, y0])
    r = np.roots(np.polyder(p))
    if not np.isreal(r).all():
        return None, None
    r = sorted(x.real for x in r)
    if p[0] > 0:
        maxim, minim = r
    else:
        minim, maxim = r
    if 0 < maxim < 1 and abs(minim - 0.5) > abs(maxim - 0.5):
        return None, None
    return minim, np.polyval(p, minim) 
Example 2
Project: ocelot   Author: ocelot-collab   File: wave.py    License: GNU General Public License v3.0 6 votes vote down vote up
def calc_phase_delay(coeff, w, w0):
    """
    expression for the phase -- coeff[0] + coeff[1]*(w - w0)/1! + coeff[2]*(w - w0)**2/2! + coeff[3]*(w - w0)**3/3!
    coeff is a list with
    coeff[0] =: measured in [rad]      --- phase
    coeff[1] =: measured in [fm s ^ 1] --- group delay
    coeff[2] =: measured in [fm s ^ 2] --- group delay dispersion (GDD)
    coeff[3] =: measured in [fm s ^ 3] --- third-order dispersion (TOD)
    ...
    """
    delta_w = w - w0
    _logger.debug('calculating phase delay')
    _logger.debug(ind_str + 'coeffs for compression = {}'.format(coeff))
    coeff_norm = [ci / (1e15) ** i / factorial(i) for i, ci in enumerate(coeff)]
    coeff_norm = list(coeff_norm)[::-1]
    _logger.debug(ind_str + 'coeffs_norm = {}'.format(coeff_norm))
    delta_phi = np.polyval(coeff_norm, delta_w)
    _logger.debug(ind_str + 'delta_phi[0] = {}'.format(delta_phi[0]))
    _logger.debug(ind_str + 'delta_phi[-1] = {}'.format(delta_phi[-1]))
    _logger.debug(ind_str + 'done')

    return delta_phi 
Example 3
Project: pyiron   Author: pyiron   File: thermo_bulk.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def get_minimum_energy_path(self, pressure=None):
        """

        Args:
            pressure:

        Returns:

        """
        if pressure is not None:
            raise NotImplemented()
        v_min_lst = []
        for c in self._coeff.T:
            v_min = np.roots(np.polyder(c, 1))
            p_der2 = np.polyder(c, 2)
            p_val2 = np.polyval(p_der2, v_min)
            v_m_lst = v_min[p_val2 > 0]
            if len(v_m_lst) > 0:
                v_min_lst.append(v_m_lst[0])
            else:
                v_min_lst.append(np.nan)
        return np.array(v_min_lst) 
Example 4
Project: pyiron   Author: pyiron   File: thermo_bulk.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def interpolate_volume(self, volumes, fit_order=None):
        """

        Args:
            volumes:
            fit_order:

        Returns:

        """
        if fit_order is not None:
            self._fit_order = fit_order
        new = self.copy()
        new.volumes = volumes
        new.energies = np.array([np.polyval(self._coeff, v) for v in volumes]).T
        return new 
Example 5
Project: pyiron   Author: pyiron   File: thermo_bulk.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def contour_pressure(self):
        """

        Returns:

        """
        try:
            import pylab as plt
        except ImportError:
            import matplotlib.pyplot as plt
        x, y = self.meshgrid()
        p_coeff = np.polyfit(self.volumes, self.pressure.T, deg=self._fit_order)
        p_grid = np.array([np.polyval(p_coeff, v) for v in self._volumes]).T
        plt.contourf(x, y, p_grid)
        plt.plot(self.get_minimum_energy_path(), self.temperatures)
        plt.xlabel("Volume [$\AA^3$]")
        plt.ylabel("Temperature [K]") 
Example 6
Project: pyiron   Author: pyiron   File: thermo_bulk.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def contour_entropy(self):
        """

        Returns:

        """
        try:
            import pylab as plt
        except ImportError:
            import matplotlib.pyplot as plt
        s_coeff = np.polyfit(self.volumes, self.entropy.T, deg=self._fit_order)
        s_grid = np.array([np.polyval(s_coeff, v) for v in self.volumes]).T
        x, y = self.meshgrid()
        plt.contourf(x, y, s_grid)
        plt.plot(self.get_minimum_energy_path(), self.temperatures)
        plt.xlabel("Volume [$\AA^3$]")
        plt.ylabel("Temperature [K]") 
Example 7
Project: allantools   Author: aewallin   File: ci.py    License: GNU Lesser General Public License v3.0 6 votes vote down vote up
def detrend(x, deg=1):
    """
    remove polynomial from data.
    used by autocorr_noise_id()

    Parameters
    ----------
    x: numpy.array
        time-series
    deg: int
        degree of polynomial to remove from x

    Returns
    -------
    x_detrended: numpy.array
        detrended time-series
    """
    t = range(len(x))
    p = np.polyfit(t, x, deg)
    residual = x - np.polyval(p, t)
    return residual

########################################################################
# Equivalent Degrees of Freedom 
Example 8
Project: pvlib-python   Author: pvlib   File: test_pvsystem.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_PVSystem_sapm_effective_irradiance(sapm_module_params, mocker):
    system = pvsystem.PVSystem(module_parameters=sapm_module_params)
    mocker.spy(pvsystem, 'sapm_effective_irradiance')

    poa_direct = 900
    poa_diffuse = 100
    airmass_absolute = 1.5
    aoi = 0
    p = (sapm_module_params['A4'], sapm_module_params['A3'],
         sapm_module_params['A2'], sapm_module_params['A1'],
         sapm_module_params['A0'])
    f1 = np.polyval(p, airmass_absolute)
    expected = f1 * (poa_direct + sapm_module_params['FD'] * poa_diffuse)
    out = system.sapm_effective_irradiance(
        poa_direct, poa_diffuse, airmass_absolute, aoi)
    pvsystem.sapm_effective_irradiance.assert_called_once_with(
        poa_direct, poa_diffuse, airmass_absolute, aoi, sapm_module_params)
    assert_allclose(out, expected, atol=0.1) 
Example 9
Project: GraphicDesignPatternByPython   Author: Relph1119   File: test_savitzky_golay.py    License: MIT License 6 votes vote down vote up
def alt_sg_coeffs(window_length, polyorder, pos):
    """This is an alternative implementation of the SG coefficients.

    It uses numpy.polyfit and numpy.polyval.  The results should be
    equivalent to those of savgol_coeffs(), but this implementation
    is slower.

    window_length should be odd.

    """
    if pos is None:
        pos = window_length // 2
    t = np.arange(window_length)
    unit = (t == pos).astype(int)
    h = np.polyval(np.polyfit(t, unit, polyorder), t)
    return h 
Example 10
Project: GraphicDesignPatternByPython   Author: Relph1119   File: test_dltisys.py    License: MIT License 6 votes vote down vote up
def test_auto(self):
        # Test dfreqresp() real part calculation.
        # 1st order low-pass filter: H(z) = 1 / (z - 0.2),
        system = TransferFunction(1, [1, -0.2], dt=0.1)
        w = [0.1, 1, 10, 100]
        w, H = dfreqresp(system, w=w)
        jw = np.exp(w * 1j)
        y = np.polyval(system.num, jw) / np.polyval(system.den, jw)

        # test real
        expected_re = y.real
        assert_almost_equal(H.real, expected_re)

        # test imag
        expected_im = y.imag
        assert_almost_equal(H.imag, expected_im) 
Example 11
Project: GraphicDesignPatternByPython   Author: Relph1119   File: test_dltisys.py    License: MIT License 6 votes vote down vote up
def test_auto(self):
        # Test bode() magnitude calculation.
        # 1st order low-pass filter: H(s) = 0.3 / (z - 0.2),
        system = TransferFunction(0.3, [1, -0.2], dt=0.1)
        w = np.array([0.1, 0.5, 1, np.pi])
        w2, mag, phase = dbode(system, w=w)
        jw = np.exp(w * 1j)
        y = np.polyval(system.num, jw) / np.polyval(system.den, jw)

        # Test mag
        expected_mag = 20.0 * np.log10(abs(y))
        assert_almost_equal(mag, expected_mag)

        # Test phase
        expected_phase = np.rad2deg(np.angle(y))
        assert_almost_equal(phase, expected_phase) 
Example 12
Project: oopt-gnpy   Author: Telecominfraproject   File: elements.py    License: BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def _nf(self, type_def, nf_model, nf_fit_coeff, gain_min, gain_flatmax, gain_target):
        # if hybrid raman, use edfa_gain_flatmax attribute, else use gain_flatmax
        #gain_flatmax = getattr(params, 'edfa_gain_flatmax', params.gain_flatmax)
        pad = max(gain_min - gain_target, 0)
        gain_target += pad
        dg = max(gain_flatmax - gain_target, 0)
        if type_def == 'variable_gain':
            g1a = gain_target - nf_model.delta_p - dg
            nf_avg = lin2db(db2lin(nf_model.nf1) + db2lin(nf_model.nf2) / db2lin(g1a))
        elif type_def == 'fixed_gain':
            nf_avg = nf_model.nf0
        elif type_def == 'openroadm':
            pin_ch = self.pin_db - lin2db(self.nch)
            # model OSNR = f(Pin)
            nf_avg = pin_ch - polyval(nf_model.nf_coef, pin_ch) + 58
        elif type_def == 'advanced_model':
            nf_avg = polyval(nf_fit_coeff, -dg)
        else:
            assert False, "Unrecognized amplifier type, this should have been checked by the JSON loader"
        return nf_avg + pad, pad 
Example 13
Project: NeuroKit   Author: neuropsychology   File: fractal_correlation.py    License: MIT License 5 votes vote down vote up
def _fractal_correlation_plot(r_vals, corr, d2):
    fit = 2 ** np.polyval(d2, np.log2(r_vals))
    plt.loglog(r_vals, corr, "bo")
    plt.loglog(r_vals, fit, "r", label=r"$D2$ = %0.3f" % d2[0])
    plt.title("Correlation Dimension")
    plt.xlabel(r"$\log_{2}$(r)")
    plt.ylabel(r"$\log_{2}$(c)")
    plt.legend()
    plt.show() 
Example 14
Project: NeuroKit   Author: neuropsychology   File: fractal_dfa.py    License: MIT License 5 votes vote down vote up
def _fractal_dfa_trends(segments, window, order=1):
    x = np.arange(window)

    coefs = np.polyfit(x[:window], segments.T, order).T

    # TODO: Could this be optimized? Something like np.polyval(x[:window], coefs)
    trends = np.array([np.polyval(coefs[j], x) for j in np.arange(len(segments))])

    return trends 
Example 15
Project: NeuroKit   Author: neuropsychology   File: fractal_dfa.py    License: MIT License 5 votes vote down vote up
def _fractal_dfa_plot(windows, fluctuations, dfa):
    fluctfit = 2 ** np.polyval(dfa, np.log2(windows))
    plt.loglog(windows, fluctuations, "bo")
    plt.loglog(windows, fluctfit, "r", label=r"$\alpha$ = %0.3f" % dfa[0])
    plt.title("DFA")
    plt.xlabel(r"$\log_{2}$(Window)")
    plt.ylabel(r"$\log_{2}$(Fluctuation)")
    plt.legend()
    plt.show()


# =============================================================================
#  Utils MDDFA
# ============================================================================= 
Example 16
Project: NeuroKit   Author: neuropsychology   File: fit_polynomial.py    License: MIT License 5 votes vote down vote up
def _fit_polynomial(y, X, order=2):
    # Generating weights and model for polynomial function with a given degree
    y_predicted = np.polyval(np.polyfit(X, y, order), X)
    return y_predicted 
Example 17
Project: NeuroKit   Author: neuropsychology   File: eventrelated_utils.py    License: MIT License 5 votes vote down vote up
def _eventrelated_rate(epoch, output={}, var="ECG_Rate"):

    # Sanitize input
    colnames = epoch.columns.values
    if len([i for i in colnames if var in i]) == 0:
        print(
            "NeuroKit warning: *_eventrelated(): input does not"
            "have an `" + var + "` column. Will skip all rate-related features."
        )
        return output

    # Get baseline
    zero = find_closest(0, epoch.index.values, return_index=True)  # Find closest to 0
    baseline = epoch[var].iloc[zero]

    signal = epoch[var].values[zero + 1 : :]
    index = epoch.index.values[zero + 1 : :]

    # Max / Min / Mean
    output[var + "_Baseline"] = baseline
    output[var + "_Max"] = np.max(signal) - baseline
    output[var + "_Min"] = np.min(signal) - baseline
    output[var + "_Mean"] = np.mean(signal) - baseline

    # Time of Max / Min
    output[var + "_Max_Time"] = index[np.argmax(signal)]
    output[var + "_Min_Time"] = index[np.argmin(signal)]

    # Modelling
    # These are experimental indices corresponding to parameters of a quadratic model
    # Instead of raw values (such as min, max etc.)
    coefs = np.polyfit(index, signal - baseline, 2)
    output[var + "_Trend_Quadratic"] = coefs[0]
    output[var + "_Trend_Linear"] = coefs[1]
    output[var + "_Trend_R2"] = fit_r2(
        y=signal - baseline, y_predicted=np.polyval(coefs, index), adjusted=False, n_parameters=3
    )

    return output 
Example 18
Project: pyberny   Author: jhrmnn   File: Math.py    License: Mozilla Public License 2.0 5 votes vote down vote up
def fit_quartic(y0, y1, g0, g1):
    """Fit constrained quartic polynomial to function values and erivatives at x = 0,1.

    Returns position and function value of minimum or None if fit fails or has
    a maximum. Quartic polynomial is constrained such that it's 2nd derivative
    is zero at just one point. This ensures that it has just one local
    extremum.  No such or two such quartic polynomials always exist. From the
    two, the one with lower minimum is chosen.
    """

    def g(y0, y1, g0, g1, c):
        a = c + 3 * (y0 - y1) + 2 * g0 + g1
        b = -2 * c - 4 * (y0 - y1) - 3 * g0 - g1
        return np.array([a, b, c, g0, y0])

    def quart_min(p):
        r = np.roots(np.polyder(p))
        is_real = np.isreal(r)
        if is_real.sum() == 1:
            minim = r[is_real][0].real
        else:
            minim = r[(r == max(-abs(r))) | (r == -max(-abs(r)))][0].real
        return minim, np.polyval(p, minim)

    # discriminant of d^2y/dx^2=0
    D = -((g0 + g1) ** 2) - 2 * g0 * g1 + 6 * (y1 - y0) * (g0 + g1) - 6 * (y1 - y0) ** 2
    if D < 1e-11:
        return None, None
    else:
        m = -5 * g0 - g1 - 6 * y0 + 6 * y1
        p1 = g(y0, y1, g0, g1, 0.5 * (m + np.sqrt(2 * D)))
        p2 = g(y0, y1, g0, g1, 0.5 * (m - np.sqrt(2 * D)))
        if p1[0] < 0 and p2[0] < 0:
            return None, None
        [minim1, minval1] = quart_min(p1)
        [minim2, minval2] = quart_min(p2)
        if minval1 < minval2:
            return minim1, minval1
        else:
            return minim2, minval2 
Example 19
Project: pyqmc   Author: WagnerGroup   File: linemin.py    License: MIT License 5 votes vote down vote up
def polyfit_relative(xfit, yfit, degree):
    p = np.polyfit(xfit, yfit, degree)
    ypred = np.polyval(p, xfit)
    resid = (ypred - yfit) ** 2
    relative_error = np.var(resid) / np.var(yfit)
    return p, relative_error 
Example 20
Project: ibllib   Author: int-brain-lab   File: test_ephysqc.py    License: MIT License 5 votes vote down vote up
def test_impeccable_dataset(self):

        fpga2bpod = np.array([11 * 1e-6, -20])  # bpod starts 20 secs before with 10 ppm drift
        fpga_trials = {
            'intervals': np.array([[0, 9.5], [10, 19.5]]),
            'stimOn_times': np.array([2, 12]),
            'goCue_times': np.array([2.0001, 12.0001]),
            'ready_tone_in': np.array([2.0001, 12.0001]),
            'stim_freeze': np.array([4., 14.]),
            'feedback_times': np.array([4.0001, 14.0001]),
            'error_tone_in': np.array([4.0001, np.nan]),
            'valve_open': np.array([np.nan, 14.0001]),
            'stimOff_times': np.array([6.0001, 15.0001]),
            'iti_in': np.array([6.0011, 15.000]),
        }

        alf_trials = {
            'goCueTrigger_times_bpod': np.polyval(fpga2bpod, fpga_trials['goCue_times'] - 0.00067),
            'response_times_bpod': np.polyval(fpga2bpod, np.array([4., 14.])),
            'intervals_bpod': np.polyval(fpga2bpod, fpga_trials['intervals']),
            # Times from session start
            'goCueTrigger_times': fpga_trials['goCue_times'] - 0.00067,
            'response_times': np.array([4., 14.]),
            'intervals': fpga_trials['intervals'],
            'stimOn_times': fpga_trials['stimOn_times'],
            'goCue_times': fpga_trials['goCue_times'],
            'feedback_times': fpga_trials['feedback_times'],
        }
        qcs, qct = ephysqc.qc_fpga_task(fpga_trials, alf_trials)
        self.assertTrue(np.all([qcs[k] for k in qcs]))
        self.assertTrue(np.all([np.all(qct[k]) for k in qct])) 
Example 21
Project: ibllib   Author: int-brain-lab   File: test_extractors.py    License: MIT License 5 votes vote down vote up
def test_sync_bpod_bonsai_poor_quality_timestamps(self):
        sync_trials_robust = raw.sync_trials_robust
        drift_pol = np.array([11 * 1e-6, -20])  # bpod starts 20 secs before with 10 ppm drift
        np.random.seed(seed=784)
        t0_full = np.cumsum(np.random.rand(50)) + .001
        t1_full = np.polyval(drift_pol, t0_full) + t0_full
        t0 = t0_full.copy()
        t1 = t1_full.copy()

        t0_, t1_ = sync_trials_robust(t0, t1)
        assert np.allclose(t1_, np.polyval(drift_pol, t0_) + t0_)

        t0_, t1_ = sync_trials_robust(t0, t1[:-1])
        assert np.allclose(t1_, np.polyval(drift_pol, t0_) + t0_)

        t0_, t1_ = sync_trials_robust(t0, t1[1:])
        assert np.allclose(t1_, np.polyval(drift_pol, t0_) + t0_)

        t0_, t1_ = sync_trials_robust(t0[1:], t1)
        assert np.allclose(t1_, np.polyval(drift_pol, t0_) + t0_)

        t0_, t1_ = sync_trials_robust(t0[:-1], t1)
        assert np.allclose(t1_, np.polyval(drift_pol, t0_) + t0_)

        t0_, t1_ = sync_trials_robust(t0, np.delete(t1, 24))
        assert np.allclose(t1_, np.polyval(drift_pol, t0_) + t0_)

        t0_, t1_ = sync_trials_robust(np.delete(t0, 12), np.delete(t1, 24))
        assert np.allclose(t1_, np.polyval(drift_pol, t0_) + t0_) 
Example 22
Project: lambda-packs   Author: ryfeus   File: _continuous_distns.py    License: MIT License 5 votes vote down vote up
def _stats(self, s):
        p = np.exp(s*s)
        mu = np.sqrt(p)
        mu2 = p*(p-1)
        g1 = np.sqrt((p-1))*(2+p)
        g2 = np.polyval([1, 2, 3, 0, -6.0], p)
        return mu, mu2, g1, g2 
Example 23
Project: lambda-packs   Author: ryfeus   File: _continuous_distns.py    License: MIT License 5 votes vote down vote up
def _stats(self):
        p = np.e
        mu = np.sqrt(p)
        mu2 = p * (p - 1)
        g1 = np.sqrt((p - 1)) * (2 + p)
        g2 = np.polyval([1, 2, 3, 0, -6.0], p)
        return mu, mu2, g1, g2 
Example 24
Project: lambda-packs   Author: ryfeus   File: _continuous_distns.py    License: MIT License 5 votes vote down vote up
def _stats(self, b, moments='mv'):
        mu, mu2, g1, g2 = None, None, None, None
        if 'm' in moments:
            mask = b > 1
            bt = np.extract(mask, b)
            mu = valarray(np.shape(b), value=np.inf)
            np.place(mu, mask, bt / (bt-1.0))
        if 'v' in moments:
            mask = b > 2
            bt = np.extract(mask, b)
            mu2 = valarray(np.shape(b), value=np.inf)
            np.place(mu2, mask, bt / (bt-2.0) / (bt-1.0)**2)
        if 's' in moments:
            mask = b > 3
            bt = np.extract(mask, b)
            g1 = valarray(np.shape(b), value=np.nan)
            vals = 2 * (bt + 1.0) * np.sqrt(bt - 2.0) / ((bt - 3.0) * np.sqrt(bt))
            np.place(g1, mask, vals)
        if 'k' in moments:
            mask = b > 4
            bt = np.extract(mask, b)
            g2 = valarray(np.shape(b), value=np.nan)
            vals = (6.0*np.polyval([1.0, 1.0, -6, -2], bt) /
                    np.polyval([1.0, -7.0, 12.0, 0.0], bt))
            np.place(g2, mask, vals)
        return mu, mu2, g1, g2 
Example 25
Project: lambda-packs   Author: ryfeus   File: _continuous_distns.py    License: MIT License 5 votes vote down vote up
def _stats(self, a):
        return (a / (a + 1.0),
                a / (a + 2.0) / (a + 1.0) ** 2,
                -2.0 * ((a - 1.0) / (a + 3.0)) * np.sqrt((a + 2.0) / a),
                6 * np.polyval([1, -1, -6, 2], a) / (a * (a + 3.0) * (a + 4))) 
Example 26
Project: auto-alt-text-lambda-api   Author: abhisuri97   File: test_polynomial.py    License: MIT License 5 votes vote down vote up
def test_polyfit(self):
        c = np.array([3., 2., 1.])
        x = np.linspace(0, 2, 7)
        y = np.polyval(c, x)
        err = [1, -1, 1, -1, 1, -1, 1]
        weights = np.arange(8, 1, -1)**2/7.0

        # check 1D case
        m, cov = np.polyfit(x, y+err, 2, cov=True)
        est = [3.8571, 0.2857, 1.619]
        assert_almost_equal(est, m, decimal=4)
        val0 = [[2.9388, -5.8776, 1.6327],
                [-5.8776, 12.7347, -4.2449],
                [1.6327, -4.2449, 2.3220]]
        assert_almost_equal(val0, cov, decimal=4)

        m2, cov2 = np.polyfit(x, y+err, 2, w=weights, cov=True)
        assert_almost_equal([4.8927, -1.0177, 1.7768], m2, decimal=4)
        val = [[8.7929, -10.0103, 0.9756],
               [-10.0103, 13.6134, -1.8178],
               [0.9756, -1.8178, 0.6674]]
        assert_almost_equal(val, cov2, decimal=4)

        # check 2D (n,1) case
        y = y[:, np.newaxis]
        c = c[:, np.newaxis]
        assert_almost_equal(c, np.polyfit(x, y, 2))
        # check 2D (n,2) case
        yy = np.concatenate((y, y), axis=1)
        cc = np.concatenate((c, c), axis=1)
        assert_almost_equal(cc, np.polyfit(x, yy, 2))

        m, cov = np.polyfit(x, yy + np.array(err)[:, np.newaxis], 2, cov=True)
        assert_almost_equal(est, m[:, 0], decimal=4)
        assert_almost_equal(est, m[:, 1], decimal=4)
        assert_almost_equal(val0, cov[:, :, 0], decimal=4)
        assert_almost_equal(val0, cov[:, :, 1], decimal=4) 
Example 27
Project: sfa-numpy   Author: spatialaudio   File: angular.py    License: MIT License 5 votes vote down vote up
def Legendre_matrix(N, ctheta):
    r"""Matrix of weighted Legendre Polynominals.

    Computes a matrix of weighted Legendre polynominals up to order N for
    the given angles

    .. math::

        L_n(\theta) = \frac{2n+1}{4 \pi} P_n(\theta)

    Parameters
    ----------
    N : int
        Maximum order.
    ctheta : (Q,) array_like
        Angles.

    Returns
    -------
    Lmn : ((N+1), Q) numpy.ndarray
        Matrix containing Legendre polynominals.
    """
    if ctheta.ndim == 0:
        M = 1
    else:
        M = len(ctheta)
    Lmn = np.zeros([N+1, M], dtype=complex)
    for n in range(N+1):
        Lmn[n, :] = (2*n+1)/(4*np.pi) * np.polyval(special.legendre(n), ctheta)

    return Lmn 
Example 28
Project: sun-position   Author: s-bear   File: sunposition.py    License: MIT License 5 votes vote down vote up
def heliocentric_longitude(jme):
        """Compute the Earth Heliocentric Longitude (L) in degrees given the Julian Ephemeris Millennium"""
        #L5, ..., L0
        Li = [sum(a*np.cos(b + c*jme) for a,b,c in abcs) for abcs in reversed(_sp._EHL_)]
        L = np.polyval(Li, jme) / 1e8
        L = np.rad2deg(L) % 360
        return L 
Example 29
Project: sun-position   Author: s-bear   File: sunposition.py    License: MIT License 5 votes vote down vote up
def heliocentric_latitude(jme):
        """Compute the Earth Heliocentric Latitude (B) in degrees given the Julian Ephemeris Millennium"""
        Bi = [sum(a*np.cos(b + c*jme) for a,b,c in abcs) for abcs in reversed(_sp._EHB_)]
        B = np.polyval(Bi, jme) / 1e8
        B = np.rad2deg(B) % 360
        return B 
Example 30
Project: sun-position   Author: s-bear   File: sunposition.py    License: MIT License 5 votes vote down vote up
def heliocentric_radius(jme):
        """Compute the Earth Heliocentric Radius (R) in astronimical units given the Julian Ephemeris Millennium"""
        Ri = [sum(a*np.cos(b + c*jme) for a,b,c in abcs) for abcs in reversed(_sp._EHR_)]
        R = np.polyval(Ri, jme) / 1e8
        return R