TODO:

incorporate 'quantum' effects of the system.
test nonlinear effects by comparing to known calculations.

First, let's compare the results from the linear and nonlinear time simulations given zero nonlinear coefficient.

Currently, Time_Sim_nonlin.py simulates a nonlinear system using purely classical dynamics. Quantum effects will be incorporated in future versions. This is done by using a non-Hermitian Hamiltonian approach. We simulate both quadrants, even though they will yield equivalent results. The derivation of the 'classical' equations of motion from the fully quantum open system model is describeds in detail in: http://aleph.physik.uni-kl.de/~korsch/papers/JPhysA_43_075306.pdf

More explicitly, we use the $-i*A$ as the equivalent linear non-Hermitian matrix for the $a$ quadrature c-number. We add the nonlinear interaction terms to this Hamiltonian. The non-Hermitian Hamiltonian equations of motion are

$$ \frac{d a}{dt} = i \frac{\partial H}{\partial a^*} \hspace{1in} \frac{d a^*}{dt} =- i \frac{\partial H^*}{\partial a}. $$

Notice we will NOT observe squeezing because the simulation is purely classical. The two quadrants are related by a complex conjugate. In future tests, we will compute the squeezing spectrum for examples in Milburn's book to see if they match our computer model.

Testing

The simplest test is to compare to the linear system in the case when we set the nonlinear coefficient to zero. We do this below.



In [34]:

    
from sympy.physics.quantum import *
from sympy.physics.quantum.boson import *
from sympy.physics.quantum.operatorordering import *



In [35]:

    
import Potapov_Code.Roots as Roots
import Potapov_Code.Potapov as Potapov
import Potapov_Code.Time_Delay_Network as Time_Delay_Network
import Potapov_Code.Time_Sims_nonlin as Time_Sims_nonlin
import Potapov_Code.functions as functions
import Potapov_Code.tests as tests
import Potapov_Code.Hamiltonian as Hamiltonian

import pickle
import numpy as np
import sympy as sp
import numpy.linalg as la
from scipy.integrate import ode
import scipy.constants as consts

import matplotlib.pyplot as plt
%matplotlib inline



In [36]:

    
from sympy import init_printing
init_printing()

Let's try a very simple example:



In [37]:

    
i = 0 ## port



In [38]:

    
displaced_freq =  5e14



In [39]:

    
## Make base system

Ex = Time_Delay_Network.Example1(max_freq=15.,center_freq = 0.)
Ex.run_Potapov()
modes = functions.spatial_modes(Ex.roots,Ex.M1,Ex.E,Ex.delays)
M = len(Ex.roots)









    



(-0.000755759666397-0.000892705709126j)



In [40]:

    
Ex.roots = map(lambda z: z + 1j *(displaced_freq ),Ex.roots)  ## <-- visible frequency ~ 500 THz



In [41]:

    
Ex.roots









    Out[41]:





[(-0.74381183771403225+500000000000000j)]



In [42]:

    
A_d,B_d,C_d,D_d = Ex.get_Potapov_ABCD(doubled=True)
A,B,C,D = Ex.get_Potapov_ABCD(doubled=False)



In [43]:

    
from imp import reload
reload(Hamiltonian)









    Out[43]:





<module 'Potapov_Code.Hamiltonian' from 'Potapov_Code/Hamiltonian.pyc'>



In [44]:

    
ham = Hamiltonian.Hamiltonian(Ex.roots,modes,Ex.delays,Omega=-1j*A,)

ham.nonlin_coeff= 1.

#length_nonlin=0.1*consts.c

ham.chi_nonlinearities = []
ham.make_chi_nonlinearity(delay_indices=[0],start_nonlin=0,
                           length_nonlin=0.1*consts.c,chi_order=3,refraction_index_func=lambda z: 1.)
H = ham.make_H()



In [45]:

    
ham.omegas









    Out[45]:





$$\left [ 7.95774715459e+13\right ]$$



In [46]:

    
import scipy.constants as consts



In [47]:

    
consts.hbar









    Out[47]:





$$1.05457180014e-34$$



In [48]:

    
phase_matching_weights[combination,pm_arr]









    Out[48]:





(2.1907224561234778e-14+2.9758549575832552e-14j)



In [49]:

    
[ham.E_field_weights[i] for i in combination]









    Out[49]:





$$\left [ 6.47765997173e-06, \quad 6.47765997173e-06, \quad 6.47765997173e-06, \quad 6.47765997173e-06\right ]$$



In [50]:

    
chi.chi_function(*chi_args)









    Out[50]:





$$1.0$$



In [51]:

    
filtering_phase_weights=False
eps=1e-5

ham.Hamiltonian_dict_nonlin = {}
for chi in ham.chi_nonlinearities:
    weight_keys = ham.make_weight_keys(chi)

    phase_matching_weights = ham.make_phase_matching_weights(
        weight_keys,chi,filtering_phase_weights=filtering_phase_weights,eps=eps)

    for combination,pm_arr in phase_matching_weights:
        omegas_to_use = map(lambda i: ham.omegas[i],combination)
        omegas_with_sign = [omega * pm for omega,pm
                            in zip(omegas_to_use,pm_arr)]
        pols = map(lambda i: ham.polarizations[i],chi.delay_indices)
        chi_args = omegas_with_sign + pols
        ham.Hamiltonian_dict_nonlin.setdefault( (combination,pm_arr),0.)
        ham.Hamiltonian_dict_nonlin[(combination,pm_arr)] += (
            chi.chi_function(*chi_args) *
            phase_matching_weights[combination,pm_arr] *
            np.prod([ham.E_field_weights[i] for i in combination]) )
        
ham.Hamiltonian_dict_nonlin









    Out[51]:





{((0, 0, 0, 0),
  (-1, -1, -1, -1)): (3.8570915606755913e-35+5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (-1, -1, -1, 1)): (3.8570915606755913e-35+5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (-1, -1, 1, -1)): (3.8570915606755913e-35+5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (-1, -1, 1, 1)): (3.8570915606755913e-35+5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (-1, 1, -1, -1)): (3.8570915606755913e-35+5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (-1, 1, -1, 1)): (3.8570915606755913e-35+5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (-1, 1, 1, -1)): (3.8570915606755913e-35+5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (-1, 1, 1, 1)): (3.8570915606755913e-35+5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (1, -1, -1, -1)): (3.8570915606755913e-35-5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (1, -1, -1, 1)): (3.8570915606755913e-35-5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (1, -1, 1, -1)): (3.8570915606755913e-35-5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (1, -1, 1, 1)): (3.8570915606755913e-35-5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (1, 1, -1, -1)): (3.8570915606755913e-35-5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (1, 1, -1, 1)): (3.8570915606755913e-35-5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (1, 1, 1, -1)): (3.8570915606755913e-35-5.2394336902903591e-35j),
 ((0, 0, 0, 0),
  (1, 1, 1, 1)): (3.8570915606755913e-35-5.2394336902903591e-35j)}



In [52]:

    
normal_order(H.expand())









    Out[52]:





$$500000000000000.0 {{a_0}^\dagger} {a_0} + 0.743811837714032 i {{a_0}^\dagger} {a_0} + 1.54283662427024 \cdot 10^{-34} {{a_0}^\dagger} \left({a_0}\right)^{3} + 1.04788673805807 \cdot 10^{-34} i {{a_0}^\dagger} \left({a_0}\right)^{3} + 2.31425493640536 \cdot 10^{-34} \left({{a_0}^\dagger}\right)^{2} \left({a_0}\right)^{2} + 1.54283662427024 \cdot 10^{-34} \left({{a_0}^\dagger}\right)^{3} {a_0} - 1.04788673805807 \cdot 10^{-34} i \left({{a_0}^\dagger}\right)^{3} {a_0} + 3.85709156067559 \cdot 10^{-35} \left({{a_0}^\dagger}\right)^{4} - 5.23943369029036 \cdot 10^{-35} i \left({{a_0}^\dagger}\right)^{4} + 3.85709156067559 \cdot 10^{-35} \left({a_0}\right)^{4} + 5.23943369029036 \cdot 10^{-35} i \left({a_0}\right)^{4}$$



In [53]:

    
Ex.roots = map(lambda z: z - 1j *(displaced_freq  ),Ex.roots)  ## <-- visible frequency ~ 500 THz

For the Kerr model, we can get detuning by just shifting the frequencies, recomputing the A matrix, and using that to reconstruct the linear Hamiltonian. More generally, we will need to rotate $a \to a e^{-i \Delta t},a^* \to a e^{i \Delta t}$ for detuning $\Delta$. This will be implemented in Potapov_Code.Hamiltonian in the future.

Alternatively, we could modify the inputs $I \to I e^{-i \Delta t}, I^* = I^* e^{i \Delta t}$.



In [54]:

    
detuning = 10.



In [55]:

    
Ex.roots = map(lambda z: z - 1j *(detuning  ),Ex.roots)  ## <-- visible frequency ~ 500 THz



In [56]:

    
Ex.roots









    Out[56]:





[(-0.74381183771403225-10j)]



In [57]:

    
A_d,B_d,C_d,D_d = Ex.get_Potapov_ABCD(doubled=True)
A,B,C,D = Ex.get_Potapov_ABCD(doubled=False)



In [192]:

    
ham.Omega = -1j*A



In [273]:

    
H = ham.make_H()

## Right now, the frequency-matching part is still under construction.
## I am adding this part by hand for now.

ham.H += 22.2*ham.Dagger(ham.a[0])**2*(ham.a[0])**2



In [274]:

    
(normal_order((ham.H).expand())).expand()









    Out[274]:





$$- 10.0 {{a_0}^\dagger} {a_0} + 0.743811837714032 i {{a_0}^\dagger} {a_0} + 1.54283662427024 \cdot 10^{-34} {{a_0}^\dagger} \left({a_0}\right)^{3} + 1.04788673805807 \cdot 10^{-34} i {{a_0}^\dagger} \left({a_0}\right)^{3} + 22.2 \left({{a_0}^\dagger}\right)^{2} \left({a_0}\right)^{2} + 1.54283662427024 \cdot 10^{-34} \left({{a_0}^\dagger}\right)^{3} {a_0} - 1.04788673805807 \cdot 10^{-34} i \left({{a_0}^\dagger}\right)^{3} {a_0} + 3.85709156067559 \cdot 10^{-35} \left({{a_0}^\dagger}\right)^{4} - 5.23943369029036 \cdot 10^{-35} i \left({{a_0}^\dagger}\right)^{4} + 3.85709156067559 \cdot 10^{-35} \left({a_0}\right)^{4} + 5.23943369029036 \cdot 10^{-35} i \left({a_0}\right)^{4}$$



In [275]:

    
## Equations of motion and initial conditions



In [276]:

    
eq_mot = ham.make_eq_motion()
y0 = np.asmatrix([0.]*2*M).T

constant inputs, different input amplitudes.



In [242]:

    
a_in_const = {}
f_const = {}
f_const_lin = {}

input_amplitudes = [1e-1,3e-1,6e-1,1.]
for amp in input_amplitudes:
    a_in_const[amp] = lambda t: np.asmatrix([amp]*np.shape(D_d)[-1]).T  ## make a sample input function
    f_const[amp] = Time_Sims_nonlin.make_f(eq_mot,B_d,a_in_const[amp])
    f_const_lin[amp] = Time_Sims_nonlin.make_f_lin(A_d,B_d,a_in_const[amp])



In [243]:

    
Y_const_lin = {}
Y_const_nonlin = {}
for amp in input_amplitudes:
    Y_const_lin[amp] = Time_Sims_nonlin.run_ODE(f_const_lin[amp], a_in_const[amp], C_d, D_d, 2*M, T = 8., dt = 0.01)  ## select f here.
    Y_const_nonlin[amp] = Time_Sims_nonlin.run_ODE(f_const[amp], a_in_const[amp], C_d, D_d, 2*M, T = 8., dt = 0.01)  ## select f here.



In [105]:

    
plt.figure(figsize=(12,6))
for amp in input_amplitudes:
    #plt.plot( [ (y[i][0,0]).real for y in Y_const_nonlin[amp]],label = 'real part of port '+str(i))
    #plt.plot( [ (y[i][0,0]).imag for y in Y_const_nonlin[amp]],label = 'imaginary part of port '+str(i))
    plt.plot( [ abs(y[i][0,0]) for y in Y_const_nonlin[amp]],label = 'abs part of port '+str(i))
    #plt.plot( [ np.arctan2((y[i][0,0]).real,(y[i][0,0]).imag) for y in Y_const_nonlin[amp]],label = 'phase part of port '+str(i))

    plt.xlabel('time',{'fontsize': 24})
    plt.ylabel('output amplitude',{'fontsize': 24})
    plt.title('Output for nonlinear coeff = ' + str(ham.nonlin_coeff) ,{'fontsize': 24})
    plt.legend(bbox_to_anchor=(1.2, 0.8),loc='upper center',fontsize=20)



In [77]:

    
## linear model for comparison



In [78]:

    
plt.figure(figsize=(12,6))
for amp in input_amplitudes:
    #plt.plot( [ (y[i][0,0]).real for y in Y_const_nonlin[amp]],label = 'real part of port '+str(i))
    #plt.plot( [ (y[i][0,0]).imag for y in Y_const_nonlin[amp]],label = 'imaginary part of port '+str(i))
    plt.plot( [ abs(y[i][0,0]) for y in Y_const_lin[amp]],label = 'abs part of port '+str(i))
    #plt.plot( [ np.arctan2((y[i][0,0]).real,(y[i][0,0]).imag) for y in Y_const_nonlin[amp]],label = 'phase part of port '+str(i))

    plt.xlabel('time',{'fontsize': 24})
    plt.ylabel('output amplitude',{'fontsize': 24})
    plt.title('Output for nonlinear coeff = ' + str(0.) ,{'fontsize': 24})
    plt.legend(bbox_to_anchor=(1.2, 0.8),loc='upper center',fontsize=20)

Ramp inputs

Nonlinear Hamiltonian



In [277]:

    
a_in = lambda t: np.asmatrix([t/300.]*np.shape(D_d)[-1]).T  ## make a sample input function
a_in_2 = lambda t: np.asmatrix([1.-t/300.]*np.shape(D_d)[-1]).T  ## make a sample input function

f = Time_Sims_nonlin.make_f(eq_mot,B_d,a_in)
# f_lin = Time_Sims_nonlin.make_f_lin(A_d,B_d,a_in)

f_2 = Time_Sims_nonlin.make_f(eq_mot,B_d,a_in_2)
# f_lin_2 = Time_Sims_nonlin.make_f_lin(A_d,B_d,a_in_2)



In [278]:

    
# Y_lin = Time_Sims_nonlin.run_ODE(f_lin, a_in, C_d, D_d, 2*M, T = 300., dt = 0.01)  ## select f here.
Y_nonlin = Time_Sims_nonlin.run_ODE(f, a_in, C_d, D_d, 2*M, T = 300., dt = 0.01)  ## select f here.

# Y_lin_2 = Time_Sims_nonlin.run_ODE(f_lin_2, a_in_2, C_d, D_d, 2*M, T = 300., dt = 0.01)  ## select f here.
Y_nonlin_2 = Time_Sims_nonlin.run_ODE(f_2, a_in_2, C_d, D_d, 2*M, T = 300., dt = 0.01)  ## select f here.



In [279]:

    
# Y_lin_2.reverse()
Y_nonlin_2.reverse()



In [280]:

    
plt.figure(figsize=(12,3))

plt.plot( np.asarray([a_in(t)[0,0] for t in np.linspace(0,300.,30001.)]),label = 'input' )
plt.xlabel('time',{'fontsize': 24})
plt.ylabel('Input amplitude',{'fontsize': 24})
plt.title('Input for nonlinear coeff = ' + str(ham.nonlin_coeff) ,{'fontsize': 24})


plt.figure(figsize=(12,6))
for i in range(1):
    #plt.plot( [ (y[i][0,0]).real for y in Y_nonlin],label = 'real part of port '+str(i))
    #plt.plot( [ (y[i][0,0]).imag for y in Y_nonlin],label = 'imaginary part of port '+str(i))
    #plt.plot( [ abs(y[i][0,0]) for y in Y_nonlin],label = 'abs part of port '+str(i))
    plt.plot( [ np.arctan2((y[i][0,0]).real,(y[i][0,0]).imag) for y in Y_nonlin],label = 'phase part of port '+str(i))

    plt.xlabel('time',{'fontsize': 24})
    plt.ylabel('output phase',{'fontsize': 24})
    plt.title('Output for nonlinear coeff = ' + str(ham.nonlin_coeff) ,{'fontsize': 24})
    plt.legend(bbox_to_anchor=(0, 1.4),loc='upper center',fontsize=20)
    
for i in range(1):
    #plt.plot( [ (y[i][0,0]).real for y in Y_nonlin_2],label = 'real part of port '+str(i)+' (inputs right to left)')
    #plt.plot( [ (y[i][0,0]).imag for y in Y_nonlin_2],label = 'imaginary part of port '+str(i)+' (inputs right to left)')
    #plt.plot( [ abs(y[i][0,0]) for y in Y_nonlin_2],label = 'abs part of port '+str(i)+' (inputs right to left)')
    plt.plot( [ np.arctan2((y[i][0,0]).real,(y[i][0,0]).imag) for y in Y_nonlin_2],label = 'phase part of port '+str(i)+' (inputs right to left)')
    plt.xlabel('time',{'fontsize': 24})
    plt.ylabel('output phase',{'fontsize': 24})
    plt.title('Output for nonlinear coeff = ' + str(ham.nonlin_coeff) ,{'fontsize': 24})
    plt.legend(bbox_to_anchor=(0.5, -0.1),loc='upper center',fontsize=20)









    



/Users/gil/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:3: DeprecationWarning: object of type <type 'float'> cannot be safely interpreted as an integer.
  app.launch_new_instance()



In [282]:

    
plt.figure(figsize=(12,3))

plt.plot( np.asarray([a_in(t)[0,0] for t in np.linspace(0,300.,30001.)]),label = 'input' )
plt.xlabel('time',{'fontsize': 24})
plt.ylabel('Input amplitude',{'fontsize': 24})
plt.title('Input for nonlinear coeff = ' + str(ham.nonlin_coeff) ,{'fontsize': 24})


plt.figure(figsize=(12,6))
for i in range(1):
    #plt.plot( [ (y[i][0,0]).real for y in Y_nonlin],label = 'real part of port '+str(i))
    #plt.plot( [ (y[i][0,0]).imag for y in Y_nonlin],label = 'imaginary part of port '+str(i))
    #plt.plot( [ abs(y[i][0,0]) for y in Y_nonlin],label = 'abs part of port '+str(i))
    plt.plot( [ abs(y[i][0,0]) for y in Y_nonlin],label = 'phase part of port '+str(i))

    plt.xlabel('time',{'fontsize': 24})
    plt.ylabel('output phase',{'fontsize': 24})
    plt.title('Output for nonlinear coeff = ' + str(ham.nonlin_coeff) ,{'fontsize': 24})
    plt.legend(bbox_to_anchor=(0, 1.4),loc='upper center',fontsize=20)
    
for i in range(1):
    #plt.plot( [ (y[i][0,0]).real for y in Y_nonlin_2],label = 'real part of port '+str(i)+' (inputs right to left)')
    #plt.plot( [ (y[i][0,0]).imag for y in Y_nonlin_2],label = 'imaginary part of port '+str(i)+' (inputs right to left)')
    #plt.plot( [ abs(y[i][0,0]) for y in Y_nonlin_2],label = 'abs part of port '+str(i)+' (inputs right to left)')
    plt.plot( [ abs((y[i][0,0])) for y in Y_nonlin_2],label = 'phase part of port '+str(i)+' (inputs right to left)')
    plt.xlabel('time',{'fontsize': 24})
    plt.ylabel('output phase',{'fontsize': 24})
    plt.title('Output for nonlinear coeff = ' + str(ham.nonlin_coeff) ,{'fontsize': 24})
    plt.legend(bbox_to_anchor=(0.5, -0.1),loc='upper center',fontsize=20)









    



/Users/gil/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:3: DeprecationWarning: object of type <type 'float'> cannot be safely interpreted as an integer.
  app.launch_new_instance()

Linear equations of motion for comparison



In [250]:

    
plt.figure(figsize=(12,6))
for i in range(1):
    #plt.plot( [ (y[i][0,0]).real for y in Y_lin],label = 'real part of port '+str(i))
    #plt.plot( [ (y[i][0,0]).imag for y in Y_lin],label = 'imaginary part of port '+str(i))
    #plt.plot( [ abs(y[i][0,0]) for y in Y_lin],label = 'abs part of port '+str(i))
    plt.plot( [ np.arctan2((y[i][0,0]).real,(y[i][0,0]).imag) for y in Y_lin],label = 'phase part of port '+str(i))

    plt.xlabel('time',{'fontsize': 24})
    plt.ylabel('output amplitude',{'fontsize': 24})
    plt.title('Output for nonlinear coeff = ' + str(ham.nonlin_coeff) ,{'fontsize': 24})
    plt.legend(bbox_to_anchor=(0, 1.4),loc='upper center',fontsize=20)
    
for i in range(1):
    #plt.plot( [ (y[i][0,0]).real for y in Y_lin_2],label = 'real part of port '+str(i)+' (inputs right to left)')
    #plt.plot( [ (y[i][0,0]).imag for y in Y_lin_2],label = 'imaginary part of port '+str(i)+' (inputs right to left)')
    #plt.plot( [ abs(y[i][0,0]) for y in Y_lin_2],label = 'abs part of port '+str(i)+' (inputs right to left)')
    plt.plot( [ np.arctan2((y[i][0,0]).real,(y[i][0,0]).imag) for y in Y_lin_2],label = 'phase part of port '+str(i)+' (inputs right to left)')
    plt.xlabel('time',{'fontsize': 24})
    plt.ylabel('output amplitude',{'fontsize': 24})
    plt.title('Output for nonlinear coeff = ' + str(ham.nonlin_coeff) ,{'fontsize': 24})
    plt.legend(bbox_to_anchor=(0.5, -0.1),loc='upper center',fontsize=20)

constant inputs, different nonlinear coefficients



In [236]:

    
nonlin_coeffs_list = [0.,1e-1,3e-1,6e-1,1.]

Y_varying_coeffs = {}

for nonlin_coeff in nonlin_coeffs_list:
    ham.nonlin_coeff= nonlin_coeff
    H = ham.make_H()
    eq_mot = ham.make_eq_motion()
    y0 = np.asmatrix([0.]*2*M).T
    a_in = lambda t: np.asmatrix([1.]*np.shape(D_d)[-1]).T  ## make a sample input function
    f = Time_Sims_nonlin.make_f(eq_mot,B_d,a_in)
    Y_varying_coeffs[nonlin_coeff] = Time_Sims_nonlin.run_ODE(f, a_in, C_d, D_d, 2*M, T = 10., dt = 0.01)  ## select f here.



In [42]:

    
for nonlin_coeff in nonlin_coeffs_list:
    plt.figure(figsize=(12,6))
    #plt.plot( [ (y[i][0,0]).real for y in Y_varying_coeffs[nonlin_coeff]],label = 'real part of port '+str(i))
    #plt.plot( [ (y[i][0,0]).imag for y in Y_varying_coeffs[nonlin_coeff]],label = 'imaginary part of port '+str(i))
    plt.plot( [ abs(y[i][0,0]) for y in Y_varying_coeffs[nonlin_coeff]],label = 'abs part of port '+str(i))
    plt.plot( [ np.arctan2((y[i][0,0]).real,(y[i][0,0]).imag) for y in Y_varying_coeffs[nonlin_coeff]],label = 'phase part of port '+str(i))

    plt.xlabel('time',{'fontsize': 24})
    plt.ylabel('output amplitude',{'fontsize': 24})
    plt.title('Output for nonlinear coeff = ' + str(nonlin_coeff) ,{'fontsize': 24})
    plt.legend(bbox_to_anchor=(1.3, 0.9),loc='upper center',fontsize=20)

Officially using the program to move the reference frame to rotating frame



In [43]:

    
displaced_freq =  5e14
i=0



In [44]:

    
## Make base system

Ex = Time_Delay_Network.Example1(max_freq=15.,center_freq = 20.)
Ex.run_Potapov()
modes = functions.spatial_modes(Ex.roots,Ex.M1,Ex.E,Ex.delays)
M = len(Ex.roots)









    



(-0.000769942088025-0.00106256892717j)






    



Potapov_Code/Time_Delay_Network.py:556: RuntimeWarning: overflow encountered in exp
  self.T_denom = lambda z: (1.-self.r* np.exp(-z*self.tau))
Potapov_Code/Time_Delay_Network.py:556: RuntimeWarning: invalid value encountered in cdouble_scalars
  self.T_denom = lambda z: (1.-self.r* np.exp(-z*self.tau))



In [45]:

    
Ex.roots









    Out[45]:





[(-0.74381183771403225+20.943951023931955j)]



In [46]:

    
Ex.roots = map(lambda z: z + 1j *(displaced_freq ),Ex.roots)  ## <-- visible frequency ~ 500 THz



In [47]:

    
ham.chi_nonlinearities[0].delay_indices









    Out[47]:





$$\left [ 0\right ]$$



In [48]:

    
A_d,B_d,C_d,D_d = Ex.get_Potapov_ABCD(doubled=True)
A,B,C,D = Ex.get_Potapov_ABCD(doubled=False)
ham = Hamiltonian.Hamiltonian(Ex.roots,modes,Ex.delays,Omega=-1j*A)
ham.nonlin_coeff= 1.

ham.chi_nonlinearities = []
ham.make_chi_nonlinearity(delay_indices=0,start_nonlin=0,
                           length_nonlin=0.1,chi_order=3,)

H = ham.make_H()



In [49]:

    
H









    Out[49]:





$$1.0 \left(\left(3.98864187600692 \cdot 10^{-35} - 4.67808639867803 \cdot 10^{-35} i\right) {{a_0}^\dagger} {a_0} {{a_0}^\dagger} {a_0} + \left(3.98864187600692 \cdot 10^{-35} - 4.67808639867803 \cdot 10^{-35} i\right) {{a_0}^\dagger} {a_0} \left({{a_0}^\dagger}\right)^{2} + \left(3.98864187600692 \cdot 10^{-35} - 4.67808639867803 \cdot 10^{-35} i\right) {{a_0}^\dagger} \left({a_0}\right)^{2} {{a_0}^\dagger} + \left(3.98864187600692 \cdot 10^{-35} - 4.67808639867803 \cdot 10^{-35} i\right) {{a_0}^\dagger} \left({a_0}\right)^{3} + \left(3.98864187600692 \cdot 10^{-35} - 4.67808639867803 \cdot 10^{-35} i\right) \left({{a_0}^\dagger}\right)^{2} {a_0} {{a_0}^\dagger} + \left(3.98864187600692 \cdot 10^{-35} - 4.67808639867803 \cdot 10^{-35} i\right) \left({{a_0}^\dagger}\right)^{2} \left({a_0}\right)^{2} + \left(3.98864187600692 \cdot 10^{-35} - 4.67808639867803 \cdot 10^{-35} i\right) \left({{a_0}^\dagger}\right)^{3} {a_0} + \left(3.98864187600692 \cdot 10^{-35} - 4.67808639867803 \cdot 10^{-35} i\right) \left({{a_0}^\dagger}\right)^{4} + \left(3.98864187600692 \cdot 10^{-35} + 4.67808639867803 \cdot 10^{-35} i\right) {a_0} {{a_0}^\dagger} {a_0} {{a_0}^\dagger} + \left(3.98864187600692 \cdot 10^{-35} + 4.67808639867803 \cdot 10^{-35} i\right) {a_0} {{a_0}^\dagger} \left({a_0}\right)^{2} + \left(3.98864187600692 \cdot 10^{-35} + 4.67808639867803 \cdot 10^{-35} i\right) {a_0} \left({{a_0}^\dagger}\right)^{2} {a_0} + \left(3.98864187600692 \cdot 10^{-35} + 4.67808639867803 \cdot 10^{-35} i\right) {a_0} \left({{a_0}^\dagger}\right)^{3} + \left(3.98864187600692 \cdot 10^{-35} + 4.67808639867803 \cdot 10^{-35} i\right) \left({a_0}\right)^{2} {{a_0}^\dagger} {a_0} + \left(3.98864187600692 \cdot 10^{-35} + 4.67808639867803 \cdot 10^{-35} i\right) \left({a_0}\right)^{2} \left({{a_0}^\dagger}\right)^{2} + \left(3.98864187600692 \cdot 10^{-35} + 4.67808639867803 \cdot 10^{-35} i\right) \left({a_0}\right)^{3} {{a_0}^\dagger} + \left(3.98864187600692 \cdot 10^{-35} + 4.67808639867803 \cdot 10^{-35} i\right) \left({a_0}\right)^{4}\right) + \left(500000000000021.0 + 0.743811837714032 i\right) {{a_0}^\dagger} {a_0}$$



In [50]:

    
## displace with method next:



In [51]:

    
ham.move_to_rotating_frame(5e14)



In [52]:

    
H = ham.H



In [55]:

    
plt.figure(figsize=(12,6))
#plt.plot( [ (y[i][0,0]).real for y in Y_rotating_frame],label = 'real part of port '+str(i))
#plt.plot( [ (y[i][0,0]).imag for y in Y_rotating_frame,label = 'imaginary part of port '+str(i))
plt.plot( [ abs(y[i][0,0]) for y in Y_rotating_frame],label = 'abs part of port '+str(i))
#plt.plot( [ np.arctan2((y[i][0,0]).real,(y[i][0,0]).imag) for y in Y_rotating_frame],label = 'phase part of port '+str(i))

plt.xlabel('time',{'fontsize': 24})
plt.ylabel('output amplitude',{'fontsize': 24})
plt.title('Output for nonlinear coeff = ' + str(ham.nonlin_coeff) ,{'fontsize': 24})
plt.legend(bbox_to_anchor=(1.3, 0.9),loc='upper center',fontsize=20)









    Out[55]:





<matplotlib.legend.Legend at 0x10d96b910>



In [ ]:



In [ ]:



In [ ]: