In [25]:
from __future__ import division
import sympy
from sympy import *
from sympy import init_printing
init_printing(use_latex='mathjax')
from sympy.utilities.autowrap import ufuncify, autowrap
from copy import deepcopy
import numpy as np
from IPython.display import display
Let's generate some values for $\theta_\text{1}$ and $\theta_\text{2}$ that we can check by hand.
In [2]:
from sympy.abc import D, k, kappa, omega, lamda, eta, alpha, theta
E_s1 = symbols('E_s1')
E_s, E_d, E_eff, h_s, h_d= symbols('E_s E_d E_eff h_s h_d')
k, kappa, omega, D, h = symbols('k kappa omega D h')
Intermediates used to calculate $\theta_\text{I}, \theta_\text{II}$:
In [3]:
intermediates = {alpha : E_eff / E_d,
eta: sqrt(k**2 + kappa**2 / E_s + omega*1j / D),
lamda: (1 - E_eff / E_s)*(k / eta)}
for key, val in intermediates.viewitems():
display(Eq(key, val))
In [4]:
theta_I_bracket_term = (-lamda / tanh(eta * h_s) +
(sinh(k * h_s) * sinh(eta * h_s) + alpha * cosh(k*h_s) * sinh(eta * h_s) -
lamda * (cosh(k*h_s) * cosh(eta*h_s) - 2 + lamda * sinh(k*h_s)/sinh(eta *h_s)))/
(cosh(k*h_s)*sinh(eta*h_s) + alpha * sinh(k*h_s) * sinh(eta*h_s) - lamda * sinh(k * h_s) * cosh(eta * h_s))
)
In [5]:
theta_I = E_s / E_eff * theta_I_bracket_term
In [6]:
theta_I
Out[6]:
In [8]:
theta_I_python = lambdify((k, E_s, E_eff, eta, h_s, alpha, lamda), theta_I)
theta_I_c = autowrap(theta_I, language='C', backend='Cython', args=(k, E_s, E_eff, eta, h_s, alpha, lamda))
theta_I_f = autowrap(theta_I, language='F95', backend='f2py', args=(k, E_s, E_eff, eta, h_s, alpha, lamda))
In [9]:
theta_II = E_s / E_d * (E_eff + (1-lamda)*E_d/tanh(k*h_d))/(E_eff / tanh(k*h_d) + (1 - lamda) * E_d)
theta_II
Out[9]:
In [10]:
{k:0.01, E_s:3.5-0.05j, E_d: 3.5-0.05j, E_eff: 3.5 - 0.5j, }
Out[10]:
In [11]:
%timeit theta_I_python(1, 3, 2, 0.2, 4, 1.02, 1.8)
In [12]:
%timeit theta_I_c(1, 3, 2, 0.2, 4, 1.02, 1.8)
In [13]:
%timeit theta_I_f(1, 3, 2, 0.2, 4, 1.02, 1.8)
In [14]:
((E_s - theta) / (E_s + theta))
Out[14]:
In [15]:
from sympy.functions.elementary.complexes import im
In [16]:
complex_result = (E_s - theta) / (E_s + theta)
In [17]:
dielectric_I = complex_result.subs({theta:theta_I}).subs(intermediates).subs(intermediates)
In [18]:
dielectric_args = (k, omega, kappa, E_s, E_d, E_eff, h_s, D)
dielectric_I_python = lambdify(dielectric_args, dielectric_I, dummify=False, modules="numpy")
In [19]:
%timeit dielectric_I_python(k=1.1, omega=432, kappa=0.3, E_s=3-0.1j, E_d=4-1j, E_eff=4-5j, h_s=0.07, D=0.0077)
In [20]:
dielectric_I_python(k=1.1, omega=432, kappa=0.3, E_s=3-0.1j, E_d=4-1j, E_eff=4-5j, h_s=0.07, D=0.0077)
Out[20]:
In [21]:
lambdify(expr=tanh(h_s * intermediates[eta]), args=[k, kappa, omega, D, h_s, E_s], modules='numpy')(0.01, 0.3, 432, 0.0077, 0.07, 3-1j)
Out[21]:
In [22]:
lambdify(expr=(h_s * intermediates[eta]), args=[k, kappa, omega, D, h_s, E_s], modules='numpy')(0.01, 0.3, 432, 1e-6, 0.07, 3-1j)
Out[22]:
In [23]:
np.tanh(100000)
Out[23]:
In [26]:
sympy.printing.lambdarepr.lambdarepr(dielectric_I)
Out[26]:
In [27]:
_eta = intermediates[eta]
In [28]:
all_k = [1, 10, 100, 1000, 10000, 100000]
all_eta = [_eta.evalf(subs={k:_k, kappa:3500, E_s:3 - 0.001j, D: 0.005, omega:300})
for _k in all_k]
In [29]:
print(str(all_eta).replace('*I', 'j').replace(",", ",\n"))
In [30]:
_lamda = intermediates[lamda]
all_lamda = [_lamda.evalf(subs=
{k:_k, E_s:3-0.001j, D:0.005,
E_eff:3 - 100j, eta: _eta_}) for _k, _eta_ in zip(all_k, all_eta)]
In [31]:
print(str(all_lamda).replace('*I', 'j').replace(",", ",\n"))
In [32]:
all_theta_I = [theta_I.evalf(subs={k:_k, E_s:3-0.001j, D:0.005,
E_eff:3 - 100j, eta: _eta_,
lamda:_lamda_, h_s:0.1, alpha: 0.65-0.0002j}) for _k, _eta, _lamda_ in zip(all_k, all_eta, all_lamda)]
In [33]:
print(str(all_theta_I).replace('*I', 'j').replace(",", ",\n"))
In [34]:
sympy.printing.lambdarepr.lambdarepr(theta_I)
Out[34]:
In [35]:
theta_I_python = lambdify((k, E_s, D, E_eff, eta, lamda, h_s, alpha), theta_I, dummify=False, modules='numpy')
In [36]:
def copy_update_dict(dictionary, **kwargs):
copied_dictionary = deepcopy(dictionary)
copied_dictionary.update(kwargs)
return copied_dictionary
In [37]:
theta_subs = {k: 1, E_s:3-0.001j, D:0.005, E_eff:3 - 100j, eta: 2020.78311260126 + 15.182507854811j, lamda:0.0001184255261724 + 0.0164941987549306j, h_s:0.1, alpha: 0.65-0.0002j}
theta_I_kwargs = {str(key): np.complex128(val) for key, val in theta_subs.items()}
theta_I_all_kwargs = [copy_update_dict(theta_I_kwargs, k=np.float(_k), eta=np.complex128(_eta_), lamda=np.complex128(_lamda_))
for _k, _eta_, _lamda_ in zip(all_k, all_eta, all_lamda)]
all_theta_I_lambda = [theta_I_python(**kwargs) for kwargs in theta_I_all_kwargs]
In [38]:
all_theta_I_lambda
Out[38]:
So for some reason the sympy version of this is implemented better.
In [39]:
theta_subs = {k: 1, E_s:3-0.001j, D:0.005, E_eff:3 - 100j, E_d: 3 - 0.001j, eta: 2020.78311260126 + 15.182507854811j, lamda:0.0001184255261724 + 0.0164941987549306j, h_s:0.1, alpha: 0.65-0.0002j}
all_theta_II = [theta_II.evalf(subs={k: _k, E_s:3-0.001j, E_eff:3 - 100j,
E_d: 3 - 0.001j, lamda:_lamda_, h_d:0.1})
for _k, _lamda_ in zip(all_k, all_lamda)]
In [40]:
print(str(all_theta_II).replace('*I', 'j').replace(",", ",\n"))
To avoid numerical errors in calculating $\theta_I$ (associated with numpy's problems with large values of hyperbolic trigonometry functions for large values of $\theta$ (see Eric Moore's github issue). We have $\theta_I$ is
$$\frac{\epsilon_{\text{s}}}{\epsilon_{\text{eff}}} \left(- \frac{\lambda}{\tanh{\left (\eta h_{\text{s}} \right )}} + \frac{\alpha \sinh{\left (\eta h_{\text{s}} \right )} \cosh{\left (h_{\text{s}} k \right )} - \lambda \left(\frac{\lambda \sinh{\left (h_{\text{s}} k \right )}}{\sinh{\left (\eta h_{\text{s}} \right )}} + \cosh{\left (\eta h_{\text{s}} \right )} \cosh{\left (h_{\text{s}} k \right )} - 2\right) + \sinh{\left (\eta h_{\text{s}} \right )} \sinh{\left (h_{\text{s}} k \right )}}{\alpha \sinh{\left (\eta h_{\text{s}} \right )} \sinh{\left (h_{\text{s}} k \right )} - \lambda \sinh{\left (h_{\text{s}} k \right )} \cosh{\left (\eta h_{\text{s}} \right )} + \sinh{\left (\eta h_{\text{s}} \right )} \cosh{\left (h_{\text{s}} k \right )}}\right)$$
In [42]:
theta_I_numexpr = lambdify((k, E_s, D, E_eff, eta, lamda, h_s, alpha), theta_I, dummify=False, modules='numexpr')
theta_subs = {k: 1, E_s:3-0.001j, D:0.005, E_eff:3 - 100j, eta: 2020.78311260126 + 15.182507854811j, lamda:0.0001184255261724 + 0.0164941987549306j, h_s:0.1, alpha: 0.65-0.0002j}
theta_I_kwargs = {str(key): np.complex128(val) for key, val in theta_subs.items()}
theta_I_all_kwargs = [copy_update_dict(theta_I_kwargs, k=np.float(_k), eta=np.complex128(_eta_), lamda=np.complex128(_lamda_))
for _k, _eta_, _lamda_ in zip(all_k, all_eta, all_lamda)]
all_theta_I_numexpr = [theta_I_numexpr(**kwargs) for kwargs in theta_I_all_kwargs]
In the limit,
$$e^{k h_{\text{s}}} \gg 1$$we can approximate,
$$\coth \eta h_\text{s} = \tanh \eta h_\text{s} = 1$$,
$$\sinh k h_\text{s} = \cosh k h_\text{s} = \exp(k h_\text{s})/2$$,
$$\sinh \eta h_\text{s} = \cosh \eta h_\text{s} = \exp(\eta h_\text{s})/2$$
In [55]:
large_k = {sinh(k*h_s): exp(k*h_s)/2, cosh(k*h_s): exp(k*h_s)/2,
sinh(eta*h_s): exp(eta*h_s)/2, cosh(eta*h_s): exp(eta*h_s)/2,
coth(eta*h_s): 1, tanh(eta*h_s): 1}
theta_I.subs(large_k)
Out[55]:
In [50]:
theta_full = theta_I.subs(intermediates).subs(intermediates)
def theta_k(_k):
return {k: _k, E_s:3-0.001j, D:0.005, E_d:4.65, E_eff:3 - 100j, h_s:0.1, kappa:3500, E_s:3 - 0.001j, D: 0.005, omega:300}
simplified_pre_sub = (E_s/E_eff*(-lamda + (1+alpha-lamda*(1 - 8*exp(-h_s*(k+eta)) + 4*lamda*exp(-2*eta*h_s)))/(1 + alpha - lamda)))
simplified = simplified_pre_sub.subs(intermediates).subs(intermediates)
In [51]:
print(theta_full.subs(theta_k(60)).evalf())
print(simplified.subs(theta_k(60)).evalf())
print(theta_I.subs(large_k).subs(intermediates).subs(intermediates).subs(theta_k(60)).evalf())
In [52]:
eta.subs(intermediates).subs(intermediates).subs(theta_k(80)).evalf()
Out[52]:
In [53]:
def theta_k2(_k, _D, _kappa):
return {k: _k, E_s:3-0.001j, D:_D, E_d:4.65, E_eff:3 - 100j, h_s:0.1, kappa:_kappa, E_s:3 - 0.001j, D: 0.005, omega:300}
In [48]:
print(theta_full.subs(theta_k2(60, 0.005/1e6, 3500*1e6)).evalf())
print(simplified.subs(theta_k2(60, 0.005/1e6, 3500*1e6)).evalf())