In [1]:
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()
import plotly.plotly as py
import plotly.graph_objs as go
import random
# import multiprocessing
# pool = multiprocessing.Pool()
In [2]:
## 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 [3]:
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
In [4]:
phi1, phi2 = sp.symbols('phi_1 phi_2')
In [5]:
## We need to find where ex1(phi1,phi2) + ex2(phi1,omega3) = 0.
ex1 = (k_symb(omega1,pol='e')+k_symb(omega2,pol='e')).expand().subs({omega1:(phi1 + phi2)/2, omega2: (phi1-phi2)/2})
ex2 = -(k_symb(omega3,pol='e')+k_symb(omega4,pol='e')).expand().subs(omega4,-phi1-omega3)
In [10]:
diff_func_4wv_1 = ufuncify([phi1,phi2], ex1)
diff_func_4wv_2 = ufuncify([phi1,omega3], ex2)
In [11]:
phi1_min = 30.
phi1_max = 34.
phi1_range = np.arange(phi1_min,phi1_max,0.01)
phi2_min = -13
phi2_max = -9
phi2_range = np.arange(phi2_min,phi2_max,0.01)
omega3_min = -26.
omega3_max = -16.
omega3_range = np.arange(omega3_min,omega3_max,0.01) ## 5e-5
print len(phi1_range),len(phi2_range),len(omega3_range)
In [13]:
eps = 1e-4
In [14]:
def eps_multiply_digitize(y,eps):
return map(lambda el: int(el/eps), y)
def make_dict_values_to_lists_of_inputs(values,inputs):
D = {}
for k, v in zip(values,inputs):
D.setdefault(k, []).append(v)
return D
In [15]:
f_phi12_omega3 = lambda phi1,phi2,omega3 : (diff_func_4wv_1(phi1_range[phi1],phi2_range[phi2])
- diff_func_4wv_2(phi1_range[phi1],omega3_range[omega3]) )
In [16]:
@timeit
def make_matching_dict_hash(phi1_range,phi2_range,omega3_range,eps=eps):
phi2_indices = range(len(phi2_range))
omega3_indices = range(len(omega3_range))
matching_dict = {}
for phi1_index,phi1 in enumerate(phi1_range):
y1 = diff_func_4wv_1(phi1,phi2_range)
y2 = diff_func_4wv_2(phi1,omega3_range)
y1_rounded = eps_multiply_digitize(y1,eps)
y1_rounded_up = [ind + 1 for ind in y1_rounded]
y2_rounded = eps_multiply_digitize(y2,eps)
y2_rounded_up = [ind + 1 for ind in y2_rounded]
D1 = make_dict_values_to_lists_of_inputs(y1_rounded+y1_rounded_up,2*phi2_indices)
D2 = make_dict_values_to_lists_of_inputs(y2_rounded+y2_rounded_up,2*omega3_indices)
inter = set(D1.keys()) & set(D2.keys())
for el in inter:
for ind1 in D1[el]:
for ind2 in D2[el]:
err = y1[ind1]-y2[ind2]
if abs(err) < eps:
matching_dict[phi1_index,ind1,ind2] = err
return matching_dict
In [17]:
#### e.g.
In [18]:
matching_dict_hash = make_matching_dict_hash(phi1_range,phi2_range,omega3_range,eps=eps)
In [21]:
len(matching_dict_hash)
Out[21]: