def solve_for_nearest( px,py,rx,ry ):
  dpx = polyder(px)
  dpy = polyder(py)
  cp = polymul( dpx, px ) + polymul( dpy, py )
  cp = polyadd( cp, -rx*dpx )      
  cp = polyadd( cp, -ry*dpy )
  t = roots(cp)
  t = real(t[isreal(t)])
  t = t[ (t>=0) * (t<=1) ]

  ##tt = linspace(0,1,100)
  ##from pylab import plot
  ##plot( polyval(px,tt), polyval(py,tt), 'k', hold = 0 )
  ##plot( [rx],[ry], 'r.' )
  ##plot( polyval(px,t[isreal(t)*(real(t)>=0)*(real(t)<=1)]),
  ##      polyval(py,t[isreal(t)*(real(t)>=0)*(real(t)<=1)]), 'o' )
  ##pdb.set_trace()

  if len(t):
    if len(t) == 1:
      return t[0]
    else:
      ux = polyval( px, t )
      uy = polyval( py, t )
      d = hypot( ux - rx, uy - ry )
      return t[ d==d.min() ][0]
  else:
    t = array([0.0,1.0])
    ux = polyval( px, t )
    uy = polyval( py, t )
    d = hypot( ux - rx, uy - ry )
    if d[0] < d[1]:
      return 0.0
    else:
      return 1.0

    def concat(self, other):
        """Create new reception parameters by combining this set of reception
        parameters with another.
        """
        # Combine the filters
        if self.filter is None:
            new_filter = other.filter
        elif other.filter is None:
            new_filter = self.filter
        elif (isinstance(self.filter, LinearFilter) and
                isinstance(other.filter, LinearFilter)):
            # Combine linear filters by multiplying their numerators and
            # denominators.
            new_filter = LinearFilter(
                np.polymul(self.filter.num, other.filter.num),
                np.polymul(self.filter.den, other.filter.den)
            )
        else:
            raise NotImplementedError(
                "Cannot combine filters of type {} and {}".format(
                    type(self.filter), type(other.filter)))

        # Combine the learning rules
        if self.learning_rule is not None and other.learning_rule is not None:
            raise NotImplementedError(
                "Cannot combine learning rules {} and {}".format(
                    self.learning_rule, other.learning_rule))

        new_learning_rule = self.learning_rule or other.learning_rule

        # Create the new reception parameters
        return ReceptionParameters(new_filter, other.width, new_learning_rule)

    def __truediv__(self, other):
        """Divide two LTI objects."""

        if isinstance(other, (int, float, complex)):
            other = _convertToTransferFunction(
                other, inputs=self.inputs,
                outputs=self.inputs)
        else:
            other = _convertToTransferFunction(other)

        if (self.inputs > 1 or self.outputs > 1 or
                other.inputs > 1 or other.outputs > 1):
            raise NotImplementedError(
                "TransferFunction.__truediv__ is currently \
                implemented only for SISO systems.")

        # Figure out the sampling time to use
        if (self.dt is None and other.dt is not None):
            dt = other.dt       # use dt from second argument
        elif (other.dt is None and self.dt is not None)\
                or (self.dt == other.dt):
            dt = self.dt        # use dt from first argument
        else:
            raise ValueError("Systems have different sampling times")

        num = polymul(self.num[0][0], other.den[0][0])
        den = polymul(self.den[0][0], other.num[0][0])

        return TransferFunction(num, den, dt)

 def decascade(self, ncascade=1):
     '''
     Reduces cascades of low order filters into smaller cascades of high order filters.
     
     ``ncascade`` is the number of cascaded filters to use, which should be
     a divisor of the original number.
     
     Note that higher order filters are often numerically unstable.
     '''
     n, m, p = self.filt_b.shape
     if p%ncascade!=0:
         raise ValueError('Number of cascades must be a divisor of original number of cascaded filters.')
     b = np.zeros((n, (m-1)*(p/ncascade)+1, ncascade))
     a = np.zeros((n, (m-1)*(p/ncascade)+1, ncascade))
     for i in range(n):
         for k in range(ncascade):
             bp = np.ones(1)
             ap = np.ones(1)
             for j in range(k*(p/ncascade), (k+1)*(p/ncascade)):
                 bp = np.polymul(bp, self.filt_b[i, ::-1, j])
                 ap = np.polymul(ap, self.filt_a[i, ::-1, j])
             bp = bp[::-1]
             ap = ap[::-1]
             a0 = ap[0]
             ap /= a0
             bp /= a0
             b[i, :len(bp), k] = bp
             a[i, :len(ap), k] = ap
     self.filt_b = np.array(b, order='F')
     self.filt_a = np.array(a, order='F')
     self.filt_state = np.zeros((b.shape[0], b.shape[1], b.shape[2]), order='F')

    def feedback(self, other=1, sign=-1):
        """Feedback interconnection between two LTI objects."""
        other = _convertToTransferFunction(other)

        if (self.inputs > 1 or self.outputs > 1 or
                other.inputs > 1 or other.outputs > 1):
            # TODO: MIMO feedback
            raise NotImplementedError("TransferFunction.feedback is currently \
only implemented for SISO functions.")

        # Figure out the sampling time to use
        if (self.dt is None and other.dt is not None):
            dt = other.dt       # use dt from second argument
        elif (other.dt is None and self.dt is not None) \
                or (self.dt == other.dt):
            dt = self.dt        # use dt from first argument
        else:
            raise ValueError("Systems have different sampling times")

        num1 = self.num[0][0]
        den1 = self.den[0][0]
        num2 = other.num[0][0]
        den2 = other.den[0][0]

        num = polymul(num1, den2)
        den = polyadd(polymul(den2, den1), -sign * polymul(num2, num1))

        return TransferFunction(num, den, dt)

def _combine_reception_params(in_reception_parameters,
                              out_reception_parameters):
    """Combine reception parameters to join two signals into one, e.g., for
    optimising out a passthrough Node.
    """
    # Construct the new reception parameters
    # Combine the filters
    filter_in = in_reception_parameters.filter
    filter_out = out_reception_parameters.filter

    if (filter_in is None or filter_out is None):
        # If either filter is None then just use the filter from the other
        # connection
        new_filter = filter_in or filter_out
    elif (isinstance(filter_in, nengo.LinearFilter) and
            isinstance(filter_out, nengo.LinearFilter)):
        # Both filters are linear filters, so multiply the numerators and
        # denominators together to get a new linear filter.
        new_num = np.polymul(filter_in.num, filter_out.num)
        new_den = np.polymul(filter_in.den, filter_out.den)

        new_filter = nengo.LinearFilter(new_num, new_den)
    else:
        raise NotImplementedError

    # Take the size in from the second reception parameter, construct the new
    # reception parameters.
    return ReceptionParameters(new_filter, out_reception_parameters.width)

def A_weighting(fs):
    """Design of an A-weighting filter.
    b, a = A_weighting(fs) designs a digital A-weighting filter for
    sampling frequency `fs`. Usage: y = scipy.signal.lfilter(b, a, x).
    Warning: `fs` should normally be higher than 20 kHz. For example,
    fs = 48000 yields a class 1-compliant filter.
    References:
       [1] IEC/CD 1672: Electroacoustics-Sound Level Meters, Nov. 1996.
    """
    # Definition of analog A-weighting filter according to IEC/CD 1672.
    f1 = 20.598997
    f2 = 107.65265
    f3 = 737.86223
    f4 = 12194.217
    A1000 = 1.9997

    NUMs = [(2*pi * f4)**2 * (10**(A1000/20)), 0, 0, 0, 0]
    DENs = polymul([1, 4*pi * f4, (2*pi * f4)**2],
                   [1, 4*pi * f1, (2*pi * f1)**2])
    DENs = polymul(polymul(DENs, [1, 2*pi * f3]),
                                 [1, 2*pi * f2])

    # Use the bilinear transformation to get the digital filter.
    # (Octave, MATLAB, and PyLab disagree about Fs vs 1/Fs)
    return bilinear(NUMs, DENs, fs)

def Closed_loop(Kz, Kp, Gz, Gp):
    """
    Return zero and pole polynomial for a closed loop function.

    Parameters
    ----------
    Kz & Gz : list
        Polynomial constants in the numerator.
    Kz & Gz : list
        Polynomial constants in the denominator.

    Returns
    -------
    Zeros_poly : list
        List of zero polynomial for closed loop function.

    Poles_poly : list
        List of pole polynomial for closed loop function.

   """

    # calculating the product of the two polynomials in the numerator
    # and denominator of transfer function GK
    Z_GK = numpy.polymul(Kz, Gz)
    P_GK = numpy.polymul(Kp, Gp)

    # calculating the polynomial of closed loop
    # sensitivity function s = 1/(1+GK)
    Zeros_poly = Z_GK
    Poles_poly = numpy.polyadd(Z_GK, P_GK)
    return Zeros_poly, Poles_poly

def compute_join_length( px, py, tlow = 0.0, thigh = 1.0 ):
  from scipy.integrate import quad
  xp  = polyder( px, 1 )
  yp  = polyder( py, 1 )
  xp2 = polymul( xp, xp )
  yp2 = polymul( yp, yp )
  p   = polyadd( xp2, yp2 )
  integrand = lambda t: sqrt( polyval( p, t ) )
  return quad(integrand, tlow, thigh) [0]

def make_C_weighting(fs):

	den = p1
	den = np.polymul(den, p1)
	den = np.polymul(den, p4)
	den = np.polymul(den, p4)
	num = np.poly1d([(2.*np.pi*f4)**2, 0, 0])

	B, A = signal.bilinear(num, den, fs)
	return B, A

def compute_join_curvature( px, py ):
  from scipy.integrate import quad
  xp  = polyder( px, 1 )
  xpp = polyder( px, 2 )
  yp  = polyder( py, 1 )
  ypp = polyder( py, 2 )
  pn = polyadd( polymul( xp, ypp ), polymul( yp, xpp )) #numerator
  pd = polyadd( polymul( xp, xp ) , polymul( yp, yp ) ) #denominator
  integrand = lambda t:  fabs(polyval( pn, t )/( polyval( pd, t )**(1.5)) ) 
  return quad(integrand, 0, 1) [0]

def residue(b,a,tol=1e-3,rtype='avg'):
    """Compute partial-fraction expansion of b(s) / a(s).

    If M = len(b) and N = len(a)

            b(s)     b[0] s**(M-1) + b[1] s**(M-2) + ... + b[M-1]
    H(s) = ------ = ----------------------------------------------
            a(s)     a[0] s**(N-1) + a[1] s**(N-2) + ... + a[N-1]

             r[0]       r[1]             r[-1]
         = -------- + -------- + ... + --------- + k(s)
           (s-p[0])   (s-p[1])         (s-p[-1])

    If there are any repeated roots (closer than tol), then the partial
    fraction expansion has terms like

            r[i]      r[i+1]              r[i+n-1]
          -------- + ----------- + ... + -----------
          (s-p[i])  (s-p[i])**2          (s-p[i])**n

    See also:  invres, poly, polyval, unique_roots
    """

    b,a = map(asarray,(b,a))
    k,b = polydiv(b,a)
    p = roots(a)
    r = p*0.0
    pout, mult = unique_roots(p,tol=tol,rtype=rtype)
    p = []
    for n in range(len(pout)):
        p.extend([pout[n]]*mult[n])
    p = asarray(p)
    # Compute the residue from the general formula
    indx = 0
    for n in range(len(pout)):
        bn = b.copy()
        pn = []
        for l in range(len(pout)):
            if l != n:
                pn.extend([pout[l]]*mult[l])
        an = atleast_1d(poly(pn))
        # bn(s) / an(s) is (s-po[n])**Nn * b(s) / a(s) where Nn is
        # multiplicity of pole at po[n]
        sig = mult[n]
        for m in range(sig,0,-1):
            if sig > m:
                # compute next derivative of bn(s) / an(s)
                term1 = polymul(polyder(bn,1),an)
                term2 = polymul(bn,polyder(an,1))
                bn = polysub(term1,term2)
                an = polymul(an,an)
            r[indx+m-1] = polyval(bn,pout[n]) / polyval(an,pout[n]) \
                          / factorial(sig-m)
        indx += sig
    return r, p, k

def mls_polynomial_coefficients(rho, degree):
    """Determine the coefficients for a MLS polynomial smoother

    Parameters
    ----------
    rho : {float}
        Spectral radius of the matrix in question
    degree : {int}
        Degree of polynomial coefficients to generate

    Returns
    -------
    Tuple of arrays (coeffs,roots) containing the
    coefficients for the (symmetric) polynomial smoother and
    the roots of polynomial prolongation smoother.

    The coefficients of the polynomial are in descending order

    References
    ----------
    .. [1] Parallel multigrid smoothing: polynomial versus Gauss--Seidel
       M. F. Adams, M. Brezina, J. J. Hu, and R. S. Tuminaro
       J. Comp. Phys., 188 (2003), pp. 593--610

    Examples
    --------
    >>> from pyamg.relaxation.chebyshev import mls_polynomial_coefficients
    >>> mls = mls_polynomial_coefficients(2.0, 2)
    >>> print mls[0] # coefficients
    [   6.4  -48.   144.  -220.   180.   -75.8   14.5]
    >>> print mls[1] # roots
    [ 1.4472136  0.5527864]
    """

    # std_roots = np.cos(np.pi * (np.arange(degree) + 0.5)/ degree)
    # print std_roots

    roots = rho/2.0 * (1.0 - np.cos(2*np.pi*(np.arange(degree,
                                       dtype='float64') + 1)/(2.0*degree+1.0)))
    # print roots
    roots = 1.0/roots

    # S_coeffs = list(-np.poly(roots)[1:][::-1])

    S = np.poly(roots)[::-1]  # monomial coefficients of S error propagator

    SSA_max = rho/((2.0*degree+1.0)**2)  # upper bound spectral radius of S^2A
    S_hat = np.polymul(S, S)  # monomial coefficients of \hat{S} propagator
    S_hat = np.hstack(((-1.0/SSA_max)*S_hat, [1]))

    # coeff for combined error propagator \hat{S}S
    coeffs = np.polymul(S_hat, S)
    coeffs = -coeffs[:-1]             # coeff for smoother

    return (coeffs, roots)

def _addSISO(num1, den1, num2, den2):
    """Return num/den = num1/den1 + num2/den2.

    Each numerator and denominator is a list of polynomial coefficients.

    """

    num = polyadd(polymul(num1, den2), polymul(num2, den1))
    den = polymul(den1, den2)

    return num, den

def Rule_7(w_start, w_end):
    # Rule 7 determining the phase of GGm at -180deg
    # this is solved visually from a plot

    w = np.logspace(w_start, w_end, 1000)

    Pz = np.polymul(G()[0], Gm()[0])
    Pp = np.polymul(G()[1], Gm()[1])
    [w, h] = scs.freqs(Pz, Pp, w)

    plt.semilogx(w, (180 / np.pi) * (phase(h) + w * Time_Delay()))
    plt.show()

def Closed_loop(Kz, Kp, Gz, Gp):
    """Kz & Gz is the polynomial constants in the numerator
    Kp & Gp is the polynomial constants in the denominator"""

    # calculating the product of the two polynomials in the numerator and denominator of transfer function GK
    Z_GK = numpy.polymul(Kz, Gz)
    P_GK = numpy.polymul(Kp, Gp)

    #calculating the polynomial of closed loop sensitivity function s = 1/(1+GK)
    Zeros_poly = Z_GK
    Poles_poly = numpy.polyadd(Z_GK, P_GK)
    return Zeros_poly, Poles_poly

def series(sys1, sys2):


    """Series connection of two systems.
    """
    if not isinstance(sys1, signal.lti):
        sys1 = signal.lti(*sys1)
    if not isinstance(sys2, signal.lti):
        sys2 = signal.lti(*sys2)
    num = np.polymul(sys1.num, sys2.num)
    den = np.polymul(sys1.den, sys2.den)
    sys = signal.lti(num, den)
    return sys

 def __init__(self, c=4):
     c_ops = self._c_ops
     if not c in c_ops:
         raise ValueError("c should be one of " + str(c_ops))
     index = c / 2 - 1
     p_1 = (1,)
     for _i in xrange(self._coefs_1_orders[index]):
         p_1 = np.polymul(p_1, self._coefs_1)
     p_2 = self._coefs_2s[index]
     p = np.polymul(p_1, p_2)
     self.poly = np.poly1d(p)
     self._diff_order = -1
     self.set_diff_order(0)

def Closed_loop (Kz,Kp,Gz,Gp):
    """ Kz and Gz are the polynomial constants in the numerator 
    Kp and Gp are the polynomial constants in the denominator""" 
    
    
    
    # calculating the product of the two polynomials in the numerator and denominator of transfer function GK
    Z_GK         =np.polymul(Kz,Gz)
    P_GK         =np.polymul(Kp,Gp)    
    # calculating the polynomial of closed loop function T=(GK/1+GK)
    Zeros_poly   =Z_GK
    Poles_poly   =np.polyadd(Z_GK,P_GK)   
    return      Zeros_poly,Poles_poly

def A_weighting(fs):
    """
    Design of an A-weighting filter.

    Designs a digital A-weighting filter for
    sampling frequency `fs`. Usage: y = lfilter(b, a, x).
    Warning: fs should normally be higher than 20 kHz. For example,
    fs = 48000 yields a class 1-compliant filter.

    fs : float
        Sampling frequency

    Since this uses the bilinear transform, frequency response around fs/2 will
    be inaccurate at lower sampling rates.  A-weighting is undefined above
    20 kHz, though.

    Example:

    from scipy.signal import freqz
    import matplotlib.pyplot as plt
    fs = 200000  # change to 48000 to see truncation
    b, a = A_weighting(fs)
    f = np.logspace(np.log10(10), np.log10(fs/2), 1000)
    w = 2*pi * f / fs
    w, h = freqz(b, a, w)
    plt.semilogx(w*fs/(2*pi), 20*np.log10(abs(h)))
    plt.grid(True, color='0.7', linestyle='-', which='both', axis='both')
    plt.axis([10, 100e3, -50, 20])

    References:
       [1] IEC/CD 1672: Electroacoustics-Sound Level Meters, Nov. 1996.
    """

    # Definition of analog A-weighting filter according to IEC/CD 1672.
    f1 = 20.598997
    f2 = 107.65265
    f3 = 737.86223
    f4 = 12194.217
    A1000 = 1.9997

    NUMs = [(2*pi * f4)**2 * (10**(A1000/20)), 0, 0, 0, 0]

    DENs = polymul([1, 4*pi * f4, (2*pi * f4)**2],
                   [1, 4*pi * f1, (2*pi * f1)**2])
    DENs = polymul(DENs, [1, 2*pi * f3])
    DENs = polymul(DENs, [1, 2*pi * f2])

    # Analog confirmed to match https://en.wikipedia.org/wiki/A-weighting#A_2

    # Use the bilinear transformation to get the digital filter.
    return bilinear(NUMs, DENs, fs)