#!/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