#!/usr/bin/python
# 13.1

import os, sys
sys.path += [os.getcwdu() + os.sep + ".." + os.sep + "python"]

import scipy as sp
from polyfuncs import polmul, poladd

def spec(A, dA, B, dB, rho):
    AA = rho * sp.convolve(A, A[::-1])
    BB = sp.convolve(B, B[::-1])
    diff1 = dA - dB
    dBB = 2 * dB
    for i in range(diff1):
        BB, dBB = polmul(BB, dBB, sp.array([0, 1]), 1)
    rbbr, drbbr = poladd(AA, 2*dA, BB, dBB)
    rbbr = rbbr.squeeze()
    rts = sp.roots(rbbr)
    rtsin = rts[dA:2*dA+1]
    b = 1
    for i in range(dA):
        b = sp.convolve(b, sp.array([1, -rtsin[i]]))
    bbr = sp.convolve(b, b[::-1])
    r = rbbr[0] / bbr[0]
    return r,b,rbbr

if __name__ == "__main__":
    A = sp.convolve([-0.5, 1], [-0.9, 1])
    dA = 2
    B = 0.5*sp.array([-0.9, 1])
    dB = 1
    rho = 1

    r, beta1, sigma = spec(A, dA, B, dB, rho)
    print "sigma", sigma
    print "beta1", beta1
    print "r", r