-
Notifications
You must be signed in to change notification settings - Fork 13
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #14 from DanielKotik/develop
Python package `optbeam`
- Loading branch information
Showing
37 changed files
with
678 additions
and
490 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,11 @@ | ||
| [bumpversion] | ||
| current_version = 1.4.2 | ||
| commit = True | ||
| tag = True | ||
| sign_tags = True | ||
|
|
||
| [bumpversion:file:setup.py] | ||
|
|
||
| [bumpversion:file:optbeam/__init__.py] | ||
|
|
||
| [bumpversion:file:CITATION.cff] |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Empty file.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,136 @@ | ||
|
|
||
| import math | ||
| import sys | ||
|
|
||
| try: | ||
| import cython | ||
| except ModuleNotFoundError: | ||
| cython = None | ||
|
|
||
| if not cython or not cython.compiled: | ||
| from math import (exp as _exp, | ||
| sqrt as _sqrt) | ||
| from cmath import exp as _cexp | ||
|
|
||
| if cython and not cython.compiled: | ||
| print("\nPlease consider compiling `%s.py` via Cython: " | ||
| "`$ cythonize -3 -i %s.py`\n" % (__name__, __name__)) | ||
|
|
||
| from .helpers import _complex_quad, Beam2d | ||
|
|
||
|
|
||
| class Beam2dCartesian(Beam2d): | ||
| """...""" | ||
|
|
||
| def profile(self, r): | ||
| """Field amplitude function. | ||
| Plane wave decomposition: calculate field amplitude at light source | ||
| position if not coinciding with beam waist. | ||
| """ | ||
| self._ry = r.y | ||
| self._rz = r.z | ||
|
|
||
| if not self.called: | ||
| print("Calculating inital field configuration. " | ||
| "This will take some time...") | ||
| self.called = True | ||
|
|
||
| try: | ||
| (result, | ||
| real_tol, | ||
| imag_tol) = _complex_quad(self if cython | ||
| and cython.compiled else self._integrand, | ||
| -self._k, self._k) | ||
| except Exception as e: | ||
| print(type(e).__name__ + ":", e) | ||
| sys.exit() | ||
|
|
||
| return result | ||
|
|
||
| def spectrum(self, k_y): | ||
| """Spectrum amplitude distribution function, f.""" | ||
| raise NotImplementedError | ||
|
|
||
| def _phase(self, k_y, x, y): | ||
| """Phase function.""" | ||
| return x*_sqrt(self._k**2 - k_y**2) + k_y*y | ||
|
|
||
| def _integrand(self, k_y): | ||
| """Integrand function.""" | ||
| return self.spectrum(k_y) * _cexp(1j*self._phase(k_y, self.x, self._ry)) | ||
|
|
||
|
|
||
| # ----------------------------------------------------------------------------- | ||
| # specific beam implementations based on Beam2dCartesian base class | ||
| # ----------------------------------------------------------------------------- | ||
| class Gauss2d(Beam2dCartesian): | ||
| """2d Gauss beam.""" | ||
|
|
||
| def __init__(self, x, params, called=False): | ||
| """...""" | ||
| self._W_y = params['W_y'] | ||
| self._norm = 2 * _sqrt(math.pi) / self._W_y | ||
| super().__init__(x, params, called) | ||
|
|
||
| def profile(self, r): | ||
| """...""" | ||
| # beam profile distribution (field amplitude) at the waist of the beam | ||
| if self.x == 0: | ||
| return self._norm * _exp(-(r.y / self._W_y)**2) | ||
| else: | ||
| return super().profile(r) | ||
|
|
||
| def spectrum(self, k_y): | ||
| """Spectrum amplitude function, f.""" | ||
| return self._f_Gauss(k_y, self._W_y) | ||
|
|
||
| def _f_Gauss(self, k_y, W_y): | ||
| """Gaussian spectrum amplitude.""" | ||
| return _exp(-(k_y*W_y/2)**2) | ||
|
|
||
|
|
||
| class IncAiry2d(Beam2dCartesian): | ||
| """2d incomplete Airy beam.""" | ||
|
|
||
| def __init__(self, x, params, called=False): | ||
| """...""" | ||
| self._W_y = params['W_y'] | ||
| self._M = params['M'] | ||
| self._W = params['W'] | ||
| super().__init__(x, params, called) | ||
|
|
||
| def profile(self, r): | ||
| """...""" | ||
| if self.x == 0: | ||
| # adjust integration boundaries | ||
| self._a = self._M-self._W | ||
| self._b = self._M+self._W | ||
|
|
||
| return super().profile(r) | ||
|
|
||
| def spectrum(self, k_y): | ||
| """...""" | ||
| return self._f_Airy(k_y, self._W_y, self._M, self._W) | ||
|
|
||
| def _f_Airy(self, k_y, W_y, M, W): | ||
| """Airy spectrum amplitude.""" | ||
| return W_y * _cexp(1.0j*(-1/3)*(k_y*W_y)**3) \ | ||
| * self._heaviside(W_y * k_y - (M - W)) \ | ||
| * self._heaviside((M + W) - W_y * k_y) | ||
|
|
||
| def _heaviside(self, x): | ||
| """Theta (Heaviside step) function.""" | ||
| return 0 if x < 0 else 1 | ||
|
|
||
| def _integrand(self, k_y): | ||
| """...""" | ||
| if self.x == 0: | ||
| xi = k_y | ||
| return _cexp(1.0j*(-(xi**3)/3 + (xi * self._ry/self._W_y))) | ||
| else: | ||
| # next line needs _integrand declared cpdef without nogil attribute, | ||
| # and will execute slower than repeating the super class integration | ||
| # routine function here | ||
| #return super(IncAiry2d, self)._integrand(k_y) | ||
| return self.spectrum(k_y) * _cexp(1j*self._phase(k_y, self.x, self._ry)) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,95 @@ | ||
|
|
||
| import sys | ||
|
|
||
| try: | ||
| import cython | ||
| except ModuleNotFoundError: | ||
| cython = None | ||
|
|
||
| if cython: | ||
| if cython.compiled: | ||
| from scipy import LowLevelCallable | ||
| else: | ||
| print("\nPlease consider compiling `%s.py` via Cython: " | ||
| "`$ cythonize -3 -i %s.py`\n" % (__name__, __name__)) | ||
|
|
||
| from scipy.integrate import quad | ||
| from types import MappingProxyType | ||
|
|
||
|
|
||
| def _real_1d_func(x, func): | ||
| """Return real part of a 1d function.""" | ||
| return func(x).real | ||
|
|
||
|
|
||
| def _imag_1d_func(x, func): | ||
| """Return imag part of a 1d function.""" | ||
| return func(x).imag | ||
|
|
||
|
|
||
| def _imag_1d_func_c(n, arr, func_ptr): | ||
| """Return imag part of a 1d function. | ||
| Cython implementation. | ||
| """ | ||
| # pure python formulation of: | ||
| # return (<Beam2dCartesian>func_ptr)(arr[0]).imag | ||
| return cython.cast(Beam2d, func_ptr)._integrand(arr[0]).imag | ||
|
|
||
|
|
||
| def _real_1d_func_c(n, arr, func_ptr): | ||
| """Return real part of a 1d function. | ||
| Cython implementation. | ||
| """ | ||
| # pure python formulation of: | ||
| # return (<Beam2dCartesian>func_ptr)(arr[0]).real | ||
| return cython.cast(Beam2d, func_ptr)._integrand(arr[0]).real | ||
|
|
||
|
|
||
| def _complex_quad(func, a, b, kwargs={}): | ||
| """Integrate real and imaginary part of the given function.""" | ||
| if cython and cython.compiled: | ||
| # pure python formulation of: cdef void *f_ptr = <void*>func | ||
| f_ptr = cython.declare(cython.p_void, cython.cast(cython.p_void, func)) | ||
|
|
||
| func_capsule = PyCapsule_New(f_ptr, cython.NULL, cython.NULL) | ||
|
|
||
| current_module = sys.modules[__name__] | ||
|
|
||
| ll_real_1d_func_c = LowLevelCallable.from_cython(current_module, | ||
| '_real_1d_func_c', | ||
| func_capsule) | ||
| ll_imag_1d_func_c = LowLevelCallable.from_cython(current_module, | ||
| '_imag_1d_func_c', | ||
| func_capsule) | ||
| real, real_tol = quad(ll_real_1d_func_c, a, b, **kwargs) | ||
| imag, imag_tol = quad(ll_imag_1d_func_c, a, b, **kwargs) | ||
| else: | ||
| real, real_tol = quad(_real_1d_func, a, b, (func,), **kwargs) | ||
| imag, imag_tol = quad(_imag_1d_func, a, b, (func,), **kwargs) | ||
|
|
||
| return real + 1j*imag, real_tol, imag_tol | ||
|
|
||
|
|
||
| class Beam2d: | ||
| """Abstract base class.""" | ||
|
|
||
| def __init__(self, x, params, called=False): | ||
| """...""" | ||
| self.x = x | ||
| self._k = params['k'] | ||
| self._params = MappingProxyType(params) # read-only view of a dict | ||
| self.called = called | ||
|
|
||
| @property | ||
| def params(self): | ||
| """Beam specific parameters. | ||
| This is a read-only property. | ||
| """ | ||
| return self._params | ||
|
|
||
| def _integrand(self, x): | ||
| """Integrand function over one coordinate x.""" | ||
| raise NotImplementedError |
Empty file.
Oops, something went wrong.