Further Sesame changes

- Refactor Sesame PDD into two files, update tests and documentation/docstrings
- Remove changes to material_system causing recursion/stack overflow errors
- add type annotations and docstrings. Black formatting
- fix DA dark current overflow and relevant test
- Improve automatic mesh for Sesame PDD
- update documentation, add tests
- make labelling of SRVs consistent between DA and Sesame.
- remove GitHub reference for Sesame installation, refer to new PyPI version
- Update outdated examples
- Fix issue where plotting error causing examples in GitHub actions to fail
- Change build for S4 to install setuptools explicitly (necessary in Python 3.12)
- Sesame now scans from long to short wavelengths in QE calculations to improve changes of convergence, and skips wavelengths where absorption = 0
- Fix issue where Fortran PDD did not work if all voltages < 0
- Upgrade from devpy to spin (renamed), and to latest version of meson in the build system

.github/workflows/build_deploy_wheels.yml

@@ -5,7 +5,7 @@ on:
     branches:
       - main
       - develop
-      - update_build_312
+
     tags:
       - "**"
 

.github/workflows/test_unit_and_examples.yml

@@ -18,7 +18,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest, windows-latest]
-        python-version: ["3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.7", "3.8", "3.9", "3.10", "3.11", "3.12"]
 
     steps:
       - name: Checkout
@@ -53,13 +53,14 @@ jobs:
 
       - name: Install Python dependecies
         run: |
-          pip install pytest meson-python ninja cython numpy git+https://github.com/scientific-python/devpy@v0.1
-          python3 -m devpy install-dependencies -test-dep
+          pip install numpy --config-settings=setup-args="-Dallow-noblas=true"
+          pip install pytest meson-python ninja cython spin
+          python3 -m spin install-dependencies -test-dep
 
       - name: Install S4
         if: matrix.os != 'windows-latest'
         run: |
-          pip install wheel
+          pip install wheel setuptools
           git clone https://github.com/phoebe-p/S4
           cd S4
           make S4_pyext
@@ -68,14 +69,14 @@ jobs:
 
       - name: Build solcore
         run: |
-          python -m devpy build -- -Dwith_pdd=true -Dinstall_test=true
+          python -m spin build -- -Dwith_pdd=true -Dinstall_test=true
 
       - name: Unit and functional tests (MacOS and Linux)
         if: matrix.os != 'windows-latest'
         env:
           SOLCORE_SPICE: ngspice
         run: |
-          python -m devpy test -- -r a -v --cov=solcore/ --ignore=solcore/tests/test_examples.py -n "auto" 
+          python -m spin test -- -r a -v --cov=solcore/ --ignore=solcore/tests/test_examples.py -n "auto" 
 
       - name: Unit and functional tests (Windows)
         if: matrix.os == 'windows-latest'
@@ -88,7 +89,7 @@ jobs:
           CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
         run: |
           python -m pip install codecov
-          python -m devpy codecov 
+          python -m spin codecov 
 
 
   test_examples:
@@ -133,13 +134,13 @@ jobs:
 
       - name: Install Python dependecies
         run: |
-          pip install pytest meson-python ninja cython numpy git+https://github.com/scientific-python/devpy@v0.1
-          python3 -m devpy install-dependencies -test-dep
+          pip install pytest meson-python ninja cython numpy spin
+          python3 -m spin install-dependencies -test-dep
 
       - name: Install S4
         if: matrix.os != 'windows-latest'
         run: |
-          pip install wheel
+          pip install wheel setuptools
           git clone https://github.com/phoebe-p/S4
           cd S4
           make S4_pyext
@@ -148,16 +149,16 @@ jobs:
 
       - name: Build solcore
         run: |
-          python -m devpy build -- -Dwith_pdd=true -Dinstall_test=true
+          python -m spin build -- -Dwith_pdd=true -Dinstall_test=true
 
       - name: Unit and functional tests (MacOS and Linux)
         if: matrix.os != 'windows-latest'
         env:
           SOLCORE_SPICE: ngspice
         run: |
-          python -m devpy test -- -r a -v solcore/tests/test_examples.py -n "auto"
+          python -m spin test -- -r a -v solcore/tests/test_examples.py -n "auto"
 
 #      - name: Unit and functional tests (Windows)
 #        if: matrix.os == 'windows-latest'
 #        run: |
-#          python -m devpy test -- -r a -v solcore/tests/test_examples.py
+#          python -m spin test -- -r a -v solcore/tests/test_examples.py

.devpy/cmds.py → .spin/cmds.py

@@ -1,6 +1,6 @@
 import click
-from devpy import util
-from devpy.cmds.meson import _get_site_packages
+from spin import util
+from spin.cmds.meson import _get_site_packages
 
 
 @click.command()

docs/source/Examples/example_3J_with_DA_solver.rst

@@ -6,128 +6,156 @@ Example of a 3J solar cell calculated with the DA solver
 .. image:: DA_qe.png
    :width: 40%
 
-- Required extra files, available in `Solcore's Github repository (Examples folder) <https://github.com/dalonsoa/solcore5>`_:
-
-    - MgF-ZnS_AR.csv
-    - in01gaas.csv
-    - Ge-Palik.csv
-
 .. code-block:: Python
 
     import numpy as np
     import matplotlib.pyplot as plt
 
     from solcore import siUnits, material, si
-    from solcore.interpolate import interp1d
     from solcore.solar_cell import SolarCell
     from solcore.structure import Junction, Layer
     from solcore.solar_cell_solver import solar_cell_solver
+    from solcore.light_source import LightSource
 
-    all_materials = []
-
+    # Define materials for the anti-reflection coating:
+    MgF2 = material("MgF2")()
+    ZnS = material("ZnScub")()
 
-    def this_dir_file(f):
-        from pathlib import Path
-        return str(Path(__file__).parent / "data" /  f)
+    ARC_layers = [Layer(si("100nm"), material=MgF2),
+                  Layer(si("50nm"), material=ZnS)]
 
+    # TOP CELL - InGaP
+    # Now we build the top cell, which requires the n and p sides of GaInP and a window
+    # layer. We also add some extra parameters needed for the calculation using the
+    # depletion approximation: the minority carriers diffusion lengths and the doping.
 
-    # We need to build the solar cell layer by layer.
-    # We start from the AR coating. In this case, we load it from an an external file
-    refl_nm = np.loadtxt(this_dir_file("MgF-ZnS_AR.csv"), unpack=True, delimiter=",")
-    ref = interp1d(x=siUnits(refl_nm[0], "nm"), y=refl_nm[1], bounds_error=False, fill_value=0)
-
-    # TOP CELL - GaInP
-    # Now we build the top cell, which requires the n and p sides of GaInP and a window layer.
-    # We also load the absorption coefficient from an external file. We also add some extra parameters needed for the
-    # calculation such as the minority carriers diffusion lengths
     AlInP = material("AlInP")
     InGaP = material("GaInP")
     window_material = AlInP(Al=0.52)
-    top_cell_n_material = InGaP(In=0.49, Nd=siUnits(2e18, "cm-3"), hole_diffusion_length=si("200nm"))
-    top_cell_p_material = InGaP(In=0.49, Na=siUnits(1e17, "cm-3"), electron_diffusion_length=si("1um"))
-
-    all_materials.append(window_material)
-    all_materials.append(top_cell_n_material)
-    all_materials.append(top_cell_p_material)
 
-    # MID CELL  - InGaAs
-    # We add manually the absorption coefficient of InGaAs since the one contained in the database doesn't cover
-    # enough range, keeping in mind that the data has to be provided as a function that takes wavelengths (m) as input and
-    # returns absorption (1/m)
-    InGaAs = material("InGaAs")
-    InGaAs_alpha = np.loadtxt(this_dir_file("in01gaas.csv"), unpack=True, delimiter=",")
-    InGaAs.alpha = interp1d(x=1240e-9 / InGaAs_alpha[0][::-1], y=InGaAs_alpha[1][::-1], bounds_error=False, fill_value=0)
+    top_cell_n_material = InGaP(In=0.49, Nd=siUnits(2e18, "cm-3"),
+                                hole_diffusion_length=si("200nm"))
+    top_cell_p_material = InGaP(In=0.49, Na=siUnits(1e17, "cm-3"),
+                                electron_diffusion_length=si("2um"))
+    # MID CELL  - GaAs
+    GaAs = material("GaAs")
 
-    mid_cell_n_material = InGaAs(In=0.01, Nd=siUnits(3e18, "cm-3"), hole_diffusion_length=si("500nm"))
-    mid_cell_p_material = InGaAs(In=0.01, Na=siUnits(1e17, "cm-3"), electron_diffusion_length=si("5um"))
-
-    all_materials.append(mid_cell_n_material)
-    all_materials.append(mid_cell_p_material)
+    mid_cell_n_material = GaAs(In=0.01, Nd=siUnits(3e18, "cm-3"),
+                               hole_diffusion_length=si("500nm"))
+    mid_cell_p_material = GaAs(In=0.01, Na=siUnits(1e17, "cm-3"),
+                               electron_diffusion_length=si("5um"))
 
     # BOTTOM CELL - Ge
-    # We add manually the absorption coefficient of Ge since the one contained in the database doesn't cover
-    # enough range.
     Ge = material("Ge")
-    Ge_alpha = np.loadtxt(this_dir_file("Ge-Palik.csv"), unpack=True, delimiter=",")
-    Ge.alpha = interp1d(x=1240e-9 / Ge_alpha[0][::-1], y=Ge_alpha[1][::-1], bounds_error=False, fill_value=0)
-
-    bot_cell_n_material = Ge(Nd=siUnits(2e18, "cm-3"), hole_diffusion_length=si("800nm"))
-    bot_cell_p_material = Ge(Na=siUnits(1e17, "cm-3"), electron_diffusion_length=si("50um"))
-
-    all_materials.append(bot_cell_n_material)
-    all_materials.append(bot_cell_p_material)
-
-    # We add some other properties to the materials, assumed the same in all cases, for simplicity.
-    # If different, we should have added them above in the definition of the materials.
-    for mat in all_materials:
-        mat.hole_mobility = 5e-2
-        mat.electron_mobility = 3.4e-3
-        mat.hole_mobility = 3.4e-3
-        mat.electron_mobility = 5e-2
-        mat.relative_permittivity = 9
-
-    # And, finally, we put everything together, adding also the surface recombination velocities. We also add some shading
-    # due to the metallisation of the cell = 8%, and indicate it has an area of 0.7x0.7 mm2 (converted to m2)
+
+    bot_cell_n_material = Ge(Nd=siUnits(2e18, "cm-3"),
+                             hole_diffusion_length=si("800nm"))
+    bot_cell_p_material = Ge(Na=siUnits(1e17, "cm-3"),
+                             electron_diffusion_length=si("50um"))
+
+    # Now that the layers are configured, we can now assemble the triple junction solar
+    # cell. Note that we also specify a metal shading of 2% and a cell area of $1cm^{2}$.
+    # SolCore calculates the EQE for all three junctions and light-IV showing the relative
+    # contribution of each sub-cell. We set "kind = 'DA'" to use the depletion
+    # approximation. We can also set the surface recombination velocities, where sn
+    # refers to the surface recombination velocity at the n-type surface, and sp refers
+    # to the SRV on the p-type side.
+
     solar_cell = SolarCell(
-        [
+            ARC_layers +
+            [
             Junction([Layer(si("25nm"), material=window_material, role='window'),
                       Layer(si("100nm"), material=top_cell_n_material, role='emitter'),
-                      Layer(si("600nm"), material=top_cell_p_material, role='base'),
-                      ], sn=1, sp=1, kind='DA'),
+                      Layer(si("400nm"), material=top_cell_p_material, role='base'),
+                      ], sn=si("1e5cm s-1"), sp=si("1e5cm s-1"), kind='DA'),
             Junction([Layer(si("200nm"), material=mid_cell_n_material, role='emitter'),
                       Layer(si("3000nm"), material=mid_cell_p_material, role='base'),
-                      ], sn=1, sp=1, kind='DA'),
+                      ], sn=si("1e5cm s-1"), sp=si("1e5cm s-1"), kind='DA'),
             Junction([Layer(si("400nm"), material=bot_cell_n_material, role='emitter'),
                       Layer(si("100um"), material=bot_cell_p_material, role='base'),
-                      ], sn=1, sp=1, kind='DA'),
-        ], reflectivity=ref, shading=0.08, cell_area=0.7 * 0.7 / 1e4)
-
-    wl = np.linspace(300, 1800, 700) * 1e-9
-    solar_cell_solver(solar_cell, 'qe', user_options={'wavelength': wl})
+                      ], sn=si("1e5cm s-1"), sp=si("1e5cm s-1"), kind='DA')
+                ],
+            shading=0.02, cell_area=1 * 1 / 1e4)
+
+    # Choose wavelength range (in m):
+    wl = np.linspace(280, 1850, 700) * 1e-9
+
+    # Calculate the EQE for the solar cell:
+    solar_cell_solver(solar_cell, 'qe', user_options={'wavelength': wl,
+                                                      'da_mode': 'green',
+                                                      'optics_method': 'TMM'
+                                                      })
+    # we pass options to use the TMM optical method to calculate realistic R, A and T
+    # values with interference in the ARC (and semiconductor) layers. We can also choose
+    # which solver mode to use for the depletion approximation. The default is 'green',
+    # which uses the (faster) Green's function method. The other method is 'bvp'.
+
+    # Plot the EQE and absorption of the individual junctions. Note that we can access
+    # the properties of the first junction (ignoring any other layers, such as the ARC,
+    # which are not part of a junction) using solar_cell(0), and the second junction using
+    # solar_cell(1), etc.
 
     plt.figure(1)
-    plt.plot(wl * 1e9, solar_cell[0].eqe(wl) * 100, 'b', label='GaInP')
-    plt.plot(wl * 1e9, solar_cell[1].eqe(wl) * 100, 'g', label='InGaAs')
-    plt.plot(wl * 1e9, solar_cell[2].eqe(wl) * 100, 'r', label='Ge')
-
+    plt.plot(wl * 1e9, solar_cell(0).eqe(wl) * 100, 'b', label='GaInP QE')
+    plt.plot(wl * 1e9, solar_cell(1).eqe(wl) * 100, 'g', label='GaAs QE')
+    plt.plot(wl * 1e9, solar_cell(2).eqe(wl) * 100, 'r', label='Ge QE')
+    plt.fill_between(wl * 1e9, solar_cell(0).layer_absorption * 100, 0, alpha=0.3,
+             label='GaInP Abs.', color='b')
+    plt.fill_between(wl * 1e9, solar_cell(1).layer_absorption * 100, 0, alpha=0.3,
+             label='GaAs Abs.', color='g')
+    plt.fill_between(wl * 1e9, solar_cell(2).layer_absorption * 100, 0, alpha=0.3,
+             label='Ge Abs.', color='r')
+
+    plt.plot(wl*1e9, 100*(1-solar_cell.reflected), '--k', label="100 - Reflectivity")
     plt.legend()
     plt.ylim(0, 100)
     plt.ylabel('EQE (%)')
     plt.xlabel('Wavelength (nm)')
+    plt.show()
+
+    # Set up the AM0 (space) solar spectrum for the light I-V calculation:
+    am0 = LightSource(source_type='standard',version='AM0',x=wl,
+                      output_units='photon_flux_per_m')
+
+
+    # Set up the voltage range for the overall cell (at which the total I-V will be
+    # calculated) as well as the internal voltages which are used to calculate the results
+    # for the individual junctions. The range of the internal_voltages should generally
+    # be wider than that for the voltages.
 
-    V = np.linspace(0, 3, 300)
-    solar_cell_solver(solar_cell, 'iv', user_options={'voltages': V, 'light_iv': True, 'wavelength': wl})
+    # this is an n-p cell, so we need to scan negative voltages
+
+    V = np.linspace(-3, 0, 300)
+    internal_voltages = np.linspace(-4, 2, 400)
+
+    # Calculate the current-voltage relationship under illumination:
+
+    solar_cell_solver(solar_cell, 'iv', user_options={'light_source': am0,
+                                                      'voltages': V,
+                                                      'internal_voltages': internal_voltages,
+                                                      'light_iv': True,
+                                                      'wavelength': wl,
+                                                      'optics_method': 'TMM',
+                                                      'mpp': True,
+                                                      })
+
+    # We pass the same options as for solving the EQE, but also set 'light_iv' and 'mpp' to
+    # True to indicate we want the IV curve under illumination and to find the maximum
+    # power point (MPP). We also pass the AM0 light source and voltages created above.
 
     plt.figure(2)
-    plt.plot(V, solar_cell.iv['IV'][1], 'k', linewidth=3, label='Total')
-    plt.plot(V, -solar_cell[0].iv(V), 'b', label='GaInP')
-    plt.plot(V, -solar_cell[1].iv(V), 'g', label='InGaAs')
-    plt.plot(V, -solar_cell[2].iv(V), 'r', label='Ge')
+    plt.plot(abs(V), -solar_cell.iv['IV'][1]/10, 'k', linewidth=3, label='3J cell')
+    plt.plot(abs(V), solar_cell(0).iv(V)/10, 'b', label='InGaP sub-cell')
+    plt.plot(abs(V), solar_cell(1).iv(V)/10, 'g', label='GaAs sub-cell')
+    plt.plot(abs(V), solar_cell(2).iv(V)/10, 'r', label='Ge sub-cell')
+    plt.text(0.5,30,f'Jsc= {abs(solar_cell.iv.Isc/10):.2f} mA.cm' + r'$^{-2}$')
+    plt.text(0.5,28,f'Voc= {abs(solar_cell.iv.Voc):.2f} V')
+    plt.text(0.5,26,f'FF= {solar_cell.iv.FF*100:.2f} %')
+    plt.text(0.5,24,f'Eta= {solar_cell.iv.Eta*100:.2f} %')
 
     plt.legend()
-    plt.ylim(0, 200)
+    plt.ylim(0, 33)
     plt.xlim(0, 3)
-    plt.ylabel('Current (A/m$^2$)')
+    plt.ylabel('Current (mA/cm$^2$)')
     plt.xlabel('Voltage (V)')
-
     plt.show()