Refractive index from:
http://refractiveindex.info/?shelf=main&book=LiNbO3&page=Zelmon-o
There are three frequencies, $\omega_1,\omega_2,\omega_3$. By the frequency matching condition, we impose $\omega_3 = -\omega_1 - \omega_2$.
The user must be watchful of the range for which the index of refraction equation is valid. In the case blow this corresponds to $\omega \in \pm [6,60]$ in units of $10^{13} Hz$.
In the code below, fixed_omega1 is the value dedicated to a fixed value of $\omega_1$. There is a single degree of freedom remaining since $\omega_3$ depends on $\omega_1,\omega_2$. We initialize a list of tuples oms with the values that will correspond to the initial values of $\omega_1$ and $\omega_2$. We must make sure that all three values are inside the admissible range. When the values are in the 'corners' of the ranges, Newton's method convergence slows.
Let's start with 3-wave mixing interactions.
The conditions we need are: \begin{align} &T \left|\omega_1+\omega_2+\omega_3 \right| \le \epsilon_1\\ &\left| L\left(k_1+k_2+k_3\right) \right|= \le \epsilon_2 \implies \left| \frac{L}{c} \left( n(\omega_1)\omega_1+n(\omega_2)\omega_2+n(\omega_3)\omega_3 \right)) \right| \le \epsilon_2. \end{align} In practice we use small cutoffs $\epsilon_1,\epsilon_2$ instead of zeros. These are dimensionless numbers. The $T$ is the duration of the simulation or decay time of the system modes (whichever is smaller), and $\epsilon_1$ maps to the fractional error in the result. $L$ is the legnth of the cavity, and $\epsilon_2$ is chosen so that $L \Delta k / 2$ is sufficiently small so that that $\text{sinc}^2(L \Delta_k /2)$ is small ($L \Delta k / 2$ should be on the order of magnitude of 10).
One way to solve a set of equations is to use Newton's method a few times. We do this here to simulateneously solve the phase and frequency matching conditions for $\chi^{(2)}$ nonlinearity. The main idea is to trace along the solutions, incrementing $\omega_1$ at each step and finding the new $\omega_2$ and $\omega_3$ at each step.
Another point of consideration when considering the cutoffs $\epsilon_1$ and $\epsilon_2$ is the number of modes being dropped. If this number if very large the approximation is not guaranteed to work well.
Later I also interpolate the dispersion relations in the appropriate region.
In [156]:
import sympy as sp
import numpy as np
import scipy.constants
from sympy.utilities.autowrap import ufuncify
import time
from scipy import interpolate
import matplotlib.pyplot as plt
%matplotlib inline
from sympy import init_printing
init_printing()
In [157]:
def plot_arr(arr):
fig = plt.figure(figsize=(15,15))
ax = fig.add_subplot(111)
cax = ax.matshow(np.asmatrix(arr), interpolation='nearest')
fig.colorbar(cax)
plt.show()
In [158]:
## from https://www.andreas-jung.com/contents/a-python-decorator-for-measuring-the-execution-time-of-methods
def timeit(method):
def timed(*args, **kw):
ts = time.time()
result = method(*args, **kw)
te = time.time()
print '%r %2.2f sec' % \
(method.__name__, te-ts)
return result
return timed
In [159]:
lambd,omega,omega1,omega2,omega3,omega4 = sp.symbols('lambda omega omega_1 omega_2 omega_3 omega_4')
l2 = lambd **2
def n_symb(pol='o'):
s = 1.
if pol == 'o':
s += 2.6734 * l2 / (l2 - 0.01764)
s += 1.2290 * l2 / (l2 - 0.05914)
s += 12.614 * l2 / (l2 - 474.6)
else:
s += 2.9804 * l2 / (l2 - 0.02047)
s += 0.5981 * l2 / (l2 - 0.0666)
s += 8.9543 * l2 / (l2 - 416.08)
return sp.sqrt(s)
def k_symb(symbol=omega,pol='o'):
'''k is accurate for omega inputs between 6-60.'''
return ((n_symb(pol=pol) * symbol )
.subs(lambd,scipy.constants.c / (symbol*1e7))) ## / scipy.constants.c
def make_newton_update_expression(dispersion_difference_function,symbols):
expression_diff = [sp.diff(dispersion_difference_function,sym) for sym in symbols]
return [- dispersion_difference_function / der for der in expression_diff]
dispersion_difference_function = k_symb()
update_expression = make_newton_update_expression(dispersion_difference_function.expand(),symbols=[omega])
In [160]:
dispersion_difference_function.expand()
Out[160]:
In [161]:
k_o_vs_omega = ufuncify([omega],k_symb())
k_e_vs_omega = ufuncify([omega],k_symb(pol='e'))
In [172]:
##### Different granularities for the range of omegas
omegas1 = np.asarray([6.+ i*5 for i in range(int(1e1))])
omegas2 = np.asarray([6.+ i*5e-1 for i in range(int(1e2))])
omegas3 = np.asarray([6.+ i*5e-2 for i in range(int(1e3))])
omegas4 = np.asarray([6.+ i*5e-3 for i in range(int(1e4))])
omegas5 = np.asarray([6.+ i*5e-4 for i in range(int(1e5))])
In [173]:
plt.plot(omegas3,k_o_vs_omega(omegas3))
plt.plot(omegas3,k_e_vs_omega(omegas3))
plt.plot(omegas3,omegas3*2.3)
plt.title('Linear relationship, ordinary dispersion, extraordinary dispersion')
Out[173]:
In [ ]:
######
In [174]:
dispersion_difference_function = k_symb(omega1,pol='o')+k_symb(omega2,pol='o')+k_symb(omega3,pol='e')
symbols=[omega1,omega2]
update_expression = make_newton_update_expression(dispersion_difference_function.expand(),
symbols=[omega1,omega2,omega3])
update_expression = [ex.subs(omega3,-omega1-omega2) for ex in update_expression]
In [175]:
newton_update_func_0 = ufuncify( symbols , update_expression[0])
newton_update_func_1 = ufuncify( symbols , update_expression[1])
newton_update_func_2 = ufuncify( symbols , update_expression[2])
omegas_dict = {}
ind_range = range(2,5)
for i in ind_range:
omegas_dict[i] = np.asarray([6.+ j*5*pow(10,-i) for j in range(1+pow(10,i))])
def trace_omegas(omegas,arg2 = 24.):
fixed_omega_1_init=omegas[0]
omegas_delta = omegas[1]-omegas[0]
oms = [(fixed_omega_1_init, arg2)]
converged_points = []
for i in range(int(1e5)):
arg2 += newton_update_func_1(*(oms[-1]))
oms.append((omegas[0],arg2))
if sum([ (oms[-1][i] - oms[-2][i])**2 for i in range(2)]) < omegas_delta * 1e-15:
break
for arg1 in omegas:
for i in range(50):
for j in range(10):
arg2 += newton_update_func_1(*(oms[-1]))
oms.append((arg1,arg2))
if sum([ (oms[-1][k] - oms[-2][k])**2 for k in range(2)]) < omegas_delta * 1e-15:
break
converged_points.append(arg2)
return np.asarray(converged_points)
In [176]:
@timeit
def timed_trace_omegas(*args,**kw):
return trace_omegas(*args, **kw)
@timeit
def interpolate_omegas(omegas_few,omegas_many,omegas_few_values):
f = interpolate.UnivariateSpline(omegas_few, omegas_few_values,k=2,s=1e-10)
return f(omegas_many)
omegas_2_converged = {}
for i in ind_range:
print "i = " + str(i)
omegas_2_converged[i] = timed_trace_omegas(omegas_dict[i])
omega_values_interpolated = interpolate_omegas(omegas_dict[3],
omegas_dict[4],omegas_2_converged[3] )
In [178]:
def trace_omegas_BENCHMARK(omegas,arg2 = 24.):
'''
Trace through triples of frequencies constrained by
:math:`omega1+omega2+omega3 = 0`.
At each step of incrementing omega1, find the new omega2
and omega3 by applying Newton's method to omega2
and the constraint to obtain omega3.
At each step apply 200 Newton's steps to converge.
Args:
omegas (list of floats):
The frequencies in the range desired
fixed_omega_3 (optional [float]): initial value of omega3.
Returns:
oms (list of floats): A history of the Newton's steps
converged_points (list of tuples):
'''
fixed_omega1=omegas[0]
omegas_delta = omegas[1]-omegas[0]
oms = [(fixed_omega1, arg2)] ## k is accurate for omega inputs between 6-60.
converged_points = []
for i in range(int(1e5)):
arg2 += newton_update_func_1(*(oms[-1]))
oms.append((omegas[0],arg2))
if sum([ (oms[-1][i] - oms[-2][i])**2 for i in range(2)]) < omegas_delta * 1e-10 :
break
for arg1 in omegas:
for i in range(200):
arg2 += newton_update_func_1(*(oms[-1]))
oms.append((arg1,arg2))
converged_points.append((arg1,arg2))
return oms,converged_points
In [179]:
def trace_omegas(omegas,arg2 = 24.):
'''
Trace through triples of frequencies constrained by
:math:`omega1+omega2+omega3 = 0`.
At each step of incrementing omega1, find the new omega2
and omega3 by applying Newton's method to omega2
and the constraint to obtain omega3.
Trace through the omegas, when finding the next
step with Newton's method occasionally check the error
and allow terminating early.
Args:
omegas (list of floats):
The frequencies in the range desired
fixed_omega_3 (optional [float]): initial value of omega3.
Returns:
oms (list of floats): A history of the Newton's steps
converged_points (list of tuples):
'''
fixed_omega1=omegas[0]
omegas_delta = omegas[1]-omegas[0]
oms = [(fixed_omega1, arg2)] ## k is accurate for omega inputs between 6-60.
converged_points = []
for i in range(int(1e5)):
arg2 += newton_update_func_1(*(oms[-1]))
oms.append((omegas[0],arg2))
if sum([ (oms[-1][i] - oms[-2][i])**2 for i in range(2)]) < omegas_delta * 1e-10 :
break
for arg1 in omegas:
for i in range(50):
for j in range(10):
arg2 += newton_update_func_1(*(oms[-1]))
oms.append((arg1,arg2))
if sum([ (oms[-1][k] - oms[-2][k])**2 for k in range(2)]) < omegas_delta * 1e-10:
break
converged_points.append((arg1,arg2))
return oms,converged_points
In [180]:
converged_dict = {}
converged_dict_BENCHMARK = {}
In [181]:
for ind,omegas in zip([3,4,5],[omegas3,omegas4,omegas5]):
t0 = time.clock()
oms,converged_points = trace_omegas(omegas)
t1 = time.clock()
print "time to converge: ", t1-t0
k_os = k_o_func(omegas)
k_es = k_e_func(omegas)
converged_omegas = [(arg1,arg2,30-arg1-arg2) for arg1,arg2 in converged_points]
converged_dict[ind] = omegas,converged_omegas
In [182]:
for ind,omegas in zip([3,4,5],[omegas3,omegas4,omegas5]):
t0 = time.clock()
oms,converged_points = trace_omegas_BENCHMARK(omegas)
t1 = time.clock()
print "time to converge: ", t1-t0
k_os = k_o_func(omegas)
k_es = k_e_func(omegas)
converged_omegas_BENCHMARK = [(arg1,arg2,30-arg1-arg2) for arg1,arg2 in converged_points]
converged_dict_BENCHMARK[ind] = omegas,converged_omegas_BENCHMARK
In [183]:
plt.figure(figsize=(20,10))
for ind,(omegas,converged_omegas) in converged_dict.iteritems():
#plt.figure(ind)
plt.plot(omegas[::pow(10,ind-3)],converged_omegas[::pow(10,ind-3)])
In [74]:
# omegas_converged = {}
In [75]:
# for i in range(3,6):
# omegas,converged_omegas = converged_dict[i]
# omegas_converged[i] = [tup[1] for tup in converged_omegas]
Note: Using UnivariateSpline with k=2 and a small value of s seems to offer a good balance of speed and precision.
In [157]:
@timeit
def interpolate_omegas(omegas_few,omegas_many,omegas_few_values):
'''
A method for interpolation from omegas_few to omegas_many.
Args:
omegas_few (list of floats):
points to use to interpolate
omegas_many (list of floats):
points to interpolate to. Should have omegas_few included.
omegas_few_values (list of floats):
The values of the omegas_few that have already been computed
Returns:
A list of length len(omegas_many) interpolating the values of omegas_few.
'''
f = interpolate.UnivariateSpline(omegas_few, omegas_few_values,k=2,s=1e-10)
return f(omegas_many)
In [158]:
omega_values_interpolated = interpolate_omegas(np.asarray(omegas3),np.asarray(omegas5[1:99900]),np.asarray(omegas_converged[3]))
In [159]:
plt.plot(omegas5[1:99900],omega_values_interpolated)
plt.plot(omegas3,omegas_converged[3])
plt.title('interpolated versus original values')
Out[159]:
In [161]:
print 'total abs error', sum(abs(omega_values_interpolated - np.asarray(omegas_converged[5][1:99900])))
print 'max error', max(abs(omega_values_interpolated - np.asarray(omegas_converged[5][1:99900])))
In [68]:
omegas_dict = {}
In [69]:
for i in range(3,6):
omegas_dict[i] = [6.+ j*5*pow(10,-i) for j in range(1+pow(10,i))]
In [70]:
print np.asarray(omegas_dict[5])
In [162]:
print np.asarray(omegas_dict[3])
In [168]:
x = omegas_dict[3]
y = omegas_dict[3]
In [169]:
len(x)
Out[169]:
In [73]:
z = np.zeros((len(x),len(y)))
In [74]:
for i in xrange(len(x)):
for j in xrange(len(y)):
z[i,j] = x[i]+y[j]
In [75]:
f = interpolate.interp2d(x,y,z,kind='q')
In [163]:
x_new = omegas_dict[4]
y_new = omegas_dict[4]
In [167]:
len(x_new)
Out[167]:
In [164]:
@timeit
def get_values(x_new,y_new):
return f(x_new,y_new)
In [165]:
z = get_values(x_new,y_new)