Fix to deal with unphysical results in incoherent calculations for incident n > 1

.github/workflows/build_deploy_wheels.yml

@@ -19,7 +19,7 @@ jobs:
 
       - name: Deploy wheels - Install dependencies
         run: |
-          python -m pip install --upgrade setuptools wheel pip twine numpy
+          python -m pip install --upgrade setuptools==65.5.1 wheel pip twine numpy
       - name: Deploy wheels - Build manylinux2014 binary wheels - py3.x - x86_64
         uses: RalfG/python-wheels-manylinux-build@v0.3.1-manylinux2014_x86_64
         with:
@@ -72,7 +72,7 @@ jobs:
 
       - name: Install python dependencies
         run: |
-          python -m pip install --upgrade setuptools wheel pip twine numpy==1.21.4
+          python -m pip install --upgrade setuptools==65.5.1 wheel pip twine numpy==1.21.4
 
       - name: Build Solcore Windows
         if: matrix.os == 'windows-latest'

.github/workflows/test_unit_and_examples.yml

@@ -43,7 +43,7 @@ jobs:
 
     - name: Install python dependencies
       run: |
-        python -m pip install --upgrade setuptools wheel pip twine
+        python -m pip install --upgrade setuptools==65.5.1 wheel pip twine
         python -m pip install -e .[dev]
 
     - name: Install S4 in Linux
@@ -127,7 +127,7 @@ jobs:
 
     - name: Install python dependencies
       run: |
-        python -m pip install --upgrade setuptools wheel pip twine
+        python -m pip install --upgrade setuptools==65.5.1 wheel pip twine
         python -m pip install -e .[dev]
 
     - name: Install Solcore Linux and MacOS

solcore/absorption_calculator/tmm_core_vec.py

@@ -173,6 +173,16 @@ def interface_R(polarization, n_i, n_f, th_i, th_f):
     Fraction of light intensity reflected at an interface.
     """
     r = interface_r(polarization, n_i, n_f, th_i, th_f)
+
+    # If the outgoing angle is pi/2, that means the light is totally internally reflected
+    # and we can set r = R = 1. If not, can get unphysical results for r.
+
+    # Note that while T actually CAN be > 1 (when you have incidence from an absorbing
+    # medium), I don't think R can ever be > 1.
+
+    r[th_f.real > (np.pi / 2) - 1e-6] = 1
+    r[th_i.real > (np.pi / 2) - 1e-6] = 1
+
     return R_from_r(r)
 
 
@@ -181,6 +191,14 @@ def interface_T(polarization, n_i, n_f, th_i, th_f):
     Fraction of light intensity transmitted at an interface.
     """
     t = interface_t(polarization, n_i, n_f, th_i, th_f)
+
+    # If the incoming angle is pi/2, that means (most likely) that the light was previously
+    # totally internally reflected. That means the light will never reach this interface and
+    # we can safely set t = T = 0; otherwise we get numerical issues which give unphysically large
+    # values of T because in T_from_t we divide by cos(th_i) which is ~ 0.
+
+    t[th_i.real > (np.pi / 2) - 1e-6] = 0
+
     return T_from_t(polarization, t, n_i, n_f, th_i, th_f)
 
 
@@ -334,12 +352,6 @@ def coh_tmm(pol, n_list, d_list, th_0, lam_vac):
             'pol': pol, 'n_list': n_list, 'd_list': d_list, 'th_0': th_0,
             'lam_vac': lam_vac}
 
-    # return {'r': r, 't': t, 'R': R, 'T': T, 'power_entering': power_entering,
-    #         'kz_list': kz_list, 'th_list': th_list,
-    #         'pol': pol, 'n_list': n_list, 'd_list': d_list, 'th_0': th_0,
-    #         'lam_vac': lam_vac}
-
-
 def coh_tmm_reverse(pol, n_list, d_list, th_0, lam_vac):
     """
     Reverses the order of the stack then runs coh_tmm.
@@ -921,7 +933,6 @@ def inc_tmm(pol, n_list, d_list, c_list, th_0, lam_vac):
         VW = np.matmul(L_list[i], VW_list[i+1].T[:, :, None])
         VW_list[i, :, :] = VW.transpose()
 
-
     # stackFB_list[n]=[F,B] means that F is light traveling forward towards n'th
     # stack and B is light traveling backwards towards n'th stack.
     # Reminder: inc_from_stack[i] = j means that the i'th stack comes after the
@@ -962,6 +973,12 @@ def inc_tmm(pol, n_list, d_list, c_list, th_0, lam_vac):
         stackFB_list_ans = np.stack(stackFB_list).transpose(2, 0, 1)
     else:
         stackFB_list_ans = []
+
+    # despite checking in interface_T and interface_R, still sometimes end up with
+    # unphysical R or T values of incident from medium with n > 1
+    R[R > 1] = 1
+    T[T < 0] = 0
+
     #('VWlist', VW_list.transpose(2, 0, 1))
     ans = {'T': T, 'R': R, 'VW_list': VW_list.transpose(2, 0, 1),
            'coh_tmm_data_list': coh_tmm_data_list,
@@ -1056,6 +1073,11 @@ def inc_position_resolved(layer, dist, inc_tmm_data, coherency_list, alphas, zer
     This function is vectorized. Analogous to position_resolved, but
     for layers (incoherent or coherent) in (partly) incoherent stacks.
     This is a new function, not from Steven Byrnes' tmm package.
+    It assumes that in incoherent layers, we can assume the absorption has
+    a Beer-Lambert profile (this is not really correct; actually, the absorption
+    profile will depend on the coherence length, as discussed in the
+    documentation for the tmm package). This is an approximation in order
+    to be able to generate absorption profiles for partly coherent stacks.
     Starting with output of inc_tmm(), calculate the Poynting vector
     and absorbed energy density a distance "dist" into layer number "layer"
     """
@@ -1071,25 +1093,32 @@ def inc_position_resolved(layer, dist, inc_tmm_data, coherency_list, alphas, zer
             A_layer = fn.run(dist[layer == l])
 
         else:
-            A_layer = beer_lambert(alphas[l] * 1e9, fraction_reaching[i], dist[layer == l] * 1e-9)
+            A_layer = beer_lambert(alphas[l] * 1e9, fraction_reaching[i], dist[layer == l] * 1e-9, A_per_layer[l])
 
         A_layer[fraction_reaching[i] < zero_threshold, :] = 0
         A_local[:, layer == l] = np.real(A_layer)
 
     return A_local
 
 
-def beer_lambert(alphas, fraction, dist):
+def beer_lambert(alphas, fraction, dist, A_total):
     """
     Calculates absorption profile according to the Beer-Lambert law given alphas (in m-1)
     and a vector of distance into the layer (in m) and the fraction of incident light
     reaching the front of the layer. This is used to calculate the absorption profile in
     incoherent layers within (partly) incoherent stacks. Vectorized over wavelengths.
+
+    Note that the Beer-Lambert law is used to generate the absorption profile in incoherent
+    layers as an approximation. The absorption profile in incoherent layers will actually
+    be dependent on the coherence length, as discussed in the documentation for the tmm package.
     """
 
     expn = np.exp(- alphas[:, None] * dist[None,:])
 
-    output = fraction[:, None]*alphas[:, None]*expn
+    A_integrated = fraction*(1-np.exp(-alphas*max(dist)))
 
-    return output/1e9
+    # scale to total absorption in layer (see docstring)
+    scale = np.divide(A_total, A_integrated, out=np.zeros_like(A_total), where=A_integrated!=0)
+    output = scale[:,None]*fraction[:, None]*alphas[:, None]*expn
 
+    return output/1e9

tests/test_tmm_core_vec.py

@@ -1,7 +1,6 @@
 import numpy as np
 from pytest import approx, raises
 
-
 def test_make_2x2_array():
     from solcore.absorption_calculator.tmm_core_vec import make_2x2_array
     z = np.zeros(3)
@@ -91,24 +90,24 @@ def test_power_entering_from_r():
 
 def test_interface_R():
     from solcore.absorption_calculator.tmm_core_vec import interface_R
-    th1 = 0.2
-    n1 = 2
-    n2 = 3
+    th1 = np.array([0.2, np.pi/2])
+    n1 = np.array([2, 2])
+    n2 = np.array([3, 3])
     th2 = np.arcsin((n1/n2)*np.sin(th1))
 
-    assert(interface_R('s', n1, n2, th1, th2) == approx(0.042193712680886314))
-    assert(interface_R('p', n1, n2, th1, th2) == approx(0.037860088421452116))
+    assert(interface_R('s', n1, n2, th1, th2) == approx([0.042193712680886314,1]))
+    assert(interface_R('p', n1, n2, th1, th2) == approx([0.037860088421452116,1]))
 
 
 def test_interface_T():
     from solcore.absorption_calculator.tmm_core_vec import interface_T
-    th1 = 0.2
-    n1 = 2
-    n2 = 3
+    th1 = np.array([0.2, np.pi/2])
+    n1 = np.array([2, 2])
+    n2 = np.array([3, 3])
     th2 = np.arcsin((n1/n2)*np.sin(th1))
 
-    assert(interface_T('s', n1, n2, th1, th2) == approx(0.9578062873191139))
-    assert(interface_T('p', n1, n2, th1, th2) == approx(0.962139911578548))
+    assert(interface_T('s', n1, n2, th1, th2) == approx([0.9578062873191139, 0]))
+    assert(interface_T('p', n1, n2, th1, th2) == approx([0.962139911578548, 0]))
 
 ### Test for coh_tmm
 def test_coh_tmm_exceptions():
@@ -1334,9 +1333,11 @@ def test_beer_lambert():
 
     alphas = np.linspace(0, 1, 5)
     fraction = np.linspace(0.2, 1, 5)
+    A_total = np.array([0, 9.99999989e-09, 2.99999992e-08, 5.99999978e-08,
+                    9.99999950e-08])
     dist = np.linspace(0, 100e-9, 4)
 
-    assert beer_lambert(alphas, fraction, dist)*1e9 == approx(np.array([[0, 0, 0, 0],
+    assert beer_lambert(alphas, fraction, dist, A_total)*1e9 == approx(np.array([[0, 0, 0, 0],
        [0.1, 0.1, 0.1, 0.1],
        [0.3, 0.3, 0.29999999, 0.29999999],
        [0.6, 0.59999999, 0.59999997, 0.59999996],