Python numpy.sinc() Examples

The following are 30 code examples of numpy.sinc(). 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 also want to check out all available functions/classes of the module numpy , or try the search function .
Example #1
Source File: TFMethods.py    From ASP with GNU General Public License v3.0 6 votes vote down vote up
def sincinterp(x):
        """
        Sinc interpolation for computation of fractional transformations.
        As appears in :
        -https://github.com/audiolabs/frft/
        ----------
        Args:
            f       : (array) Complex valued input array
            a       : (float) Alpha factor
        Returns:
            ret     : (array) Real valued synthesised data
        """
        N = len(x)
        y = np.zeros(2 * N - 1, dtype=x.dtype)
        y[:2 * N:2] = x
        xint = fftconvolve( y[:2 * N], np.sinc(np.arange(-(2 * N - 3), (2 * N - 2)).T / 2),)
        return xint[2 * N - 3: -2 * N + 3] 
Example #2
Source File: fourier_test.py    From odl with Mozilla Public License 2.0 6 votes vote down vote up
def test_fourier_trafo_charfun_1d():
    # Characteristic function of [0, 1], its Fourier transform is
    # given by exp(-1j * y / 2) * sinc(y/2)
    def char_interval(x):
        return (x >= 0) & (x <= 1)

    def char_interval_ft(x):
        return np.exp(-1j * x / 2) * sinc(x / 2) / np.sqrt(2 * np.pi)

    # Base version
    discr = odl.uniform_discr(-2, 2, 40, impl='numpy')
    dft_base = FourierTransform(discr)

    # Complex version, should be as good
    discr = odl.uniform_discr(-2, 2, 40, impl='numpy', dtype='complex64')
    dft_complex = FourierTransform(discr)

    # Without shift
    discr = odl.uniform_discr(-2, 2, 40, impl='numpy', dtype='complex64')
    dft_complex_shift = FourierTransform(discr, shift=False)

    for dft in [dft_base, dft_complex, dft_complex_shift]:
        func_true_ft = dft.range.element(char_interval_ft)
        func_dft = dft(char_interval)
        assert (func_dft - func_true_ft).norm() < 5e-6 
Example #3
Source File: test_wrapper_bohamiann.py    From RoBO with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def setUp(self):
        X_task_1 = np.random.rand(10, 2)
        y_task_1 = np.sinc(X_task_1 * 10 - 5).sum(axis=1)
        X_task_1 = np.concatenate((X_task_1, np.zeros([10, 1])), axis=1)

        X_task_2 = np.random.rand(10, 2)
        y_task_2 = np.sinc(X_task_2 * 2 - 4).sum(axis=1)
        X_task_2 = np.concatenate((X_task_2, np.ones([10, 1])), axis=1)

        X_task_3 = np.random.rand(10, 2)
        y_task_3 = np.sinc(X_task_3 * 8 - 6).sum(axis=1)
        X_task_3 = np.concatenate((X_task_3, 2 * np.ones([10, 1])), axis=1)

        self.X = np.concatenate((X_task_1, X_task_2, X_task_3), axis=0)
        self.y = np.concatenate((y_task_1, y_task_2, y_task_3), axis=0)
        self.model = WrapperBohamiann()
        self.model.train(self.X, self.y) 
Example #4
Source File: test_incumbent_estimation.py    From RoBO with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def test_projected_incumbent_estimation(self):
        X = np.random.randn(20, 10)
        y = np.sinc(X).sum(axis=1)

        class DemoModel(object):
            def train(self, X, y):
                self.X = X
                self.y = y

            def predict(self, X):
                return self.y, np.ones(self.y.shape[0])

        model = DemoModel()
        model.train(X, y)
        inc, inc_val = incumbent_estimation.projected_incumbent_estimation(model, X, proj_value=1)
        b = np.argmin(y)

        assert inc[-1] == 1
        assert np.all(inc[:-1] == X[b])
        assert inc_val == y[b] 
Example #5
Source File: fourier_test.py    From odl with Mozilla Public License 2.0 6 votes vote down vote up
def test_fourier_trafo_freq_shifted_charfun_1d():
    # Frequency-shifted characteristic function: mult. with
    # exp(-1j * b * x) corresponds to shifting the FT by b.
    def fshift_char_interval(x):
        return np.exp(-1j * x * np.pi) * ((x >= -0.5) & (x <= 0.5))

    def fshift_char_interval_ft(x):
        return sinc((x + np.pi) / 2) / np.sqrt(2 * np.pi)

    # Number of points is very important here (aliasing)
    discr = odl.uniform_discr(-2, 2, 400, impl='numpy',
                              dtype='complex64')
    dft = FourierTransform(discr)
    func_true_ft = dft.range.element(fshift_char_interval_ft)
    func_dft = dft(fshift_char_interval)
    assert (func_dft - func_true_ft).norm() < 0.001 
Example #6
Source File: fourier_test.py    From odl with Mozilla Public License 2.0 6 votes vote down vote up
def test_dft_with_known_pairs_2d():

    # Frequency-shifted product of characteristic functions
    def fshift_char_rect(x):
        # Characteristic function of the cuboid
        # [-1, 1] x [1, 2]
        return (x[0] >= -1) & (x[0] <= 1) & (x[1] >= 1) & (x[1] <= 2)

    def fshift_char_rect_ft(x):
        # FT is a product of shifted and frequency-shifted sinc functions
        # 1st comp.: 2 * sinc(y)
        # 2nd comp.: exp(-1j * y * 3/2) * sinc(y/2)
        # Overall factor: (2 * pi)^(-1)
        return (
            2 * sinc(x[0])
            * np.exp(-1j * x[1] * 3 / 2) * sinc(x[1] / 2)
            / (2 * np.pi)
        )

    discr = odl.uniform_discr([-2] * 2, [2] * 2, (100, 400), impl='numpy',
                              dtype='complex64')
    dft = FourierTransform(discr)
    func_true_ft = dft.range.element(fshift_char_rect_ft)
    func_dft = dft(fshift_char_rect)
    assert (func_dft - func_true_ft).norm() < 0.001 
Example #7
Source File: interpolation.py    From scarlet with MIT License 6 votes vote down vote up
def lanczos(dx, a=3):
    """Lanczos kernel

    Parameters
    ----------
    dx: float
        amount to shift image
    a: int
        Lanczos window size parameter

    Returns
    -------
    result: array-like
        1D Lanczos kernel
    """
    if np.abs(dx) > 1:
        raise ValueError("The fractional shift dx must be between -1 and 1")
    window = np.arange(-a + 1, a + 1) + np.floor(dx)
    y = np.sinc(dx - window) * np.sinc((dx - window) / a)
    return y, window.astype(int) 
Example #8
Source File: fourier_test.py    From odl with Mozilla Public License 2.0 6 votes vote down vote up
def test_fourier_trafo_scaling():
    # Test if the FT scales correctly

    # Characteristic function of [0, 1], its Fourier transform is
    # given by exp(-1j * y / 2) * sinc(y/2)
    def char_interval(x):
        return (x >= 0) & (x <= 1)

    def char_interval_ft(x):
        return np.exp(-1j * x / 2) * sinc(x / 2) / np.sqrt(2 * np.pi)

    discr = odl.uniform_discr(-2, 2, 40, impl='numpy', dtype='complex128')
    dft = FourierTransform(discr)

    for factor in (2, 1j, -2.5j, 1 - 4j):
        func_true_ft = factor * dft.range.element(char_interval_ft)
        func_dft = dft(factor * discr.element(char_interval))
        assert (func_dft - func_true_ft).norm() < 1e-6 
Example #9
Source File: beamformer.py    From setk with Apache License 2.0 6 votes vote down vote up
def diffuse_covar(num_bins, dist_mat, sr=16000, c=340, diag_eps=0.1):
    """
    Compute covarance matrices of the spherically isotropic noise field
        Gamma(omega)_{ij} = sinc(omega * tau_{ij}) = sinc(2 * pi * f * tau_{ij})
    Arguments:
        num_bins: number of the FFT points
        dist_mat: distance matrix
        sr: sample rate
        c: sound of the speed
    Return:
        covar: covariance matrix, F x N x N
    """
    N, _ = dist_mat.shape
    eps = np.eye(N) * diag_eps
    covar = np.zeros([num_bins, N, N])
    for f in range(num_bins):
        omega = np.pi * f * sr / (num_bins - 1)
        covar[f] = np.sinc(dist_mat * omega / c) + eps
    return covar 
Example #10
Source File: quaternion.py    From RelativePose with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def expmap_to_quaternion(e):
    """
    Convert axis-angle rotations (aka exponential maps) to quaternions.
    Stable formula from "Practical Parameterization of Rotations Using the Exponential Map".
    Expects a tensor of shape (*, 3), where * denotes any number of dimensions.
    Returns a tensor of shape (*, 4).
    """
    assert e.shape[-1] == 3
    
    original_shape = list(e.shape)
    original_shape[-1] = 4
    e = e.reshape(-1, 3)

    theta = np.linalg.norm(e, axis=1).reshape(-1, 1)
    w = np.cos(0.5*theta).reshape(-1, 1)
    xyz = 0.5*np.sinc(0.5*theta/np.pi)*e
    return np.concatenate((w, xyz), axis=1).reshape(original_shape) 
Example #11
Source File: test_propagation.py    From hcipy with MIT License 6 votes vote down vote up
def test_fraunhofer_propagation_rectangular():
	for num_pix in [512, 1024]:
		pupil_grid = make_pupil_grid(num_pix)
		focal_grid = make_focal_grid(16, 8)

		for size in [[1,1], [0.75,1], [0.75,0.75]]:
			aperture = evaluate_supersampled(rectangular_aperture(size), pupil_grid, 8)

			for focal_length in [1, 1.3]:
				prop = FraunhoferPropagator(pupil_grid, focal_grid, focal_length=focal_length)

				for wavelength in [1, 0.8]:
					wf = Wavefront(aperture, wavelength)
					img = prop(wf).electric_field
					img /= img[np.argmax(np.abs(img))]

					k_x, k_y = np.array(size) / wavelength / focal_length
					reference = (np.sinc(k_x * focal_grid.x) * np.sinc(k_y * focal_grid.y))

					if num_pix == 512:
						assert np.abs(img - reference).max() < 5e-5
					elif num_pix == 1024:
						assert np.abs(img - reference).max() < 2e-5
					else:
						assert False # This should never happen. 
Example #12
Source File: powder.py    From xrayutilities with GNU General Public License v2.0 6 votes vote down vote up
def conv_receiver_slit(self):
        """
        compute the rectangular convolution for the receiver slit or SiPSD
        pixel size

        Returns
        -------
        array-like
            the convolver
        """
        me = self.get_function_name()  # the name of the convolver, as a string
        # The receiver slit convolution is a top-hat of angular half-width
        # a=(slit_width/2)/diffractometer_radius
        # which has Fourier transform of sin(a omega)/(a omega)
        # NOTE! numpy's sinc(x) is sin(pi x)/(pi x), not sin(x)/x
        if self.param_dicts[me].get("slit_width", None) is None:
            return None

        epsr = (self.param_dicts["conv_receiver_slit"]["slit_width"] /
                self.param_dicts["conv_global"]["diffractometer_radius"])
        return self.general_tophat(me, epsr) 
Example #13
Source File: catalog.py    From nbodykit with GNU General Public License v3.0 6 votes vote down vote up
def CompensateTSC(w, v):
    """
    Return the Fourier-space kernel that accounts for the convolution of
    the gridded field with the TSC window function in configuration space.

    .. note::
        see equation 18 (with p=3) of
        `Jing et al 2005 <https://arxiv.org/abs/astro-ph/0409240>`_

    Parameters
    ----------
    w : list of arrays
        the list of "circular" coordinate arrays, ranging from
        :math:`[-\pi, \pi)`.
    v : array_like
        the field array
    """
    for i in range(3):
        wi = w[i]
        tmp = (numpy.sinc(0.5 * wi / numpy.pi) ) ** 3
        v = v / tmp
    return v 
Example #14
Source File: catalog.py    From nbodykit with GNU General Public License v3.0 6 votes vote down vote up
def CompensatePCS(w, v):
    """
    Return the Fourier-space kernel that accounts for the convolution of
    the gridded field with the PCS window function in configuration space.

    .. note::
        see equation 18 (with p=4) of
        `Jing et al 2005 <https://arxiv.org/abs/astro-ph/0409240>`_

    Parameters
    ----------
    w : list of arrays
        the list of "circular" coordinate arrays, ranging from
        :math:`[-\pi, \pi)`.
    v : array_like
        the field array
    """
    for i in range(3):
        wi = w[i]
        tmp = (numpy.sinc(0.5 * wi / numpy.pi) ) ** 4
        v = v / tmp
    return v 
Example #15
Source File: catalog.py    From nbodykit with GNU General Public License v3.0 6 votes vote down vote up
def CompensateCIC(w, v):
    """
    Return the Fourier-space kernel that accounts for the convolution of
    the gridded field with the CIC window function in configuration space

    .. note::
        see equation 18 (with p=2) of
        `Jing et al 2005 <https://arxiv.org/abs/astro-ph/0409240>`_

    Parameters
    ----------
    w : list of arrays
        the list of "circular" coordinate arrays, ranging from
        :math:`[-\pi, \pi)`.
    v : array_like
        the field array
    """
    for i in range(3):
        wi = w[i]
        tmp = (numpy.sinc(0.5 * wi / numpy.pi) ) ** 2
        tmp[wi == 0.] = 1.
        v = v / tmp
    return v 
Example #16
Source File: rotor_parameterisation.py    From clifford with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def val_exp(B_val):
    """
    Fast implementation of the translation and rotation specific exp function
    """
    t_val = imt_func(B_val, no.value)

    phiP_val = B_val - mult_with_ninf(t_val)
    phi = np.sqrt(-float(gmt_func(phiP_val, phiP_val)[0]))
    P_val = phiP_val / phi

    P_n_val = gmt_func(P_val, I3.value)
    t_nor_val = gmt_func(imt_func(t_val, P_n_val), P_n_val)
    t_par_val = t_val - t_nor_val

    coef_val = np.sin(phi) * P_val
    coef_val[0] += np.cos(phi)

    R_val = coef_val + gmt_func(coef_val, mult_with_ninf(t_nor_val)) + \
        np.sinc(phi/np.pi) * mult_with_ninf(t_par_val)
    return R_val 
Example #17
Source File: transforms.py    From spectral_connectivity with GNU General Public License v3.0 6 votes vote down vote up
def _get_taper_eigenvalues(tapers, half_bandwidth, time_index):
    '''Finds the eigenvalues of the original spectral concentration
    problem using the autocorr sequence technique from Percival and Walden,
    1993 pg 390

    Parameters
    ----------
    tapers : array, shape (n_tapers, n_time_samples_per_window)
    half_bandwidth : float
    time_index : array, (n_time_samples_per_window,)

    Returns
    -------
    eigenvalues : array, shape (n_tapers,)

    '''

    ideal_filter = 4 * half_bandwidth * np.sinc(
        2 * half_bandwidth * time_index)
    ideal_filter[0] = 2 * half_bandwidth
    n_time_samples_per_window = len(time_index)
    return np.dot(
        _auto_correlation(tapers)[:, :n_time_samples_per_window],
        ideal_filter) 
Example #18
Source File: utilities.py    From pyroomacoustics with MIT License 6 votes vote down vote up
def fractional_delay(t0):
    """
    Creates a fractional delay filter using a windowed sinc function.
    The length of the filter is fixed by the module wide constant
    `frac_delay_length` (default 81).

    Parameters
    ----------
    t0: float
        The delay in fraction of sample. Typically between -1 and 1.

    Returns
    -------
    numpy array
        A fractional delay filter with specified delay.
    """

    N = constants.get("frac_delay_length")

    return np.hanning(N) * np.sinc(np.arange(N) - (N - 1) / 2 - t0) 
Example #19
Source File: prony.py    From parametric_modeling with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
def main():
    """Test driver"""
    # From pp. 149-150
    x = np.ones(21)
    p = q = 1
    print('x: {}\np: {}\nq: {}'.format(x,p,q))
    b,a,err = prony(x, p, q)
    print('a: {}\nb: {}\nerr: {}'.format(a,b,err))
 
    # From pp. 152-153
    # Note that these results don't match the book, but they do match the
    # MATLAB version. So I'm either setting things up wrong or this is an
    # errata in the book.
    p = q = 5
    nd = 5
    n = np.arange(11)
    i = np.sinc((n-nd)/2)/2
    b,a,err = prony(i, p, q)
    print('a: {}\nb: {}\nerr: {}'.format(a,b,err)) 
Example #20
Source File: Filter.py    From urh with GNU General Public License v3.0 6 votes vote down vote up
def design_windowed_sinc_lpf(fc, bw):
        N = Filter.get_filter_length_from_bandwidth(bw)

        # Compute sinc filter impulse response
        h = np.sinc(2 * fc * (np.arange(N) - (N - 1) / 2.))

        # We use blackman window function
        w = np.blackman(N)

        # Multiply sinc filter with window function
        h = h * w

        # Normalize to get unity gain
        h_unity = h / np.sum(h)

        return h_unity 
Example #21
Source File: fourier_test.py    From odl with Mozilla Public License 2.0 5 votes vote down vote up
def sinc(x):
    # numpy.sinc scales by pi, we don't want that
    return np.sinc(x / np.pi)


# ---- DiscreteFourierTransform ---- # 
Example #22
Source File: test_limits.py    From numdifftools with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def test_sinx_div_x(self):

        def fun(x):
            return np.sin(x) / x

        for path in ['radial', 'spiral']:
            lim_f = Limit(fun, path=path, full_output=True)

            x = np.arange(-10, 10) / np.pi
            lim_f0, err = lim_f(x * np.pi)
            assert_array_almost_equal(lim_f0, np.sinc(x))
            assert np.all(err.error_estimate < 1.0e-14) 
Example #23
Source File: ft_utils.py    From odl with Mozilla Public License 2.0 5 votes vote down vote up
def _interp_kernel_ft(norm_freqs, interp):
    """Scaled FT of a one-dimensional interpolation kernel.

    For normalized frequencies ``-1/2 <= xi <= 1/2``, this
    function returns::

        sinc(pi * xi)**k / sqrt(2 * pi)

    where ``k=1`` for 'nearest' and ``k=2`` for 'linear' interpolation.

    Parameters
    ----------
    norm_freqs : `numpy.ndarray`
        Normalized frequencies between -1/2 and 1/2
    interp : {'nearest', 'linear'}
        Type of interpolation kernel

    Returns
    -------
    ker_ft : `numpy.ndarray`
        Values of the kernel FT at the given frequencies
    """
    # Numpy's sinc(x) is equal to the 'math' sinc(pi * x)
    ker_ft = np.sinc(norm_freqs)
    interp_ = str(interp).lower()
    if interp_ == 'nearest':
        pass
    elif interp_ == 'linear':
        ker_ft *= ker_ft
    else:
        raise ValueError("`interp` '{}' not understood".format(interp))

    ker_ft /= np.sqrt(2 * np.pi)
    return ker_ft 
Example #24
Source File: test_quantity_non_ufuncs.py    From Carnets with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def test_sinc(self):
        q = [0., 3690., -270., 690.] * u.deg
        out = np.sinc(q)
        expected = np.sinc(q.to_value(u.radian)) * u.one
        assert isinstance(out, u.Quantity)
        assert np.all(out == expected)
        with pytest.raises(u.UnitsError):
            np.sinc(1.*u.one) 
Example #25
Source File: utils.py    From pystoi with MIT License 5 votes vote down vote up
def _resample_window_oct(p, q):
    """Port of Octave code to Python"""

    gcd = np.gcd(p, q)
    if gcd > 1:
        p /= gcd
        q /= gcd

    # Properties of the antialiasing filter
    log10_rejection = -3.0
    stopband_cutoff_f = 1. / (2 * max(p, q))
    roll_off_width = stopband_cutoff_f / 10

    # Determine filter length
    rejection_dB = -20 * log10_rejection
    L = np.ceil((rejection_dB - 8) / (28.714 * roll_off_width))

    # Ideal sinc filter
    t = np.arange(-L, L + 1)
    ideal_filter = 2 * p * stopband_cutoff_f \
        * np.sinc(2 * stopband_cutoff_f * t)

    # Determine parameter of Kaiser window
    if (rejection_dB >= 21) and (rejection_dB <= 50):
        beta = 0.5842 * (rejection_dB - 21)**0.4 \
            + 0.07886 * (rejection_dB - 21)
    elif rejection_dB > 50:
        beta = 0.1102 * (rejection_dB - 8.7)
    else:
        beta = 0.0

    # Apodize ideal filter response
    h = np.kaiser(2 * L + 1, beta) * ideal_filter

    return h 
Example #26
Source File: interpolation.py    From scarlet with MIT License 5 votes vote down vote up
def sinc2D(y, x):
    """
    2-D sinc function based on the product of 2 1-D sincs

    Parameters
    ----------
        x, y: arrays
            Coordinates where to evaluate the 2-D sinc
    Returns
    -------
    result: array
        2-D sinc evaluated in x and y
    """
    return np.dot(np.sinc(y), np.sinc(x)) 
Example #27
Source File: waveform.py    From QuLab with MIT License 5 votes vote down vote up
def sinc(bw):
    width = 100 / bw
    return Waveform(bounds=(-0.5 * width, 0.5 * width, +np.inf),
                    seq=(_zero, _basic_wave(SINC, bw), _zero)) 
Example #28
Source File: utilities.py    From pyroomacoustics with MIT License 5 votes vote down vote up
def low_pass_dirac(t0, alpha, Fs, N):
    """
    Creates a vector containing a lowpass Dirac of duration T sampled at Fs
    with delay t0 and attenuation alpha.

    If t0 and alpha are 2D column vectors of the same size, then the function
    returns a matrix with each line corresponding to pair of t0/alpha values.
    """

    return alpha * np.sinc(np.arange(N) - Fs * t0) 
Example #29
Source File: Spectrum.py    From PyRAT with Mozilla Public License 2.0 5 votes vote down vote up
def spec_correction(array, alpha=1.00, fix=True, cutoff=0.05, func=(None, None)):
        if array.ndim == 1:
            array = array[np.newaxis, ...]
        corr = np.empty_like(array)
        cent = [None] * len(array)
        for k, arr in enumerate(array):
            spec = arr / np.mean(arr)
            spec = filters.uniform_filter(spec, len(array) // 16, mode='wrap')
            offset = np.argmax(spec)
            spec = np.roll(spec, -offset)
            up_bound = 0
            while up_bound < len(spec) and spec[up_bound] > cutoff * spec[0]:
                up_bound += 1
            lo_bound = -1
            while lo_bound > -len(spec) and spec[lo_bound] > cutoff * spec[0]:
                lo_bound -= 1

            if func[0] == "Hamming" or func[0] is None:
                w = func[1] - (1.0 - func[1]) * np.cos(
                    2 * np.pi * np.arange(up_bound - lo_bound) / (up_bound - lo_bound - 1))
            elif func[0] == "Lanczos":
                w = np.sinc(2 * np.arange(up_bound - lo_bound) / (up_bound - lo_bound - 1) - 1)

            if fix is True:
                corr[k, ...] = 1.0 / spec
            else:
                corr[k, ...] = 1.0

            corr[k, up_bound - 1:lo_bound + 1] = 0.0

            corr[k, lo_bound:] *= w[0:-lo_bound]
            corr[k, :up_bound] *= w[-up_bound:]
            corr[k, ...] = np.roll(corr[k, ...], offset)
            corr[k, ...] /= np.mean(corr[k, ...])

            if offset > len(spec) // 2:
                offset -= len(spec)
            cent[k] = (up_bound + lo_bound) // 2 + offset
        return np.squeeze(corr), cent 
Example #30
Source File: rotor_parameterisation.py    From clifford with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def extractRotorComponents(R):
    """
    Extracts the translation and rotation information from a TR rotor

    From :cite:`wareham-interpolation`.
    """
    phi = np.arccos(R[()])             # scalar
    phi2 = phi * phi                  # scalar
    # Notice: np.sinc(pi * x)/(pi x)
    phi_sinc = np.sinc(phi/np.pi)             # scalar
    phiP = ((R(2)*ninf)|ep)/(phi_sinc)
    t_normal_n = -((phiP * R(4))/(phi2 * phi_sinc))
    t_perpendicular_n = -(phiP * (phiP * R(2))(2))/(phi2 * phi_sinc)
    return phiP, t_normal_n, t_perpendicular_n