From 54235d1be2fbc7c2129d2e50383a652200b882bb Mon Sep 17 00:00:00 2001 From: Daniel Kotik Date: Thu, 27 Feb 2020 15:24:33 +0100 Subject: [PATCH] Use params dictionary --- Airy_2d/Airy2d.py | 39 +++++++++++++++++++++++++-------------- 1 file changed, 25 insertions(+), 14 deletions(-) diff --git a/Airy_2d/Airy2d.py b/Airy_2d/Airy2d.py index 537d00a..17004d7 100644 --- a/Airy_2d/Airy2d.py +++ b/Airy_2d/Airy2d.py @@ -137,6 +137,8 @@ def main(args): shift = source_shift + r_w chi_rad = math.radians(chi_deg) + params = dict(W_y=w_0, k=k1, M=M, W=W) + # -------------------------------------------------------------------------- # placement of the dielectric interface within the computational cell # -------------------------------------------------------------------------- @@ -172,12 +174,16 @@ def Delta_x(alpha): # -------------------------------------------------------------------------- # beam profile distribution (field amplitude) at the waist of the beam # -------------------------------------------------------------------------- - def Gauss(r, W_y=w_0): + def Gauss(r, params): """Gauss profile.""" + W_y = params['W_y'] + return math.exp(-(r.y / W_y)**2) - def Ai_inc(r, W_y=w_0, M=M, W=W): + def Ai_inc(r, params): """Incomplete Airy function.""" + W_y, M, W = params['W_y'], params['M'], params['W'] + (result, real_tol, imag_tol) = complex_quad(lambda xi: @@ -186,8 +192,8 @@ def Ai_inc(r, W_y=w_0, M=M, W=W): return result if test_output: - print("w_0:", w_0) - print("Airy function 1:", Ai_inc(mp.Vector3(1, -0.3, 1), w_0, 0, 4)) + print("w_0:", params['W_y']) + print("Airy function 1:", Ai_inc(mp.Vector3(1, -0.3, 1), params)) # -------------------------------------------------------------------------- # spectrum amplitude distribution @@ -196,24 +202,28 @@ def Heaviside(x): """Theta (Heaviside step) function.""" return 0 if x < 0 else 1 - def f_Gauss(k_y, W_y=w_0): + def f_Gauss(k_y, params): """Gaussian spectrum amplitude.""" + W_y = params['W_y'] + return math.exp(-(k_y*W_y/2)**2) - def f_Airy(k_y, W_y=w_0, M=M, W=W): + def f_Airy(k_y, params): """Airy spectrum amplitude.""" + W_y, M, W = params['W_y'], params['M'], params['W'] + return W_y*cmath.exp(1.0j*(-1/3)*(k_y*W_y)**3) \ * Heaviside(W_y*k_y - (M-W)) * Heaviside((M+W) - (W_y*k_y)) if test_output: - print("Airy spectrum:", f_Airy(0.2, w_0, 0, 4)) + print("Airy spectrum:", f_Airy(0.2, params)) # -------------------------------------------------------------------------- # plane wave decomposition # (purpose: calculate field amplitude at light source position if not # coinciding with beam waist) # -------------------------------------------------------------------------- - def psi(r, f, x): + def psi(r, f, x, params): """Field amplitude function.""" try: getattr(psi, "called") @@ -230,7 +240,8 @@ def phase(k_y, x, y): (result, real_tol, imag_tol) = complex_quad(lambda k_y: - f(k_y) * cmath.exp(1j*phase(k_y, x, r.y)), + f(k_y, params) * + cmath.exp(1j*phase(k_y, x, r.y)), -k1, k1, limit=100) except Exception as e: print(type(e).__name__ + ":", e) @@ -246,7 +257,7 @@ def phase(k_y, x, y): r = mp.Vector3(0, y, z) print() - print("psi :", psi(r, f_Airy, x)) + print("psi :", psi(r, f_Airy, x, params)) sys.exit() # -------------------------------------------------------------------------- @@ -275,10 +286,10 @@ def phase(k_y, x, y): component=mp.Ez if s_pol else mp.Ey, size=mp.Vector3(0, 9, 0), center=mp.Vector3(source_shift, 0, 0), - #amp_func=lambda r: Gauss(r, w_0) - #amp_func=lambda r: Ai_inc(r, w_0, M, W) - #amp_func=lambda r: psi(r, lambda k_y: f_Gauss(k_y, w_0), shift) - amp_func=lambda r: psi(r, lambda k_y: f_Airy(k_y, w_0, M, W), shift) + #amp_func=lambda r: Gauss(r, params) + #amp_func=lambda r: Ai_inc(r, params) + #amp_func=lambda r: psi(r, f_Gauss, shift, params) + amp_func=lambda r: psi(r, f_Airy, shift, params) ) ]