Skip to content

Incoherent Solver

A 2×2 isotropic solver that supports partly-incoherent stacks. Real coatings often sit on a thick substrate (e.g. a 1 mm glass slide) that is far thicker than the source's coherence length, so light reflected from its two surfaces does not interfere in practice. A fully-coherent solver would produce dense, unphysical Fabry–Perot ripples; this solver lets you mark each interior layer as coherent ('c') or incoherent ('i') via a c_list.

IncoherentIsotropicFilmSolver subclasses IsotropicFilmSolver and shares its constructor and UX, with two differences:

  • it takes an extra c_list argument (one 'c'/'i' code per interior layer), and
  • simulate() returns real power coefficients (Rs, Rp, Ts, Tp) rather than complex amplitudes, because incoherent calculations discard phase by design.

The two semi-infinite media are always treated as incoherent. The solver is fully differentiable through layer thicknesses, so it can be used for inverse design of stacks that include thick substrates.

Isotropic only

Incoherent handling is currently available only for the 2×2 isotropic solver. Anisotropic incoherent TMM is tracked as future work.

difftmm.IncoherentIsotropicFilmSolver 

IncoherentIsotropicFilmSolver(mat_in, mat_out, mat_ls, c_list, thickness_ls=None, thickness_min=0.0, thickness_max=1000.0, batch_size=1, sigmoid_param=False, device=torch.device('cuda'))

Bases: IsotropicFilmSolver

Multi-layer coating solver with partly-incoherent layer support.

Same UX as :class:IsotropicFilmSolver but additionally accepts a c_list argument that marks each interior layer as coherent ('c') or incoherent ('i'). Returns real power coefficients (Rs, Rp, Ts, Tp) rather than complex amplitudes, because incoherent calculations discard phase by design.

Initialize the incoherent isotropic film solver.

Parameters:

Name Type Description Default
mat_in

Incident medium. Float/complex scalar, str material name, or Material.

 required
mat_out

Exit medium. Float/complex scalar, str material name, or Material.

 required
mat_ls

Refractive indices of interior layers (length N). Each entry is a float/complex scalar, str material name, or Material.

 required
c_list

list of 'c'/'i' codes, length N. One per interior layer.

 required
thickness_ls

thicknesses of interior layers in um, length N. If None, random.

 None
thickness_min

minimum thickness in um.

 0.0
thickness_max

maximum thickness in um. Defaults to 1000 (um) so substrate-sized layers fit without special configuration.

 1000.0
batch_size

number of film stacks in the batch.

 1
sigmoid_param

if True, use sigmoid parameterization.

 False
device

torch device.

 device('cuda')

Raises:

Type Description
ValueError

if c_list length does not match the number of interior layers, or contains anything other than 'c' / 'i'.
Source code in difftmm-src/difftmm/film_solver_incoherent.py 
 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
 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
def __init__(
    self,
    mat_in,
    mat_out,
    mat_ls,
    c_list,
    thickness_ls=None,
    thickness_min=0.0,
    thickness_max=1000.0,   # in um; default sized for thick substrates
    batch_size=1,
    sigmoid_param=False,
    device=torch.device("cuda"),
):
    """Initialize the incoherent isotropic film solver.

    Args:
        mat_in: Incident medium. Float/complex scalar, str material name, or Material.
        mat_out: Exit medium. Float/complex scalar, str material name, or Material.
        mat_ls: Refractive indices of interior layers (length N). Each entry is a
            float/complex scalar, str material name, or Material.
        c_list: list of 'c'/'i' codes, length N. One per interior layer.
        thickness_ls: thicknesses of interior layers in um, length N. If None,
                      random.
        thickness_min: minimum thickness in um.
        thickness_max: maximum thickness in um. Defaults to 1000 (um) so
            substrate-sized layers fit without special configuration.
        batch_size: number of film stacks in the batch.
        sigmoid_param: if True, use sigmoid parameterization.
        device: torch device.

    Raises:
        ValueError: if ``c_list`` length does not match the number of
            interior layers, or contains anything other than ``'c'`` / ``'i'``.
    """
    # Validate c_list before delegating to the parent — fail fast with a
    # helpful error rather than letting the parent allocate state we then
    # discard.
    n_layers_provided = (
        mat_ls.shape[0] if torch.is_tensor(mat_ls) else len(mat_ls)
    )
    if len(c_list) != n_layers_provided:
        raise ValueError(
            f"c_list length {len(c_list)} does not match number of interior "
            f"layers {n_layers_provided}."
        )
    for code in c_list:
        if code not in ("c", "i"):
            raise ValueError("c_list entries must be 'c' or 'i'.")
    self.c_list = list(c_list)

    super().__init__(
        mat_in=mat_in,
        mat_out=mat_out,
        mat_ls=mat_ls,
        thickness_ls=thickness_ls,
        thickness_min=thickness_min,
        thickness_max=thickness_max,
        batch_size=batch_size,
        sigmoid_param=sigmoid_param,
        device=device,
    )

coh_stack_power_RT_isotropic staticmethod 

coh_stack_power_RT_isotropic(n_layers_1d, d_1d, wv_1d, n_in, n_out, theta_1d)

Power-domain (Rs, Rp, Ts, Tp) for a coherent isotropic stack.

Thin wrapper around create_jones_matrix_isotropic that converts complex amplitudes to real power coefficients. Used as a building block for the incoherent TMM solver, where each coherent stack contributes its forward/backward (R, T) to the intensity transfer matrix.

Layer convention: layers are ordered from the incident medium (n_in) to the exit medium (n_out), i.e. the same physical order used by tmm_numpy.coh_tmm. Internally, layers are reversed before being passed to create_jones_matrix_isotropic, which stores them in the reverse physical order due to its characteristic-matrix accumulation convention.

Parameters:

Name Type Description Default
n_layers_1d

refractive index of each interior layer in physical order (n_in-side first), shape (batch, n_layer). Complex.

 required
d_1d

thickness of each interior layer in um, same physical ordering as n_layers_1d, shape (batch, n_layer). Real.

 required
wv_1d

wavelengths in um, shape (batch, n_wls). Real.

 required
n_in

scalar incident refractive index (top medium).

 required
n_out

scalar exit refractive index (bottom medium).

 required
theta_1d

incident angles in radians, shape (batch, n_angles). Real, in [0, pi/2].

 required

Returns:

Type Description

Rs, Rp, Ts, Tp: real tensors, each shape (batch, n_wls, n_angles), in [0, 1].
Note

This helper is public-API: it is exposed via difftmm.__init__ in Task 9. Other modules in this package use it as a building block for the incoherent TMM solver, but external users may also call it directly when they need power coefficients without phase.

Source code in difftmm-src/difftmm/film_solver_incoherent.py 
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
137
138
139
140
141
142
143
144
145
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
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
@staticmethod
def coh_stack_power_RT_isotropic(
    n_layers_1d,
    d_1d,
    wv_1d,
    n_in,
    n_out,
    theta_1d,
):
    """Power-domain (Rs, Rp, Ts, Tp) for a coherent isotropic stack.

    Thin wrapper around ``create_jones_matrix_isotropic`` that converts
    complex amplitudes to real power coefficients. Used as a building block
    for the incoherent TMM solver, where each coherent stack contributes its
    forward/backward (R, T) to the intensity transfer matrix.

    Layer convention: layers are ordered from the incident medium (n_in) to
    the exit medium (n_out), i.e. the same physical order used by
    ``tmm_numpy.coh_tmm``.  Internally, layers are reversed before being
    passed to ``create_jones_matrix_isotropic``, which stores them in the
    reverse physical order due to its characteristic-matrix accumulation
    convention.

    Args:
        n_layers_1d: refractive index of each interior layer in physical order
            (n_in-side first), shape (batch, n_layer). Complex.
        d_1d: thickness of each interior layer in um, same physical ordering
            as n_layers_1d, shape (batch, n_layer). Real.
        wv_1d: wavelengths in um, shape (batch, n_wls). Real.
        n_in: scalar incident refractive index (top medium).
        n_out: scalar exit refractive index (bottom medium).
        theta_1d: incident angles in radians, shape (batch, n_angles). Real, in [0, pi/2].

    Returns:
        Rs, Rp, Ts, Tp: real tensors, each shape (batch, n_wls, n_angles), in [0, 1].

    Note:
        This helper is public-API: it is exposed via ``difftmm.__init__`` in
        Task 9. Other modules in this package use it as a building block for
        the incoherent TMM solver, but external users may also call it
        directly when they need power coefficients without phase.
    """
    # Layer convention: callers pass layers in physical order (n_in-side first),
    # matching the convention used by tmm_numpy.coh_tmm. However,
    # create_jones_matrix_isotropic empirically interprets its n_layers_1d
    # argument as ordered from the n_out-side first -- verified by direct
    # comparison with tmm_numpy.coh_tmm in test_coh_stack_power_RT_matches_tmm_numpy.
    # We therefore flip along the layer axis before delegating. If the upstream
    # convention is ever clarified or changed, update this flip and the test
    # will tell us immediately.
    n_layers_rev = torch.flip(n_layers_1d, dims=[1])
    d_rev = torch.flip(d_1d, dims=[1])

    ts, tp, rs, rp = create_jones_matrix_isotropic(
        n_layers_rev, d_rev, wv_1d, n_in, n_out, theta_1d
    )

    # Reflectance is |r|^2 for both polarizations.
    Rs = (rs.real ** 2 + rs.imag ** 2)
    Rp = (rp.real ** 2 + rp.imag ** 2)

    # Transmittance: |t|^2 must be scaled by an intensity-correction factor that
    # depends on how the amplitudes ts/tp are normalized inside
    # _compute_isotropic_tmm. That solver builds the characteristic matrix using
    # two *different* effective indices:
    #     n_eff_s = n * cos(theta)   for s-pol
    #     n_eff_p = n / cos(theta)   for p-pol
    # So ts/tp are normalized to those effective indices, and the power-correction
    # is the ratio of effective indices -- NOT the |E|^2-based ratio
    #     (n_out * conj(cos)) / (n_in * conj(cos))
    # used by tmm_numpy.T_from_t. Concretely:
    #     s-pol: T = |t|^2 * Re(n_out * cos_th_out) / Re(n_in * cos_th_in)
    #     p-pol: T = |t|^2 * Re(n_out / cos_th_out) / Re(n_in / cos_th_in)
    # Energy conservation (R+T = 1 for real lossless stacks) and agreement with
    # tmm_numpy.coh_tmm to ~1e-7 (see test_coh_stack_power_RT_matches_tmm_numpy)
    # confirm this is the consistent power formula for this amplitude normalization.
    device = n_layers_1d.device
    dtype = torch.complex64

    n_in_t = torch.tensor(n_in, dtype=dtype, device=device)
    n_out_t = torch.tensor(n_out, dtype=dtype, device=device)
    cos_th_in = torch.cos(theta_1d.to(dtype)).unsqueeze(1)  # (batch, 1, angles)
    sin_th_in = torch.sin(theta_1d.to(dtype)).unsqueeze(1)
    sin_th_out = n_in_t * sin_th_in / n_out_t
    cos_th_out = torch.sqrt(1 - sin_th_out ** 2)

    # s-pol: ratio of (n * cos_theta) quantities.
    num_s = (n_out_t * cos_th_out).real   # n_out * cos(th_out)
    den_s = (n_in_t * cos_th_in).real    # n_in  * cos(th_in)
    Ts = (ts.real ** 2 + ts.imag ** 2) * (num_s / den_s)

    # p-pol: ratio of (n / cos_theta) quantities.
    num_p = (n_out_t / cos_th_out).real   # n_out / cos(th_out)
    den_p = (n_in_t / cos_th_in).real    # n_in  / cos(th_in)
    Tp = (tp.real ** 2 + tp.imag ** 2) * (num_p / den_p)

    return Rs, Rp, Ts, Tp

group_layers_by_coherence staticmethod 

group_layers_by_coherence(c_list: Sequence[str]) -> Dict[str, object]

Group a layer-coherence list into coherent stacks and incoherent layers.

A "stack" is a maximal run of consecutive coherent layers, bracketed on each side by an incoherent layer (which is required to be present because the first and last layers are semi-infinite and must be 'i').

Parameters:

Name Type Description Default
c_list Sequence[str]

Per-layer coherence flags. Each entry is 'i' (incoherent) or 'c' (coherent). First and last entries must be 'i' because those layers are semi-infinite.

 required

Returns:

Type Description
Dict[str, object]

Dict with keys: - num_inc_layers (int): number of incoherent layers. - num_stacks (int): number of coherent stacks. - inc_alllayer_indices (List[int]): for each incoherent layer i, its index in the original layer list. - stack_alllayer_indices (List[List[int]]): for each stack s, the original indices of layers in [bracketing_inc, coh_layers..., bracketing_inc]. - stack_after_inc (List[Optional[int]]): for each incoherent layer i, the stack-index of the stack immediately after it, or None if the next layer is also incoherent (or there is no next layer). - inc_after_stack (List[int]): for each stack s, the incoherent-layer index that immediately precedes the stack.

Raises:

Type Description
ValueError

if the first or last entry is not 'i', or any entry is neither 'i' nor 'c'.
Source code in difftmm-src/difftmm/film_solver_incoherent.py 
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
252
253
254
255
256
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
@staticmethod
def group_layers_by_coherence(c_list: Sequence[str]) -> Dict[str, object]:
    """Group a layer-coherence list into coherent stacks and incoherent layers.

    A "stack" is a maximal run of consecutive coherent layers, bracketed on
    each side by an incoherent layer (which is required to be present because
    the first and last layers are semi-infinite and must be 'i').

    Args:
        c_list: Per-layer coherence flags. Each entry is 'i' (incoherent) or
            'c' (coherent). First and last entries must be 'i' because those
            layers are semi-infinite.

    Returns:
        Dict with keys:
            - num_inc_layers (int): number of incoherent layers.
            - num_stacks (int): number of coherent stacks.
            - inc_alllayer_indices (List[int]): for each incoherent layer i,
              its index in the original layer list.
            - stack_alllayer_indices (List[List[int]]): for each stack s,
              the original indices of layers in [bracketing_inc, coh_layers..., bracketing_inc].
            - stack_after_inc (List[Optional[int]]): for each incoherent layer i,
              the stack-index of the stack immediately after it, or None if the
              next layer is also incoherent (or there is no next layer).
            - inc_after_stack (List[int]): for each stack s, the incoherent-layer
              index that immediately precedes the stack.

    Raises:
        ValueError: if the first or last entry is not 'i', or any entry is
            neither 'i' nor 'c'.
    """
    if len(c_list) < 2 or c_list[0] != "i" or c_list[-1] != "i":
        raise ValueError("c_list must start and end with 'i' (semi-infinite layers).")
    for code in c_list:
        if code not in ("i", "c"):
            raise ValueError("c_list entries must be 'i' or 'c'.")

    inc_alllayer_indices: List[int] = []
    stack_alllayer_indices: List[List[int]] = []
    stack_after_inc: List[Optional[int]] = []
    inc_after_stack: List[int] = []

    inc_index = -1  # incremented when we visit an 'i' layer
    in_stack = False
    current_stack: List[int] = []

    for layer_idx, code in enumerate(c_list):
        if code == "i":
            inc_index += 1
            inc_alllayer_indices.append(layer_idx)
            if in_stack:
                # Close out the stack with this incoherent layer as its right bracket.
                current_stack.append(layer_idx)
                stack_alllayer_indices.append(current_stack)
                # Whoever opened this stack was the previous incoherent layer,
                # which is at incoherent index (inc_index - 1).
                inc_after_stack.append(inc_index - 1)
                current_stack = []
                in_stack = False
                # This 'i' has no stack following it *yet*; will be patched if the next 'c' opens one.
                stack_after_inc.append(None)
            else:
                stack_after_inc.append(None)
        else:  # 'c'
            if not in_stack:
                # Open a stack: left bracket is the previous incoherent layer.
                in_stack = True
                current_stack = [inc_alllayer_indices[-1], layer_idx]
                # The most recent incoherent layer is followed by a stack whose
                # index will be len(stack_alllayer_indices) once the stack closes.
                stack_after_inc[-1] = len(stack_alllayer_indices)
            else:
                current_stack.append(layer_idx)

    return {
        "num_inc_layers": len(inc_alllayer_indices),
        "num_stacks": len(stack_alllayer_indices),
        "inc_alllayer_indices": inc_alllayer_indices,
        "stack_alllayer_indices": stack_alllayer_indices,
        "stack_after_inc": stack_after_inc,
        "inc_after_stack": inc_after_stack,
    }

interface_power_RT staticmethod 

interface_power_RT(n_i, n_f, cos_i, cos_f)

Fresnel power reflectance/transmittance at a single interface.

All inputs are complex tensors broadcastable to a common shape. Returns (Rs, Rp, Ts, Tp) as real tensors of that shape.

Math

r_s = (n_i cos_i - n_f cos_f) / (n_i cos_i + n_f cos_f) r_p = (n_f cos_i - n_i cos_f) / (n_f cos_i + n_i cos_f) t_s = 2 n_i cos_i / (n_i cos_i + n_f cos_f) t_p = 2 n_i cos_i / (n_f cos_i + n_i cos_f) R = |r|^2 T_s = |t_s|^2 * Re(n_f cos_f) / Re(n_i cos_i) T_p = |t_p|^2 * Re(n_f conj(cos_f)) / Re(n_i conj(cos_i))

Source code in difftmm-src/difftmm/film_solver_incoherent.py 
290
291
292
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
324
325
326
327
@staticmethod
def interface_power_RT(n_i, n_f, cos_i, cos_f):
    """Fresnel power reflectance/transmittance at a single interface.

    All inputs are complex tensors broadcastable to a common shape.
    Returns (Rs, Rp, Ts, Tp) as real tensors of that shape.

    Math:
        r_s = (n_i cos_i - n_f cos_f) / (n_i cos_i + n_f cos_f)
        r_p = (n_f cos_i - n_i cos_f) / (n_f cos_i + n_i cos_f)
        t_s = 2 n_i cos_i / (n_i cos_i + n_f cos_f)
        t_p = 2 n_i cos_i / (n_f cos_i + n_i cos_f)
        R = |r|^2
        T_s = |t_s|^2 * Re(n_f cos_f) / Re(n_i cos_i)
        T_p = |t_p|^2 * Re(n_f conj(cos_f)) / Re(n_i conj(cos_i))
    """
    n_i_cos_i = n_i * cos_i
    n_f_cos_f = n_f * cos_f
    n_f_cos_i = n_f * cos_i
    n_i_cos_f = n_i * cos_f

    r_s = (n_i_cos_i - n_f_cos_f) / (n_i_cos_i + n_f_cos_f)
    r_p = (n_f_cos_i - n_i_cos_f) / (n_f_cos_i + n_i_cos_f)
    t_s = 2 * n_i_cos_i / (n_i_cos_i + n_f_cos_f)
    t_p = 2 * n_i_cos_i / (n_f_cos_i + n_i_cos_f)

    Rs = r_s.real ** 2 + r_s.imag ** 2
    Rp = r_p.real ** 2 + r_p.imag ** 2

    ts_sq = t_s.real ** 2 + t_s.imag ** 2
    tp_sq = t_p.real ** 2 + t_p.imag ** 2

    Ts = ts_sq * (n_f_cos_f.real / n_i_cos_i.real)
    Tp = tp_sq * (
        (n_f * torch.conj(cos_f)).real / (n_i * torch.conj(cos_i)).real
    )

    return Rs, Rp, Ts, Tp

create_intensity_RT_isotropic staticmethod 

create_intensity_RT_isotropic(n_layers_1d, d_1d, wv_1d, n_in, n_out, theta_1d, c_list)

Partly-incoherent intensity TMM for isotropic multi-layer films.

Each interior layer is marked coherent ('c') or incoherent ('i') via c_list. The two semi-infinite media (n_in, n_out) are always incoherent, so c_list describes only the interior layers and the full coherence sequence is ['i'] + c_list + ['i'].

The algorithm
  1. Group consecutive coherent layers into stacks.
  2. For each stack, run the existing coherent 2x2 solver in both directions to get (R_fwd, T_fwd, R_bwd, T_bwd).
  3. Compute single-pass absorption P_i for each interior incoherent layer.
  4. Build real 2x2 intensity transfer matrices L_i and multiply.
  5. Read off total R and T.

Parameters:

Name Type Description Default
n_layers_1d

refractive indices of interior layers, shape (batch, n_layer). Complex.

 required
d_1d

thicknesses of interior layers in um, shape (batch, n_layer). Real.

 required
wv_1d

wavelengths in um, shape (batch, n_wls). Real.

 required
n_in

scalar incident refractive index (top medium).

 required
n_out

scalar exit refractive index (bottom medium).

 required
theta_1d

incident angles in radians, shape (batch, n_angles). Real, in [0, pi/2].

 required
c_list

list of 'c'/'i' codes, length n_layer (interior layers only).

 required

Returns:

Type Description

Rs, Rp, Ts, Tp: real tensors, each shape (batch, n_wls, n_angles), in [0, 1].
Note

Inside any coherent sub-stack, the refractive indices of the bracketing incoherent layers must be identical across the batch dimension. Heterogeneous per-batch indices in stack-bracketing layers raise ValueError. Per-batch interior coherent-layer indices are still supported.

Source code in difftmm-src/difftmm/film_solver_incoherent.py 
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
365
366
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
417
418
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
@staticmethod
def create_intensity_RT_isotropic(
    n_layers_1d,
    d_1d,
    wv_1d,
    n_in,
    n_out,
    theta_1d,
    c_list,
):
    """Partly-incoherent intensity TMM for isotropic multi-layer films.

    Each interior layer is marked coherent ('c') or incoherent ('i') via
    ``c_list``. The two semi-infinite media (n_in, n_out) are always
    incoherent, so ``c_list`` describes only the *interior* layers and the
    full coherence sequence is ``['i'] + c_list + ['i']``.

    The algorithm:
      1. Group consecutive coherent layers into stacks.
      2. For each stack, run the existing coherent 2x2 solver in both
         directions to get (R_fwd, T_fwd, R_bwd, T_bwd).
      3. Compute single-pass absorption P_i for each interior incoherent
         layer.
      4. Build real 2x2 intensity transfer matrices L_i and multiply.
      5. Read off total R and T.

    Args:
        n_layers_1d: refractive indices of interior layers, shape (batch, n_layer). Complex.
        d_1d: thicknesses of interior layers in um, shape (batch, n_layer). Real.
        wv_1d: wavelengths in um, shape (batch, n_wls). Real.
        n_in: scalar incident refractive index (top medium).
        n_out: scalar exit refractive index (bottom medium).
        theta_1d: incident angles in radians, shape (batch, n_angles). Real, in [0, pi/2].
        c_list: list of 'c'/'i' codes, length n_layer (interior layers only).

    Returns:
        Rs, Rp, Ts, Tp: real tensors, each shape (batch, n_wls, n_angles), in [0, 1].

    Note:
        Inside any coherent sub-stack, the refractive indices of the
        bracketing incoherent layers must be identical across the batch
        dimension. Heterogeneous per-batch indices in stack-bracketing
        layers raise ValueError. Per-batch interior coherent-layer indices
        are still supported.
    """
    if len(c_list) != n_layers_1d.shape[1]:
        raise ValueError(
            f"c_list length {len(c_list)} does not match interior-layer count "
            f"{n_layers_1d.shape[1]}."
        )

    device = n_layers_1d.device
    real_dtype = torch.float32
    complex_dtype = torch.complex64

    batchsize = n_layers_1d.shape[0]
    num_wv = wv_1d.shape[1]
    num_angles = theta_1d.shape[1]

    # Full coherence sequence including semi-infinite endpoints.
    full_c_list = ["i"] + list(c_list) + ["i"]
    groups = IncoherentIsotropicFilmSolver.group_layers_by_coherence(full_c_list)

    num_inc = groups["num_inc_layers"]
    inc_alllayer_indices = groups["inc_alllayer_indices"]
    stack_alllayer_indices = groups["stack_alllayer_indices"]
    stack_after_inc = groups["stack_after_inc"]

    # Build full per-layer index tensors with the endpoints in place.
    n_in_col = torch.full(
        (batchsize, 1), n_in, dtype=complex_dtype, device=device
    )
    n_out_col = torch.full(
        (batchsize, 1), n_out, dtype=complex_dtype, device=device
    )
    n_full = torch.cat([n_in_col, n_layers_1d.to(complex_dtype), n_out_col], dim=1)

    d_in_col = torch.full((batchsize, 1), float("inf"), dtype=real_dtype, device=device)
    d_out_col = torch.full((batchsize, 1), float("inf"), dtype=real_dtype, device=device)
    d_full = torch.cat([d_in_col, d_1d.to(real_dtype), d_out_col], dim=1)

    # Snell's law to get cos(theta) in each layer.
    sin_th_in = torch.sin(theta_1d).to(complex_dtype)  # (batch, angles)
    n_full_b = n_full.unsqueeze(-1)           # (batch, n_full, 1)
    n_in_b = n_in_col.unsqueeze(-1)           # (batch, 1, 1)
    sin_th_in_b = sin_th_in.unsqueeze(1)      # (batch, 1, angles)
    sin_th_layers = n_in_b * sin_th_in_b / n_full_b   # (batch, n_full, angles)
    cos_th_layers = torch.sqrt(1 - sin_th_layers ** 2)  # complex

    # Single-pass absorption P_i for each *interior* incoherent layer.
    # wv_b shape: (batch, 1, n_wv, 1) for broadcasting over n_int_inc and angles.
    wv_b = wv_1d.unsqueeze(-1).unsqueeze(1)  # (batch, 1, n_wv, 1)
    interior_inc_alllayer = [
        inc_alllayer_indices[i] for i in range(1, num_inc - 1)
    ]
    if len(interior_inc_alllayer) > 0:
        idx = torch.tensor(interior_inc_alllayer, dtype=torch.long, device=device)
        n_inc_interior = n_full.index_select(1, idx)          # (batch, n_int_inc)
        d_inc_interior = d_full.index_select(1, idx)          # (batch, n_int_inc)
        cos_inc_interior = cos_th_layers.index_select(1, idx) # (batch, n_int_inc, angles)
        n_inc_b = n_inc_interior.unsqueeze(-1).unsqueeze(-1)  # (batch, n_int_inc, 1, 1)
        d_inc_b = d_inc_interior.unsqueeze(-1).unsqueeze(-1)  # (batch, n_int_inc, 1, 1)
        cos_inc_b = cos_inc_interior.unsqueeze(2)             # (batch, n_int_inc, 1, angles)
        imag_part = (n_inc_b * cos_inc_b).imag
        P_interior = torch.exp(-4 * torch.pi * d_inc_b.real * imag_part / wv_b.real)
        # shape: (batch, n_int_inc, n_wv, angles)
        P_interior = torch.clamp(P_interior, min=1e-30)
    else:
        P_interior = None

    Rs_total, Rp_total, Ts_total, Tp_total = IncoherentIsotropicFilmSolver._inc_total_RT_for_all_pols(
        n_full=n_full,
        d_full=d_full,
        cos_th_layers=cos_th_layers,
        wv_1d=wv_1d,
        theta_1d=theta_1d,
        n_in=n_in,
        n_out=n_out,
        groups=groups,
        P_interior=P_interior,
    )

    return Rs_total, Rp_total, Ts_total, Tp_total

simulate 

simulate(theta, wvln)

Compute (Rs, Rp, Ts, Tp) for the configured stack.

Parameters:

Name Type Description Default
theta

angles in radians. 1D (n_angles,) or 2D (batch, n_angles).

 required
wvln

wavelengths in um. Scalar, list, or 1D tensor.

 required

Returns:

Type Description

Rs, Rp, Ts, Tp: real tensors of shape (batch, n_wvlns, n_angles).
Source code in difftmm-src/difftmm/film_solver_incoherent.py 
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
def simulate(self, theta, wvln):
    """Compute (Rs, Rp, Ts, Tp) for the configured stack.

    Args:
        theta: angles in radians. 1D ``(n_angles,)`` or 2D ``(batch, n_angles)``.
        wvln: wavelengths in um. Scalar, list, or 1D tensor.

    Returns:
        Rs, Rp, Ts, Tp: real tensors of shape ``(batch, n_wvlns, n_angles)``.
    """
    theta, wv_batch, d_batch, n_in_t, n_out_t, n_layers_t = (
        self._prepare_simulate_inputs(theta, wvln)
    )
    # The incoherent kernel treats refractive indices as constant across
    # wavelengths (no wavelength dim in n_layers_1d / n_in / n_out). For
    # dispersive Materials, the first-wavelength value is used; for
    # non-dispersive materials this is exact.
    n_in_val = n_in_t[0, 0].item()
    n_out_val = n_out_t[0, 0].item()
    n_layers_const = n_layers_t[:, 0, :]
    return IncoherentIsotropicFilmSolver.create_intensity_RT_isotropic(
        n_layers_const, d_batch, wv_batch,
        n_in_val, n_out_val, theta, self.c_list,
    )