6.2.3. WPM 3D with storing intermediate data

3D drawings are very nice, but sometimes we have not enough memory to store data. In this case we can use the following functions at Scalar_fields_XY scheme: WPM and BPM.

It is quite easy to execute them at vacuum, as only a XY initial field is required. However it is also possible to use a 3d refraction index structure. In order to not storing it, a function with the XYZ definition is required. However, the data are only extracted for each z position, and then no problems with storage are produced.

In my compute I have got the following results: - 128 x 128 x 128 = 600 ms - 256 x 256 x 256 = 5 seconds - 512 x 512 x 512 = 44 seconds (I could not resolve it in XYZ). - 1024 x 1024 x 1024 = 7 minutes (I could not resolve it in XYZ).

6.2.3.1. Modules

from diffractio import np, plt, sp, um, mm, degrees
from diffractio.scalar_fields_XYZ import Scalar_field_XYZ

from diffractio.scalar_masks_XY import Scalar_mask_XY
from diffractio.scalar_sources_XY import Scalar_source_XY
from diffractio.scalar_masks_XYZ import Scalar_mask_XYZ

from diffractio.utils_math import get_k, ndgrid, nearest
from diffractio.scalar_fields_XY import WPM_schmidt_kernel

import time
import sys

num_z = 256

x = np.linspace(-50 * um, 50 * um, 256)
y = np.linspace(-50 * um, 50 * um, 256)
z = np.linspace(0, 70 * um, num_z)
z_obs = z[-1]
wavelength = 1 * um

t0 = Scalar_mask_XY(x, y, wavelength)
t0.circle(r0=(0, 0), radius=45 * um, angle=0)

u0 = Scalar_source_XY(x, y, wavelength)
u0.gauss_beam(r0=(0, 0), w0=45 * um, z0=0)

u_iter = u0 * t0

def fn(x, y, z, wavelegth):
    u = Scalar_mask_XYZ(x, y, z, wavelength)
    u.sphere(r0=(0, 0, 25),
             radius=20 * um,
             refractive_index=2,
             angles=(0, 0, 0))
    index = u.n
    del u
    return np.squeeze(index)

# Example of structure at a given k-position in z distance

refractive_index = fn(x, y, np.array([
    44,
]), wavelength)

%%time
u_iter.WPM(fn, z_ini=z[0], z_end= z[-1], dz=(z[-1]-z[0])/num_z,
                            has_edges=True, matrix=False, verbose=True)

Time = 5.99 s, time/loop = 23.41 ms
CPU times: user 6 s, sys: 24.9 ms, total: 6.03 s
Wall time: 6.01 s

6.2.3.2. Comparison with 3D

We perform the same experiment, but storing data. Finally, a comparison between both procedures is performed.

u_3d = Scalar_mask_XYZ(x, y, z, wavelength)
u_3d.n = fn(x, y, z, wavelength)

u_3d.incident_field(u0=u0 * t0)

%%time
u_3d.clear_field()
u_3d.WPM(verbose=True, has_edges=True)

Time = 7.17 s, time/loop = 28.0 ms
CPU times: user 7.24 s, sys: 135 ms, total: 7.38 s
Wall time: 7.35 s

u_3d.draw_XZ(y0=0, logarithm=True)
plt.axis('scaled')

<Figure size 432x288 with 0 Axes>

u_3d.draw_XY(z0=50, logarithm=True, normalize='', has_colorbar='vertical')

6.2.3.3. Differences between both algorithms

I_3d_final = np.abs(u_3d.u[:, :, -1])**2
I_iter = np.abs(u_iter.u)**2

u_diffs = Scalar_mask_XY(x, y, wavelength)
u_diffs.u = (I_iter - I_3d_final) / I_iter.max()

u_diffs.draw(has_colorbar='vertical', logarithm=False)
print("Differences = {:2.2f} %".format((100 * np.abs(
    (I_iter - I_3d_final) / I_iter.max()).mean())))

Differences = 3.06 %