

In [1]:

    
from scipy.sparse.linalg import eigs
import scipy.sparse as sparse
import numpy as np
import numpy.linalg as la
import sympy as sp

import matplotlib.pyplot as plt
%matplotlib inline

import Potapov_Code.Time_Delay_Network as networks
from Potapov_Code.functions import spatial_modes
from decimal import Decimal

from sympy import init_printing
init_printing() 

from sympy.utilities.autowrap import ufuncify
import time



In [2]:

    
import Potapov_Code.Time_Delay_Network



In [3]:

    
from Potapov_Code.functions import gcd_lst



In [4]:

    
X1 = Potapov_Code.Time_Delay_Network.Example3()



In [6]:

    
def _find_commensurate(delays):
    '''
    Find the 'gcd' but for Decimal numbers.

    Args:
        delays(list of Demicals): numbers whose gcd will be found.

    Returns:
        Decimal gcd.
    '''
    mult = min([d.as_tuple().exponent for d in delays])
    power = 10**-mult
    delays = map(lambda x: x*power,delays)
    int_gcd = gcd_lst(delays)
    return int_gcd/power

Decimal_delays = map(lambda x: Decimal(str(x)),X1.delays)
Decimal_gcd = _find_commensurate(Decimal_delays)

xs = [sp.symbols('x_'+str(i)) for i in range(4) ]
E_sym = sp.Matrix(np.zeros_like(X1.M1))
for i,delay in enumerate(Decimal_delays):
    E_sym[i,i] = xs[i]
M1_sym = sp.Matrix(X1.M1)
num, den = (E_sym - M1_sym).det().as_numer_denom()

Alternatively, we can Taylor expand about the $z$ variable only. This gives just Newton's method. To first order we get

$$ \Delta z = - \frac{f(z,T+\Delta T)}{\frac{\partial f}{\partial z (z,T+\Delta T}}. $$

This can be used as an iterative procedure with the substitution $z \mapsto z + \Delta_z$. In addition, dependence on $\Delta T $ can be incorporated by evaluation $\Delta T \equiv \Delta T(z + \Delta z)$.



In [7]:

    
num, den









    Out[7]:





$$\left ( 1.0 x_{0}^{2} x_{1}^{2} x_{2} x_{3} + 0.32 x_{0}^{2} x_{1}^{2} - 0.72 x_{0} x_{1} x_{2} x_{3} - 0.8352 x_{0} x_{1} + 0.1296 x_{2} x_{3} + 0.2592, \quad 1.0 x_{0} x_{1} - 0.36\right )$$



In [54]:

    
z, z_Delta = sp.symbols('z dz',complex=True)
Ts = [sp.symbols('T_'+str(i),real=True) for i in range(4)]
x,y = sp.symbols('x y', real = True)



In [16]:

    
num2 = num.subs({x: sp.exp(-z*T) for x,T in zip(xs,Ts)})



In [48]:

    
num2









    Out[48]:





$$0.2592 + 0.1296 e^{- T_{2} z} e^{- T_{3} z} - 0.8352 e^{- T_{0} z} e^{- T_{1} z} - 0.72 e^{- T_{0} z} e^{- T_{1} z} e^{- T_{2} z} e^{- T_{3} z} + 0.32 e^{- 2 T_{0} z} e^{- 2 T_{1} z} + 1.0 e^{- 2 T_{0} z} e^{- 2 T_{1} z} e^{- T_{2} z} e^{- T_{3} z}$$



In [17]:

    
num3 = num2.subs({sp.exp(-z*T): sp.cos(-x*T)*(1j*sp.sin(y*T)+sp.cos(y*T)) for T in Ts})



In [49]:

    
num3









    Out[49]:





$$1.0 \left(1.0 i \sin{\left (T_{0} y \right )} + \cos{\left (T_{0} y \right )}\right)^{2} \left(1.0 i \sin{\left (T_{1} y \right )} + \cos{\left (T_{1} y \right )}\right)^{2} \left(1.0 i \sin{\left (T_{2} y \right )} + \cos{\left (T_{2} y \right )}\right) \left(1.0 i \sin{\left (T_{3} y \right )} + \cos{\left (T_{3} y \right )}\right) \cos^{2}{\left (T_{0} x \right )} \cos^{2}{\left (T_{1} x \right )} \cos{\left (T_{2} x \right )} \cos{\left (T_{3} x \right )} + 0.32 \left(1.0 i \sin{\left (T_{0} y \right )} + \cos{\left (T_{0} y \right )}\right)^{2} \left(1.0 i \sin{\left (T_{1} y \right )} + \cos{\left (T_{1} y \right )}\right)^{2} \cos^{2}{\left (T_{0} x \right )} \cos^{2}{\left (T_{1} x \right )} - 0.72 \left(1.0 i \sin{\left (T_{0} y \right )} + \cos{\left (T_{0} y \right )}\right) \left(1.0 i \sin{\left (T_{1} y \right )} + \cos{\left (T_{1} y \right )}\right) \left(1.0 i \sin{\left (T_{2} y \right )} + \cos{\left (T_{2} y \right )}\right) \left(1.0 i \sin{\left (T_{3} y \right )} + \cos{\left (T_{3} y \right )}\right) \cos{\left (T_{0} x \right )} \cos{\left (T_{1} x \right )} \cos{\left (T_{2} x \right )} \cos{\left (T_{3} x \right )} - 0.8352 \left(1.0 i \sin{\left (T_{0} y \right )} + \cos{\left (T_{0} y \right )}\right) \left(1.0 i \sin{\left (T_{1} y \right )} + \cos{\left (T_{1} y \right )}\right) \cos{\left (T_{0} x \right )} \cos{\left (T_{1} x \right )} + 0.1296 \left(1.0 i \sin{\left (T_{2} y \right )} + \cos{\left (T_{2} y \right )}\right) \left(1.0 i \sin{\left (T_{3} y \right )} + \cos{\left (T_{3} y \right )}\right) \cos{\left (T_{2} x \right )} \cos{\left (T_{3} x \right )} + 0.2592$$



In [18]:

    
num_real,num_imag = num3.expand().as_real_imag()



In [57]:



In [44]:

    
[x,y]+[Ts]









    Out[44]:





$$\left [ x, \quad y, \quad \left [ T_{0}, \quad T_{1}, \quad T_{2}, \quad T_{3}\right ]\right ]$$



In [122]:

    
f_r =  ufuncify( [x,y]+Ts, num_real)
f_i =  ufuncify( [x,y]+Ts, num_imag)



In [47]:

    
f_r(1,2,3,4,5,6)+f_i(1,2,3,4,5,6)*1j









    Out[47]:





(0.14730830564958425-0.40346430265728966j)



In [52]:

    
## testing
num3.subs({x:1,y:2,Ts[0]:3,Ts[1]:4,Ts[2]:5,Ts[3]:6}).evalf()









    Out[52]:





$$0.147308305649584 - 0.40346430265729 i$$



In [89]:

    
#Consolidating
def get_real_imag_func(expression):
    '''
    Takes an expression in terms of Ts and sp.exp(-z*T) for T in Ts. 
    Here :math:`z = x + i y` is a complex number.
    '''
    D = {sp.exp(-z*T): sp.cos(-x*T)*(1j*sp.sin(y*T)+sp.cos(y*T)) for T in Ts}
    expression2 = expression.subs(D)
    num_real,num_imag = expression2.expand().as_real_imag()
    f_r =  ufuncify( [x,y]+Ts, num_real)
    f_i =  ufuncify( [x,y]+Ts, num_imag)
    return lambda x,y,Ts: f_r(x,y,*Ts)+f_i(x,y,*Ts)*1j



In [90]:

    
F = get_real_imag_func(num2)



In [91]:

    
##  testing
F(1,2,(3,4,5,6))









    Out[91]:





(0.14730830564958425-0.40346430265728966j)



In [104]:

    
get_real_imag_func(sp.diff(num2,z))(1,2,(3,4,5,6))









    Out[104]:





(-0.33226636779973417+1.7256217430980794j)



In [117]:

    
def blah(args):
    for i in range(500):
        F(*args)



In [161]:

    
ls = [[0.]*5]*5



In [163]:

    
ls[1][2] = 4



In [164]:

    
ls









    Out[164]:





$$\left [ \left [ 0.0, \quad 0.0, \quad 4, \quad 0.0, \quad 0.0\right ], \quad \left [ 0.0, \quad 0.0, \quad 4, \quad 0.0, \quad 0.0\right ], \quad \left [ 0.0, \quad 0.0, \quad 4, \quad 0.0, \quad 0.0\right ], \quad \left [ 0.0, \quad 0.0, \quad 4, \quad 0.0, \quad 0.0\right ], \quad \left [ 0.0, \quad 0.0, \quad 4, \quad 0.0, \quad 0.0\right ]\right ]$$



In [165]:

    
ls2 = np.zeros((5,5))



In [166]:

    
ls2









    Out[166]:





array([[ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.]])



In [169]:

    
ls2[2,3] = 4



In [181]:

    
ls2









    Out[181]:





array([[ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  4.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.]])



In [174]:

    
ls2[2]









    Out[174]:





array([ 0.,  0.,  0.,  4.,  0.])



In [175]:

    
delays = [1,4,5,6,7]



In [179]:

    
map(sum,zip(ls2[2],delays))









    Out[179]:





$$\left [ 1.0, \quad 4.0, \quad 5.0, \quad 10.0, \quad 7.0\right ]$$



In [183]:

    
round(1.111111,2)









    Out[183]:





$$1.11$$



In [ ]: