Commit 1fd06c6a authored by Leon Merten Lohse's avatar Leon Merten Lohse
Browse files

use consistent scaling and implement CS propagator

parent cb594a78
%% Cell type:code id: tags:
``` python
import sys
import numpy as np
import xraylib as xrl
import matplotlib.pyplot as plt
from pathlib import Path
# path to fresnel (adjust if necessary)
sys.path.append(f'{Path.home()}/fresnel/')
import fresnel.solver as solver
```
%% Cell type:code id: tags:
``` python
from platform import python_version
print(python_version())
```
%% Output
3.8.7
%% Cell type:code id: tags:
``` python
%matplotlib inline
```
%% Cell type:code id: tags:
``` python
import fresnel.solver as solver
```
%% Cell type:code id: tags:
``` python
#xrl.GetCompoundDataNISTByName("Kapton Polyimide Film")
```
%% Cell type:code id: tags:
``` python
# simulation box
energy = 12. # keV
wavelength = 1.2398 / (energy * 1e3) * 1e-3 # mm
nx = 1010
nz = 530
wg_radius = 25e-6
xtot = 10 * wg_radius # mm
ztot = 0.4 # mm
# units
dx = xtot / nx / wavelength
dz = ztot / nz / wavelength
print(f"dx: {dx}, dz: {dz} wavelengths")
```
%% Output
dx: 2.395787247703638, dz: 7304.89092884732 wavelengths
%% Cell type:code id: tags:
``` python
# materials and geometry
density_Ge = xrl.ElementDensity(xrl.SymbolToAtomicNumber('Ge'))
density_Polyimide = 1.42
density_Si = xrl.ElementDensity(xrl.SymbolToAtomicNumber('Si'))
n_Ge = xrl.Refractive_Index('Ge', energy, density_Ge)
n_Si = xrl.Refractive_Index('Si', energy, density_Si)
n_PI = xrl.Refractive_Index('Kapton Polyimide Film', energy, density_Polyimide)
xx = np.linspace(-xtot/2, xtot/2, nx)
core = (np.abs(xx) < wg_radius)
cladding = ~core
refractive_index = n_Ge * cladding + 1 * core
```
%% Cell type:code id: tags:
``` python
boundary = (0, 0,)
u0 = np.exp(-xx**2 / (2 * (2 * wg_radius)**2)).astype(np.complex128)
```
%% Cell type:code id: tags:
``` python
propagator = solver.Propagator2d(refractive_index, u0, dz, dx)
```
%% Cell type:code id: tags:
``` python
field = np.zeros((nz, nx), dtype=np.complex128)
field[0,...] = u0
for iz in range(1, nz):
boundary = (0,0)
field[iz, ...] = propagator.step(refractive_index, boundary)
```
%% Cell type:code id: tags:
``` python
result = field.transpose()
fig, ax = plt.subplots()
im = ax.imshow(np.abs(result[:,:])**2, aspect='auto', clim=(0,3), cmap='jet', extent=[0, ztot, -xtot/2 *1e6, xtot/2 * 1e6])
ax.set_ylim(-2e6 * wg_radius, 2e6 * wg_radius )
fig.colorbar(im)
ax.set_xlabel(r'$z$ (mm)')
ax.set_ylabel(r'$x$ (nm)')
```
%% Output
Text(0, 0.5, '$x$ (nm)')
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
......
%% Cell type:code id: tags:
``` python
import sys
import numpy as np
import xraylib as xrl
import matplotlib.pyplot as plt
from pathlib import Path
# path to fresnel (adjust if necessary)
sys.path.append(f'{Path.home()}/fresnel/')
import fresnel.solver as solver
```
%% Cell type:code id: tags:
``` python
%matplotlib inline
```
%% Cell type:code id: tags:
``` python
import fresnel.solver as solver
```
%% Cell type:code id: tags:
``` python
energy = 12. # keV
# geometry
wavelength = 1.2398 / (energy * 1e3) * 1e-3 # mm
nx = 1024
ny = 512
nz = 1000
ztot = 0.4 # mm
wg_radius = 25e-6
xtot = 10 * wg_radius # mm
ytot = 10 * wg_radius
# units
dx = xtot / nx / wavelength
dy = ytot / ny / wavelength
dz = ztot / nz / wavelength
print(f"dx: {dx}, dy: {dy}, dz: {dz} wavelengths")
```
%% Output
dx: 2.36303234392644, dy: 4.72606468785288, dz: 3871.592192289079 wavelengths
%% Cell type:code id: tags:
``` python
density_Ge = xrl.ElementDensity(xrl.SymbolToAtomicNumber('Ge'))
density_Polyimide = 1.42
density_Si = xrl.ElementDensity(xrl.SymbolToAtomicNumber('Si'))
n_Ge = xrl.Refractive_Index('Ge', energy, density_Ge)
n_Si = xrl.Refractive_Index('Si', energy, density_Si)
n_PI = xrl.Refractive_Index('Kapton Polyimide Film', energy, density_Polyimide)
print('Ge', n_Ge -1)
print('Si', n_Si -1)
print('PI', n_PI -1)
xx = np.linspace(-xtot/2, xtot/2, nx).reshape([1,-1])
yy = np.linspace(-ytot/2, ytot/2, ny).reshape([-1,1])
core = (np.abs(xx) < wg_radius) * (np.abs(yy) < wg_radius)
cladding = ~core
refractive_index = n_Ge * cladding + 1 * core
```
%% Output
Ge (-6.426292114225518e-06+7.120128111534574e-07j)
Si (-3.3845180545943876e-06+3.810012288290064e-08j)
PI (-2.1023289580313076e-06+2.2059592416020462e-09j)
%% Cell type:code id: tags:
``` python
u0 = np.exp(-(xx**2 + yy**2)/ (2 * (2 * wg_radius)**2)).astype(np.complex128)
boundary = (0, 0, 0, 0)
propagator = solver.Propagator3d(refractive_index, u0, dz, dy, dx)
```
%% Cell type:code id: tags:
``` python
field = np.zeros((nz, ny, nx), dtype=np.complex128)
field[0,...] = u0
for iz in range(1, nz):
boundary = (0, 0, 0, 0)
field[iz, ...] = propagator.step(refractive_index, boundary)
if iz % 100 == 0:
print(iz)
```
%% Output
100
200
300
400
500
600
700
800
900
%% Cell type:code id: tags:
``` python
result = field[:,ny//2, :].transpose()
fig, ax = plt.subplots()
im = ax.imshow(np.abs(result)**2, aspect='auto', clim=(0,4), cmap='cubehelix', extent=[0, ztot, -xtot/2 *1e6, xtot/2 * 1e6])
ax.set_ylim(-2e6 * wg_radius, 2e6 * wg_radius )
fig.colorbar(im)
ax.set_xlabel(r'$z$ (mm)')
ax.set_ylabel(r'$x$ (nm)')
result = field[:,:, nx//2].transpose()
fig, ax = plt.subplots()
im = ax.imshow(np.abs(result)**2, aspect='auto', clim=(0,4), cmap='cubehelix', extent=[0, ztot, -xtot/2 *1e6, xtot/2 * 1e6])
ax.set_ylim(-2e6 * wg_radius, 2e6 * wg_radius )
fig.colorbar(im)
ax.set_xlabel(r'$z$ (um)')
ax.set_ylabel(r'$y$ (nm)')
```
%% Output
Text(0, 0.5, '$y$ (nm)')
%% Cell type:code id: tags:
``` python
```
%% Cell type:code id: tags:
``` python
```
......
This diff is collapsed.
......@@ -7,12 +7,12 @@ _k = 2 * np.pi # all lengths are in units of the wavelength
class Solver2d():
"""
Solver for 2d PDEs of the form
u_z = A * u_xx + F * u
Az(x) * u_z = -Axx(x) * u_xx + F(x,z) * u
"""
dtype = np.complex128
def __init__(self, A0, F0, u0, dz, dx):
def __init__(self, Az, Axx, F0, u0, dz, dx):
self._nx = u0.shape[-1]
self._ones = np.ones((self._nx,), dtype = self.dtype)
......@@ -23,77 +23,167 @@ class Solver2d():
self.dx = dx
self.dz = dz
self._compute_coefficients(A0, F0)
self.rz = self._compute_rz(Az)
self.rxx = self._compute_rxx(Axx)
self.f = self._compute_f(F0)
def _compute_rz(self, Az):
return Az * self._ones
def _compute_rxx(self, Axx):
return -Axx * self.dz / 2. / self.dx**2 * self._ones
def _compute_f(self, F):
return F * self.dz / 2. * self._ones
def step(self, F, boundary):
up = self.u
fp = self.f
def _compute_coefficients(self, A, F):
self.f = self._compute_f(F)
u = self._boundary_slice
self.rx = A * self.dz / ( 2. * self.dx**2 ) * self._ones
self.rc = F * self.dz / 2. * self._ones
# set boundary values
u[0] = boundary[0]
u[-1] = boundary[1]
self.u = _solver.step1d_A0F(self.rz, self.rxx, fp, self.f, up, u)
def step(self, A, F, boundary):
return self.u
up = self.u
rxp = self.rx
rcp = self.rc
self._compute_coefficients(A, F)
class Propagator2d(Solver2d):
def __init__(self, n0, u0, dz, dx):
Az = 2j * _k
Axx = 1
F0 = self._compute_potential(n0)
super().__init__(Az, Axx, F0, u0, dz, dx)
def _compute_potential(self,n):
return _k**2 * (1 - n*n)
def step(self, n, boundary):
F = self._compute_potential(n)
return super().step(F, boundary)
class Solver2dfull():
"""
Solver for 2d PDEs of the form
Az(x) * u_z = Axx(x) * u_xx + Ax(x) * u_x + F(x,z) * u
"""
dtype = np.complex128
def __init__(self, Az, Axx, Ax, F0, u0, dz, dx):
self._nx = u0.shape[-1]
self._ones = np.ones((self._nx,), dtype = self.dtype)
self._boundary_slice = np.zeros_like(self._ones)
self.u = u0 * self._ones
self.dx = dx
self.dz = dz
self.rz = self._compute_rz(Az)
self.rxx = self._compute_rxx(Axx)
self.rx = self._compute_rx(Ax)
self.f = self._compute_f(F0)
def _compute_rz(self, Az):
return Az * self._ones
def _compute_rxx(self, Axx):
return -Axx * self.dz / 2. / self.dx**2 * self._ones
def _compute_rx(self, Ax):
return -Ax * self.dz / 4. / self.dx * self._ones
def _compute_f(self, F):
return F * self.dz / 2. * self._ones
def step(self, F, boundary):
up = self.u
fp = self.f
self.f = self._compute_f(F)
u = self._boundary_slice
rx = self.rx
rc = self.rc
# set boundary values
u[0] = boundary[0]
u[-1] = boundary[1]
self.u = _solver.step_AF(up, rxp, rcp, u, rx, rc)
self.u = _solver.step1d_AAF(self.rz, self.rxx, self.rx, fp, self.f, up, u)
return self.u
class Propagator2d(Solver2d):
class PropagatorCS(Solver2dfull):
def __init__(self, n0, u0, dz, dx):
A, F = self._compute_helmholtz_coefficients(n0)
def __init__(self, n0, u0, dz, dx):
super().__init__(A, F, u0, dz, dx)
nx = u0.shape[-1]
self._x = np.linspace(-nx * dx / 2, nx * dx / 2, nx)
Az = 2j * _k * self._x
Axx = self._x
Ax = 1
F0 = self._compute_potential(n0)
super().__init__(Az, Axx, Ax, F0, u0, dz, dx)
def _compute_helmholtz_coefficients(self,n):
A = 1j / (2. * _k)
F = 1j * _k * (n*n - 1) / 2
def _compute_potential(self,n):
return _k**2 * (1 - n*n)
return A, F
def step(self, n, boundary):
A, F = self._compute_helmholtz_coefficients(n)
return super().step(A, F, boundary)
F = self._compute_potential(n) * self._x
return super().step(F, boundary)
class Solver3d():
"""
Solver for equations of the form
u_z = A * u_xx + C * u_yy + F * u
Abbreviations:
'rx' : A * dz / dx**2
'rc' : C * dz / dy**2
'rf : F * dz / 2
u_z = A * u_xx + A * u_yy + C * u
"""
dtype = np.complex128
def __init__(self, A0, F0, u0, dz, dy, dx):
def __init__(self, Az, Axx, Ayy, F0, u