Skip to content

odak.wave

odak.wave

Provides necessary definitions for merging geometric optics with wave theory and classical approaches in the wave theory as well. See "Introduction to Fourier Optcs" from Joseph Goodman for the theoratical explanation.

adaptive_sampling_angular_spectrum(field, k, distance, dx, wavelength)

A definition to calculate adaptive sampling angular spectrum based beam propagation. For more Zhang, Wenhui, Hao Zhang, and Guofan Jin. "Adaptive-sampling angular spectrum method with full utilization of space-bandwidth product." Optics Letters 45.16 (2020): 4416-4419.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
def adaptive_sampling_angular_spectrum(field, k, distance, dx, wavelength):
    """
    A definition to calculate adaptive sampling angular spectrum based beam propagation. For more Zhang, Wenhui, Hao Zhang, and Guofan Jin. "Adaptive-sampling angular spectrum method with full utilization of space-bandwidth product." Optics Letters 45.16 (2020): 4416-4419.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    iflag = -1
    eps = 10**(-12)
    nv, nu = field.shape
    l = nu*dx
    x = np.linspace(-l/2, l/2, nu)
    y = np.linspace(-l/2, l/2, nv)
    X, Y = np.meshgrid(x, y)
    fx = np.linspace(-1./2./dx, 1./2./dx, nu)
    fy = np.linspace(-1./2./dx, 1./2./dx, nv)
    FX, FY = np.meshgrid(fx, fy)
    forig = 1./2./dx
    fc2 = 1./2*(nu/wavelength/np.abs(distance))**0.5
    ss = np.abs(fc2)/forig
    zc = nu*dx**2/wavelength
    K = nu/2/np.amax(np.abs(fx))
    m = 2
    nnu2 = m*nu
    nnv2 = m*nv
    fxn = np.linspace(-1./2./dx, 1./2./dx, nnu2)
    fyn = np.linspace(-1./2./dx, 1./2./dx, nnv2)
    if np.abs(distance) > zc*2:
        fxn = fxn*ss
        fyn = fyn*ss
    FXN, FYN = np.meshgrid(fxn, fyn)
    Hn = np.exp(1j*k*distance*(1-(FXN*wavelength)**2-(FYN*wavelength)**2)**0.5)
    FX = FXN/np.amax(FXN)*np.pi
    FY = FYN/np.amax(FYN)*np.pi
    t_2 = nufft2(field, FX*ss, FY*ss, size=[nnv2, nnu2], sign=iflag, eps=eps)
    FX = FX/np.amax(FX)*np.pi
    FY = FY/np.amax(FY)*np.pi
    result = nuifft2(Hn*t_2, FX*ss, FY*ss, size=[nv, nu], sign=-iflag, eps=eps)
    return result

add_phase(field, new_phase)

Definition for adding a phase to a given complex field.

Parameters:

  • field – 
           Complex field.
  • new_phase – 
           Complex phase.

Returns:

  • new_field ( complex64 ) –

    Complex field.

Source code in odak/wave/__init__.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
def add_phase(field, new_phase):
    """
    Definition for adding a phase to a given complex field.

    Parameters
    ----------
    field        : np.complex64
                   Complex field.
    new_phase    : np.complex64
                   Complex phase.

    Returns
    -------
    new_field    : np.complex64
                   Complex field.
    """
    phase = calculate_phase(field)
    amplitude = calculate_amplitude(field)
    new_field = amplitude*np.cos(phase+new_phase) + \
        1j*amplitude*np.sin(phase+new_phase)
    return new_field

add_random_phase(field)

Definition for adding a random phase to a given complex field.

Parameters:

  • field – 
           Complex field.

Returns:

  • new_field ( complex64 ) –

    Complex field.

Source code in odak/wave/__init__.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
def add_random_phase(field):
    """
    Definition for adding a random phase to a given complex field.

    Parameters
    ----------
    field        : np.complex64
                   Complex field.

    Returns
    -------
    new_field    : np.complex64
                   Complex field.
    """
    random_phase = np.pi*np.random.random(field.shape)
    new_field = add_phase(field, random_phase)
    return new_field

adjust_phase_only_slm_range(native_range, working_wavelength, native_wavelength)

Definition for calculating the phase range of the Spatial Light Modulator (SLM) for a given wavelength. Here you prove maximum angle as the lower bound is typically zero. If the lower bound isn't zero in angles, you can use this very same definition for calculating lower angular bound as well.

Parameters:

  • native_range – 
                 Native range of the phase only SLM in radians (i.e. two pi).
  • working_wavelength (float) – 
                 Wavelength of the illumination source or some working wavelength.
  • native_wavelength – 
                 Wavelength which the SLM is designed for.

Returns:

  • new_range ( float ) –

    Calculated phase range in radians.

Source code in odak/wave/__init__.py 
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
def adjust_phase_only_slm_range(native_range, working_wavelength, native_wavelength):
    """
    Definition for calculating the phase range of the Spatial Light Modulator (SLM) for a given wavelength. Here you prove maximum angle as the lower bound is typically zero. If the lower bound isn't zero in angles, you can use this very same definition for calculating lower angular bound as well.

    Parameters
    ----------
    native_range       : float
                         Native range of the phase only SLM in radians (i.e. two pi).
    working_wavelength : float
                         Wavelength of the illumination source or some working wavelength.
    native_wavelength  : float
                         Wavelength which the SLM is designed for.

    Returns
    -------
    new_range          : float
                         Calculated phase range in radians.
    """
    new_range = native_range/working_wavelength*native_wavelength
    return new_range

angular_spectrum(field, k, distance, dx, wavelength)

A definition to calculate angular spectrum based beam propagation.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
def angular_spectrum(field, k, distance, dx, wavelength):
    """
    A definition to calculate angular spectrum based beam propagation.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    nv, nu = field.shape
    fx = np.linspace(-1. / 2. /dx, 1. /2. /dx, nu)
    fy = np.linspace(-1. / 2. /dx, 1. /2. /dx, nv)
    FX, FY = np.meshgrid(fx, fy)
    H = np.exp(1j * k * distance * (1 - (FX * wavelength) ** 2 - (FY * wavelength) ** 2) ** 0.5)
    U1 = np.fft.fftshift(np.fft.fft2(field))
    U2 = H * U1
    result = np.fft.ifft2(np.fft.ifftshift(U2))
    return result

band_extended_angular_spectrum(field, k, distance, dx, wavelength)

A definition to calculate bandextended angular spectrum based beam propagation. For more Zhang, Wenhui, Hao Zhang, and Guofan Jin. "Band-extended angular spectrum method for accurate diffraction calculation in a wide propagation range." Optics Letters 45.6 (2020): 1543-1546.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
def band_extended_angular_spectrum(field, k, distance, dx, wavelength):
    """
    A definition to calculate bandextended angular spectrum based beam propagation. For more Zhang, Wenhui, Hao Zhang, and Guofan Jin. "Band-extended angular spectrum method for accurate diffraction calculation in a wide propagation range." Optics Letters 45.6 (2020): 1543-1546.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    iflag = -1
    eps = 10**(-12)
    nv, nu = field.shape
    l = nu*dx
    x = np.linspace(-l/2, l/2, nu)
    y = np.linspace(-l/2, l/2, nv)
    X, Y = np.meshgrid(x, y)
    Z = X**2+Y**2
    fx = np.linspace(-1./2./dx, 1./2./dx, nu)
    fy = np.linspace(-1./2./dx, 1./2./dx, nv)
    FX, FY = np.meshgrid(fx, fy)
    K = nu/2/np.amax(fx)
    fcn = 1./2*(nu/wavelength/np.abs(distance))**0.5
    ss = np.abs(fcn)/np.amax(np.abs(fx))
    zc = nu*dx**2/wavelength
    if np.abs(distance) < zc:
        fxn = fx
        fyn = fy
    else:
        fxn = fx*ss
        fyn = fy*ss
    FXN, FYN = np.meshgrid(fxn, fyn)
    Hn = np.exp(1j*k*distance*(1-(FXN*wavelength)**2-(FYN*wavelength)**2)**0.5)
    X = X/np.amax(X)*np.pi
    Y = Y/np.amax(Y)*np.pi
    t_asmNUFT = nufft2(field, X*ss, Y*ss, sign=iflag, eps=eps)
    result = nuifft2(Hn*t_asmNUFT, X*ss, Y*ss, sign=-iflag, eps=eps)
    return result

band_limited_angular_spectrum(field, k, distance, dx, wavelength)

A definition to calculate bandlimited angular spectrum based beam propagation. For more Matsushima, Kyoji, and Tomoyoshi Shimobaba. "Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields." Optics express 17.22 (2009): 19662-19673.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
def band_limited_angular_spectrum(field, k, distance, dx, wavelength):
    """
    A definition to calculate bandlimited angular spectrum based beam propagation. For more Matsushima, Kyoji, and Tomoyoshi Shimobaba. "Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields." Optics express 17.22 (2009): 19662-19673.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    nv, nu = field.shape
    fx = np.linspace(-1. / 2. /dx, 1. /2. /dx, nu)
    fy = np.linspace(-1. / 2. /dx, 1. /2. /dx, nv)
    FX, FY = np.meshgrid(fx, fy)
    H_exp = np.exp(1j * k * distance * (1 - (FX * wavelength) ** 2 - (FY * wavelength) ** 2) ** 0.5)

    x = dx * float(nu)
    y = dx * float(nv)
    fx_max = 1 / np.sqrt((2 * distance * (1 / x))**2 + 1) / wavelength
    fy_max = 1 / np.sqrt((2 * distance * (1 / y))**2 + 1) / wavelength
    H_filter = ((np.abs(FX) < fx_max) & (np.abs(FY) < fy_max))
    H = H_filter * H_exp

    U1 = np.fft.fftshift(np.fft.fft2(field))
    U2 = H * U1
    result = np.fft.ifft2(np.fft.ifftshift(U2))
    return result

calculate_intensity(field)

Definition to calculate intensity of a single or multiple given electric field(s).

Parameters:

  • field – 
           Electric fields or an electric field.

Returns:

  • intensity ( float ) –

    Intensity or intensities of electric field(s).

Source code in odak/wave/__init__.py 
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
def calculate_intensity(field):
    """
    Definition to calculate intensity of a single or multiple given electric field(s).

    Parameters
    ----------
    field        : ndarray.complex or complex
                   Electric fields or an electric field.

    Returns
    -------
    intensity    : float
                   Intensity or intensities of electric field(s).
    """
    intensity = np.abs(field)**2
    return intensity

distance_between_two_points(point1, point2)

Definition to calculate distance between two given points.

Parameters:

  • point1 – 
          First point in X,Y,Z.
  • point2 – 
          Second point in X,Y,Z.

Returns:

  • distance ( float ) –

    Distance in between given two points.

Source code in odak/tools/vector.py 
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
def distance_between_two_points(point1, point2):
    """
    Definition to calculate distance between two given points.

    Parameters
    ----------
    point1      : list
                  First point in X,Y,Z.
    point2      : list
                  Second point in X,Y,Z.

    Returns
    ----------
    distance    : float
                  Distance in between given two points.
    """
    point1 = np.asarray(point1)
    point2 = np.asarray(point2)
    if len(point1.shape) == 1 and len(point2.shape) == 1:
        distance = np.sqrt(np.sum((point1-point2)**2))
    elif len(point1.shape) == 2 or len(point2.shape) == 2:
        distance = np.sqrt(np.sum((point1-point2)**2, axis=1))
    return distance

double_convergence(nx, ny, k, r, dx)

A definition to generate initial phase for a Gerchberg-Saxton method. For more details consult Sun, Peng, et al. "Holographic near-eye display system based on double-convergence light Gerchberg-Saxton algorithm." Optics express 26.8 (2018): 10140-10151.

Parameters:

  • nx – 
         Size of the output along X.
  • ny – 
         Size of the output along Y.
  • k – 
         See odak.wave.wavenumber for more.
  • r – 
         The distance between location of a light source and an image plane.
  • dx – 
         Pixel pitch.

Returns:

  • function ( ndarray ) –

    Generated phase pattern for a Gerchberg-Saxton method.

Source code in odak/wave/lens.py 
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def double_convergence(nx, ny, k, r, dx):
    """
    A definition to generate initial phase for a Gerchberg-Saxton method. For more details consult Sun, Peng, et al. "Holographic near-eye display system based on double-convergence light Gerchberg-Saxton algorithm." Optics express 26.8 (2018): 10140-10151.

    Parameters
    ----------
    nx         : int
                 Size of the output along X.
    ny         : int
                 Size of the output along Y.
    k          : odak.wave.wavenumber
                 See odak.wave.wavenumber for more.
    r          : float
                 The distance between location of a light source and an image plane.
    dx         : float
                 Pixel pitch.

    Returns
    -------
    function   : ndarray
                 Generated phase pattern for a Gerchberg-Saxton method.
    """
    size = [ny, nx]
    x = np.linspace(-size[0]*dx/2, size[0]*dx/2, size[0])
    y = np.linspace(-size[1]*dx/2, size[1]*dx/2, size[1])
    X, Y = np.meshgrid(x, y)
    Z = X**2+Y**2
    w = np.exp(1j*k*Z/r)
    return w

electric_field_per_plane_wave(amplitude, opd, k, phase=0, w=0, t=0)

Definition to return state of a plane wave at a particular distance and time.

Parameters:

  • amplitude – 
           Amplitude of a wave.
  • opd – 
           Optical path difference in mm.
  • k – 
           Wave number of a wave, see odak.wave.parameters.wavenumber for more.
  • phase – 
           Initial phase of a wave.
  • w – 
           Rotation speed of a wave, see odak.wave.parameters.rotationspeed for more.
  • t – 
           Time in seconds.

Returns:

  • field ( complex ) –

    A complex number that provides the resultant field in the complex form A*e^(j(wt+phi)).

Source code in odak/wave/vector.py 
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
def electric_field_per_plane_wave(amplitude, opd, k, phase=0, w=0, t=0):
    """
    Definition to return state of a plane wave at a particular distance and time.

    Parameters
    ----------
    amplitude    : float
                   Amplitude of a wave.
    opd          : float
                   Optical path difference in mm.
    k            : float
                   Wave number of a wave, see odak.wave.parameters.wavenumber for more.
    phase        : float
                   Initial phase of a wave.
    w            : float
                   Rotation speed of a wave, see odak.wave.parameters.rotationspeed for more.
    t            : float
                   Time in seconds.

    Returns
    -------
    field        : complex
                   A complex number that provides the resultant field in the complex form A*e^(j(wt+phi)).
    """
    field = amplitude*np.exp(1j*(-w*t+opd*k+phase))/opd**2
    return field

fraunhofer(field, k, distance, dx, wavelength)

A definition to calculate Fraunhofer based beam propagation.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
def fraunhofer(field, k, distance, dx, wavelength):
    """
    A definition to calculate Fraunhofer based beam propagation.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    nv, nu = field.shape
    l = nu*dx
    l2 = wavelength*distance/dx
    dx2 = wavelength*distance/l
    fx = np.linspace(-l2/2., l2/2., nu)
    fy = np.linspace(-l2/2., l2/2., nv)
    FX, FY = np.meshgrid(fx, fy)
    FZ = FX**2+FY**2
    c = np.exp(1j*k*distance)/(1j*wavelength*distance) * \
        np.exp(1j*k/(2*distance)*FZ)
    result = c*np.fft.ifftshift(np.fft.fft2(np.fft.fftshift(field)))*dx**2
    return result

fraunhofer_equal_size_adjust(field, distance, dx, wavelength)

A definition to match the physical size of the original field with the propagated field.

Parameters:

  • field – 
               Complex field (MxN).
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • new_field ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
def fraunhofer_equal_size_adjust(field, distance, dx, wavelength):
    """
    A definition to match the physical size of the original field with the propagated field.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    new_field        : np.complex
                       Final complex field (MxN).
    """
    nv, nu = field.shape
    l1 = nu*dx
    l2 = wavelength*distance/dx
    m = l1/l2
    px = int(m*nu)
    py = int(m*nv)
    nx = int(field.shape[0]/2-px/2)
    ny = int(field.shape[1]/2-py/2)
    new_field = np.copy(field[nx:nx+px, ny:ny+py])
    return new_field

fraunhofer_inverse(field, k, distance, dx, wavelength)

A definition to calculate Inverse Fraunhofer based beam propagation.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
def fraunhofer_inverse(field, k, distance, dx, wavelength):
    """
    A definition to calculate Inverse Fraunhofer based beam propagation.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    distance = np.abs(distance)
    nv, nu = field.shape
    l = nu*dx
    l2 = wavelength*distance/dx
    dx2 = wavelength*distance/l
    fx = np.linspace(-l2/2., l2/2., nu)
    fy = np.linspace(-l2/2., l2/2., nv)
    FX, FY = np.meshgrid(fx, fy)
    FZ = FX**2+FY**2
    c = np.exp(1j*k*distance)/(1j*wavelength*distance) * \
        np.exp(1j*k/(2*distance)*FZ)
    result = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(field/dx**2/c)))
    return result

generate_complex_field(amplitude, phase)

Definition to generate a complex field with a given amplitude and phase.

Parameters:

  • amplitude – 
                Amplitude of the field.
  • phase – 
                Phase of the field.

Returns:

  • field ( ndarray ) –

    Complex field.

Source code in odak/wave/__init__.py 
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
def generate_complex_field(amplitude, phase):
    """
    Definition to generate a complex field with a given amplitude and phase.

    Parameters
    ----------
    amplitude         : ndarray
                        Amplitude of the field.
    phase             : ndarray
                        Phase of the field.

    Returns
    -------
    field             : ndarray
                        Complex field.
    """
    field = amplitude*np.cos(phase)+1j*amplitude*np.sin(phase)
    return field

gerchberg_saxton(field, n_iterations, distance, dx, wavelength, slm_range=6.28, propagation_type='IR Fresnel', initial_phase=None)

Definition to compute a hologram using an iterative method called Gerchberg-Saxton phase retrieval algorithm. For more on the method, see: Gerchberg, Ralph W. "A practical algorithm for the determination of phase from image and diffraction plane pictures." Optik 35 (1972): 237-246.

Parameters:

  • field – 
               Complex field (MxN).
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.
  • slm_range – 
               Typically this is equal to two pi. See odak.wave.adjust_phase_only_slm_range() for more.
  • propagation_type (str, default: 'IR Fresnel' ) – 
               Type of the propagation (IR Fresnel, TR Fresnel, Fraunhofer).
  • initial_phase – 
               Phase to be added to the initial value.

Returns:

  • hologram ( complex ) –

    Calculated complex hologram.

  • reconstruction ( complex ) –

    Calculated reconstruction using calculated hologram.

Source code in odak/wave/classical.py 
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
def gerchberg_saxton(field, n_iterations, distance, dx, wavelength, slm_range=6.28, propagation_type='IR Fresnel', initial_phase=None):
    """
    Definition to compute a hologram using an iterative method called Gerchberg-Saxton phase retrieval algorithm. For more on the method, see: Gerchberg, Ralph W. "A practical algorithm for the determination of phase from image and diffraction plane pictures." Optik 35 (1972): 237-246.

    Parameters
    ----------
    field            : np.complex64
                       Complex field (MxN).
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.
    slm_range        : float
                       Typically this is equal to two pi. See odak.wave.adjust_phase_only_slm_range() for more.
    propagation_type : str
                       Type of the propagation (IR Fresnel, TR Fresnel, Fraunhofer).
    initial_phase    : np.complex64
                       Phase to be added to the initial value.

    Returns
    -------
    hologram         : np.complex
                       Calculated complex hologram.
    reconstruction   : np.complex
                       Calculated reconstruction using calculated hologram. 
    """
    k = wavenumber(wavelength)
    target = calculate_amplitude(field)
    hologram = generate_complex_field(np.ones(field.shape), 0)
    hologram = zero_pad(hologram)
    if type(initial_phase) == type(None):
        hologram = add_random_phase(hologram)
    else:
        initial_phase = zero_pad(initial_phase)
        hologram = add_phase(hologram, initial_phase)
    center = [int(hologram.shape[0]/2.), int(hologram.shape[1]/2.)]
    orig_shape = [int(field.shape[0]/2.), int(field.shape[1]/2.)]
    for i in tqdm(range(n_iterations), leave=False):
        reconstruction = propagate_beam(
            hologram, k, distance, dx, wavelength, propagation_type)
        new_target = calculate_amplitude(reconstruction)
        new_target[
            center[0]-orig_shape[0]:center[0]+orig_shape[0],
            center[1]-orig_shape[1]:center[1]+orig_shape[1]
        ] = target
        reconstruction = generate_complex_field(
            new_target, calculate_phase(reconstruction))
        hologram = propagate_beam(
            reconstruction, k, -distance, dx, wavelength, propagation_type)
        hologram = generate_complex_field(1, calculate_phase(hologram))
        hologram = hologram[
            center[0]-orig_shape[0]:center[0]+orig_shape[0],
            center[1]-orig_shape[1]:center[1]+orig_shape[1],
        ]
        hologram = zero_pad(hologram)
    reconstruction = propagate_beam(
        hologram, k, distance, dx, wavelength, propagation_type)
    hologram = hologram[
        center[0]-orig_shape[0]:center[0]+orig_shape[0],
        center[1]-orig_shape[1]:center[1]+orig_shape[1]
    ]
    reconstruction = reconstruction[
        center[0]-orig_shape[0]:center[0]+orig_shape[0],
        center[1]-orig_shape[1]:center[1]+orig_shape[1]
    ]
    return hologram, reconstruction

gerchberg_saxton_3d(fields, n_iterations, distances, dx, wavelength, slm_range=6.28, propagation_type='IR Fresnel', initial_phase=None, target_type='no constraint', coefficients=None)

Definition to compute a multi plane hologram using an iterative method called Gerchberg-Saxton phase retrieval algorithm. For more on the method, see: Zhou, Pengcheng, et al. "30.4: Multi‐plane holographic display with a uniform 3D Gerchberg‐Saxton algorithm." SID Symposium Digest of Technical Papers. Vol. 46. No. 1. 2015.

Parameters:

  • fields – 
               Complex fields (MxN).
  • distances – 
               Propagation distances.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.
  • slm_range – 
               Typically this is equal to two pi. See odak.wave.adjust_phase_only_slm_range() for more.
  • propagation_type (str, default: 'IR Fresnel' ) – 
               Type of the propagation (IR Fresnel, TR Fresnel, Fraunhofer).
  • initial_phase – 
               Phase to be added to the initial value.
  • target_type – 
               Target type. `No constraint` targets the input target as is. `Double constraint` follows the idea in this paper, which claims to suppress speckle: Chang, Chenliang, et al. "Speckle-suppressed phase-only holographic three-dimensional display based on double-constraint Gerchberg–Saxton algorithm." Applied optics 54.23 (2015): 6994-7001.

Returns:

  • hologram ( complex ) –

    Calculated complex hologram.

Source code in odak/wave/classical.py 
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
def gerchberg_saxton_3d(fields, n_iterations, distances, dx, wavelength, slm_range=6.28, propagation_type='IR Fresnel', initial_phase=None, target_type='no constraint', coefficients=None):
    """
    Definition to compute a multi plane hologram using an iterative method called Gerchberg-Saxton phase retrieval algorithm. For more on the method, see: Zhou, Pengcheng, et al. "30.4: Multi‐plane holographic display with a uniform 3D Gerchberg‐Saxton algorithm." SID Symposium Digest of Technical Papers. Vol. 46. No. 1. 2015.

    Parameters
    ----------
    fields           : np.complex64
                       Complex fields (MxN).
    distances        : list
                       Propagation distances.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.
    slm_range        : float
                       Typically this is equal to two pi. See odak.wave.adjust_phase_only_slm_range() for more.
    propagation_type : str
                       Type of the propagation (IR Fresnel, TR Fresnel, Fraunhofer).
    initial_phase    : np.complex64
                       Phase to be added to the initial value.
    target_type      : str
                       Target type. `No constraint` targets the input target as is. `Double constraint` follows the idea in this paper, which claims to suppress speckle: Chang, Chenliang, et al. "Speckle-suppressed phase-only holographic three-dimensional display based on double-constraint Gerchberg–Saxton algorithm." Applied optics 54.23 (2015): 6994-7001. 

    Returns
    -------
    hologram         : np.complex
                       Calculated complex hologram.
    """
    k = wavenumber(wavelength)
    targets = calculate_amplitude(np.asarray(fields)).astype(np.float64)
    hologram = generate_complex_field(np.ones(targets[0].shape), 0)
    hologram = zero_pad(hologram)
    if type(initial_phase) == type(None):
        hologram = add_random_phase(hologram)
    else:
        initial_phase = zero_pad(initial_phase)
        hologram = add_phase(hologram, initial_phase)
    center = [int(hologram.shape[0]/2.), int(hologram.shape[1]/2.)]
    orig_shape = [int(fields[0].shape[0]/2.), int(fields[0].shape[1]/2.)]
    holograms = np.zeros(
        (len(distances), hologram.shape[0], hologram.shape[1]), dtype=np.complex64)
    for i in tqdm(range(n_iterations), leave=False):
        for distance_id in tqdm(range(len(distances)), leave=False):
            distance = distances[distance_id]
            reconstruction = propagate_beam(
                hologram, k, distance, dx, wavelength, propagation_type)
            if target_type == 'double constraint':
                if type(coefficients) == type(None):
                    raise Exception(
                        "Provide coeeficients of alpha,beta and gamma for double constraint.")
                alpha = coefficients[0]
                beta = coefficients[1]
                gamma = coefficients[2]
                target_current = 2*alpha * \
                    np.copy(targets[distance_id])-beta * \
                    calculate_amplitude(reconstruction)
                target_current[target_current == 0] = gamma * \
                    np.abs(reconstruction[target_current == 0])
            elif target_type == 'no constraint':
                target_current = np.abs(targets[distance_id])
            new_target = calculate_amplitude(reconstruction)
            new_target[
                center[0]-orig_shape[0]:center[0]+orig_shape[0],
                center[1]-orig_shape[1]:center[1]+orig_shape[1]
            ] = target_current
            reconstruction = generate_complex_field(
                new_target, calculate_phase(reconstruction))
            hologram_layer = propagate_beam(
                reconstruction, k, -distance, dx, wavelength, propagation_type)
            hologram_layer = generate_complex_field(
                1., calculate_phase(hologram_layer))
            hologram_layer = hologram_layer[
                center[0]-orig_shape[0]:center[0]+orig_shape[0],
                center[1]-orig_shape[1]:center[1]+orig_shape[1]
            ]
            hologram_layer = zero_pad(hologram_layer)
            holograms[distance_id] = hologram_layer
        hologram = np.sum(holograms, axis=0)
    hologram = hologram[
        center[0]-orig_shape[0]:center[0]+orig_shape[0],
        center[1]-orig_shape[1]:center[1]+orig_shape[1]
    ]
    return hologram

impulse_response_fresnel(field, k, distance, dx, wavelength)

A definition to calculate impulse response based Fresnel approximation for beam propagation.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
def impulse_response_fresnel(field, k, distance, dx, wavelength):
    """
    A definition to calculate impulse response based Fresnel approximation for beam propagation.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).

    """
    nv, nu = field.shape
    x = np.linspace(-nu / 2 * dx, nu / 2 * dx, nu)
    y = np.linspace(-nv / 2 * dx, nv / 2 * dx, nv)
    X, Y = np.meshgrid(x, y)
    h = 1. / (1j * wavelength * distance) * np.exp(1j * k / (2 * distance) * (X ** 2 + Y ** 2))
    H = np.fft.fft2(np.fft.fftshift(h))*dx**2
    U1 = np.fft.fft2(np.fft.fftshift(field))
    U2 = H * U1
    result = np.fft.ifftshift(np.fft.ifft2(U2)) /(dx**2)

    return result

linear_grating(nx, ny, every=2, add=3.14, axis='x')

A definition to generate a linear grating.

Parameters:

  • nx – 
         Size of the output along X.
  • ny – 
         Size of the output along Y.
  • every – 
         Add the add value at every given number.
  • add – 
         Angle to be added.
  • axis – 
         Axis eiter X,Y or both.

Returns:

  • field ( ndarray ) –

    Linear grating term.

Source code in odak/wave/lens.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
def linear_grating(nx, ny, every=2, add=3.14, axis='x'):
    """
    A definition to generate a linear grating.

    Parameters
    ----------
    nx         : int
                 Size of the output along X.
    ny         : int
                 Size of the output along Y.
    every      : int
                 Add the add value at every given number.
    add        : float
                 Angle to be added.
    axis       : string
                 Axis eiter X,Y or both.

    Returns
    -------
    field      : ndarray
                 Linear grating term.
    """
    grating = np.zeros((nx, ny), dtype=np.complex64)
    if axis == 'x':
        grating[::every, :] = np.exp(1j*add)
    if axis == 'y':
        grating[:, ::every] = np.exp(1j*add)
    if axis == 'xy':
        checker = np.indices((nx, ny)).sum(axis=0) % every
        checker += 1
        checker = checker % 2
        grating = np.exp(1j*checker*add)
    return grating

nufft2(field, fx, fy, size=None, sign=1, eps=10 ** -12)

A definition to take 2D Non-Uniform Fast Fourier Transform (NUFFT).

Parameters:

  • field – 
          Input field.
  • fx – 
          Frequencies along x axis.
  • fy – 
          Frequencies along y axis.
  • size – 
          Size.
  • sign – 
          Sign of the exponential used in NUFFT kernel.
  • eps – 
          Accuracy of NUFFT.

Returns:

  • result ( ndarray ) –

    Inverse NUFFT of the input field.

Source code in odak/tools/matrix.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
def nufft2(field, fx, fy, size=None, sign=1, eps=10**(-12)):
    """
    A definition to take 2D Non-Uniform Fast Fourier Transform (NUFFT).

    Parameters
    ----------
    field       : ndarray
                  Input field.
    fx          : ndarray
                  Frequencies along x axis.
    fy          : ndarray
                  Frequencies along y axis.
    size        : list
                  Size.
    sign        : float
                  Sign of the exponential used in NUFFT kernel.
    eps         : float
                  Accuracy of NUFFT.

    Returns
    ----------
    result      : ndarray
                  Inverse NUFFT of the input field.
    """
    try:
        import finufft
    except:
        print('odak.tools.nufft2 requires finufft to be installed: pip install finufft')
    image = np.copy(field).astype(np.complex128)
    result = finufft.nufft2d2(
        fx.flatten(), fy.flatten(), image, eps=eps, isign=sign)
    if type(size) == type(None):
        result = result.reshape(field.shape)
    else:
        result = result.reshape(size)
    return result

nuifft2(field, fx, fy, size=None, sign=1, eps=10 ** -12)

A definition to take 2D Adjoint Non-Uniform Fast Fourier Transform (NUFFT).

Parameters:

  • field – 
          Input field.
  • fx – 
          Frequencies along x axis.
  • fy – 
          Frequencies along y axis.
  • size – 
          Shape of the NUFFT calculated for an input field.
  • sign – 
          Sign of the exponential used in NUFFT kernel.
  • eps – 
          Accuracy of NUFFT.

Returns:

  • result ( ndarray ) –

    NUFFT of the input field.

Source code in odak/tools/matrix.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
def nuifft2(field, fx, fy, size=None, sign=1, eps=10**(-12)):
    """
    A definition to take 2D Adjoint Non-Uniform Fast Fourier Transform (NUFFT).

    Parameters
    ----------
    field       : ndarray
                  Input field.
    fx          : ndarray
                  Frequencies along x axis.
    fy          : ndarray
                  Frequencies along y axis.
    size        : list or ndarray
                  Shape of the NUFFT calculated for an input field.
    sign        : float
                  Sign of the exponential used in NUFFT kernel.
    eps         : float
                  Accuracy of NUFFT.

    Returns
    ----------
    result      : ndarray
                  NUFFT of the input field.
    """
    try:
        import finufft
    except:
        print('odak.tools.nuifft2 requires finufft to be installed: pip install finufft')
    image = np.copy(field).astype(np.complex128)
    if type(size) == type(None):
        result = finufft.nufft2d1(
            fx.flatten(),
            fy.flatten(),
            image.flatten(),
            image.shape,
            eps=eps,
            isign=sign
        )
    else:
        result = finufft.nufft2d1(
            fx.flatten(),
            fy.flatten(),
            image.flatten(),
            (size[0], size[1]),
            eps=eps,
            isign=sign
        )
    result = np.asarray(result)
    return result

prism_phase_function(nx, ny, k, angle, dx=0.001, axis='x')

A definition to generate 2D phase function that represents a prism. See Goodman's Introduction to Fourier Optics book for more.

Parameters:

  • nx – 
         Size of the output along X.
  • ny – 
         Size of the output along Y.
  • k – 
         See odak.wave.wavenumber for more.
  • angle – 
         Tilt angle of the prism in degrees.
  • dx – 
         Pixel pitch.
  • axis – 
         Axis of the prism.

Returns:

  • prism ( ndarray ) –

    Generated phase function for a prism.

Source code in odak/wave/lens.py 
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
def prism_phase_function(nx, ny, k, angle, dx=0.001, axis='x'):
    """
    A definition to generate 2D phase function that represents a prism. See Goodman's Introduction to Fourier Optics book for more.

    Parameters
    ----------
    nx         : int
                 Size of the output along X.
    ny         : int
                 Size of the output along Y.
    k          : odak.wave.wavenumber
                 See odak.wave.wavenumber for more.
    angle      : float
                 Tilt angle of the prism in degrees.
    dx         : float
                 Pixel pitch.
    axis       : str
                 Axis of the prism.

    Returns
    -------
    prism      : ndarray
                 Generated phase function for a prism.
    """
    angle = np.radians(angle)
    size = [ny, nx]
    x = np.linspace(-size[0]*dx/2, size[0]*dx/2, size[0])
    y = np.linspace(-size[1]*dx/2, size[1]*dx/2, size[1])
    X, Y = np.meshgrid(x, y)
    if axis == 'y':
        prism = np.exp(-1j*k*np.sin(angle)*Y)
    elif axis == 'x':
        prism = np.exp(-1j*k*np.sin(angle)*X)
    return prism

produce_phase_only_slm_pattern(hologram, slm_range, filename=None, bits=8, default_range=6.28, illumination=None)

Definition for producing a pattern for a phase only Spatial Light Modulator (SLM) using a given field.

Parameters:

  • hologram – 
                 Input holographic field.
  • slm_range – 
                 Range of the phase only SLM in radians for a working wavelength (i.e. two pi). See odak.wave.adjust_phase_only_slm_range() for more.
  • filename – 
                 Optional variable, if provided the patterns will be save to given location.
  • bits – 
                 Quantization bits.
  • default_range – 
                 Default range of phase only SLM.
  • illumination – 
                 Spatial illumination distribution.

Returns:

  • pattern ( complex64 ) –

    Adjusted phase only pattern.

  • hologram_digital ( int ) –

    Digital representation of the hologram.

Source code in odak/wave/__init__.py 
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
def produce_phase_only_slm_pattern(hologram, slm_range, filename=None, bits=8, default_range=6.28, illumination=None):
    """
    Definition for producing a pattern for a phase only Spatial Light Modulator (SLM) using a given field.

    Parameters
    ----------
    hologram           : np.complex64
                         Input holographic field.
    slm_range          : float
                         Range of the phase only SLM in radians for a working wavelength (i.e. two pi). See odak.wave.adjust_phase_only_slm_range() for more.
    filename           : str
                         Optional variable, if provided the patterns will be save to given location.
    bits               : int
                         Quantization bits.
    default_range      : float 
                         Default range of phase only SLM.
    illumination       : np.ndarray
                         Spatial illumination distribution.

    Returns
    -------
    pattern            : np.complex64
                         Adjusted phase only pattern.
    hologram_digital   : np.int
                         Digital representation of the hologram.
    """
    #hologram_phase   = calculate_phase(hologram) % default_range
    hologram_phase = calculate_phase(hologram)
    hologram_phase = hologram_phase % slm_range
    hologram_phase /= slm_range
    hologram_phase *= 2**bits
    hologram_phase = hologram_phase.astype(np.int32)
    hologram_digital = np.copy(hologram_phase)
    if type(filename) != type(None):
        save_image(
            filename,
            hologram_phase,
            cmin=0,
            cmax=2**bits
        )
    hologram_phase = hologram_phase.astype(np.float64)
    hologram_phase *= slm_range/2**bits
    if type(illumination) == type(None):
        A = 1.
    else:
        A = illumination
    return A*np.cos(hologram_phase)+A*1j*np.sin(hologram_phase), hologram_digital

propagate_beam(field, k, distance, dx, wavelength, propagation_type='IR Fresnel')

Definitions for Fresnel Impulse Response (IR), Angular Spectrum (AS), Bandlimited Angular Spectrum (BAS), Fresnel Transfer Function (TF), Fraunhofer diffraction in accordence with "Computational Fourier Optics" by David Vuelz. For more on Bandlimited Fresnel impulse response also known as Bandlimited Angular Spectrum method see "Band-limited Angular Spectrum Method for Numerical Simulation of Free-Space Propagation in Far and Near Fields".

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.
  • propagation_type (str, default: 'IR Fresnel' ) – 
               Type of the propagation (IR Fresnel, Angular Spectrum, Bandlimited Angular Spectrum, TR Fresnel, Fraunhofer).

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
def propagate_beam(field, k, distance, dx, wavelength, propagation_type='IR Fresnel'):
    """
    Definitions for Fresnel Impulse Response (IR), Angular Spectrum (AS), Bandlimited Angular Spectrum (BAS), Fresnel Transfer Function (TF), Fraunhofer diffraction in accordence with "Computational Fourier Optics" by David Vuelz. For more on Bandlimited Fresnel impulse response also known as Bandlimited Angular Spectrum method see "Band-limited Angular Spectrum Method for Numerical Simulation of Free-Space Propagation in Far and Near Fields".

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.
    propagation_type : str
                       Type of the propagation (IR Fresnel, Angular Spectrum, Bandlimited Angular Spectrum, TR Fresnel, Fraunhofer).

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    if propagation_type == 'Rayleigh-Sommerfeld':
        result = rayleigh_sommerfeld(field, k, distance, dx, wavelength)
    elif propagation_type == 'Angular Spectrum':
        result = angular_spectrum(field, k, distance, dx, wavelength)
    elif propagation_type == 'Impulse Response Fresnel':
        result = impulse_response_fresnel(field, k, distance, dx, wavelength)
    elif propagation_type == 'Bandlimited Angular Spectrum':
        result = band_limited_angular_spectrum(
            field, k, distance, dx, wavelength)
    elif propagation_type == 'Bandextended Angular Spectrum':
        result = band_extended_angular_spectrum(
            field, k, distance, dx, wavelength)
    elif propagation_type == 'Adaptive Sampling Angular Spectrum':
        result = adaptive_sampling_angular_spectrum(
            field, k, distance, dx, wavelength)
    elif propagation_type == 'Transfer Function Fresnel':
        result = transfer_function_fresnel(field, k, distance, dx, wavelength)
    elif propagation_type == 'Fraunhofer':
        result = fraunhofer(field, k, distance, dx, wavelength)
    elif propagation_type == 'Fraunhofer Inverse':
        result = fraunhofer_inverse(field, k, distance, dx, wavelength)
    else:
        raise Exception("Unknown propagation type selected.")
    return result

propagate_field(points0, points1, field0, wave_number, direction=1)

Definition to propagate a field from points to an another points in space: propagate a given array of spherical sources to given set of points in space.

Parameters:

  • points0 – 
            Start points (i.e. odak.tools.grid_sample).
  • points1 – 
            End points (ie. odak.tools.grid_sample).
  • field0 – 
            Field for given starting points.
  • wave_number – 
            Wave number of a wave, see odak.wave.wavenumber for more.
  • direction – 
            For propagating in forward direction set as 1, otherwise -1.

Returns:

  • field1 ( ndarray ) –

    Field for given end points.

Source code in odak/wave/vector.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
def propagate_field(points0, points1, field0, wave_number, direction=1):
    """
    Definition to propagate a field from points to an another points in space: propagate a given array of spherical sources to given set of points in space.

    Parameters
    ----------
    points0       : ndarray
                    Start points (i.e. odak.tools.grid_sample).
    points1       : ndarray
                    End points (ie. odak.tools.grid_sample).
    field0        : ndarray
                    Field for given starting points.
    wave_number   : float
                    Wave number of a wave, see odak.wave.wavenumber for more.
    direction     : float
                    For propagating in forward direction set as 1, otherwise -1.

    Returns
    -------
    field1        : ndarray
                    Field for given end points.
    """
    field1 = np.zeros(points1.shape[0], dtype=np.complex64)
    for point_id in range(points0.shape[0]):
        point = points0[point_id]
        distances = distance_between_two_points(
            point,
            points1
        )
        field1 += electric_field_per_plane_wave(
            calculate_amplitude(field0[point_id]),
            distances*direction,
            wave_number,
            phase=calculate_phase(field0[point_id])
        )
    return field1

propagate_plane_waves(field, opd, k, w=0, t=0)

Definition to propagate a field representing a plane wave at a particular distance and time.

Parameters:

  • field – 
           Complex field.
  • opd – 
           Optical path difference in mm.
  • k – 
           Wave number of a wave, see odak.wave.parameters.wavenumber for more.
  • w – 
           Rotation speed of a wave, see odak.wave.parameters.rotationspeed for more.
  • t – 
           Time in seconds.

Returns:

  • new_field ( complex ) –

    A complex number that provides the resultant field in the complex form A*e^(j(wt+phi)).

Source code in odak/wave/vector.py 
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def propagate_plane_waves(field, opd, k, w=0, t=0):
    """
    Definition to propagate a field representing a plane wave at a particular distance and time.

    Parameters
    ----------
    field        : complex
                   Complex field.
    opd          : float
                   Optical path difference in mm.
    k            : float
                   Wave number of a wave, see odak.wave.parameters.wavenumber for more.
    w            : float
                   Rotation speed of a wave, see odak.wave.parameters.rotationspeed for more.
    t            : float
                   Time in seconds.

    Returns
    -------
    new_field     : complex
                    A complex number that provides the resultant field in the complex form A*e^(j(wt+phi)).
    """
    new_field = field*np.exp(1j*(-w*t+opd*k))/opd**2
    return new_field

quadratic_phase_function(nx, ny, k, focal=0.4, dx=0.001, offset=[0, 0])

A definition to generate 2D quadratic phase function, which is typically use to represent lenses.

Parameters:

  • nx – 
         Size of the output along X.
  • ny – 
         Size of the output along Y.
  • k – 
         See odak.wave.wavenumber for more.
  • focal – 
         Focal length of the quadratic phase function.
  • dx – 
         Pixel pitch.
  • offset – 
         Deviation from the center along X and Y axes.

Returns:

  • function ( ndarray ) –

    Generated quadratic phase function.

Source code in odak/wave/lens.py 
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
def quadratic_phase_function(nx, ny, k, focal=0.4, dx=0.001, offset=[0, 0]):
    """ 
    A definition to generate 2D quadratic phase function, which is typically use to represent lenses.

    Parameters
    ----------
    nx         : int
                 Size of the output along X.
    ny         : int
                 Size of the output along Y.
    k          : odak.wave.wavenumber
                 See odak.wave.wavenumber for more.
    focal      : float
                 Focal length of the quadratic phase function.
    dx         : float
                 Pixel pitch.
    offset     : list
                 Deviation from the center along X and Y axes.

    Returns
    -------
    function   : ndarray
                 Generated quadratic phase function.
    """
    size = [nx, ny]
    x = np.linspace(-size[0]*dx/2, size[0]*dx/2, size[0])-offset[1]*dx
    y = np.linspace(-size[1]*dx/2, size[1]*dx/2, size[1])-offset[0]*dx
    X, Y = np.meshgrid(x, y)
    Z = X**2+Y**2
    qwf = np.exp(1j*k*0.5*np.sin(Z/focal))
    return qwf

rayleigh_resolution(diameter, focal=None, wavelength=0.0005)

Definition to calculate rayleigh resolution limit of a lens with a certain focal length and an aperture. Lens is assumed to be focusing a plane wave at a focal distance.

Parameter

diameter : float Diameter of a lens. focal : float Focal length of a lens, when focal length is provided, spatial resolution is provided at the focal plane. When focal length isn't provided angular resolution is provided. wavelength : float Wavelength of light.

Returns:

  • resolution ( float ) –

    Resolvable angular or spatial spot size, see focal in parameters to know what to expect.

Source code in odak/wave/__init__.py 
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
def rayleigh_resolution(diameter, focal=None, wavelength=0.0005):
    """
    Definition to calculate rayleigh resolution limit of a lens with a certain focal length and an aperture. Lens is assumed to be focusing a plane wave at a focal distance.

    Parameter
    ---------
    diameter    : float
                  Diameter of a lens.
    focal       : float
                  Focal length of a lens, when focal length is provided, spatial resolution is provided at the focal plane. When focal length isn't provided angular resolution is provided.
    wavelength  : float
                  Wavelength of light.

    Returns
    --------
    resolution  : float
                  Resolvable angular or spatial spot size, see focal in parameters to know what to expect.

    """
    resolution = 1.22*wavelength/diameter
    if type(focal) != type(None):
        resolution *= focal
    return resolution

rayleigh_sommerfeld(field, k, distance, dx, wavelength)

Definition to compute beam propagation using Rayleigh-Sommerfeld's diffraction formula (Huygens-Fresnel Principle). For more see Section 3.5.2 in Goodman, Joseph W. Introduction to Fourier optics. Roberts and Company Publishers, 2005.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
def rayleigh_sommerfeld(field, k, distance, dx, wavelength):
    """
    Definition to compute beam propagation using Rayleigh-Sommerfeld's diffraction formula (Huygens-Fresnel Principle). For more see Section 3.5.2 in Goodman, Joseph W. Introduction to Fourier optics. Roberts and Company Publishers, 2005.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    nv, nu = field.shape
    x = np.linspace(-nv * dx / 2, nv * dx / 2, nv)
    y = np.linspace(-nu * dx / 2, nu * dx / 2, nu)
    X, Y = np.meshgrid(x, y)
    Z = X ** 2 + Y ** 2
    result = np.zeros(field.shape, dtype=np.complex64)
    direction = int(distance/np.abs(distance))
    for i in range(nu):
        for j in range(nv):
            if field[i, j] != 0:
                r01 = np.sqrt(distance ** 2 + (X - X[i, j]) ** 2 + (Y - Y[i, j]) ** 2) * direction
                cosnr01 = np.cos(distance / r01)
                result += field[i, j] * np.exp(1j * k * r01) / r01 * cosnr01
    result *= 1. / (1j * wavelength)
    return result

rotationspeed(wavelength, c=3 * 10 ** 11)

Definition for calculating rotation speed of a wave (w in A*e^(j(wt+phi))).

Parameters:

  • wavelength – 
           Wavelength of a wave in mm.
  • c – 
           Speed of wave in mm/seconds. Default is the speed of light in the void!

Returns:

  • w ( float ) –

    Rotation speed.

Source code in odak/wave/__init__.py 
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
def rotationspeed(wavelength, c=3*10**11):
    """
    Definition for calculating rotation speed of a wave (w in A*e^(j(wt+phi))).

    Parameters
    ----------
    wavelength   : float
                   Wavelength of a wave in mm.
    c            : float
                   Speed of wave in mm/seconds. Default is the speed of light in the void!

    Returns
    -------
    w            : float
                   Rotation speed.

    """
    f = c*wavelength
    w = 2*np.pi*f
    return w

set_amplitude(field, amplitude)

Definition to keep phase as is and change the amplitude of a given field.

Parameters:

  • field – 
           Complex field.
  • amplitude – 
           Amplitudes.

Returns:

  • new_field ( complex64 ) –

    Complex field.

Source code in odak/wave/__init__.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
def set_amplitude(field, amplitude):
    """
    Definition to keep phase as is and change the amplitude of a given field.

    Parameters
    ----------
    field        : np.complex64
                   Complex field.
    amplitude    : np.array or np.complex64
                   Amplitudes.

    Returns
    -------
    new_field    : np.complex64
                   Complex field.
    """
    amplitude = calculate_amplitude(amplitude)
    phase = calculate_phase(field)
    new_field = amplitude*np.cos(phase)+1j*amplitude*np.sin(phase)
    return new_field

transfer_function_fresnel(field, k, distance, dx, wavelength)

A definition to calculate convolution based Fresnel approximation for beam propagation.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
def transfer_function_fresnel(field, k, distance, dx, wavelength):
    """
    A definition to calculate convolution based Fresnel approximation for beam propagation.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).

    """
    nv, nu = field.shape
    L = nu*dx
    fx = np.linspace(-1. / 2. /dx, 1. /2. /dx, nu)
    fy = np.linspace(-1. / 2. /dx, 1. /2. /dx, nv)
    FX, FY = np.meshgrid(fx, fy)
    H = np.exp(-1j * distance * (k - np.pi * wavelength * (FX**2 + FY**2) ))
    U1 = np.fft.fft2(np.fft.fftshift(field)) * ((1/L)**2)
    U2 = np.fft.fftshift(H)*U1
    result = np.fft.ifftshift(np.fft.ifft2(U2)) / ((1/L)**2)
    return result

wavenumber(wavelength)

Definition for calculating the wavenumber of a plane wave.

Parameters:

  • wavelength – 
           Wavelength of a wave in mm.

Returns:

  • k ( float ) –

    Wave number for a given wavelength.

Source code in odak/wave/__init__.py 
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
def wavenumber(wavelength):
    """
    Definition for calculating the wavenumber of a plane wave.

    Parameters
    ----------
    wavelength   : float
                   Wavelength of a wave in mm.

    Returns
    -------
    k            : float
                   Wave number for a given wavelength.
    """
    k = 2*np.pi/wavelength
    return k

zero_pad(field, size=None, method='center')

Definition to zero pad a MxN array to 2Mx2N array.

Parameters:

  • field – 
                Input field MxN array.
  • size – 
                Size to be zeropadded.
  • method – 
                Zeropad either by placing the content to center or to the left.

Returns:

  • field_zero_padded ( ndarray ) –

    Zeropadded version of the input field.

Source code in odak/tools/matrix.py 
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
def zero_pad(field, size=None, method='center'):
    """
    Definition to zero pad a MxN array to 2Mx2N array.

    Parameters
    ----------
    field             : ndarray
                        Input field MxN array.
    size              : list
                        Size to be zeropadded.
    method            : str
                        Zeropad either by placing the content to center or to the left.

    Returns
    ----------
    field_zero_padded : ndarray
                        Zeropadded version of the input field.
    """
    if type(size) == type(None):
        hx = int(np.ceil(field.shape[0])/2)
        hy = int(np.ceil(field.shape[1])/2)
    else:
        hx = int(np.ceil((size[0]-field.shape[0])/2))
        hy = int(np.ceil((size[1]-field.shape[1])/2))
    if method == 'center':
        field_zero_padded = np.pad(
            field, ([hx, hx], [hy, hy]), constant_values=(0, 0))
    elif method == 'left aligned':
        field_zero_padded = np.pad(
            field, ([0, 2*hx], [0, 2*hy]), constant_values=(0, 0))
    if type(size) != type(None):
        field_zero_padded = field_zero_padded[0:size[0], 0:size[1]]
    return field_zero_padded

adaptive_sampling_angular_spectrum(field, k, distance, dx, wavelength)

A definition to calculate adaptive sampling angular spectrum based beam propagation. For more Zhang, Wenhui, Hao Zhang, and Guofan Jin. "Adaptive-sampling angular spectrum method with full utilization of space-bandwidth product." Optics Letters 45.16 (2020): 4416-4419.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
def adaptive_sampling_angular_spectrum(field, k, distance, dx, wavelength):
    """
    A definition to calculate adaptive sampling angular spectrum based beam propagation. For more Zhang, Wenhui, Hao Zhang, and Guofan Jin. "Adaptive-sampling angular spectrum method with full utilization of space-bandwidth product." Optics Letters 45.16 (2020): 4416-4419.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    iflag = -1
    eps = 10**(-12)
    nv, nu = field.shape
    l = nu*dx
    x = np.linspace(-l/2, l/2, nu)
    y = np.linspace(-l/2, l/2, nv)
    X, Y = np.meshgrid(x, y)
    fx = np.linspace(-1./2./dx, 1./2./dx, nu)
    fy = np.linspace(-1./2./dx, 1./2./dx, nv)
    FX, FY = np.meshgrid(fx, fy)
    forig = 1./2./dx
    fc2 = 1./2*(nu/wavelength/np.abs(distance))**0.5
    ss = np.abs(fc2)/forig
    zc = nu*dx**2/wavelength
    K = nu/2/np.amax(np.abs(fx))
    m = 2
    nnu2 = m*nu
    nnv2 = m*nv
    fxn = np.linspace(-1./2./dx, 1./2./dx, nnu2)
    fyn = np.linspace(-1./2./dx, 1./2./dx, nnv2)
    if np.abs(distance) > zc*2:
        fxn = fxn*ss
        fyn = fyn*ss
    FXN, FYN = np.meshgrid(fxn, fyn)
    Hn = np.exp(1j*k*distance*(1-(FXN*wavelength)**2-(FYN*wavelength)**2)**0.5)
    FX = FXN/np.amax(FXN)*np.pi
    FY = FYN/np.amax(FYN)*np.pi
    t_2 = nufft2(field, FX*ss, FY*ss, size=[nnv2, nnu2], sign=iflag, eps=eps)
    FX = FX/np.amax(FX)*np.pi
    FY = FY/np.amax(FY)*np.pi
    result = nuifft2(Hn*t_2, FX*ss, FY*ss, size=[nv, nu], sign=-iflag, eps=eps)
    return result

angular_spectrum(field, k, distance, dx, wavelength)

A definition to calculate angular spectrum based beam propagation.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
def angular_spectrum(field, k, distance, dx, wavelength):
    """
    A definition to calculate angular spectrum based beam propagation.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    nv, nu = field.shape
    fx = np.linspace(-1. / 2. /dx, 1. /2. /dx, nu)
    fy = np.linspace(-1. / 2. /dx, 1. /2. /dx, nv)
    FX, FY = np.meshgrid(fx, fy)
    H = np.exp(1j * k * distance * (1 - (FX * wavelength) ** 2 - (FY * wavelength) ** 2) ** 0.5)
    U1 = np.fft.fftshift(np.fft.fft2(field))
    U2 = H * U1
    result = np.fft.ifft2(np.fft.ifftshift(U2))
    return result

band_extended_angular_spectrum(field, k, distance, dx, wavelength)

A definition to calculate bandextended angular spectrum based beam propagation. For more Zhang, Wenhui, Hao Zhang, and Guofan Jin. "Band-extended angular spectrum method for accurate diffraction calculation in a wide propagation range." Optics Letters 45.6 (2020): 1543-1546.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
def band_extended_angular_spectrum(field, k, distance, dx, wavelength):
    """
    A definition to calculate bandextended angular spectrum based beam propagation. For more Zhang, Wenhui, Hao Zhang, and Guofan Jin. "Band-extended angular spectrum method for accurate diffraction calculation in a wide propagation range." Optics Letters 45.6 (2020): 1543-1546.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    iflag = -1
    eps = 10**(-12)
    nv, nu = field.shape
    l = nu*dx
    x = np.linspace(-l/2, l/2, nu)
    y = np.linspace(-l/2, l/2, nv)
    X, Y = np.meshgrid(x, y)
    Z = X**2+Y**2
    fx = np.linspace(-1./2./dx, 1./2./dx, nu)
    fy = np.linspace(-1./2./dx, 1./2./dx, nv)
    FX, FY = np.meshgrid(fx, fy)
    K = nu/2/np.amax(fx)
    fcn = 1./2*(nu/wavelength/np.abs(distance))**0.5
    ss = np.abs(fcn)/np.amax(np.abs(fx))
    zc = nu*dx**2/wavelength
    if np.abs(distance) < zc:
        fxn = fx
        fyn = fy
    else:
        fxn = fx*ss
        fyn = fy*ss
    FXN, FYN = np.meshgrid(fxn, fyn)
    Hn = np.exp(1j*k*distance*(1-(FXN*wavelength)**2-(FYN*wavelength)**2)**0.5)
    X = X/np.amax(X)*np.pi
    Y = Y/np.amax(Y)*np.pi
    t_asmNUFT = nufft2(field, X*ss, Y*ss, sign=iflag, eps=eps)
    result = nuifft2(Hn*t_asmNUFT, X*ss, Y*ss, sign=-iflag, eps=eps)
    return result

band_limited_angular_spectrum(field, k, distance, dx, wavelength)

A definition to calculate bandlimited angular spectrum based beam propagation. For more Matsushima, Kyoji, and Tomoyoshi Shimobaba. "Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields." Optics express 17.22 (2009): 19662-19673.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
def band_limited_angular_spectrum(field, k, distance, dx, wavelength):
    """
    A definition to calculate bandlimited angular spectrum based beam propagation. For more Matsushima, Kyoji, and Tomoyoshi Shimobaba. "Band-limited angular spectrum method for numerical simulation of free-space propagation in far and near fields." Optics express 17.22 (2009): 19662-19673.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    nv, nu = field.shape
    fx = np.linspace(-1. / 2. /dx, 1. /2. /dx, nu)
    fy = np.linspace(-1. / 2. /dx, 1. /2. /dx, nv)
    FX, FY = np.meshgrid(fx, fy)
    H_exp = np.exp(1j * k * distance * (1 - (FX * wavelength) ** 2 - (FY * wavelength) ** 2) ** 0.5)

    x = dx * float(nu)
    y = dx * float(nv)
    fx_max = 1 / np.sqrt((2 * distance * (1 / x))**2 + 1) / wavelength
    fy_max = 1 / np.sqrt((2 * distance * (1 / y))**2 + 1) / wavelength
    H_filter = ((np.abs(FX) < fx_max) & (np.abs(FY) < fy_max))
    H = H_filter * H_exp

    U1 = np.fft.fftshift(np.fft.fft2(field))
    U2 = H * U1
    result = np.fft.ifft2(np.fft.ifftshift(U2))
    return result

fraunhofer(field, k, distance, dx, wavelength)

A definition to calculate Fraunhofer based beam propagation.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
def fraunhofer(field, k, distance, dx, wavelength):
    """
    A definition to calculate Fraunhofer based beam propagation.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    nv, nu = field.shape
    l = nu*dx
    l2 = wavelength*distance/dx
    dx2 = wavelength*distance/l
    fx = np.linspace(-l2/2., l2/2., nu)
    fy = np.linspace(-l2/2., l2/2., nv)
    FX, FY = np.meshgrid(fx, fy)
    FZ = FX**2+FY**2
    c = np.exp(1j*k*distance)/(1j*wavelength*distance) * \
        np.exp(1j*k/(2*distance)*FZ)
    result = c*np.fft.ifftshift(np.fft.fft2(np.fft.fftshift(field)))*dx**2
    return result

fraunhofer_equal_size_adjust(field, distance, dx, wavelength)

A definition to match the physical size of the original field with the propagated field.

Parameters:

  • field – 
               Complex field (MxN).
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • new_field ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
def fraunhofer_equal_size_adjust(field, distance, dx, wavelength):
    """
    A definition to match the physical size of the original field with the propagated field.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    new_field        : np.complex
                       Final complex field (MxN).
    """
    nv, nu = field.shape
    l1 = nu*dx
    l2 = wavelength*distance/dx
    m = l1/l2
    px = int(m*nu)
    py = int(m*nv)
    nx = int(field.shape[0]/2-px/2)
    ny = int(field.shape[1]/2-py/2)
    new_field = np.copy(field[nx:nx+px, ny:ny+py])
    return new_field

fraunhofer_inverse(field, k, distance, dx, wavelength)

A definition to calculate Inverse Fraunhofer based beam propagation.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
def fraunhofer_inverse(field, k, distance, dx, wavelength):
    """
    A definition to calculate Inverse Fraunhofer based beam propagation.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    distance = np.abs(distance)
    nv, nu = field.shape
    l = nu*dx
    l2 = wavelength*distance/dx
    dx2 = wavelength*distance/l
    fx = np.linspace(-l2/2., l2/2., nu)
    fy = np.linspace(-l2/2., l2/2., nv)
    FX, FY = np.meshgrid(fx, fy)
    FZ = FX**2+FY**2
    c = np.exp(1j*k*distance)/(1j*wavelength*distance) * \
        np.exp(1j*k/(2*distance)*FZ)
    result = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(field/dx**2/c)))
    return result

gerchberg_saxton(field, n_iterations, distance, dx, wavelength, slm_range=6.28, propagation_type='IR Fresnel', initial_phase=None)

Definition to compute a hologram using an iterative method called Gerchberg-Saxton phase retrieval algorithm. For more on the method, see: Gerchberg, Ralph W. "A practical algorithm for the determination of phase from image and diffraction plane pictures." Optik 35 (1972): 237-246.

Parameters:

  • field – 
               Complex field (MxN).
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.
  • slm_range – 
               Typically this is equal to two pi. See odak.wave.adjust_phase_only_slm_range() for more.
  • propagation_type (str, default: 'IR Fresnel' ) – 
               Type of the propagation (IR Fresnel, TR Fresnel, Fraunhofer).
  • initial_phase – 
               Phase to be added to the initial value.

Returns:

  • hologram ( complex ) –

    Calculated complex hologram.

  • reconstruction ( complex ) –

    Calculated reconstruction using calculated hologram.

Source code in odak/wave/classical.py 
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
def gerchberg_saxton(field, n_iterations, distance, dx, wavelength, slm_range=6.28, propagation_type='IR Fresnel', initial_phase=None):
    """
    Definition to compute a hologram using an iterative method called Gerchberg-Saxton phase retrieval algorithm. For more on the method, see: Gerchberg, Ralph W. "A practical algorithm for the determination of phase from image and diffraction plane pictures." Optik 35 (1972): 237-246.

    Parameters
    ----------
    field            : np.complex64
                       Complex field (MxN).
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.
    slm_range        : float
                       Typically this is equal to two pi. See odak.wave.adjust_phase_only_slm_range() for more.
    propagation_type : str
                       Type of the propagation (IR Fresnel, TR Fresnel, Fraunhofer).
    initial_phase    : np.complex64
                       Phase to be added to the initial value.

    Returns
    -------
    hologram         : np.complex
                       Calculated complex hologram.
    reconstruction   : np.complex
                       Calculated reconstruction using calculated hologram. 
    """
    k = wavenumber(wavelength)
    target = calculate_amplitude(field)
    hologram = generate_complex_field(np.ones(field.shape), 0)
    hologram = zero_pad(hologram)
    if type(initial_phase) == type(None):
        hologram = add_random_phase(hologram)
    else:
        initial_phase = zero_pad(initial_phase)
        hologram = add_phase(hologram, initial_phase)
    center = [int(hologram.shape[0]/2.), int(hologram.shape[1]/2.)]
    orig_shape = [int(field.shape[0]/2.), int(field.shape[1]/2.)]
    for i in tqdm(range(n_iterations), leave=False):
        reconstruction = propagate_beam(
            hologram, k, distance, dx, wavelength, propagation_type)
        new_target = calculate_amplitude(reconstruction)
        new_target[
            center[0]-orig_shape[0]:center[0]+orig_shape[0],
            center[1]-orig_shape[1]:center[1]+orig_shape[1]
        ] = target
        reconstruction = generate_complex_field(
            new_target, calculate_phase(reconstruction))
        hologram = propagate_beam(
            reconstruction, k, -distance, dx, wavelength, propagation_type)
        hologram = generate_complex_field(1, calculate_phase(hologram))
        hologram = hologram[
            center[0]-orig_shape[0]:center[0]+orig_shape[0],
            center[1]-orig_shape[1]:center[1]+orig_shape[1],
        ]
        hologram = zero_pad(hologram)
    reconstruction = propagate_beam(
        hologram, k, distance, dx, wavelength, propagation_type)
    hologram = hologram[
        center[0]-orig_shape[0]:center[0]+orig_shape[0],
        center[1]-orig_shape[1]:center[1]+orig_shape[1]
    ]
    reconstruction = reconstruction[
        center[0]-orig_shape[0]:center[0]+orig_shape[0],
        center[1]-orig_shape[1]:center[1]+orig_shape[1]
    ]
    return hologram, reconstruction

gerchberg_saxton_3d(fields, n_iterations, distances, dx, wavelength, slm_range=6.28, propagation_type='IR Fresnel', initial_phase=None, target_type='no constraint', coefficients=None)

Definition to compute a multi plane hologram using an iterative method called Gerchberg-Saxton phase retrieval algorithm. For more on the method, see: Zhou, Pengcheng, et al. "30.4: Multi‐plane holographic display with a uniform 3D Gerchberg‐Saxton algorithm." SID Symposium Digest of Technical Papers. Vol. 46. No. 1. 2015.

Parameters:

  • fields – 
               Complex fields (MxN).
  • distances – 
               Propagation distances.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.
  • slm_range – 
               Typically this is equal to two pi. See odak.wave.adjust_phase_only_slm_range() for more.
  • propagation_type (str, default: 'IR Fresnel' ) – 
               Type of the propagation (IR Fresnel, TR Fresnel, Fraunhofer).
  • initial_phase – 
               Phase to be added to the initial value.
  • target_type – 
               Target type. `No constraint` targets the input target as is. `Double constraint` follows the idea in this paper, which claims to suppress speckle: Chang, Chenliang, et al. "Speckle-suppressed phase-only holographic three-dimensional display based on double-constraint Gerchberg–Saxton algorithm." Applied optics 54.23 (2015): 6994-7001.

Returns:

  • hologram ( complex ) –

    Calculated complex hologram.

Source code in odak/wave/classical.py 
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
def gerchberg_saxton_3d(fields, n_iterations, distances, dx, wavelength, slm_range=6.28, propagation_type='IR Fresnel', initial_phase=None, target_type='no constraint', coefficients=None):
    """
    Definition to compute a multi plane hologram using an iterative method called Gerchberg-Saxton phase retrieval algorithm. For more on the method, see: Zhou, Pengcheng, et al. "30.4: Multi‐plane holographic display with a uniform 3D Gerchberg‐Saxton algorithm." SID Symposium Digest of Technical Papers. Vol. 46. No. 1. 2015.

    Parameters
    ----------
    fields           : np.complex64
                       Complex fields (MxN).
    distances        : list
                       Propagation distances.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.
    slm_range        : float
                       Typically this is equal to two pi. See odak.wave.adjust_phase_only_slm_range() for more.
    propagation_type : str
                       Type of the propagation (IR Fresnel, TR Fresnel, Fraunhofer).
    initial_phase    : np.complex64
                       Phase to be added to the initial value.
    target_type      : str
                       Target type. `No constraint` targets the input target as is. `Double constraint` follows the idea in this paper, which claims to suppress speckle: Chang, Chenliang, et al. "Speckle-suppressed phase-only holographic three-dimensional display based on double-constraint Gerchberg–Saxton algorithm." Applied optics 54.23 (2015): 6994-7001. 

    Returns
    -------
    hologram         : np.complex
                       Calculated complex hologram.
    """
    k = wavenumber(wavelength)
    targets = calculate_amplitude(np.asarray(fields)).astype(np.float64)
    hologram = generate_complex_field(np.ones(targets[0].shape), 0)
    hologram = zero_pad(hologram)
    if type(initial_phase) == type(None):
        hologram = add_random_phase(hologram)
    else:
        initial_phase = zero_pad(initial_phase)
        hologram = add_phase(hologram, initial_phase)
    center = [int(hologram.shape[0]/2.), int(hologram.shape[1]/2.)]
    orig_shape = [int(fields[0].shape[0]/2.), int(fields[0].shape[1]/2.)]
    holograms = np.zeros(
        (len(distances), hologram.shape[0], hologram.shape[1]), dtype=np.complex64)
    for i in tqdm(range(n_iterations), leave=False):
        for distance_id in tqdm(range(len(distances)), leave=False):
            distance = distances[distance_id]
            reconstruction = propagate_beam(
                hologram, k, distance, dx, wavelength, propagation_type)
            if target_type == 'double constraint':
                if type(coefficients) == type(None):
                    raise Exception(
                        "Provide coeeficients of alpha,beta and gamma for double constraint.")
                alpha = coefficients[0]
                beta = coefficients[1]
                gamma = coefficients[2]
                target_current = 2*alpha * \
                    np.copy(targets[distance_id])-beta * \
                    calculate_amplitude(reconstruction)
                target_current[target_current == 0] = gamma * \
                    np.abs(reconstruction[target_current == 0])
            elif target_type == 'no constraint':
                target_current = np.abs(targets[distance_id])
            new_target = calculate_amplitude(reconstruction)
            new_target[
                center[0]-orig_shape[0]:center[0]+orig_shape[0],
                center[1]-orig_shape[1]:center[1]+orig_shape[1]
            ] = target_current
            reconstruction = generate_complex_field(
                new_target, calculate_phase(reconstruction))
            hologram_layer = propagate_beam(
                reconstruction, k, -distance, dx, wavelength, propagation_type)
            hologram_layer = generate_complex_field(
                1., calculate_phase(hologram_layer))
            hologram_layer = hologram_layer[
                center[0]-orig_shape[0]:center[0]+orig_shape[0],
                center[1]-orig_shape[1]:center[1]+orig_shape[1]
            ]
            hologram_layer = zero_pad(hologram_layer)
            holograms[distance_id] = hologram_layer
        hologram = np.sum(holograms, axis=0)
    hologram = hologram[
        center[0]-orig_shape[0]:center[0]+orig_shape[0],
        center[1]-orig_shape[1]:center[1]+orig_shape[1]
    ]
    return hologram

impulse_response_fresnel(field, k, distance, dx, wavelength)

A definition to calculate impulse response based Fresnel approximation for beam propagation.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
def impulse_response_fresnel(field, k, distance, dx, wavelength):
    """
    A definition to calculate impulse response based Fresnel approximation for beam propagation.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).

    """
    nv, nu = field.shape
    x = np.linspace(-nu / 2 * dx, nu / 2 * dx, nu)
    y = np.linspace(-nv / 2 * dx, nv / 2 * dx, nv)
    X, Y = np.meshgrid(x, y)
    h = 1. / (1j * wavelength * distance) * np.exp(1j * k / (2 * distance) * (X ** 2 + Y ** 2))
    H = np.fft.fft2(np.fft.fftshift(h))*dx**2
    U1 = np.fft.fft2(np.fft.fftshift(field))
    U2 = H * U1
    result = np.fft.ifftshift(np.fft.ifft2(U2)) /(dx**2)

    return result

propagate_beam(field, k, distance, dx, wavelength, propagation_type='IR Fresnel')

Definitions for Fresnel Impulse Response (IR), Angular Spectrum (AS), Bandlimited Angular Spectrum (BAS), Fresnel Transfer Function (TF), Fraunhofer diffraction in accordence with "Computational Fourier Optics" by David Vuelz. For more on Bandlimited Fresnel impulse response also known as Bandlimited Angular Spectrum method see "Band-limited Angular Spectrum Method for Numerical Simulation of Free-Space Propagation in Far and Near Fields".

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.
  • propagation_type (str, default: 'IR Fresnel' ) – 
               Type of the propagation (IR Fresnel, Angular Spectrum, Bandlimited Angular Spectrum, TR Fresnel, Fraunhofer).

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
def propagate_beam(field, k, distance, dx, wavelength, propagation_type='IR Fresnel'):
    """
    Definitions for Fresnel Impulse Response (IR), Angular Spectrum (AS), Bandlimited Angular Spectrum (BAS), Fresnel Transfer Function (TF), Fraunhofer diffraction in accordence with "Computational Fourier Optics" by David Vuelz. For more on Bandlimited Fresnel impulse response also known as Bandlimited Angular Spectrum method see "Band-limited Angular Spectrum Method for Numerical Simulation of Free-Space Propagation in Far and Near Fields".

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.
    propagation_type : str
                       Type of the propagation (IR Fresnel, Angular Spectrum, Bandlimited Angular Spectrum, TR Fresnel, Fraunhofer).

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    if propagation_type == 'Rayleigh-Sommerfeld':
        result = rayleigh_sommerfeld(field, k, distance, dx, wavelength)
    elif propagation_type == 'Angular Spectrum':
        result = angular_spectrum(field, k, distance, dx, wavelength)
    elif propagation_type == 'Impulse Response Fresnel':
        result = impulse_response_fresnel(field, k, distance, dx, wavelength)
    elif propagation_type == 'Bandlimited Angular Spectrum':
        result = band_limited_angular_spectrum(
            field, k, distance, dx, wavelength)
    elif propagation_type == 'Bandextended Angular Spectrum':
        result = band_extended_angular_spectrum(
            field, k, distance, dx, wavelength)
    elif propagation_type == 'Adaptive Sampling Angular Spectrum':
        result = adaptive_sampling_angular_spectrum(
            field, k, distance, dx, wavelength)
    elif propagation_type == 'Transfer Function Fresnel':
        result = transfer_function_fresnel(field, k, distance, dx, wavelength)
    elif propagation_type == 'Fraunhofer':
        result = fraunhofer(field, k, distance, dx, wavelength)
    elif propagation_type == 'Fraunhofer Inverse':
        result = fraunhofer_inverse(field, k, distance, dx, wavelength)
    else:
        raise Exception("Unknown propagation type selected.")
    return result

rayleigh_sommerfeld(field, k, distance, dx, wavelength)

Definition to compute beam propagation using Rayleigh-Sommerfeld's diffraction formula (Huygens-Fresnel Principle). For more see Section 3.5.2 in Goodman, Joseph W. Introduction to Fourier optics. Roberts and Company Publishers, 2005.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
def rayleigh_sommerfeld(field, k, distance, dx, wavelength):
    """
    Definition to compute beam propagation using Rayleigh-Sommerfeld's diffraction formula (Huygens-Fresnel Principle). For more see Section 3.5.2 in Goodman, Joseph W. Introduction to Fourier optics. Roberts and Company Publishers, 2005.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).
    """
    nv, nu = field.shape
    x = np.linspace(-nv * dx / 2, nv * dx / 2, nv)
    y = np.linspace(-nu * dx / 2, nu * dx / 2, nu)
    X, Y = np.meshgrid(x, y)
    Z = X ** 2 + Y ** 2
    result = np.zeros(field.shape, dtype=np.complex64)
    direction = int(distance/np.abs(distance))
    for i in range(nu):
        for j in range(nv):
            if field[i, j] != 0:
                r01 = np.sqrt(distance ** 2 + (X - X[i, j]) ** 2 + (Y - Y[i, j]) ** 2) * direction
                cosnr01 = np.cos(distance / r01)
                result += field[i, j] * np.exp(1j * k * r01) / r01 * cosnr01
    result *= 1. / (1j * wavelength)
    return result

transfer_function_fresnel(field, k, distance, dx, wavelength)

A definition to calculate convolution based Fresnel approximation for beam propagation.

Parameters:

  • field – 
               Complex field (MxN).
  • k – 
               Wave number of a wave, see odak.wave.wavenumber for more.
  • distance – 
               Propagation distance.
  • dx – 
               Size of one single pixel in the field grid (in meters).
  • wavelength – 
               Wavelength of the electric field.

Returns:

  • result ( complex ) –

    Final complex field (MxN).

Source code in odak/wave/classical.py 
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
def transfer_function_fresnel(field, k, distance, dx, wavelength):
    """
    A definition to calculate convolution based Fresnel approximation for beam propagation.

    Parameters
    ----------
    field            : np.complex
                       Complex field (MxN).
    k                : odak.wave.wavenumber
                       Wave number of a wave, see odak.wave.wavenumber for more.
    distance         : float
                       Propagation distance.
    dx               : float
                       Size of one single pixel in the field grid (in meters).
    wavelength       : float
                       Wavelength of the electric field.

    Returns
    -------
    result           : np.complex
                       Final complex field (MxN).

    """
    nv, nu = field.shape
    L = nu*dx
    fx = np.linspace(-1. / 2. /dx, 1. /2. /dx, nu)
    fy = np.linspace(-1. / 2. /dx, 1. /2. /dx, nv)
    FX, FY = np.meshgrid(fx, fy)
    H = np.exp(-1j * distance * (k - np.pi * wavelength * (FX**2 + FY**2) ))
    U1 = np.fft.fft2(np.fft.fftshift(field)) * ((1/L)**2)
    U2 = np.fft.fftshift(H)*U1
    result = np.fft.ifftshift(np.fft.ifft2(U2)) / ((1/L)**2)
    return result

double_convergence(nx, ny, k, r, dx)

A definition to generate initial phase for a Gerchberg-Saxton method. For more details consult Sun, Peng, et al. "Holographic near-eye display system based on double-convergence light Gerchberg-Saxton algorithm." Optics express 26.8 (2018): 10140-10151.

Parameters:

  • nx – 
         Size of the output along X.
  • ny – 
         Size of the output along Y.
  • k – 
         See odak.wave.wavenumber for more.
  • r – 
         The distance between location of a light source and an image plane.
  • dx – 
         Pixel pitch.

Returns:

  • function ( ndarray ) –

    Generated phase pattern for a Gerchberg-Saxton method.

Source code in odak/wave/lens.py 
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def double_convergence(nx, ny, k, r, dx):
    """
    A definition to generate initial phase for a Gerchberg-Saxton method. For more details consult Sun, Peng, et al. "Holographic near-eye display system based on double-convergence light Gerchberg-Saxton algorithm." Optics express 26.8 (2018): 10140-10151.

    Parameters
    ----------
    nx         : int
                 Size of the output along X.
    ny         : int
                 Size of the output along Y.
    k          : odak.wave.wavenumber
                 See odak.wave.wavenumber for more.
    r          : float
                 The distance between location of a light source and an image plane.
    dx         : float
                 Pixel pitch.

    Returns
    -------
    function   : ndarray
                 Generated phase pattern for a Gerchberg-Saxton method.
    """
    size = [ny, nx]
    x = np.linspace(-size[0]*dx/2, size[0]*dx/2, size[0])
    y = np.linspace(-size[1]*dx/2, size[1]*dx/2, size[1])
    X, Y = np.meshgrid(x, y)
    Z = X**2+Y**2
    w = np.exp(1j*k*Z/r)
    return w

linear_grating(nx, ny, every=2, add=3.14, axis='x')

A definition to generate a linear grating.

Parameters:

  • nx – 
         Size of the output along X.
  • ny – 
         Size of the output along Y.
  • every – 
         Add the add value at every given number.
  • add – 
         Angle to be added.
  • axis – 
         Axis eiter X,Y or both.

Returns:

  • field ( ndarray ) –

    Linear grating term.

Source code in odak/wave/lens.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
def linear_grating(nx, ny, every=2, add=3.14, axis='x'):
    """
    A definition to generate a linear grating.

    Parameters
    ----------
    nx         : int
                 Size of the output along X.
    ny         : int
                 Size of the output along Y.
    every      : int
                 Add the add value at every given number.
    add        : float
                 Angle to be added.
    axis       : string
                 Axis eiter X,Y or both.

    Returns
    -------
    field      : ndarray
                 Linear grating term.
    """
    grating = np.zeros((nx, ny), dtype=np.complex64)
    if axis == 'x':
        grating[::every, :] = np.exp(1j*add)
    if axis == 'y':
        grating[:, ::every] = np.exp(1j*add)
    if axis == 'xy':
        checker = np.indices((nx, ny)).sum(axis=0) % every
        checker += 1
        checker = checker % 2
        grating = np.exp(1j*checker*add)
    return grating

prism_phase_function(nx, ny, k, angle, dx=0.001, axis='x')

A definition to generate 2D phase function that represents a prism. See Goodman's Introduction to Fourier Optics book for more.

Parameters:

  • nx – 
         Size of the output along X.
  • ny – 
         Size of the output along Y.
  • k – 
         See odak.wave.wavenumber for more.
  • angle – 
         Tilt angle of the prism in degrees.
  • dx – 
         Pixel pitch.
  • axis – 
         Axis of the prism.

Returns:

  • prism ( ndarray ) –

    Generated phase function for a prism.

Source code in odak/wave/lens.py 
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
def prism_phase_function(nx, ny, k, angle, dx=0.001, axis='x'):
    """
    A definition to generate 2D phase function that represents a prism. See Goodman's Introduction to Fourier Optics book for more.

    Parameters
    ----------
    nx         : int
                 Size of the output along X.
    ny         : int
                 Size of the output along Y.
    k          : odak.wave.wavenumber
                 See odak.wave.wavenumber for more.
    angle      : float
                 Tilt angle of the prism in degrees.
    dx         : float
                 Pixel pitch.
    axis       : str
                 Axis of the prism.

    Returns
    -------
    prism      : ndarray
                 Generated phase function for a prism.
    """
    angle = np.radians(angle)
    size = [ny, nx]
    x = np.linspace(-size[0]*dx/2, size[0]*dx/2, size[0])
    y = np.linspace(-size[1]*dx/2, size[1]*dx/2, size[1])
    X, Y = np.meshgrid(x, y)
    if axis == 'y':
        prism = np.exp(-1j*k*np.sin(angle)*Y)
    elif axis == 'x':
        prism = np.exp(-1j*k*np.sin(angle)*X)
    return prism

quadratic_phase_function(nx, ny, k, focal=0.4, dx=0.001, offset=[0, 0])

A definition to generate 2D quadratic phase function, which is typically use to represent lenses.

Parameters:

  • nx – 
         Size of the output along X.
  • ny – 
         Size of the output along Y.
  • k – 
         See odak.wave.wavenumber for more.
  • focal – 
         Focal length of the quadratic phase function.
  • dx – 
         Pixel pitch.
  • offset – 
         Deviation from the center along X and Y axes.

Returns:

  • function ( ndarray ) –

    Generated quadratic phase function.

Source code in odak/wave/lens.py 
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
def quadratic_phase_function(nx, ny, k, focal=0.4, dx=0.001, offset=[0, 0]):
    """ 
    A definition to generate 2D quadratic phase function, which is typically use to represent lenses.

    Parameters
    ----------
    nx         : int
                 Size of the output along X.
    ny         : int
                 Size of the output along Y.
    k          : odak.wave.wavenumber
                 See odak.wave.wavenumber for more.
    focal      : float
                 Focal length of the quadratic phase function.
    dx         : float
                 Pixel pitch.
    offset     : list
                 Deviation from the center along X and Y axes.

    Returns
    -------
    function   : ndarray
                 Generated quadratic phase function.
    """
    size = [nx, ny]
    x = np.linspace(-size[0]*dx/2, size[0]*dx/2, size[0])-offset[1]*dx
    y = np.linspace(-size[1]*dx/2, size[1]*dx/2, size[1])-offset[0]*dx
    X, Y = np.meshgrid(x, y)
    Z = X**2+Y**2
    qwf = np.exp(1j*k*0.5*np.sin(Z/focal))
    return qwf

calculate_amplitude(field)

Definition to calculate amplitude of a single or multiple given electric field(s).

Parameters:

  • field – 
           Electric fields or an electric field.

Returns:

  • amplitude ( float ) –

    Amplitude or amplitudes of electric field(s).

Source code in odak/wave/utils.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
def calculate_amplitude(field):
    """ 
    Definition to calculate amplitude of a single or multiple given electric field(s).

    Parameters
    ----------
    field        : ndarray.complex or complex
                   Electric fields or an electric field.

    Returns
    -------
    amplitude    : float
                   Amplitude or amplitudes of electric field(s).
    """
    amplitude = np.abs(field)
    return amplitude

calculate_phase(field, deg=False)

Definition to calculate phase of a single or multiple given electric field(s).

Parameters:

  • field – 
           Electric fields or an electric field.
  • deg – 
           If set True, the angles will be returned in degrees.

Returns:

  • phase ( float ) –

    Phase or phases of electric field(s) in radians.

Source code in odak/wave/utils.py 
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def calculate_phase(field, deg=False):
    """ 
    Definition to calculate phase of a single or multiple given electric field(s).

    Parameters
    ----------
    field        : ndarray.complex or complex
                   Electric fields or an electric field.
    deg          : bool
                   If set True, the angles will be returned in degrees.

    Returns
    -------
    phase        : float
                   Phase or phases of electric field(s) in radians.
    """
    phase = np.angle(field)
    if deg == True:
        phase *= 180./np.pi
    return phase

electric_field_per_plane_wave(amplitude, opd, k, phase=0, w=0, t=0)

Definition to return state of a plane wave at a particular distance and time.

Parameters:

  • amplitude – 
           Amplitude of a wave.
  • opd – 
           Optical path difference in mm.
  • k – 
           Wave number of a wave, see odak.wave.parameters.wavenumber for more.
  • phase – 
           Initial phase of a wave.
  • w – 
           Rotation speed of a wave, see odak.wave.parameters.rotationspeed for more.
  • t – 
           Time in seconds.

Returns:

  • field ( complex ) –

    A complex number that provides the resultant field in the complex form A*e^(j(wt+phi)).

Source code in odak/wave/vector.py 
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
def electric_field_per_plane_wave(amplitude, opd, k, phase=0, w=0, t=0):
    """
    Definition to return state of a plane wave at a particular distance and time.

    Parameters
    ----------
    amplitude    : float
                   Amplitude of a wave.
    opd          : float
                   Optical path difference in mm.
    k            : float
                   Wave number of a wave, see odak.wave.parameters.wavenumber for more.
    phase        : float
                   Initial phase of a wave.
    w            : float
                   Rotation speed of a wave, see odak.wave.parameters.rotationspeed for more.
    t            : float
                   Time in seconds.

    Returns
    -------
    field        : complex
                   A complex number that provides the resultant field in the complex form A*e^(j(wt+phi)).
    """
    field = amplitude*np.exp(1j*(-w*t+opd*k+phase))/opd**2
    return field

propagate_field(points0, points1, field0, wave_number, direction=1)

Definition to propagate a field from points to an another points in space: propagate a given array of spherical sources to given set of points in space.

Parameters:

  • points0 – 
            Start points (i.e. odak.tools.grid_sample).
  • points1 – 
            End points (ie. odak.tools.grid_sample).
  • field0 – 
            Field for given starting points.
  • wave_number – 
            Wave number of a wave, see odak.wave.wavenumber for more.
  • direction – 
            For propagating in forward direction set as 1, otherwise -1.

Returns:

  • field1 ( ndarray ) –

    Field for given end points.

Source code in odak/wave/vector.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
def propagate_field(points0, points1, field0, wave_number, direction=1):
    """
    Definition to propagate a field from points to an another points in space: propagate a given array of spherical sources to given set of points in space.

    Parameters
    ----------
    points0       : ndarray
                    Start points (i.e. odak.tools.grid_sample).
    points1       : ndarray
                    End points (ie. odak.tools.grid_sample).
    field0        : ndarray
                    Field for given starting points.
    wave_number   : float
                    Wave number of a wave, see odak.wave.wavenumber for more.
    direction     : float
                    For propagating in forward direction set as 1, otherwise -1.

    Returns
    -------
    field1        : ndarray
                    Field for given end points.
    """
    field1 = np.zeros(points1.shape[0], dtype=np.complex64)
    for point_id in range(points0.shape[0]):
        point = points0[point_id]
        distances = distance_between_two_points(
            point,
            points1
        )
        field1 += electric_field_per_plane_wave(
            calculate_amplitude(field0[point_id]),
            distances*direction,
            wave_number,
            phase=calculate_phase(field0[point_id])
        )
    return field1

propagate_plane_waves(field, opd, k, w=0, t=0)

Definition to propagate a field representing a plane wave at a particular distance and time.

Parameters:

  • field – 
           Complex field.
  • opd – 
           Optical path difference in mm.
  • k – 
           Wave number of a wave, see odak.wave.parameters.wavenumber for more.
  • w – 
           Rotation speed of a wave, see odak.wave.parameters.rotationspeed for more.
  • t – 
           Time in seconds.

Returns:

  • new_field ( complex ) –

    A complex number that provides the resultant field in the complex form A*e^(j(wt+phi)).

Source code in odak/wave/vector.py 
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def propagate_plane_waves(field, opd, k, w=0, t=0):
    """
    Definition to propagate a field representing a plane wave at a particular distance and time.

    Parameters
    ----------
    field        : complex
                   Complex field.
    opd          : float
                   Optical path difference in mm.
    k            : float
                   Wave number of a wave, see odak.wave.parameters.wavenumber for more.
    w            : float
                   Rotation speed of a wave, see odak.wave.parameters.rotationspeed for more.
    t            : float
                   Time in seconds.

    Returns
    -------
    new_field     : complex
                    A complex number that provides the resultant field in the complex form A*e^(j(wt+phi)).
    """
    new_field = field*np.exp(1j*(-w*t+opd*k))/opd**2
    return new_field