<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 TODO
Step1:  Let's try a very simple example
Step2:  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}$.
Step3:  constant inputs, different input amplitudes.
Step4:  Ramp inputs
Nonlinear Hamiltonian
Step5:  Linear equations of motion for comparison
Step6:  constant inputs, different nonlinear coefficients
Step7:  Officially using the program to move the reference frame to rotating frame
<ASSISTANT_TASK:>
Python Code:
from sympy.physics.quantum import *
from sympy.physics.quantum.boson import *
from sympy.physics.quantum.operatorordering import *
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
from sympy import init_printing
init_printing() 
Explanation: 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.
End of explanation
i = 0 ## port
displaced_freq =  5e14
## 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)
Ex.roots = map(lambda z: z + 1j *(displaced_freq ),Ex.roots)  ## <-- visible frequency ~ 500 THz
Ex.roots
A_d,B_d,C_d,D_d = Ex.get_Potapov_ABCD(doubled=True)
A,B,C,D = Ex.get_Potapov_ABCD(doubled=False)
from imp import reload
reload(Hamiltonian)
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()
ham.omegas
import scipy.constants as consts
consts.hbar
phase_matching_weights[combination,pm_arr]
[ham.E_field_weights[i] for i in combination]
chi.chi_function(*chi_args)
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
normal_order(H.expand())
Ex.roots = map(lambda z: z - 1j *(displaced_freq  ),Ex.roots)  ## <-- visible frequency ~ 500 THz
Explanation: Let's try a very simple example:
End of explanation
detuning = 10.
Ex.roots = map(lambda z: z - 1j *(detuning  ),Ex.roots)  ## <-- visible frequency ~ 500 THz
Ex.roots
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.Omega = -1j*A
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
(normal_order((ham.H).expand())).expand()
## Equations of motion and initial conditions
eq_mot = ham.make_eq_motion()
y0 = np.asmatrix([0.]*2*M).T
Explanation: 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}$.
End of explanation
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])
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.
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)
## linear model for comparison
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)
Explanation: constant inputs, different input amplitudes.
End of explanation
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)
# 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.
# Y_lin_2.reverse()
Y_nonlin_2.reverse()
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)
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)
Explanation: Ramp inputs
Nonlinear Hamiltonian
End of explanation
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)
Explanation: Linear equations of motion for comparison
End of explanation
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.
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)
Explanation: constant inputs, different nonlinear coefficients
End of explanation
displaced_freq =  5e14
i=0
## 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)
Ex.roots
Ex.roots = map(lambda z: z + 1j *(displaced_freq ),Ex.roots)  ## <-- visible frequency ~ 500 THz
ham.chi_nonlinearities[0].delay_indices
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()
H
## displace with method next:
ham.move_to_rotating_frame(5e14)
H = ham.H
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)
Explanation: Officially using the program to move the reference frame to rotating frame
End of explanation
<END_TASK>