diff --git a/.devpy/cmds.py b/.devpy/cmds.py index fe800010..38f8153b 100644 --- a/.devpy/cmds.py +++ b/.devpy/cmds.py @@ -1,6 +1,6 @@ import click from devpy import util -from devpy.cmds.util import get_site_packages, set_pythonpath, run +from devpy.cmds.meson import _get_site_packages @click.command() @@ -25,19 +25,19 @@ def install_dependencies(test_dep=False, doc_dep=False): config = util.get_config() default_dependencies = config["project"]['dependencies'] print("Installing dependencies", default_dependencies) - run( + util.run( ["pip", "install"] + list(default_dependencies), ) if test_dep: test_dependencies = config["project.optional-dependencies"]['test'] print("Installing test-dependencies", test_dependencies) - run( + util.run( ["pip", "install"] + list(test_dependencies), ) if doc_dep: doc_dependencies = config["project.optional-dependencies"]['doc'] print("Installing doc-dependencies", doc_dependencies) - run( + util.run( ["pip", "install"] + list(doc_dependencies), ) @@ -60,9 +60,9 @@ def codecov(build_dir, codecov_args): codecov parameters """ - site_path = get_site_packages(build_dir) + site_path = _get_site_packages() - run( + util.run( ["codecov"] + list(codecov_args), cwd=site_path, - ) \ No newline at end of file + ) diff --git a/.github/workflows/test_unit_and_examples.yml b/.github/workflows/test_unit_and_examples.yml index 5cdf0951..33c7524d 100644 --- a/.github/workflows/test_unit_and_examples.yml +++ b/.github/workflows/test_unit_and_examples.yml @@ -53,7 +53,7 @@ jobs: - name: Install Python dependecies run: | - pip install pytest meson-python ninja cython numpy git+https://github.com/scientific-python/devpy + pip install pytest meson-python ninja cython numpy git+https://github.com/scientific-python/devpy@v0.1 python3 -m devpy install-dependencies -test-dep - name: Install S4 @@ -125,7 +125,7 @@ jobs: - name: Install Python dependecies run: | - pip install pytest meson-python ninja cython numpy git+https://github.com/scientific-python/devpy + pip install pytest meson-python ninja cython numpy git+https://github.com/scientific-python/devpy@v0.1 python3 -m devpy install-dependencies -test-dep - name: Install S4 diff --git a/docs/source/Solvers/solving_solar_cells.rst b/docs/source/Solvers/solving_solar_cells.rst index ce9066d4..898a50d8 100755 --- a/docs/source/Solvers/solving_solar_cells.rst +++ b/docs/source/Solvers/solving_solar_cells.rst @@ -97,7 +97,7 @@ The options available as well as the default values are: If the Boltzmann approximation should be used in the detailed balance solver. Any other choice will result in using the full Plank equation, which will be slower, in general. - :doc:`Depletion approximation solver options ` - - **da_mode** = 'bvp' + - **da_mode** = 'green' Selects the numerical approximation method for the drift-diffusion equation in the depletion approximation solver. Possible values are “bvp” for numerical solution using the `solve_bvp` method of the `scipy.integrate` module or 'green' for a semi-analytic solution using Green's functions. The latter is expected to be faster. - :doc:`Poisson-drift diffusion solver options ` diff --git a/examples/MJ_solar_cell_external_reflectance.py b/examples/MJ_solar_cell_external_reflectance.py new file mode 100755 index 00000000..ab43515b --- /dev/null +++ b/examples/MJ_solar_cell_external_reflectance.py @@ -0,0 +1,125 @@ +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 + +all_materials = [] + +def this_dir_file(f): + return "data/" + f + +# 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) + +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) + +# BOTTOM CELL - Ge +Ge = material("Ge") + +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 + +# 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) +solar_cell = SolarCell( + [ + 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'), + 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'), + 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, + 'da_mode': 'green'}) + +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.legend() +plt.ylim(0, 100) +plt.ylabel('EQE (%)') +plt.xlabel('Wavelength (nm)') +plt.show() + +V = np.linspace(0, 3, 300) +solar_cell_solver(solar_cell, 'iv', user_options={'voltages': V, + 'light_iv': True, + 'wavelength': wl, + 'da_mode': 'green' + }) + +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.legend() +plt.ylim(0, 230) +plt.xlim(0, 3) +plt.ylabel('Current (A/m$^2$)') +plt.xlabel('Voltage (V)') + +plt.show() diff --git a/examples/MJ_solar_cell_using_DA.py b/examples/MJ_solar_cell_using_DA.py index fbcff9d8..948e0cdd 100755 --- a/examples/MJ_solar_cell_using_DA.py +++ b/examples/MJ_solar_cell_using_DA.py @@ -2,45 +2,48 @@ 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): - return "data/" + f +ARC_layers = [Layer(si("100nm"), material=MgF2), + Layer(si("50nm"), material=ZnS)] -# 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 - 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. -# 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")) +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")) + +# For convenience we will set the carrier transport properties to default values for +# all the materials in one place. To do that we add the materials to a list called +# all_materials: + +all_materials = [] 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) +# 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")) +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")) all_materials.append(mid_cell_n_material) all_materials.append(mid_cell_p_material) @@ -48,63 +51,110 @@ def this_dir_file(f): # BOTTOM CELL - Ge Ge = material("Ge") -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")) +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. +# We now set the electronic properties of these materials to some default values. +# Although in practice there will be differences, for this example we set them to +# reasonable generic values for simplicity. + 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 -# 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) +# 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. + 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'), + Layer(si("400nm"), material=top_cell_p_material, role='base'), ], sn=1, sp=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'), 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=1, sp=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': 'bvp', + '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'. 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 + len(ARC_layers)].eqe(wl) * 100, 'b', label='GaInP QE') +plt.plot(wl * 1e9, solar_cell[1 + len(ARC_layers)].eqe(wl) * 100, 'g', label='GaAs QE') +plt.plot(wl * 1e9, solar_cell[2 + len(ARC_layers)].eqe(wl) * 100, 'r', label='Ge QE') +plt.plot(wl * 1e9, solar_cell[0 + len(ARC_layers)].layer_absorption * 100, '--b', + label='GaInP Abs.') +plt.plot(wl * 1e9, solar_cell[1 + len(ARC_layers)].layer_absorption * 100, '--g', + label='GaAs Abs.') +plt.plot(wl * 1e9, solar_cell[2 + len(ARC_layers)].layer_absorption * 100, '--r', + label='Ge Abs.') + +plt.plot(wl*1e9, solar_cell.reflected*100, '--k', label="Reflected") plt.legend() plt.ylim(0, 100) plt.ylabel('EQE (%)') plt.xlabel('Wavelength (nm)') plt.show() +# Set up the AM0 (space) solar spectrum +am0 = LightSource(source_type='standard',version='AM0',x=wl, + output_units='photon_flux_per_m') + V = np.linspace(0, 3, 300) -solar_cell_solver(solar_cell, 'iv', user_options={'voltages': V, 'light_iv': True, 'wavelength': wl}) + +# Calculate the current-voltage relationship under illumination: + +solar_cell_solver(solar_cell, 'iv', user_options={'light_source': am0, + 'voltages': V, + 'light_iv': True, + 'wavelength': wl, + 'optics_method': 'TMM', + 'mpp': True, + # 'da_mode': 'bvp' + }) +# 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(V, solar_cell.iv['IV'][1]/10, 'k', linewidth=3, label='3J cell') +plt.plot(V, -solar_cell[0 + len(ARC_layers)].iv(V)/10, 'b', label='InGaP sub-cell') +plt.plot(V, -solar_cell[1 + len(ARC_layers)].iv(V)/10, 'g', label='GaAs sub-cell') +plt.plot(V, -solar_cell[2 + len(ARC_layers)].iv(V)/10, 'r', label='Ge sub-cell') +plt.text(0.5,30,f'Jsc= {solar_cell.iv.Isc/10:.2f} mA.cm' + r'$^{-2}$') +plt.text(0.5,28,f'Voc= {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, 230) +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() +plt.show() \ No newline at end of file diff --git a/examples/analytic_depletion_approximation.py b/examples/analytic_depletion_approximation.py new file mode 100644 index 00000000..a13a3c42 --- /dev/null +++ b/examples/analytic_depletion_approximation.py @@ -0,0 +1,461 @@ +import numpy as np +from solcore import material +from solcore.constants import q, kb +import matplotlib.pyplot as plt +from scipy.integrate import solve_bvp, quad_vec +from functools import partial +from solcore.interpolate import interp1d +import warnings +from pytest import approx + +GaAs = material("GaAs")() + +# make n and p the same +# width = 850e-9 +width_top = 300e-9 +width_bottom = 1000000e-9 + + +def get_J_sc_diffusion_vs_WL(xa, xb, g, D, L, y0, S, wl, ph, side="top"): + zz = np.linspace(xa, xb, 1001, endpoint=False) + # excluding the last point - depending on the mesh/floating point errors, + # sometimes this is actually in the next layer + + gg = g(zz) * ph # generation rate * photon flux + + out = np.zeros_like(wl) + sol_success = np.zeros_like(wl) # keep track of whether solve_bvp is converging + solutions = [] + + for i in range(len(wl)): + if np.all( + gg[:, i] == 0 + ): # no reason to solve anything if no generation at this wavelength + out[i] = 0 + sol_success[i] = 0 + + else: + + def A(x): + return np.interp(x, zz, gg[:, i]) / D + y0 / L**2 + + # generation and n0/p0 term in differential equation + # (eq. 6.15 & 6.20 in Jenny Nelson, Physics of Solar Cells) + + def fun(x, y): + # differential equation (eq. 6.15 & 6.20 in Jenny Nelson, + # Physics of Solar Cells) + # solve_bvp solves equation of form: + # dy / dx = f(x, y, p), a <= x <= b + # in this case y = [n or p, dn/dx or dp/dx] + # y[0] = carrier concentration (n or p) + # y[1] = carrier concentration gradient (dn/dx or dp/dx) + out1 = y[1] # by definition! dy/dx = dy/dx + + out2 = y[0] / L**2 - A(x) + # actually solving the differential equation (6.15 & 6.20) + + return np.vstack((out1, out2)) + + # boundary conditions for solve_bvp: + if side == "top": + + def bc(ya, yb): + left = ya[1] - S / D * (ya[0] - y0) + # eq. 6.18 - b.c. at front of junction (surface recombination) + + right = yb[0] - y0 + # eq. 6.17 - b.c. edge of depletion region, top half of junction + # added - y0 (generally very small so makes almost no difference) + return np.array([left, right]) + + else: + + def bc(ya, yb): + left = ya[0] - y0 + print("left", left, y0) + # eq. 6.21 - b.c. edge of depletion region, bottom half of junction + # added - y0 (generally very small so makes almost no difference) + + right = yb[1] + S / D * (yb[0] - y0) + # eq. 6.22 - b.c. at back of junction (surface recombination) + # changed sign! Current is going the other way + return np.array([left, right]) + + guess = y0 * np.ones((2, zz.size)) + guess[1] = np.zeros_like(guess[0]) + + solution = solve_bvp(fun, bc, zz, guess, max_nodes=2 * zz.shape[0]) + # increase max_nodes to avoid "too many mesh points" message + + sol_success[i] = solution.status + + if side == "top": + out[i] = solution.y[1][-1] + # current at edge of depletion region (top half of junction), eq. 6.33 + else: + out[i] = solution.y[1][0] + # current at edge of depletion region (bottom half of junction), eq 6.38 + + solutions.append(solution) + + # give a warning f any of the solution statuses are not 0 using warnings.warn: + if np.any(sol_success != 0): + warnings.warn( + "Depletion approximation (DA) EQE calculation: " + "solve_bvp did not converge as expected for some wavelengths", + RuntimeWarning, + ) + + return out, solutions + + +def get_J_sc_diffusion_green(xa, xb, g, D, L, y0, S, ph, side="top"): + """Computes the derivative of the minority carrier concentration at the edge of the + junction by approximating the convolution integral resulting from applying the + Green's function method to the drift-diffusion equation. + + :param xa: Coordinate at the start the junction. + :param xb: Coordinate at the end the junction. + :param g: Carrier generation rate at point x (expected as function). + :param D: Diffusion constant. + :param L: Diffusion length. + :param y0: Carrier equilibrium density. + :param S: Surface recombination velocity. + :param ph: Light spectrum. + :param side: String to indicate the edge of interest. Either 'top' or 'bottom'. + + :return: The derivative of the minority carrier concentration at the edge of the + junction. + """ + + xbL = (xb - xa) / L + crvel = S / D * L + ph_over_D = ph / D + + print("xbL", xbL) + + # if L too low in comparison to junction width, avoid nan's + if xbL > 1.0e2: + print("low L") + if side == "top": + fun = partial(_conv_exp_top, xa=xa, xb=xb, g=g, L=L, phoD=ph_over_D) + else: + fun = partial(_conv_exp_bottom, xa=xa, g=g, L=L, phoD=ph_over_D) + cp = 1.0 + + else: + if side == "top": + cp = -np.cosh(xbL) - crvel * np.sinh(xbL) + + fun = partial( + _conv_green_top, xa=xa, xb=xb, g=g, L=L, phoD=ph_over_D, crvel=crvel + ) + else: + cp = np.cosh(xbL) + crvel * np.sinh(xbL) + + fun = partial( + _conv_green_bottom, xb=xb, g=g, L=L, phoD=ph_over_D, crvel=-crvel + ) + + out, err, info = quad_vec(fun, xa, xb, epsrel=1.0e-5, full_output=True) + # print(info) + return out.squeeze() / cp + + +def _conv_exp_top(x, xa, xb, g, L, phoD): + """Convolution of the carrier generation rate with the approximate Green's function + kernel at point x. To be used with the numerical integration routine to compute the + minority carrier derivative on the top edge. This kernel approximates the original + one when the diffusion length is 2 orders of magnitude higher than the junction + width by assuming that sinh(x) = cosh(x) = .5 * exp(x). + + :param x: Coordinate in the junction (variable to be integrated). + :param xa: Coordinate at the start the junction. + :param xb: Coordinate at the end the junction. + :param g: Carrier generation rate at point x (expected as function). + :param L: Diffusion length. + :param phoD: Light spectrum divided by the diffusion constant D. + """ + xc = (xa - x) / L + xv = np.array( + [ + xa + xb - x, + ] + ) + Pkern = -np.exp(xc) + Gx = g(xv) * phoD + return Pkern * Gx + + +def _conv_exp_bottom(x, xa, g, L, phoD): + """Convolution of the carrier generation rate with the approximate Green's function + kernel at point x. To be used with the numerical integration routine to compute the + minority carrier derivative on the bottom edge. This kernel approximates the + original one when the diffusion length is 2 orders of magnitude higher than + the junction width by assuming that sinh(x) = cosh(x) = .5 * exp(x). + + :param x: Coordinate in the junction (variable to be integrated). + :param xa: Coordinate at the start the junction. + :param g: Carrier generation rate at point x (expected as function). + :param L: Diffusion length. + :param phoD: Light spectrum divided by the diffusion constant D. + """ + xc = (xa - x) / L + xv = np.array( + [ + x, + ] + ) + Pkern = np.exp(xc) + Gx = g(xv) * phoD + return Pkern * Gx + + +def _conv_green_top(x, xa, xb, g, L, phoD, crvel): + """Convolution of the carrier generation rate with the Green's function kernel at + point x. To be used with the numerical integration routine to compute the minority + carrier derivative on the top edge. + + :param x: Coordinate in the junction (variable to be integrated). + :param xa: Coordinate at the start the junction. + :param xb: Coordinate at the end the junction. + :param g: Carrier generation rate at point x (expected as function). + :param L: Diffusion length. + :param phoD: Light spectrum divided by the diffusion constant D. + :param crvel: Coefficient computed as S / D * L, with S the surface recombination + velocity. + """ + xc = (xb - x) / L + xv = np.array( + [ + xa + xb - x, + ] + ) + Pkern = np.cosh(xc) + crvel * np.sinh(xc) + Gx = g(xv) * phoD + return Pkern * Gx + + +def _conv_green_bottom(x, xb, g, L, phoD, crvel): + """Convolution of the carrier generation rate with the Green's function kernel at + point x. To be used with the numerical integration routine to compute the minority + carrier derivative on the bottom edge. + + :param x: Coordinate in the junction (variable to be integrated). + :param xb: Coordinate at the end the junction. + :param g: Carrier generation rate at point x (expected as function). + :param L: Diffusion length. + :param phoD: Light spectrum divided by the diffusion constant D. + :param crvel: Coefficient computed as S / D * L, with S the surface recombination + velocity. + """ + xc = (xb - x) / L + xv = np.array( + [ + x, + ] + ) + Pkern = np.cosh(xc) - crvel * np.sinh(xc) + Gx = g(xv) * phoD + return Pkern * Gx + + +bs = 1e20 + +Dn, Ln, Sn = 0.001263, 7e-06, 1000 + +alpha = GaAs.alpha(800e-9) +alpha = 1000 + + +V = 0 + +T = 273 + +kT = kb * T + +wp = 100e-9 +R = 0 + +xp = 0 + +Na = 6e16 * 1e6 +Nd = 3e17 * 1e6 + +ni = GaAs.ni +xi = 0 + +Dp, Lp, Sp = 0.001263, 7e-06, 400 +p0 = ni**2 / Nd + +es = GaAs.permittivity + +Vbi = (kT / q) * np.log(Nd * Na / ni**2) + +wn = (-xi + np.sqrt(xi**2 + 2.0 * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / ( + 1 + Nd / Na +) +wp = (-xi + np.sqrt(xi**2 + 2.0 * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / ( + 1 + Na / Nd +) + +n0 = ni**2 / Na + +wl = np.array([400, 500]) * 1e-9 +alpha_vec = np.array([alpha, alpha * 1.1]) +zz = np.linspace(0, width_top + width_bottom, 2001) +generation = interp1d( + zz, (1 - R) * alpha_vec * np.exp(-alpha_vec * zz[:, None]), axis=0 +) + +all_sols_n = get_J_sc_diffusion_vs_WL( + 0, width_top - wp, generation, Dn, Ln, n0, Sn, wl, bs, "top" +)[1] +all_sols_p = get_J_sc_diffusion_vs_WL( + width_top + wn, + width_top + width_bottom, + generation, + Dp, + Lp, + p0, + Sp, + wl, + bs, + "bottom", +)[1] +n_sol_solcore = all_sols_n[0].y[0] +p_sol_solcore = all_sols_p[0].y[0] + +xbb = width_top - wp - (width_top - wp) / 1001.0 + +all_sols_n_green = get_J_sc_diffusion_green( + 0, xbb, generation, Dn, Ln, n0, Sn, bs, "top" +)[0] + +xbb = width_top + width_bottom - (width_top + width_bottom - width_top + wn) / 1001.0 +all_sols_p_green = get_J_sc_diffusion_green( + width_top + wn, xbb, generation, Dp, Lp, p0, Sp, bs, "bottom" +)[0] + + +def n_mathematica(x, L, alpha, xa, D, S, n0, A): + aL = alpha * L + a1 = -aL * np.exp(2 * alpha * (x + xa) + (2 * x / L)) + aL * np.exp( + 2 * alpha * (x + xa) + (2 * xa / L) + ) + + a2 = np.exp(2 * alpha * x + alpha * xa + (xa / L)) - np.exp( + alpha * x + 2 * alpha * xa + (x / L) + ) + + a3 = np.exp(((aL + 1) * (2 * x + xa)) / L) - np.exp(((aL + 1) * (x + 2 * xa)) / L) + + b1 = np.exp(((aL + 1) * (2 * x + xa)) / L) - np.exp(((aL + 1) * (x + 2 * xa)) / L) + + b2 = np.exp(alpha * x + 2 * alpha * xa + (x / L)) - np.exp( + 2 * alpha * x + alpha * xa + (xa / L) + ) + + b3 = -np.exp(2 * alpha * (x + xa) + (2 * x / L)) + np.exp( + 2 * alpha * (x + xa) + (2 * xa / L) + ) + + numerator = ( + A + * L**2 + * np.exp(-2 * alpha * (x + xa) - (x / L)) + * (D * (a1 + a2 + a3) + L * S * (b1 + b2 + b3)) + ) + denominator = ( + D + * (alpha**2 * L**2 - 1) + * (D * (np.exp(2 * xa / L) + 1) + L * S * (np.exp(2 * xa / L) - 1)) + ) + result = numerator / denominator + n0 + + return result + + +def p_mathematica(x, L, alpha, xb, D, S, p0, A): + aL = alpha * L + denominator = ( + D + * (-1 + alpha**2 * L**2) + * (D * (1 + np.exp((2 * xb) / L)) + (-1 + np.exp((2 * xb) / L)) * L * S) + ) + + a1 = np.exp((2 * xb) / L + 2 * alpha * (x + xb)) - np.exp( + (((1 + aL) * (x + 2 * xb)) / L) + ) + + a2 = np.exp((2 * x) / L + 2 * alpha * (x + xb)) - np.exp( + alpha * x + x / L + 2 * alpha * xb + ) + + a3 = ( + alpha * np.exp(2 * alpha * x + alpha * xb + xb / L) * L + - alpha * np.exp(((1 + aL) * (2 * x + xb)) / L) * L + ) + + b1 = np.exp(((1 + aL) * (2 * x + xb)) / L) - np.exp(((1 + aL) * (x + 2 * xb)) / L) + + b2 = +np.exp(alpha * x + x / L + 2 * alpha * xb) - np.exp( + 2 * alpha * x + alpha * xb + xb / L + ) + + b3 = np.exp((2 * xb) / L + 2 * alpha * (x + xb)) - np.exp( + (2 * x) / L + 2 * alpha * (x + xb) + ) + + numerator = ( + A + * np.exp(-(x / L) - 2 * alpha * (x + xb)) + * L**2 + * (D * (a1 + a2 + a3) + (b1 + b2 + b3) * L * S) + ) + + p = p0 + numerator / denominator + return p + + +x_0 = np.linspace(0, width_top - wp, 1001, endpoint=False) + +pref = (1 - R) * alpha * bs + +n_ma = n_mathematica(x_0, Ln, alpha, (width_top - wp), Dn, Sn, n0, pref) + +pref = (1 - R) * alpha * bs * np.exp(-alpha * (width_top + wn)) +x_1 = np.linspace(0, width_bottom - wn, 1001, endpoint=False) + +p_ma = p_mathematica(x_1, Lp, alpha, width_bottom - wn, Dp, Sp, p0, pref) + +plt.figure() +plt.semilogy(x_0 * 1e9, n_ma) +plt.semilogy(x_0 * 1e9, n_sol_solcore, "--") +plt.tight_layout() +plt.show() + +plt.figure() +plt.semilogy(x_1 * 1e9, p_ma) +plt.semilogy(x_1 * 1e9, p_sol_solcore, "--") +plt.tight_layout() +plt.show() + +# plt.figure() +# plt.plot(x_1*1e9, p_ma/p_sol_solcore) +# plt.show() + +print(n_ma == approx(n_sol_solcore, rel=0.02)) +print(p_ma == approx(p_sol_solcore, rel=0.02)) + +plt.figure() +plt.semilogy(x_0 * 1e9, all_sols_n[0].y[1]) +plt.show() + +plt.figure() +plt.semilogy(x_1 * 1e9, all_sols_p[0].y[1]) +plt.show() + +print(all_sols_n[0].y[1][-1], all_sols_n_green) +print(all_sols_p[0].y[1][0], all_sols_p_green) diff --git a/pyproject.toml b/pyproject.toml index 99b138a2..487d78c8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -67,7 +67,7 @@ dev = ["pre-commit"] package = 'solcore' [tool.devpy.commands] -"Build" = ["devpy.build", "devpy.test"] +"Build" = ["devpy.cmds.meson.build", "devpy.cmds.meson.test"] "Extensions" = ['.devpy/cmds.py:codecov', '.devpy/cmds.py:install_dependencies'] [tool.pytest.ini_options] diff --git a/solcore/absorption_calculator/transfer_matrix.py b/solcore/absorption_calculator/transfer_matrix.py index 2f0aad74..f70849fd 100755 --- a/solcore/absorption_calculator/transfer_matrix.py +++ b/solcore/absorption_calculator/transfer_matrix.py @@ -1,5 +1,6 @@ -""" This module serves as interface between solcore structures (junctions, layers, materials...) and the -transfer matrix package developed by Steven Byrnes and included in the PyPi repository. +""" This module serves as interface between solcore structures (junctions, layers, +materials...) and the transfer matrix package developed by Steven Byrnes and included in +the PyPi repository. """ import numpy as np @@ -17,16 +18,20 @@ def np_cache(function): - """This function was taken from https://stackoverflow.com/questions/52331944/cache-decorator-for-numpy-arrays/52332109#52332109 - - Creates a cacheable version of a function which takes a 1D numpy array as input, by using a wrapping function - which converts the array to a tuple. It returns a function which returns the same output as the input function, - but can be cached, avoiding a bottleneck when optical constants in a material are looked up repeatedly. - - If the input argument is hashable already - eg. a single float number - then the function is not chached and is returned "as is". + """This function was taken from + https://stackoverflow.com/questions/52331944/cache-decorator-for-numpy-arrays/52332109#52332109 # noqa: E501 + + Creates a cacheable version of a function which takes a 1D numpy array as input, by + using a wrapping function which converts the array to a tuple. It returns a + function which returns the same output as the input function, but can be cached, + avoiding a bottleneck when optical constants in a material are looked up repeatedly. + + If the input argument is hashable already - eg. a single float number - then the + function is not chached and is returned "as is". :param: function: the function of which a cacheable version is to be created :return: wrapper: the cacheable version of the function""" + @lru_cache() def cached_wrapper(hashable_array): array = np.array(hashable_array) @@ -47,15 +52,18 @@ def wrapper(array): class OptiStack(object): - """ Class that contains an optical structure: a sequence of layers with a thickness and a complex refractive index. + """Class that contains an optical structure: a sequence of layers with a thickness + and a complex refractive index. - It serves as an intermediate step between solcore layers and materials and the stack of thicknesses and - and n and k values necessary to run calculations involving TMM. When creating an OptiStack object, the thicknesses - of all the layers forming the Solcore structure and the optical data of the materials of the layers are extracted + It serves as an intermediate step between solcore layers and materials and the + stack of thicknesses and n and k values necessary to run calculations involving TMM. + When creating an OptiStack object, the thicknesses of all the layers forming the + Solcore structure and the optical data of the materials of the layers are extracted and arranged in such a way they can be easily and fastly read by the TMM functions. - In addition to a solcore structure with Layers, it can also take a list where each element represent a layer - written as a list and contains the layer thickness and the dielectrical model, the raw n and k data as a function + In addition to a solcore structure with Layers, it can also take a list where each + element represent a layer written as a list and contains the layer thickness and + the dielectrical model, the raw n and k data as a function of wavelengths, or a whole Device structure as the type used in the PDD model. In summary, this class acepts: @@ -69,40 +77,55 @@ class OptiStack(object): solcore.Layer, solcore.Layer ] - This allows for maximum flexibility when creating the optical model, allowing to construct the stack with - experimental data, modelled data and known material properties from the database. + This allows for maximum flexibility when creating the optical model, allowing to + construct the stack with experimental data, modelled data and known material + properties from the database. - Yet anther way of defining the layers mixes experimental data with a DielectricModel within the same layer but in - spectrally distinct regions. The syntaxis for the layer is: + Yet anther way of defining the layers mixes experimental data with a DielectricModel + within the same layer but in spectrally distinct regions. The syntaxis for the + layer is: layer = [thickness, wavelength, n, k, DielectricModel, mixing] - where mixing is a list containing three elements: [the mixing point (nm), the mixing width (nm), zero or one] - depending if the mixing function should be increasing with the wavelength or decreasing. If increasing (zero), the - Dielectric model will be used at long wavelengths and the experimental data at short wavelengths. If decreasing - (one) the oposite is done. The mixing point and mixing width control how smooth is the transition between one and - the other type of data. - - Extra layers such as he semi-infinite, air-like first and last medium, and a back highly absorbing layer are - included at runtime to fulfill the requirements of the TMM solver or to solve some of its limitations. + where mixing is a list containing three elements: + [the mixing point (nm), the mixing width (nm), zero or one] + depending if the mixing function should be increasing with the wavelength or + decreasing. If increasing (zero), the Dielectric model will be used at long + wavelengths and the experimental data at short wavelengths. If decreasing (one) + the oposite is done. The mixing point and mixing width control how smooth is the + transition between one and the other type of data. + + Extra layers such as he semi-infinite, air-like first and last medium, and a back + highly absorbing layer are included at runtime to fulfill the requirements of the + TMM solver or to solve some of its limitations. """ - def __init__(self, structure=(), no_back_reflection=False, substrate=None, incidence=None, **kwargs): - """ Class constructor. It takes a Solcore structure and extract the thickness and optical data from the - Layers and the materials. Option is given to indicate if the reflexion from the back of the structure must be - supressed, usefull for ellipsometry calculations. This is done by creating an artificial highly absorbing but - not reflecting layer just at the back. - - Alternativelly, it can take simply a list of [thickness, DielectricModel] or [thickness, wavelength, n, k] for - each layer accounting for the refractive index of the layers. The three options can be mixed for maximum - flexibility. + def __init__( + self, + structure=(), + no_back_reflection=False, + substrate=None, + incidence=None, + **kwargs, + ): + """Class constructor. It takes a Solcore structure and extract the thickness + and optical data from the Layers and the materials. Option is given to + indicate if the reflection from the back of the structure must be suppressed, + useful for ellipsometry calculations. This is done by creating an artificial + highly absorbing but not reflecting layer just at the back. + + Alternatively, it can take simply a list of [thickness, DielectricModel] or + [thickness, wavelength, n, k] for each layer accounting for the refractive + index of the layers. The three options can be mixed for maximum flexibility. :param structure: A list with one or more layers. - :param no_back_reflection: If reflection from the back must be suppressed. Default=False. - :param substrate: a semi-infinite transmission medium. Note that if no_back_reflection is set to True, - adding a substrate won't make any difference. - :param incidence: a semi-infinite transmission medium. Note that add things may happen if - this material is absorbing + :param no_back_reflection: If reflection from the back must be suppressed. + Default=False. + :param substrate: a semi-infinite transmission medium. Note that if + no_back_reflection is set to True, adding a substrate won't make any + difference. + :param incidence: a semi-infinite incidence medium. Note that add things may + happen if this material is absorbing """ self.widths = [] @@ -118,34 +141,41 @@ def __init__(self, structure=(), no_back_reflection=False, substrate=None, incid self.no_back_reflection = no_back_reflection - if 'no_back_reflexion' in kwargs: - warn('The no_back_reflexion warning is deprecated. Use no_back_reflection instead.', FutureWarning) - self.no_back_reflection = kwargs['no_back_reflexion'] + if "no_back_reflexion" in kwargs: + warn( + "The no_back_reflexion warning is deprecated. " + "Use no_back_reflection instead.", + FutureWarning, + ) + self.no_back_reflection = kwargs["no_back_reflexion"] def get_indices(self, wl): - """ Returns the complex refractive index of the stack. + """Returns the complex refractive index of the stack. :param wl: Wavelength of the light in nm. - :return: A list with the complex refractive index of each layer, including the semi-infinite front and back - layers and, opionally, the back absorbing layer used to suppress back surface relfexion. + :return: A list with the complex refractive index of each layer, including the + semi-infinite front and back layers and, opionally, the back absorbing layer + used to suppress back surface relfexion. """ out = [] - wl_m = solcore.si(wl, 'nm') + wl_m = solcore.si(wl, "nm") - if hasattr(self, 'n_sub'): - n1 = self.n_sub(wl_m) + self.k_sub(wl_m)*1.0j + if hasattr(self, "n_sub"): + n1 = self.n_sub(wl_m) + self.k_sub(wl_m) * 1.0j else: - if hasattr(wl, 'size'): + if hasattr(wl, "size"): n1 = np.ones_like(wl, dtype=complex) else: n1 = 1 # incidence medium! - if hasattr(self, 'n_sup'): - n0 = self.n_sup(wl_m) #+ self.k_sup(wl_m)*1.0j ignore complex part for now to avoid errors + if hasattr(self, "n_sup"): + n0 = self.n_sup( + wl_m + ) # + self.k_sup(wl_m)*1.0j ignore complex part for now to avoid errors else: - if hasattr(wl, 'size'): + if hasattr(wl, "size"): n0 = np.ones_like(wl, dtype=complex) else: n0 = 1 @@ -155,19 +185,22 @@ def get_indices(self, wl): # substrate irrelevant if no_back_reflection = True if self.no_back_reflection: - return [n0] + out + [self.n_data[-1](wl_m) + self._k_absorbing(wl_m) * 1.0j, n1] # look at last entry in stack, + return ( + [n0] + + out + + [self.n_data[-1](wl_m) + self._k_absorbing(wl_m) * 1.0j, n1] + ) # look at last entry in stack, # make high;y absorbing layer based on it. else: - return [n0] + out + [n1] - def get_widths(self): - """ Returns the widths of the layers of the stack. + """Returns the widths of the layers of the stack. - :return: A list with the widths each layer, including the semi-infinite front and back layers and, opionally, - the back absorbing layer used to suppress back surface relfexion, defined as 1 mm thick. + :return: A list with the widths each layer, including the semi-infinite front + and back layers and, opionally, the back absorbing layer used to suppress + back surface relfection, defined as 1 mm thick. """ if self.no_back_reflection: @@ -176,8 +209,9 @@ def get_widths(self): return [np.inf] + self.widths + [np.inf] def _k_absorbing(self, wl): - """ k value of the back highly absorbing layer. It is the maximum between the bottom layer of the stack or a - finite, small value that will absorb all light within the absorbing layer thickness. + """k value of the back highly absorbing layer. It is the maximum between the + bottom layer of the stack or a finite, small value that will absorb all light + within the absorbing layer thickness. :param wl: Wavelength of the light in nm. :return: The k value at each wavelength. @@ -186,24 +220,26 @@ def _k_absorbing(self, wl): @staticmethod def _k_dummy(wl): - """ Dummy k value to be used with the dielectric model, which produces the refractive index as a complex - number. + """Dummy k value to be used with the dielectric model, which produces the + refractive index as a complex number. :param wl: Wavelength of the light in nm. :return: The k value at each wavelength... set to zero. """ - return 0. + return 0.0 def add_layers(self, layers): - """ Generic function to add layers to the OptiStack. Internally, it calls add_solcore_layer, - add_modelled_layer or add_raw_nk_layer. + """Generic function to add layers to the OptiStack. Internally, it calls + add_solcore_layer, add_modelled_layer or add_raw_nk_layer. - :param layers: A list with the layers to add (even if it is just one layer) It can be one or more and it can - mixed, Solcore-based and modelled layers. + :param layers: A list with the layers to add (even if it is just one layer) + It can be one or more and it can be mixed, Solcore-based and modelled + layers. :return: None """ try: - # If the input is a whole device structure, we get just the layers information + # If the input is a whole device structure, we get just the layers + # information if type(layers) is dict: layers = ToStructure(layers) @@ -217,9 +253,9 @@ def add_layers(self, layers): for layer in layers: self.layers.append(layer) - if 'Layer' in str(type(layer)): + if "Layer" in str(type(layer)): self._add_solcore_layer(layer) - elif 'DielectricConstantModel' in str(type(layer[1])): + elif "DielectricConstantModel" in str(type(layer[1])): self._add_modelled_layer(layer) else: self._add_raw_nk_layer(layer) @@ -227,16 +263,20 @@ def add_layers(self, layers): self.num_layers += 1 except Exception as err: - print('Error when adding a new layer to the OptiStack. {}'.format(err)) + print("Error when adding a new layer to the OptiStack. {}".format(err)) def remove_layer(self, idx): - """ Removes layer with index idx from the OptiStack + """Removes layer with index idx from the OptiStack :param idx: Index of the layer to remove :return: None """ if idx >= self.num_layers: - print('Error when removing layers. idx must be: 0 <= idx <= {}.'.format(self.num_layers - 1)) + print( + "Error when removing layers. idx must be: 0 <= idx <= {}.".format( + self.num_layers - 1 + ) + ) return self.widths.pop(idx) @@ -246,14 +286,17 @@ def remove_layer(self, idx): self.num_layers -= 1 def swap_layers(self, idx1, idx2): - """ Swaps two layers in the OptiStack. + """Swaps two layers in the OptiStack. :param idx1: The index of one of the layers. :param idx2: The index of the other. :return: None """ if idx1 >= self.num_layers or idx2 >= self.num_layers: - print('Error when removing layers. idx must be: 0 <= idx1, idx2 <= {}.'.format(self.num_layers - 1)) + print( + "Error when removing layers. idx must be: 0 <= idx1, idx2 <= {}.". + format(self.num_layers - 1) + ) return self.widths[idx1], self.widths[idx2] = self.widths[idx2], self.widths[idx1] @@ -261,7 +304,6 @@ def swap_layers(self, idx1, idx2): self.n_data[idx1], self.n_data[idx2] = self.n_data[idx2], self.n_data[idx1] self.k_data[idx1], self.k_data[idx2] = self.k_data[idx2], self.k_data[idx1] - def set_widths(self, widths): """Changes the widths of the layers in the stack. @@ -271,28 +313,28 @@ def set_widths(self, widths): if type(widths) is np.ndarray: widths = widths.tolist() - assert len(widths) == self.num_layers, \ - 'Error: The list of widths must have as many elements (now {}) as the ' \ - 'number of layers (now {}).'.format(len(widths), self.num_layers) + assert len(widths) == self.num_layers, ( + "Error: The list of widths must have as many elements (now {}) as the " + "number of layers (now {}).".format(len(widths), self.num_layers) + ) self.widths = widths - def _add_solcore_layer(self, layer): - """ Adds a Solcore layer to the end (bottom) of the stack, extracting its thickness and n and k data. + """Adds a Solcore layer to the end (bottom) of the stack, extracting its + thickness and n and k data. :param layer: The Solcore layer :return: None """ - self.widths.append(solcore.asUnit(layer.width, 'nm')) + self.widths.append(solcore.asUnit(layer.width, "nm")) self.models.append([]) self.n_data.append(np_cache(layer.material.n)) self.k_data.append(np_cache(layer.material.k)) - def _add_modelled_layer(self, layer): - """ Adds a layer to the end (bottom) of the stack. The layer must be defined as a list containing the layer - thickness in nm and a dielectric model. + """Adds a layer to the end (bottom) of the stack. The layer must be defined as + a list containing the layer thickness in nm and a dielectric model. :param layer: The new layer to add as [thickness, DielectricModel] :return: None @@ -303,13 +345,15 @@ def _add_modelled_layer(self, layer): self.k_data.append(np_cache(self._k_dummy)) def _add_raw_nk_layer(self, layer): - """ Adds a layer to the end (bottom) of the stack. The layer must be defined as a list containing the layer - thickness in nm, the wavelength, the n and the k data as array-like objects. + """Adds a layer to the end (bottom) of the stack. The layer must be defined as + a list containing the layer thickness in nm, the wavelength, the n and the k + data as array-like objects. :param layer: The new layer to add as [thickness, wavelength, n, k] :return: None """ - # We make sure that the wavelengths are increasing, revering the arrays otherwise. + # We make sure that the wavelengths are increasing, reversing the arrays + # otherwise. if layer[1][0] > layer[1][-1]: layer[1] = layer[1][::-1] layer[2] = layer[2][::-1] @@ -319,52 +363,93 @@ def _add_raw_nk_layer(self, layer): if len(layer) >= 5: self.models.append(layer[4]) - c = solcore.si(layer[5][0], 'nm') - w = solcore.si(layer[5][1], 'nm') + c = solcore.si(layer[5][0], "nm") + w = solcore.si(layer[5][1], "nm") d = layer[5][2] # = 0 for increasing, =1 for decreasing def mix(x): - out = 1 + np.exp(-(x - c) / w) out = d + (-1) ** d * 1 / out return out - n_data = np_cache(lambda x: self.models[-1].n_and_k(x) * mix(x) + (1 - mix(x)) * interp1d( - x=solcore.si(layer[1], 'nm'), y=layer[2], fill_value=layer[2][-1])(x)) - k_data = np_cache(lambda x: interp1d(x=solcore.si(layer[1], 'nm'), y=layer[3], fill_value=layer[3][-1])(x)) + n_data = np_cache( + lambda x: self.models[-1].n_and_k(x) * mix(x) + + (1 - mix(x)) + * interp1d( + x=solcore.si(layer[1], "nm"), y=layer[2], fill_value=layer[2][-1] + )(x) + ) + k_data = np_cache( + lambda x: interp1d( + x=solcore.si(layer[1], "nm"), y=layer[3], fill_value=layer[3][-1] + )(x) + ) self.n_data.append(n_data) self.k_data.append(k_data) else: self.models.append([]) - self.n_data.append(np_cache(interp1d(x=solcore.si(layer[1], 'nm'), y=layer[2], fill_value=layer[2][-1]))) - self.k_data.append(np_cache(interp1d(x=solcore.si(layer[1], 'nm'), y=layer[3], fill_value=layer[3][-1]))) - - -def calculate_rat(structure, wavelength, angle=0, pol='u', - coherent=True, coherency_list=None, no_back_reflection=True, **kwargs): - """ Calculates the reflected, absorbed and transmitted intensity of the structure for the wavelengths and angles - defined. - - :param structure: A solcore Structure object with layers and materials or a OptiStack object. + self.n_data.append( + np_cache( + interp1d( + x=solcore.si(layer[1], "nm"), + y=layer[2], + fill_value=layer[2][-1], + ) + ) + ) + self.k_data.append( + np_cache( + interp1d( + x=solcore.si(layer[1], "nm"), + y=layer[3], + fill_value=layer[3][-1], + ) + ) + ) + + +def calculate_rat( + structure, + wavelength, + angle=0, + pol="u", + coherent=True, + coherency_list=None, + no_back_reflection=True, + **kwargs, +): + """Calculates the reflected, absorbed and transmitted intensity of the structure + for the wavelengths and angles defined. + + :param structure: A solcore Structure object with layers and materials or a + OptiStack object. :param wavelength: Wavelengths (in nm) in which calculate the data. An array. - :param angle: Angle (in degrees) of the incident light. Default: 0 (normal incidence). + :param angle: Angle (in degrees) of the incident light. + Default: 0 (normal incidence). :param pol: Polarisation of the light: 's', 'p' or 'u'. Default: 'u' (unpolarised). - :param coherent: If the light is coherent or not. If not, a coherency list must be added. - :param coherency_list: A list indicating in which layers light should be treated as coeherent ('c') and in which - incoherent ('i'). It needs as many elements as layers in the structure. - :param no_back_reflection: If reflexion from the back must be supressed. Default=True. + :param coherent: If the light is coherent or not. If not, a coherency list + must be added. + :param coherency_list: A list indicating in which layers light should be treated as + coherent ('c') and in which incoherent ('i'). It needs as many elements as + layers in the structure. + :param no_back_reflection: If reflection from the back must be supressed. + Default=True. :return: A dictionary with the R, A and T at the specified wavelengths and angle. """ num_wl = len(wavelength) - if 'no_back_reflexion' in kwargs: - warn('The no_back_reflexion argument is deprecated. Use no_back_reflection instead.', FutureWarning) - no_back_reflection = kwargs['no_back_reflexion'] + if "no_back_reflexion" in kwargs: + warn( + "The no_back_reflexion argument is deprecated. " + "Use no_back_reflection instead.", + FutureWarning, + ) + no_back_reflection = kwargs["no_back_reflexion"] - if 'OptiStack' in str(type(structure)): + if "OptiStack" in str(type(structure)): stack = structure stack.no_back_reflection = no_back_reflection else: @@ -372,173 +457,276 @@ def calculate_rat(structure, wavelength, angle=0, pol='u', if not coherent: if coherency_list is not None: - assert len(coherency_list) == stack.num_layers, \ - 'Error: The coherency list must have as many elements (now {}) as the ' \ - 'number of layers (now {}).'.format(len(coherency_list), stack.num_layers) + assert len(coherency_list) == stack.num_layers, ( + "Error: The coherency list must have as many elements (now {}) as the " + "number of layers (now {}).".format( + len(coherency_list), stack.num_layers + ) + ) if stack.no_back_reflection: - coherency_list = ['i'] + coherency_list + ['i', 'i'] + coherency_list = ["i"] + coherency_list + ["i", "i"] else: - coherency_list = ['i'] + coherency_list + ['i'] + coherency_list = ["i"] + coherency_list + ["i"] else: - raise Exception('Error: For incoherent or partly incoherent calculations you must supply the ' - 'coherency_list parameter with as many elements as the number of layers in the ' - 'structure') - - output = {'R': np.zeros(num_wl), 'A': np.zeros(num_wl), 'T': np.zeros(num_wl), 'all_p': [], 'all_s': [], 'out': [], - 'out_s': [], 'out_p': [], 'A_per_layer_s': [], 'A_per_layer_p': []} - - if pol in 'sp': + raise Exception( + "Error: For incoherent or partly incoherent calculations you must " + "supply the coherency_list parameter with as many elements as the " + "number of layers in the structure" + ) + + output = { + "R": np.zeros(num_wl), + "A": np.zeros(num_wl), + "T": np.zeros(num_wl), + "all_p": [], + "all_s": [], + "out": [], + "out_s": [], + "out_p": [], + "A_per_layer_s": [], + "A_per_layer_p": [], + } + + if pol in "sp": if coherent: - out = tmm.coh_tmm(pol, stack.get_indices(wavelength), stack.get_widths(), angle * degree, wavelength) + out = tmm.coh_tmm( + pol, + stack.get_indices(wavelength), + stack.get_widths(), + angle * degree, + wavelength, + ) A_per_layer = tmm.absorp_in_each_layer(out) - output['R'] = out['R'] - output['A'] = 1 - out['R'] - out['T'] - output['T'] = out['T'] - output['A_per_layer'] = A_per_layer - output['out'] = out + output["R"] = out["R"] + output["A"] = 1 - out["R"] - out["T"] + output["T"] = out["T"] + output["A_per_layer"] = A_per_layer + output["out"] = out else: - out = tmm.inc_tmm(pol, stack.get_indices(wavelength), stack.get_widths(), coherency_list, angle * degree, wavelength) + out = tmm.inc_tmm( + pol, + stack.get_indices(wavelength), + stack.get_widths(), + coherency_list, + angle * degree, + wavelength, + ) A_per_layer = np.array(tmm.inc_absorp_in_each_layer(out)) - output['R'] = out['R'] - output['A'] = 1 - out['R'] - out['T'] - output['T'] = out['T'] - output['A_per_layer'] = A_per_layer - output['out'] = out + output["R"] = out["R"] + output["A"] = 1 - out["R"] - out["T"] + output["T"] = out["T"] + output["A_per_layer"] = A_per_layer + output["out"] = out else: if coherent: - out_p = tmm.coh_tmm('p', stack.get_indices(wavelength), stack.get_widths(), angle * degree, wavelength) - out_s = tmm.coh_tmm('s', stack.get_indices(wavelength), stack.get_widths(), angle * degree, wavelength) + out_p = tmm.coh_tmm( + "p", + stack.get_indices(wavelength), + stack.get_widths(), + angle * degree, + wavelength, + ) + out_s = tmm.coh_tmm( + "s", + stack.get_indices(wavelength), + stack.get_widths(), + angle * degree, + wavelength, + ) A_per_layer_p = tmm.absorp_in_each_layer(out_p) A_per_layer_s = tmm.absorp_in_each_layer(out_s) - output['R'] = 0.5 * (out_p['R'] + out_s['R']) - output['T'] = 0.5 * (out_p['T'] + out_s['T']) - output['A'] = 1 - output['R'] - output['T'] - output['A_per_layer'] = 0.5*(A_per_layer_p + A_per_layer_s) + output["R"] = 0.5 * (out_p["R"] + out_s["R"]) + output["T"] = 0.5 * (out_p["T"] + out_s["T"]) + output["A"] = 1 - output["R"] - output["T"] + output["A_per_layer"] = 0.5 * (A_per_layer_p + A_per_layer_s) - output['out_s'] = out_s - output['out_p'] = out_p - output['A_per_layer_s'] = A_per_layer_s - output['A_per_layer_p'] = A_per_layer_p + output["out_s"] = out_s + output["out_p"] = out_p + output["A_per_layer_s"] = A_per_layer_s + output["A_per_layer_p"] = A_per_layer_p else: - - - out_p = tmm.inc_tmm('p', stack.get_indices(wavelength), stack.get_widths(), coherency_list, angle * degree, wavelength) - out_s = tmm.inc_tmm('s', stack.get_indices(wavelength), stack.get_widths(), coherency_list, angle * degree, wavelength) + out_p = tmm.inc_tmm( + "p", + stack.get_indices(wavelength), + stack.get_widths(), + coherency_list, + angle * degree, + wavelength, + ) + out_s = tmm.inc_tmm( + "s", + stack.get_indices(wavelength), + stack.get_widths(), + coherency_list, + angle * degree, + wavelength, + ) A_per_layer_p = np.array(tmm.inc_absorp_in_each_layer(out_p)) A_per_layer_s = np.array(tmm.inc_absorp_in_each_layer(out_s)) - output['R'] = 0.5 * (out_p['R'] + out_s['R']) - output['T'] = 0.5 * (out_p['T'] + out_s['T']) - output['A'] = 1 - output['R'] - output['T'] - output['all_p'] = out_p['power_entering_list'] - output['all_s'] = out_s['power_entering_list'] - output['A_per_layer'] = 0.5*(A_per_layer_p + A_per_layer_s) - - output['out_s'] = out_s - output['out_p'] = out_p - output['A_per_layer_s'] = A_per_layer_s - output['A_per_layer_p'] = A_per_layer_p + output["R"] = 0.5 * (out_p["R"] + out_s["R"]) + output["T"] = 0.5 * (out_p["T"] + out_s["T"]) + output["A"] = 1 - output["R"] - output["T"] + output["all_p"] = out_p["power_entering_list"] + output["all_s"] = out_s["power_entering_list"] + output["A_per_layer"] = 0.5 * (A_per_layer_p + A_per_layer_s) + output["out_s"] = out_s + output["out_p"] = out_p + output["A_per_layer_s"] = A_per_layer_s + output["A_per_layer_p"] = A_per_layer_p return output -def calculate_ellipsometry(structure, wavelength, angle, no_back_reflection=True, **kwargs): - """ Calculates the ellipsometric parameters psi and delta. It can only deal with coherent light and the whole stack - (including back surface) is considered, so caution must be taken when comparing the simulated results with - experiments where the back surface is rough or layers are thick and coherent light propagation makes no sense. +def calculate_ellipsometry( + structure, wavelength, angle, no_back_reflection=True, **kwargs +): + """Calculates the ellipsometric parameters psi and delta. It can only deal with + coherent light and the whole stack (including back surface) is considered, + so caution must be taken when comparing the simulated results with experiments + where the back surface is rough or layers are thick and coherent light propagation + makes no sense. - The optional argument no_back_reflection can be included to add an extra layer on the back absorbing all light that - reaches that position without any reflexion, to remove the reflexion from the back surface. + The optional argument no_back_reflection can be included to add an extra layer on + the back absorbing all light that reaches that position without any reflection, + to remove the reflection from the back surface. :param structure: A solcore structure with layers and materials. :param wavelength: Wavelengths (in nm) in which calculate the data. An array. - :param angle: A tupple or list with the angles (in degrees) in which to calculate the data. - :param no_back_reflection: If reflexion from the back must be suppressed. Default=True. - :return: A dictionary with psi and delta at the specified wavelengths and angles (2D arrays). + :param angle: A tupple or list with the angles (in degrees) in which to calculate + the data. + :param no_back_reflection: If reflection from the back must be suppressed. + Default=True. + :return: A dictionary with psi and delta at the specified wavelengths and angles + (2D arrays). """ - if 'no_back_reflexion' in kwargs: - warn('The no_back_reflexion argument is deprecated. Use no_back_reflection instead.', FutureWarning) - no_back_reflection = kwargs['no_back_reflexion'] + if "no_back_reflexion" in kwargs: + warn( + "The no_back_reflexion argument is deprecated. " + "Use no_back_reflection instead.", + FutureWarning, + ) + no_back_reflection = kwargs["no_back_reflexion"] num_wl = len(wavelength) num_ang = len(angle) - if 'OptiStack' in str(type(structure)): + if "OptiStack" in str(type(structure)): stack = structure stack.no_back_reflection = no_back_reflection else: - if hasattr(structure, 'substrate'): + if hasattr(structure, "substrate"): substrate = structure.substrate else: substrate = None - stack = OptiStack(structure, no_back_reflection=no_back_reflection, substrate=substrate) + stack = OptiStack( + structure, no_back_reflection=no_back_reflection, substrate=substrate + ) - output = {'psi': np.zeros((num_wl, num_ang)), 'Delta': np.zeros((num_wl, num_ang))} + output = {"psi": np.zeros((num_wl, num_ang)), "Delta": np.zeros((num_wl, num_ang))} for i, ang in enumerate(angle): - out = tmm.ellips(stack.get_indices(wavelength), stack.get_widths(), ang * degree, wavelength) - output['psi'][:, i] = out['psi'] / degree + out = tmm.ellips( + stack.get_indices(wavelength), stack.get_widths(), ang * degree, wavelength + ) + output["psi"][:, i] = out["psi"] / degree # We revere the sign of Delta in order to use Woollam sign convention - output['Delta'][:, i] = -out['Delta'] / degree + output["Delta"][:, i] = -out["Delta"] / degree - output['Delta'][:, i] = np.where(output['Delta'][:, i] > 0, -output['Delta'][:, i], output['Delta'][:, i]) - output['Delta'][:, i] = np.where(output['Delta'][:, i] < 180, 180 - output['Delta'][:, i], - output['Delta'][:, i]) + output["Delta"][:, i] = np.where( + output["Delta"][:, i] > 0, -output["Delta"][:, i], output["Delta"][:, i] + ) + output["Delta"][:, i] = np.where( + output["Delta"][:, i] < 180, + 180 - output["Delta"][:, i], + output["Delta"][:, i], + ) return output - -def calculate_absorption_profile(structure, wavelength, z_limit=None, steps_size=2, dist=None, - no_back_reflection=True, angle=0, pol = 'u', - coherent=True, coherency_list=None, zero_threshold=1e-6, RAT_out=None, **kwargs): - """ It calculates the absorbed energy density within the material. From the documentation: - - 'In principle this has units of [power]/[volume], but we can express it as a multiple of incoming light power - density on the material, which has units [power]/[area], so that absorbed energy density has units of 1/[length].' - - Integrating this absorption profile in the whole stack gives the same result that the absorption obtained with - calculate_rat as long as the spacial mesh (controlled by steps_thinest_layer) is fine enough. If the structure is - very thick and the mesh not thin enough, the calculation might diverge at short wavelengths. +def calculate_absorption_profile( + structure, + wavelength, + z_limit=None, + steps_size=2, + dist=None, + no_back_reflection=True, + angle=0, + pol="u", + coherent=True, + coherency_list=None, + zero_threshold=1e-6, + RAT_out=None, + **kwargs, +): + """It calculates the absorbed energy density within the material. + From the documentation: + + 'In principle this has units of [power]/[volume], but we can express it as a + multiple of incoming light power density on the material, which has units [ + power]/[area], so that absorbed energy density has units of 1/[length].' + + Integrating this absorption profile in the whole stack gives the same result that + the absorption obtained with calculate_rat as long as the spatial mesh (controlled + by steps_thinest_layer) is fine enough. If the structure is very thick and the + mesh not thin enough, the calculation might diverge at short wavelengths. :param structure: A solcore structure with layers and materials. :param wavelength: Wavelengths in which calculate the data (in nm). An array :param z_limit: Maximum value in the z direction - :param steps_size: if the dist is not specified, the step size in nm to use in the depth-dependent calculation + :param steps_size: if the dist is not specified, the step size in nm to use in the + depth-dependent calculation :param dist: the positions (in nm) at which to calculate depth-dependent absorption - :param no_back_reflection: whether to suppress reflections from the back interface (True) or not (False) + :param no_back_reflection: whether to suppress reflections from the back interface + (True) or not (False) :param angle: incidence of angle in degrees :param pol: polarization of incident light: 's', 'p' or 'u' (unpolarized) - :param coherent: True if all the layers are to be treated coherently, False otherwise - :param coherency_list: if coherent is False, a list of 'c' (coherent) or 'i' (incoherent) for each layer - :param zero_threshold: when the fraction of incident light absorbed in a layer is less than this value, the absorption - profile is completely set to zero for both coherent and incoherent calculations. This is applied on a wavelength-by-wavelength - basis and is intended to prevent errors where integrating a weak absorption profile in a layer over many points leads to - calculated EQE > total absorption in that layer. + :param coherent: True if all the layers are to be treated coherently, + False otherwise + :param coherency_list: if coherent is False, a list of 'c' (coherent) or + 'i' (incoherent) for each layer + :param zero_threshold: when the fraction of incident light absorbed in a layer is + less than this value, the absorption profile is completely set to zero for + both coherent and incoherent calculations. This is applied on a + wavelength-by-wavelength basis and is intended to prevent errors where + integrating a weak absorption profile in a layer over many + points leads to calculated EQE > total absorption in that layer. :param RAT_out: output from calculate_rat for the same stack & options - :return: A dictionary containing the positions (in nm) and a 2D array with the absorption in the structure as a - function of the position and the wavelength. + :return: A dictionary containing the positions (in nm) and a 2D array with the + absorption in the structure as a function of the position and the wavelength. """ - if 'no_back_reflexion' in kwargs: - warn('The no_back_reflexion warning is deprecated. Use no_back_reflection instead.', FutureWarning) - no_back_reflection = kwargs['no_back_reflexion'] + if "no_back_reflexion" in kwargs: + warn( + "The no_back_reflexion warning is deprecated. " + "Use no_back_reflection instead.", + FutureWarning, + ) + no_back_reflection = kwargs["no_back_reflexion"] if RAT_out is None: # R, A per layer, T not yet calculated, calculate now - RAT_out = calculate_rat(structure, wavelength, angle, pol=pol, coherent=coherent, coherency_list=coherency_list, - no_back_reflection=no_back_reflection) + RAT_out = calculate_rat( + structure, + wavelength, + angle, + pol=pol, + coherent=coherent, + coherency_list=coherency_list, + no_back_reflection=no_back_reflection, + ) num_wl = len(wavelength) - if 'OptiStack' in str(type(structure)): + if "OptiStack" in str(type(structure)): stack = structure stack.no_back_reflection = no_back_reflection else: @@ -551,86 +739,104 @@ def calculate_absorption_profile(structure, wavelength, z_limit=None, steps_size if not coherent: if coherency_list is not None: - if stack.no_back_reflection: - coherency_list = ['i'] + coherency_list + ['i', 'i'] + coherency_list = ["i"] + coherency_list + ["i", "i"] else: - coherency_list = ['i'] + coherency_list + ['i'] + coherency_list = ["i"] + coherency_list + ["i"] else: - raise Exception('Error: For incoherent or partly incoherent calculations you must supply the ' - 'coherency_list parameter with as many elements as the number of layers in the ' - 'structure') + raise Exception( + "Error: For incoherent or partly incoherent calculations you must " + "supply the coherency_list parameter with as many elements as the " + "number of layers in the structure" + ) - output = {'position': dist, 'absorption': np.zeros((num_wl, len(dist)))} + output = {"position": dist, "absorption": np.zeros((num_wl, len(dist)))} - if pol in 'sp': + layer, d_in_layer = tmm.find_in_structure_with_inf(stack.get_widths(), dist) + if pol in "sp": if coherent: - A_per_layer = RAT_out['A_per_layer'] - no_abs_in_layer = np.where(A_per_layer[:-1,:] < zero_threshold) - RAT_out['out']['vw_list'][no_abs_in_layer[0], no_abs_in_layer[1], :] = 0 + A_per_layer = RAT_out["A_per_layer"] + no_abs_in_layer = np.where(A_per_layer[:-1, :] < zero_threshold) + RAT_out["out"]["vw_list"][no_abs_in_layer[0], no_abs_in_layer[1], :] = 0 - layer, d_in_layer = tmm.find_in_structure_with_inf(stack.get_widths(), dist) - data = tmm.position_resolved(layer, d_in_layer, RAT_out['out']) - output['absorption'] = data['absor'] + data = tmm.position_resolved(layer, d_in_layer, RAT_out["out"]) + output["absorption"] = data["absor"] else: - layer, d_in_layer = tmm.find_in_structure_with_inf(stack.get_widths(), dist) - data = tmm.inc_position_resolved(layer, d_in_layer, RAT_out['out'], coherency_list, - 4*np.pi*np.imag(stack.get_indices(wavelength))/wavelength, zero_threshold) - output['absorption'] = data + data = tmm.inc_position_resolved( + layer, + d_in_layer, + RAT_out["out"], + coherency_list, + 4 * np.pi * np.imag(stack.get_indices(wavelength)) / wavelength, + zero_threshold, + ) + output["absorption"] = data else: - - - if coherent: - - A_per_layer_s = RAT_out['A_per_layer_s'] - A_per_layer_p = RAT_out['A_per_layer_p'] - no_abs_in_layer_s = np.where(A_per_layer_s[:-1,:] < zero_threshold) + A_per_layer_s = RAT_out["A_per_layer_s"] + A_per_layer_p = RAT_out["A_per_layer_p"] + no_abs_in_layer_s = np.where(A_per_layer_s[:-1, :] < zero_threshold) no_abs_in_layer_p = np.where(A_per_layer_p[:-1, :] < zero_threshold) - RAT_out['out_s']['vw_list'][no_abs_in_layer_s[0], no_abs_in_layer_s[1], :] = 0 - RAT_out['out_p']['vw_list'][no_abs_in_layer_p[0], no_abs_in_layer_p[1], :] = 0 + RAT_out["out_s"]["vw_list"][ + no_abs_in_layer_s[0], no_abs_in_layer_s[1], : + ] = 0 + RAT_out["out_p"]["vw_list"][ + no_abs_in_layer_p[0], no_abs_in_layer_p[1], : + ] = 0 - layer, d_in_layer = tmm.find_in_structure_with_inf(stack.get_widths(), dist) + data_s = tmm.position_resolved(layer, d_in_layer, RAT_out["out_s"]) + data_p = tmm.position_resolved(layer, d_in_layer, RAT_out["out_p"]) - data_s = tmm.position_resolved(layer, d_in_layer, RAT_out['out_s']) - data_p = tmm.position_resolved(layer, d_in_layer, RAT_out['out_p']) - - output['absorption'] = 0.5*(data_s['absor'] + data_p['absor']) + output["absorption"] = 0.5 * (data_s["absor"] + data_p["absor"]) else: - layer, d_in_layer = tmm.find_in_structure_with_inf(stack.get_widths(), dist) - data_s = tmm.inc_position_resolved(layer, d_in_layer, RAT_out['out_s'], coherency_list, - 4*np.pi*np.imag(stack.get_indices(wavelength))/wavelength, zero_threshold) - data_p = tmm.inc_position_resolved(layer, d_in_layer, RAT_out['out_p'], coherency_list, - 4*np.pi*np.imag(stack.get_indices(wavelength))/wavelength, zero_threshold) - - output['absorption'] = 0.5*(data_s + data_p) + data_s = tmm.inc_position_resolved( + layer, + d_in_layer, + RAT_out["out_s"], + coherency_list, + 4 * np.pi * np.imag(stack.get_indices(wavelength)) / wavelength, + zero_threshold, + ) + data_p = tmm.inc_position_resolved( + layer, + d_in_layer, + RAT_out["out_p"], + coherency_list, + 4 * np.pi * np.imag(stack.get_indices(wavelength)) / wavelength, + zero_threshold, + ) + + output["absorption"] = 0.5 * (data_s + data_p) return output -if __name__ == '__main__': +if __name__ == "__main__": import matplotlib.pyplot as plt from solcore import material, si from solcore.structure import Layer, Structure - GaAs = material('GaAs')(T=300) - InGaAs = material('InGaAs')(T=300, In=0.1) - - my_structure = Structure([ - Layer(si(3000, 'nm'), material=InGaAs), - Layer(si(30, 'um'), material=GaAs), + GaAs = material("GaAs")(T=300) + InGaAs = material("InGaAs")(T=300, In=0.1) - ]) + my_structure = Structure( + [ + Layer(si(3000, "nm"), material=InGaAs), + Layer(si(30, "um"), material=GaAs), + ] + ) wavelength = np.linspace(450, 1100, 300) - out = calculate_rat(my_structure, wavelength, coherent=True, no_back_reflection=False) + out = calculate_rat( + my_structure, wavelength, coherent=True, no_back_reflection=False + ) # # # plt.plot(wavelength, out['R'], 'b', label='Reflexion') # plt.plot(wavelength, out['A'], 'r', label='Absorption') @@ -641,11 +847,11 @@ def calculate_absorption_profile(structure, wavelength, z_limit=None, steps_size angles = [60, 65, 70] out = calculate_ellipsometry(my_structure, wavelength, angle=angles) # - plt.plot(wavelength, out['psi'][:, 0], 'b', label='psi') - plt.plot(wavelength, out['Delta'][:, 0], 'r', label='Delta') + plt.plot(wavelength, out["psi"][:, 0], "b", label="psi") + plt.plot(wavelength, out["Delta"][:, 0], "r", label="Delta") for i in range(1, len(angles)): - plt.plot(wavelength, out['psi'][:, i], 'b') - plt.plot(wavelength, out['Delta'][:, i], 'r') + plt.plot(wavelength, out["psi"][:, i], "b") + plt.plot(wavelength, out["Delta"][:, i], "r") # # plt.legend() diff --git a/solcore/analytic_solar_cells/depletion_approximation.py b/solcore/analytic_solar_cells/depletion_approximation.py index 69a6777e..e341f2a5 100755 --- a/solcore/analytic_solar_cells/depletion_approximation.py +++ b/solcore/analytic_solar_cells/depletion_approximation.py @@ -8,24 +8,29 @@ from solcore.state import State from solcore.light_source import LightSource +import warnings + da_options = State() -da_options.da_mode = 'bvp' +da_options.da_mode = "green" + def identify_layers(junction): # First we have to figure out if we are talking about a PN, NP, PIN or NIP junction # We search for the emitter and check if it is n-type or p-type idx = 0 - pn_or_np = 'pn' + pn_or_np = "pn" homojunction = True for layer in junction: - if layer.role.lower() != 'emitter': + if layer.role.lower() != "emitter": idx += 1 else: Na = 0 Nd = 0 - if hasattr(layer.material, 'Na'): Na = layer.material.Na - if hasattr(layer.material, 'Nd'): Nd = layer.material.Nd + if hasattr(layer.material, "Na"): + Na = layer.material.Na + if hasattr(layer.material, "Nd"): + Nd = layer.material.Nd if Na < Nd: pn_or_np = "np" nRegion = junction[idx] @@ -37,10 +42,10 @@ def identify_layers(junction): break # Now we check for an intrinsic region and, if there is, for the base. - if junction[idx + 1].role.lower() == 'intrinsic': + if junction[idx + 1].role.lower() == "intrinsic": iRegion = junction[idx + 1] - if junction[idx + 2].role.lower() == 'base': + if junction[idx + 2].role.lower() == "base": if pn_or_np == "pn": nRegion = junction[idx + 2] @@ -48,16 +53,23 @@ def identify_layers(junction): pRegion = junction[idx + 2] id_bottom = idx + 2 - homojunction = homojunction and nRegion.material.material_string == pRegion.material.material_string - homojunction = homojunction and nRegion.material.material_string == iRegion.material.material_string + homojunction = ( + homojunction + and nRegion.material.material_string == pRegion.material.material_string + ) + homojunction = ( + homojunction + and nRegion.material.material_string == iRegion.material.material_string + ) else: raise RuntimeError( - 'ERROR processing junctions: A layer following the "intrinsic" layer must be defined as ' - '"base".') + "ERROR processing junctions: A layer following the " + '"intrinsic" layer must be defined as "base".' + ) # If there is no intrinsic region, we check directly the base - elif junction[idx + 1].role.lower() == 'base': + elif junction[idx + 1].role.lower() == "base": if pn_or_np == "pn": nRegion = junction[idx + 1] @@ -67,22 +79,24 @@ def identify_layers(junction): iRegion = None id_bottom = idx + 1 - homojunction = homojunction and nRegion.material.material_string == pRegion.material.material_string + homojunction = ( + homojunction + and nRegion.material.material_string == pRegion.material.material_string + ) else: raise RuntimeError( - 'ERROR processing junctions: A layer following the "emitter" must be defined as "intrinsic"' - 'or "base".') + 'ERROR processing junctions: A layer following the "emitter" ' + 'must be defined as "intrinsic" or "base".' + ) # We assert that we are really working with an homojunction - assert homojunction, 'ERROR: The DA solver only works with homojunctions, for now.' - + assert homojunction, "ERROR: The DA solver only works with homojunctions, for now." return id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np def identify_parameters(junction, T, pRegion, nRegion, iRegion): - kbT = kb * T xp = pRegion.width @@ -98,7 +112,8 @@ def identify_parameters(junction, T, pRegion, nRegion, iRegion): else: es = nRegion.material.permittivity # equal for n and p. I hope. - # For the diffusion length, subscript n and p refer to the carriers, electrons and holes + # For the diffusion length, subscript n and p refer to the carriers, + # electrons and holes if hasattr(junction, "ln"): ln = junction.ln else: @@ -109,7 +124,8 @@ def identify_parameters(junction, T, pRegion, nRegion, iRegion): else: lp = nRegion.material.hole_diffusion_length - # For the diffusion coefficient, n and p refer to the regions, n side and p side. Yeah, it's confusing... + # For the diffusion coefficient, n and p refer to the regions, + # n side and p side. Yeah, it's confusing... if hasattr(junction, "mup"): dp = junction.mup * kbT / q else: @@ -129,32 +145,42 @@ def identify_parameters(junction, T, pRegion, nRegion, iRegion): def iv_depletion(junction, options): - """ Calculates the IV curve of a junction object using the depletion approximation as described in J. Nelson, “The Physics of Solar Cells”, Imperial College Press (2003). The junction is then updated with an "iv" function that calculates the IV curve at any voltage. + """Calculates the IV curve of a junction object using the depletion approximation as + described in J. Nelson, “The Physics of Solar Cells”, Imperial College Press (2003). + The junction is then updated with an "iv" function that calculates the IV curve at + any voltage. :param junction: A junction object. :param options: Solver options. :return: None. """ - science_reference('Depletion approximation', - 'J. Nelson, “The Physics of Solar Cells”, Imperial College Press (2003).') + science_reference( + "Depletion approximation", + "J. Nelson, “The Physics of Solar Cells”, Imperial College Press (2003).", + ) junction.voltage = options.internal_voltages T = options.T kbT = kb * T id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np = identify_layers(junction) - xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd, Na, ni, es = identify_parameters(junction, T, pRegion, nRegion, iRegion) + xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd, Na, ni, es = identify_parameters( + junction, T, pRegion, nRegion, iRegion + ) niSquared = ni**2 - Vbi = (kbT / q) * np.log(Nd * Na / niSquared) if not hasattr(junction, "Vbi") else junction.Vbi # Jenny p146 - - #Na, Nd, ni, niSquared, xi, ln, lp, xn, xp, sn, sp, dn, dp, es, id_top, id_bottom, Vbi, pn_or_np = process_junction(junction, options) + Vbi = ( + (kbT / q) * np.log(Nd * Na / niSquared) + if not hasattr(junction, "Vbi") + else junction.Vbi + ) # Jenny p146 - R_shunt = min(junction.R_shunt, 1e14) if hasattr(junction, 'R_shunt') else 1e14 + R_shunt = min(junction.R_shunt, 1e14) if hasattr(junction, "R_shunt") else 1e14 - # And now we account for the possible applied voltage, which can be, at most, equal to Vbi + # And now we account for the possible applied voltage, + # which can be, at most, equal to Vbi V = np.where(junction.voltage < Vbi - 0.001, junction.voltage, Vbi - 0.001) wn, wp = get_depletion_widths(junction, es, Vbi, V, Na, Nd, xi) @@ -178,24 +204,32 @@ def iv_depletion(junction, options): min_bot, min_top = niSquared / Na, niSquared / Nd JtopDark = get_j_dark(x_top, w_top, l_top, s_top, d_top, V, min_top, T) - JbotDark = get_j_dark(x_bottom, w_bottom, l_bottom, s_bottom, d_bottom, V, min_bot, T) + JbotDark = get_j_dark( + x_bottom, w_bottom, l_bottom, s_bottom, d_bottom, V, min_bot, T + ) - # hereby we define the subscripts to refer to the layer in which the current is generated: + # hereby we define the subscripts to refer to the layer in which + # the current is generated: if pn_or_np == "pn": JnDark, JpDark = JbotDark, JtopDark else: JpDark, JnDark = JbotDark, JtopDark - # These might not be the right lifetimes. Actually, they are not as they include all recombination processes, not - # just SRH recombination, which is what the equation in Jenny, p159 refers to. Let´ leave them, for now. - lifetime_n = ln ** 2 / dn - lifetime_p = lp ** 2 / dp # Jenny p163 - - # Here we use the full version of the SRH recombination term as calculated by Sah et al. Works for positive bias - # and moderately negative ones. - - science_reference('SRH current term.', - 'C. T. Sah, R. N. Noyce, and W. Shockley, “Carrier Generation and Recombination in P-N Junctions and P-N Junction Characteristics,” presented at the Proceedings of the IRE, 1957, vol. 45, no. 9, pp. 1228–1243.') + # These might not be the right lifetimes. Actually, they are not as + # they include all recombination processes, not just SRH recombination, + # which is what the equation in Jenny, p159 refers to. Let´s leave them, for now. + lifetime_n = ln**2 / dn + lifetime_p = lp**2 / dp # Jenny p163 + + # Here we use the full version of the SRH recombination term as calculated by + # Sah et al. Works for positive bias and moderately negative ones. + + science_reference( + "SRH current term.", + "C. T. Sah, R. N. Noyce, and W. Shockley, “Carrier Generation and " + "Recombination in P-N Junctions and P-N Junction Characteristics,” presented " + "at the Proceedings of the IRE, 1957, vol. 45, no. 9, pp. 1228–1243.", + ) Jrec = get_Jsrh(ni, V, Vbi, lifetime_p, lifetime_n, w, kbT) J_sc_top = 0 @@ -203,7 +237,6 @@ def iv_depletion(junction, options): J_sc_scr = 0 if options.light_iv: - widths = [] for layer in junction: widths.append(layer.width) @@ -212,30 +245,54 @@ def iv_depletion(junction, options): g = junction.absorbed wl = options.wavelength - wl_sp, ph = options.light_source.spectrum(output_units='photon_flux_per_m', x=wl) + wl_sp, ph = options.light_source.spectrum( + output_units="photon_flux_per_m", x=wl + ) id_v0 = np.argmin(abs(V)) # The contribution from the Emitter (top side). xa = cum_widths[id_top] xb = cum_widths[id_top + 1] - w_top[id_v0] - if options.da_mode == 'bvp': - deriv = get_J_sc_diffusion(xa, xb, g, d_top, l_top, min_top, s_top, wl, ph, side='top') + if options.da_mode == "bvp": + deriv = get_J_sc_diffusion( + xa, xb, g, d_top, l_top, min_top, s_top, wl, ph, side="top" + ) else: - xbb = xb - (xb - xa)/1001. - deriv = get_J_sc_diffusion_green(xa, xbb, g, d_top, l_top, min_top, s_top, ph, side='top') - deriv = np.trapz(deriv, wl) + xbb = xb - (xb - xa) / 1001.0 + + zz = np.linspace(xa, xb, 1001, endpoint=False) + gg = ph * g(zz) + g_vs_z = np.trapz(gg, wl, axis=1) + g_vs_z[np.isnan(g_vs_z)] = 0 + g_vs_z = interp1d(zz, g_vs_z, axis=0) + + deriv = get_J_sc_diffusion_green( + xa, xbb, g_vs_z, d_top, l_top, s_top, 1, side="top" + ) + J_sc_top = q * d_top * abs(deriv) # The contribution from the Base (bottom side). xa = cum_widths[id_bottom] + w_bottom[id_v0] xb = cum_widths[id_bottom + 1] - if options.da_mode == 'bvp': - deriv = get_J_sc_diffusion(xa, xb, g, d_bottom, l_bottom, min_bot, s_bottom, wl, ph, side='bottom') + if options.da_mode == "bvp": + deriv = get_J_sc_diffusion( + xa, xb, g, d_bottom, l_bottom, min_bot, s_bottom, wl, ph, side="bottom" + ) else: - xbb = xb - (xb - xa)/1001. - deriv = get_J_sc_diffusion_green(xa, xbb, g, d_bottom, l_bottom, min_bot, s_bottom, ph, side='bottom') - deriv = np.trapz(deriv, wl) + xbb = xb - (xb - xa) / 1001.0 + zz = np.linspace(xa, xb, 1001, endpoint=False) + gg = ph * g(zz) + g_vs_z = np.trapz(gg, wl, axis=1) + g_vs_z[np.isnan(g_vs_z)] = 0 + g_vs_z = interp1d(zz, g_vs_z, axis=0) + + deriv = get_J_sc_diffusion_green( + xa, xbb, g_vs_z, d_bottom, l_bottom, s_bottom, 1, side="bottom" + ) + # photogeneration is included in g_vs_z, so set ph=1 + J_sc_bot = q * d_bottom * abs(deriv) # The contribution from the SCR (includes the intrinsic region, if present). @@ -244,18 +301,34 @@ def iv_depletion(junction, options): J_sc_scr = q * get_J_sc_SCR(xa, xb, g, wl, ph) # And, finally, we output the currents - junction.current = Jrec + JnDark + JpDark + V / R_shunt - J_sc_top - J_sc_bot - J_sc_scr - junction.iv = interp1d(junction.voltage, junction.current, kind='linear', bounds_error=False, assume_sorted=True, - fill_value=(junction.current[0], junction.current[-1])) - junction.region_currents = State({"Jn_dif": JnDark, "Jp_dif": JpDark, "Jscr_srh": Jrec, - "J_sc_top": J_sc_top, "J_sc_bot": J_sc_bot, "J_sc_scr": J_sc_scr}) - - -def get_j_dark(x, w, l, s, d, V, minority, T): + junction.current = ( + Jrec + JnDark + JpDark + V / R_shunt - J_sc_top - J_sc_bot - J_sc_scr + ) + junction.iv = interp1d( + junction.voltage, + junction.current, + kind="linear", + bounds_error=False, + assume_sorted=True, + fill_value=(junction.current[0], junction.current[-1]), + ) + junction.region_currents = State( + { + "Jn_dif": JnDark, + "Jp_dif": JpDark, + "Jscr_srh": Jrec, + "J_sc_top": J_sc_top, + "J_sc_bot": J_sc_bot, + "J_sc_scr": J_sc_scr, + } + ) + + +def get_j_dark(x, w, L, s, d, V, minority, T): """ :param x: width of top junction :param w: depletion width in top junction - :param l: diffusion length + :param L: diffusion length :param s: surface recombination velocity :param d: diffusion coefficient :param V: voltage @@ -266,24 +339,31 @@ def get_j_dark(x, w, l, s, d, V, minority, T): """ # We calculate some fractions - harg = (x - w) / l + harg = (x - w) / L sinh_harg = np.sinh(harg) cosh_harg = np.cosh(harg) - lsod = (l * s) / d + lsod = (L * s) / d # And then we are ready to calculate the different currents # Missing the voltage dependent part of these equations. # They should be 6.34 and 6.39, not 6.62 and 6.63 - J_dark = (q * d * minority / l) * (np.exp(q * V / kb / T) - 1) * \ - ((lsod * cosh_harg + sinh_harg) / (lsod * sinh_harg + cosh_harg)) + J_dark = ( + (q * d * minority / L) + * (np.exp(q * V / kb / T) - 1) + * ((lsod * cosh_harg + sinh_harg) / (lsod * sinh_harg + cosh_harg)) + ) return J_dark def get_Jsrh(ni, V, Vbi, tp, tn, w, kbT, dEt=0): - science_reference('SRH current term.', - 'C. T. Sah, R. N. Noyce, and W. Shockley, “Carrier Generation and Recombination in P-N Junctions and P-N Junction Characteristics,” presented at the Proceedings of the IRE, 1957, vol. 45, no. 9, pp. 1228–1243.') + science_reference( + "SRH current term.", + "C. T. Sah, R. N. Noyce, and W. Shockley, “Carrier Generation and Recombination" + " in P-N Junctions and P-N Junction Characteristics,” presented at the " + "Proceedings of the IRE, 1957, vol. 45, no. 9, pp. 1228–1243.", + ) out = np.zeros(V.shape) @@ -294,7 +374,8 @@ def get_Jsrh(ni, V, Vbi, tp, tn, w, kbT, dEt=0): def forward(ni, V, Vbi, tp, tn, w, kbT, dEt=0): - """ Equation 27 of Sah's paper. Strictly speaking, it is not valid for intermediate negative bias. """ + """Equation 27 of Sah's paper. Strictly speaking, it is not valid for + intermediate negative bias.""" J0 = 2 * q * ni * w / np.sqrt(tn * tp) f = factor(V, Vbi, tp, tn, kbT, dEt) @@ -304,8 +385,8 @@ def forward(ni, V, Vbi, tp, tn, w, kbT, dEt=0): def factor(V, Vbi, tp, tn, kbT, dEt=0): - """ The integral of Eq. 27 in Sah's paper. While it is coninuum (in principle) it has to be done in two parts. - (or three) """ + """The integral of Eq. 27 in Sah's paper. While it is continuum (in principle) it + has to be done in two parts (or three)""" trap = q * dEt / kbT + np.log(tp / tn) / 2 b = np.exp(-q * V / kbT / 2) * np.cosh(trap) @@ -324,7 +405,7 @@ def factor(V, Vbi, tp, tn, kbT, dEt=0): # For b values >= 1 and <=1e7 m = (b >= 1) * (b <= 1e7) xx = b[m] - yy = np.sqrt(xx ** 2 - 1) + yy = np.sqrt(xx**2 - 1) top = np.log(abs(z2[m] - yy + xx) / abs(z2[m] + yy + xx)) bot = np.log(abs(z1[m] - yy + xx) / abs(z1[m] + yy + xx)) @@ -342,7 +423,7 @@ def factor(V, Vbi, tp, tn, kbT, dEt=0): return out -def get_J_sc_diffusion(xa, xb, g, D, L, y0, S, wl, ph, side='top'): +def get_J_sc_diffusion(xa, xb, g, D, L, y0, S, wl, ph, side="top"): """ :param xa: :param xb: @@ -358,6 +439,8 @@ def get_J_sc_diffusion(xa, xb, g, D, L, y0, S, wl, ph, side='top'): :return: out """ + # see comments in get_J_sc_diffusion_vs_WL for details on how differential + # equations are solved zz = np.linspace(xa, xb, 1001, endpoint=False) gg = g(zz) * ph @@ -365,39 +448,54 @@ def get_J_sc_diffusion(xa, xb, g, D, L, y0, S, wl, ph, side='top'): g_vs_z[np.isnan(g_vs_z)] = 0 - A = lambda x: np.interp(x, zz, g_vs_z) / D + y0 / L ** 2 + def A(x): + return np.interp(x, zz, g_vs_z) / D + y0 / L**2 def fun(x, y): out1 = y[1] - out2 = y[0] / L ** 2 - A(x) + out2 = y[0] / L**2 - A(x) return np.vstack((out1, out2)) - if side == 'top': + if side == "top": + def bc(ya, yb): left = ya[1] - S / D * (ya[0] - y0) - right = yb[0] + right = yb[0] - y0 return np.array([left, right]) + else: + def bc(ya, yb): - left = ya[0] - right = yb[1] - S / D * (yb[0] - y0) + left = ya[0] - y0 + right = yb[1] + S / D * (yb[0] - y0) return np.array([left, right]) guess = y0 * np.ones((2, zz.size)) guess[1] = np.zeros_like(guess[0]) - solution = solve_bvp(fun, bc, zz, guess) + solution = solve_bvp(fun, bc, zz, guess, max_nodes=2 * zz.shape[0]) - if side == 'top': + if side == "top": out = solution.y[1][-1] else: out = solution.y[1][0] + if solution.status != 0: + warnings.warn( + "Depletion approximation (DA) I-V calculation: " + "solve_bvp did not converge as expected", + RuntimeWarning, + ) + return out def _conv_exp_top(x, xa, xb, g, L, phoD): - """Convolution of the carrier generation rate with the approximate Green's function kernel at point x. To be used with the numerical integration routine to compute the minority carrier derivative on the top edge. This kernel approximates the original one when the diffusion length is 2 orders of magnitude higher than the junction width by assuming that sinh(x) = cosh(x) = .5 * exp(x). + """Convolution of the carrier generation rate with the approximate Green's function + kernel at point x. To be used with the numerical integration routine to compute the + minority carrier derivative on the top edge. This kernel approximates the original + one when the diffusion length is 2 orders of magnitude higher than the junction + width by assuming that sinh(x) = cosh(x) = .5 * exp(x). :param x: Coordinate in the junction (variable to be integrated). :param xa: Coordinate at the start the junction. @@ -407,14 +505,22 @@ def _conv_exp_top(x, xa, xb, g, L, phoD): :param phoD: Light spectrum divided by the diffusion constant D. """ xc = (xa - x) / L - xv = np.array([xa + xb - x, ]) + xv = np.array( + [ + xa + xb - x, + ] + ) Pkern = -np.exp(xc) Gx = g(xv) * phoD - return Pkern*Gx + return Pkern * Gx def _conv_exp_bottom(x, xa, g, L, phoD): - """Convolution of the carrier generation rate with the approximate Green's function kernel at point x. To be used with the numerical integration routine to compute the minority carrier derivative on the bottom edge. This kernel approximates the original one when the diffusion length is 2 orders of magnitude higher than the junction width by assuming that sinh(x) = cosh(x) = .5 * exp(x). + """Convolution of the carrier generation rate with the approximate Green's function + kernel at point x. To be used with the numerical integration routine to compute the + minority carrier derivative on the bottom edge. This kernel approximates the + original one when the diffusion length is 2 orders of magnitude higher than + the junction width by assuming that sinh(x) = cosh(x) = .5 * exp(x). :param x: Coordinate in the junction (variable to be integrated). :param xa: Coordinate at the start the junction. @@ -423,14 +529,20 @@ def _conv_exp_bottom(x, xa, g, L, phoD): :param phoD: Light spectrum divided by the diffusion constant D. """ xc = (xa - x) / L - xv = np.array([x, ]) + xv = np.array( + [ + x, + ] + ) Pkern = np.exp(xc) Gx = g(xv) * phoD - return Pkern*Gx + return Pkern * Gx def _conv_green_top(x, xa, xb, g, L, phoD, crvel): - """Convolution of the carrier generation rate with the Green's function kernel at point x. To be used with the numerical integration routine to compute the minority carrier derivative on the top edge. + """Convolution of the carrier generation rate with the Green's function kernel at + point x. To be used with the numerical integration routine to compute the minority + carrier derivative on the top edge. :param x: Coordinate in the junction (variable to be integrated). :param xa: Coordinate at the start the junction. @@ -438,34 +550,49 @@ def _conv_green_top(x, xa, xb, g, L, phoD, crvel): :param g: Carrier generation rate at point x (expected as function). :param L: Diffusion length. :param phoD: Light spectrum divided by the diffusion constant D. - :param crvel: Coefficient computed as S / D * L, with S the surface recombination velocity. + :param crvel: Coefficient computed as S / D * L, with S the surface recombination + velocity. """ xc = (xb - x) / L - xv = np.array([xa + xb - x, ]) + xv = np.array( + [ + xa + xb - x, + ] + ) Pkern = np.cosh(xc) + crvel * np.sinh(xc) Gx = g(xv) * phoD - return Pkern*Gx + + return Pkern * Gx def _conv_green_bottom(x, xb, g, L, phoD, crvel): - """Convolution of the carrier generation rate with the Green's function kernel at point x. To be used with the numerical integration routine to compute the minority carrier derivative on the bottom edge. + """Convolution of the carrier generation rate with the Green's function kernel at + point x. To be used with the numerical integration routine to compute the minority + carrier derivative on the bottom edge. :param x: Coordinate in the junction (variable to be integrated). :param xb: Coordinate at the end the junction. :param g: Carrier generation rate at point x (expected as function). :param L: Diffusion length. :param phoD: Light spectrum divided by the diffusion constant D. - :param crvel: Coefficient computed as S / D * L, with S the surface recombination velocity. + :param crvel: Coefficient computed as S / D * L, with S the surface recombination + velocity. """ xc = (xb - x) / L - xv = np.array([x, ]) + xv = np.array( + [ + x, + ] + ) Pkern = np.cosh(xc) - crvel * np.sinh(xc) Gx = g(xv) * phoD - return Pkern*Gx + return Pkern * Gx -def get_J_sc_diffusion_green(xa, xb, g, D, L, y0, S, ph, side='top'): - """Computes the derivative of the minority carrier concentration at the edge of the junction by approximating the convolution integral resulting from applying the Green's function method to the drift-diffusion equation. +def get_J_sc_diffusion_green(xa, xb, g, D, L, S, ph, side="top"): + """Computes the derivative of the minority carrier concentration at the edge of the + junction by approximating the convolution integral resulting from applying the + Green's function method to the drift-diffusion equation. :param xa: Coordinate at the start the junction. :param xb: Coordinate at the end the junction. @@ -477,41 +604,46 @@ def get_J_sc_diffusion_green(xa, xb, g, D, L, y0, S, ph, side='top'): :param ph: Light spectrum. :param side: String to indicate the edge of interest. Either 'top' or 'bottom'. - :return: The derivative of the minority carrier concentration at the edge of the junction. + :return: The derivative of the minority carrier concentration at the edge of the + junction. """ - science_reference('DA Green\'s function method.', - 'T. Vasileiou, J. M. Llorens, J. Buencuerpo, J. M. Ripalda, D. Izzo and L. Summerer, “Light absorption enhancement and radiation hardening for triple junction solar cell through bioinspired nanostructures,” Bioinspir. Biomim., vol. 16, no. 5, pp. 056010, 2021.') + science_reference( + "DA Green's function method.", + "T. Vasileiou, J. M. Llorens, J. Buencuerpo, J. M. Ripalda, D. Izzo and " + "L. Summerer, “Light absorption enhancement and radiation hardening for triple " + "junction solar cell through bioinspired nanostructures,” " + "Bioinspir. Biomim., vol. 16, no. 5, pp. 056010, 2021.", + ) xbL = (xb - xa) / L crvel = S / D * L ph_over_D = ph / D # if L too low in comparison to junction width, avoid nan's - if xbL > 1.e2: - if side == 'top': - cadd = -y0 / L - fun = partial(_conv_exp_top, xa=xa, xb=xb, g=g, - L=L, phoD=ph_over_D) + if xbL > 1.0e2: + if side == "top": + fun = partial(_conv_exp_top, xa=xa, xb=xb, g=g, L=L, phoD=ph_over_D) else: - cadd = y0 / L - fun = partial(_conv_exp_bottom, xa=xa, g=g, - L=L, phoD=ph_over_D) - cp = 1. + fun = partial(_conv_exp_bottom, xa=xa, g=g, L=L, phoD=ph_over_D) + cp = 1.0 + else: - if side == 'top': + if side == "top": cp = -np.cosh(xbL) - crvel * np.sinh(xbL) - cadd = (np.sinh(xbL) + crvel * np.cosh(xbL)) * y0 / L - fun = partial(_conv_green_top, xa=xa, xb=xb, g=g, - L=L, phoD=ph_over_D, crvel=crvel) + + fun = partial( + _conv_green_top, xa=xa, xb=xb, g=g, L=L, phoD=ph_over_D, crvel=crvel + ) else: - cp = np.cosh(xbL) - crvel * np.sinh(xbL) - cadd = (np.sinh(xbL) - crvel * np.cosh(xbL)) * y0 / L - fun = partial(_conv_green_bottom, xb=xb, g=g, - L=L, phoD=ph_over_D, crvel=crvel) + cp = np.cosh(xbL) + crvel * np.sinh(xbL) - out, err = quad_vec(fun, xa, xb, epsrel=1.e-5) - return (out.squeeze() + cadd) / cp + fun = partial( + _conv_green_bottom, xb=xb, g=g, L=L, phoD=ph_over_D, crvel=-crvel + ) + + out, err = quad_vec(fun, xa, xb, epsrel=1.0e-5) + return out.squeeze() / cp def get_J_sc_SCR(xa, xb, g, wl, ph): @@ -523,27 +655,37 @@ def get_J_sc_SCR(xa, xb, g, wl, ph): def qe_depletion(junction, options): - """ Calculates the QE curve of a junction object using the depletion approximation as described in J. Nelson, “The Physics of Solar Cells”, Imperial College Press (2003). The junction is then updated with an "iqe" and several "eqe" functions that calculates the QE curve at any wavelength. + """Calculates the QE curve of a junction object using the depletion approximation as + described in J. Nelson, “The Physics of Solar Cells”, Imperial College Press (2003). + The junction is then updated with an "iqe" and several "eqe" functions that + calculates the QE curve at any wavelength. :param junction: A junction object. :param options: Solver options. :return: None. """ - science_reference('Depletion approximation', - 'J. Nelson, “The Physics of Solar Cells”, Imperial College Press (2003).') - + science_reference( + "Depletion approximation", + "J. Nelson, “The Physics of Solar Cells”, Imperial College Press (2003).", + ) # First we have to figure out if we are talking about a PN, NP, PIN or NIP junction T = options.T kbT = kb * T id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np = identify_layers(junction) - xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd, Na, ni, es = identify_parameters(junction, T, pRegion, nRegion, iRegion) + xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd, Na, ni, es = identify_parameters( + junction, T, pRegion, nRegion, iRegion + ) - niSquared = ni ** 2 + niSquared = ni**2 - Vbi = (kbT / q) * np.log(Nd * Na / niSquared) if not hasattr(junction, "Vbi") else junction.Vbi # Jenny p146 + Vbi = ( + (kbT / q) * np.log(Nd * Na / niSquared) + if not hasattr(junction, "Vbi") + else junction.Vbi + ) # Jenny p146 wn, wp = get_depletion_widths(junction, es, Vbi, 0, Na, Nd, xi) @@ -569,29 +711,41 @@ def qe_depletion(junction, options): g = junction.absorbed wl = options.wavelength - wl_sp, ph = LightSource(source_type='black body', x=wl, T=6000).spectrum(output_units='photon_flux_per_m', x=wl) - + wl_sp, ph = LightSource(source_type="black body", x=wl, T=6000).spectrum( + output_units="photon_flux_per_m", x=wl + ) + # wl_sp, ph = options.light_source.spectrum(output_units='photon_flux_per_m', x=wl) + ph = 1e10 * np.ones_like(ph) # The contribution from the Emitter (top side). xa = cum_widths[id_top] xb = cum_widths[id_top + 1] - w_top - if options.da_mode == 'bvp': - deriv = get_J_sc_diffusion_vs_WL(xa, xb, g, d_top, l_top, min_top, s_top, wl, ph, side='top') + if options.da_mode == "bvp": + deriv = get_J_sc_diffusion_vs_WL( + xa, xb, g, d_top, l_top, min_top, s_top, wl, ph, side="top" + ) else: - xbb = xb - (xb - xa)/1001. - deriv = get_J_sc_diffusion_green(xa, xbb, g, d_top, l_top, min_top, s_top, ph, side='top') - j_sc_top = d_top * abs(deriv) + xbb = xb - (xb - xa) / 1001.0 + deriv = get_J_sc_diffusion_green( + xa, xbb, g, d_top, l_top, s_top, ph, side="top" + ) + j_sc_top = d_top * abs(deriv) # The contribution from the Base (bottom side). xa = cum_widths[id_bottom] + w_bottom xb = cum_widths[id_bottom + 1] - if options.da_mode == 'bvp': - deriv = get_J_sc_diffusion_vs_WL(xa, xb, g, d_bottom, l_bottom, min_bot, s_bottom, wl, ph, side='bottom') + if options.da_mode == "bvp": + deriv = get_J_sc_diffusion_vs_WL( + xa, xb, g, d_bottom, l_bottom, min_bot, s_bottom, wl, ph, side="bottom" + ) else: - xbb = xb - (xb - xa)/1001. - deriv = get_J_sc_diffusion_green(xa, xbb, g, d_bottom, l_bottom, min_bot, s_bottom, ph, side='bottom') + xbb = xb - (xb - xa) / 1001.0 + deriv = get_J_sc_diffusion_green( + xa, xbb, g, d_bottom, l_bottom, s_bottom, ph, side="bottom" + ) + j_sc_bot = d_bottom * abs(deriv) # The contribution from the SCR (includes the intrinsic region, if present). @@ -602,10 +756,15 @@ def qe_depletion(junction, options): # The total light absorbed, but not necessarily collected, is: xa = cum_widths[id_top] xb = cum_widths[id_bottom + 1] - zz = np.linspace(xa, xb, 10001) + zz = np.linspace(xa, xb, 1001) gg = g(zz) * ph current_absorbed = np.trapz(gg, zz, axis=0) + # why does this happen sometimes? + # j_sc_top[j_sc_top < 0] = 0 + # j_sc_bot[j_sc_bot < 0] = 0 + # j_sc_scr[j_sc_scr < 0] = 0 + # Now, we put everything together j_sc = j_sc_top + j_sc_bot + j_sc_scr @@ -614,22 +773,54 @@ def qe_depletion(junction, options): eqe_base = j_sc_bot / ph eqe_scr = j_sc_scr / ph - iqe = j_sc / current_absorbed - iqe[np.isnan(iqe)] = 0 # if zero current_absorbed, get NaN in previous line; want 0 IQE + iqe = np.divide(j_sc, current_absorbed, + out=np.zeros_like(j_sc), where=current_absorbed!=0) junction.iqe = interp1d(wl, iqe) - junction.eqe = interp1d(wl, eqe, kind='linear', bounds_error=False, assume_sorted=True, - fill_value=(eqe[0], eqe[-1])) - junction.eqe_emitter = interp1d(wl, eqe_emitter, kind='linear', bounds_error=False, assume_sorted=True, - fill_value=(eqe_emitter[0], eqe_emitter[-1])) - junction.eqe_base = interp1d(wl, eqe_base, kind='linear', bounds_error=False, assume_sorted=True, - fill_value=(eqe_base[0], eqe_base[-1])) - junction.eqe_scr = interp1d(wl, eqe_scr, kind='linear', bounds_error=False, assume_sorted=True, - fill_value=(eqe_scr[0], eqe_scr[-1])) + junction.eqe = interp1d( + wl, + eqe, + kind="linear", + bounds_error=False, + assume_sorted=True, + fill_value=(eqe[0], eqe[-1]), + ) + junction.eqe_emitter = interp1d( + wl, + eqe_emitter, + kind="linear", + bounds_error=False, + assume_sorted=True, + fill_value=(eqe_emitter[0], eqe_emitter[-1]), + ) + junction.eqe_base = interp1d( + wl, + eqe_base, + kind="linear", + bounds_error=False, + assume_sorted=True, + fill_value=(eqe_base[0], eqe_base[-1]), + ) + junction.eqe_scr = interp1d( + wl, + eqe_scr, + kind="linear", + bounds_error=False, + assume_sorted=True, + fill_value=(eqe_scr[0], eqe_scr[-1]), + ) + junction.qe = State( + { + "WL": wl, + "IQE": junction.iqe(wl), + "EQE": junction.eqe(wl), + "EQE_emitter": junction.eqe_emitter(wl), + "EQE_base": junction.eqe_base(wl), + "EQE_scr": junction.eqe_scr(wl), + } + ) - junction.qe = State({'WL': wl, 'IQE': junction.iqe(wl), 'EQE': junction.eqe(wl), 'EQE_emitter': junction.eqe_emitter(wl), - 'EQE_base': junction.eqe_base(wl), 'EQE_scr': junction.eqe_scr(wl)}) def get_J_sc_SCR_vs_WL(xa, xb, g, wl, ph): zz = np.linspace(xa, xb, 1001, endpoint=False) @@ -639,69 +830,125 @@ def get_J_sc_SCR_vs_WL(xa, xb, g, wl, ph): return out -def get_J_sc_diffusion_vs_WL(xa, xb, g, D, L, y0, S, wl, ph, side='top'): - zz = np.linspace(xa, xb, 1001, endpoint=False) # excluding the last point - depending on the mesh/floating point errors, sometimes this is actually in the next layer - gg = g(zz) * ph +def get_J_sc_diffusion_vs_WL(xa, xb, g, D, L, y0, S, wl, ph, side="top"): + zz = np.linspace(xa, xb, 1001, endpoint=False) + # excluding the last point - depending on the mesh/floating point errors, + # sometimes this is actually in the next layer + + gg = g(zz) * ph # generation rate * photon flux + out = np.zeros_like(wl) + sol_success = np.zeros_like(wl) # keep track of whether solve_bvp is converging for i in range(len(wl)): - - if np.all(gg[:,i] == 0): # no reason to solve anything if no generation at this wavelength + if np.all( + gg[:, i] == 0 + ): # no reason to solve anything if no generation at this wavelength out[i] = 0 + sol_success[i] = 0 else: - A = lambda x: np.interp(x, zz, gg[:, i]) / D + y0 / L ** 2 + + def A(x): + return np.interp(x, zz, gg[:, i]) / D + y0 / L**2 + + # generation and n0/p0 term in differential equation + # (eq. 6.15 & 6.20 in Jenny Nelson, Physics of Solar Cells) def fun(x, y): - out1 = y[1] - out2 = y[0] / L ** 2 - A(x) + # differential equation (eq. 6.15 & 6.20 in Jenny Nelson, + # Physics of Solar Cells) + # solve_bvp solves equation of form: + # dy / dx = f(x, y, p), a <= x <= b + # in this case y = [n or p, dn/dx or dp/dx] + # y[0] = carrier concentration (n or p) + # y[1] = carrier concentration gradient (dn/dx or dp/dx) + out1 = y[1] # by definition! dy/dx = dy/dx + + out2 = y[0] / L**2 - A(x) + # actually solving the differential equation (6.15 & 6.20) + return np.vstack((out1, out2)) - if side == 'top': + # boundary conditions for solve_bvp: + if side == "top": + def bc(ya, yb): left = ya[1] - S / D * (ya[0] - y0) - right = yb[0] + # eq. 6.18 - b.c. at front of junction (surface recombination) + + right = yb[0] - y0 + # eq. 6.17 - b.c. edge of depletion region, top half of junction + # added - y0 (generally very small so makes almost no difference) return np.array([left, right]) + else: + def bc(ya, yb): - left = ya[0] - right = yb[1] - S / D * (yb[0] - y0) + left = ya[0] - y0 + # eq. 6.21 - b.c. edge of depletion region, bottom half of junction + # added - y0 (generally very small so makes almost no difference) + + right = yb[1] + S / D * (yb[0] - y0) + # eq. 6.22 - b.c. at back of junction (surface recombination) + # changed sign! Current is going the other way return np.array([left, right]) guess = y0 * np.ones((2, zz.size)) guess[1] = np.zeros_like(guess[0]) - solution = solve_bvp(fun, bc, zz, guess) - if side == 'top': + solution = solve_bvp(fun, bc, zz, guess, max_nodes=2 * zz.shape[0]) + # increase max_nodes to avoid "too many mesh points" message + + sol_success[i] = solution.status + + if side == "top": out[i] = solution.y[1][-1] + # current at edge of depletion region (top half of junction), eq. 6.33 else: out[i] = solution.y[1][0] + # current at edge of depletion region (bottom half of junction), eq 6.38 + + # give a warning f any of the solution statuses are not 0 using warnings.warn: + if np.any(sol_success != 0): + warnings.warn( + "Depletion approximation (DA) EQE calculation: " + "solve_bvp did not converge as expected for some wavelengths", + RuntimeWarning, + ) return out def get_depletion_widths(junction, es, Vbi, V, Na, Nd, xi): - if not hasattr(junction, "wp") or not hasattr(junction, "wn"): - - if hasattr(junction, "depletion_approximation") and junction.depletion_approximation == "one-sided abrupt": + if ( + hasattr(junction, "depletion_approximation") + and junction.depletion_approximation == "one-sided abrupt" + ): print("using one-sided abrupt junction approximation for depletion width") one_sided = True else: one_sided = False - if one_sided: - science_reference("Sze abrupt junction approximation", - "Sze: The Physics of Semiconductor Devices, 2nd edition, John Wiley & Sons, Inc (2007)") + science_reference( + "Sze abrupt junction approximation", + "Sze: The Physics of Semiconductor Devices, " + "2nd edition, John Wiley & Sons, Inc (2007)", + ) wn = np.sqrt(2 * es * (Vbi - V) / (q * Nd)) wp = np.sqrt(2 * es * (Vbi - V) / (q * Na)) else: - wn = (-xi + np.sqrt(xi ** 2 + 2. * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / (1 + Nd / Na) - wp = (-xi + np.sqrt(xi ** 2 + 2. * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / (1 + Na / Nd) + wn = ( + -xi + np.sqrt(xi**2 + 2.0 * es * (Vbi - V) / q * (1 / Na + 1 / Nd)) + ) / (1 + Nd / Na) + wp = ( + -xi + np.sqrt(xi**2 + 2.0 * es * (Vbi - V) / q * (1 / Na + 1 / Nd)) + ) / (1 + Na / Nd) wn = wn if not hasattr(junction, "wn") else junction.wn wp = wp if not hasattr(junction, "wp") else junction.wp - return wn, wp \ No newline at end of file + return wn, wp diff --git a/solcore/optics/rcwa.py b/solcore/optics/rcwa.py index e922c56b..cbc1f91f 100755 --- a/solcore/optics/rcwa.py +++ b/solcore/optics/rcwa.py @@ -186,14 +186,17 @@ def calculate_absorption_rcwa(tmm_out, initial=1): def diff_absorption(z): idx = all_z.searchsorted(z) - idx = np.where(idx <= len(all_z) - 2, idx, len(all_z) - 2) - try: - z1 = all_z[idx] - z2 = all_z[idx + 1] + idx = np.where(idx <= len(all_z) - 1, idx, len(all_z) - 1) + idx = np.where(idx > 0, idx, 1) - f = (z - z1) / (z2 - z1) + try: + z1 = all_z[idx - 1] + z2 = all_z[idx] - out = f * all_abs[:, idx] + (1 - f) * all_abs[:, idx + 1] + f = np.divide(z - z1, z2 - z1, + out=np.zeros_like(z), where=np.abs(z2-z1) > 1e-12) + # this is to avoid divide by zero errors (|f| gets very large) when z1 = z2 + out = (1-f) * all_abs[:, idx - 1] + f * all_abs[:, idx] except IndexError: out = all_abs[:, idx] diff --git a/solcore/optics/tmm.py b/solcore/optics/tmm.py index 4176a979..4cb71f69 100755 --- a/solcore/optics/tmm.py +++ b/solcore/optics/tmm.py @@ -191,14 +191,17 @@ def calculate_absorption_tmm(tmm_out, initial=1): def diff_absorption(z): idx = all_z.searchsorted(z) - idx = np.where(idx <= len(all_z) - 2, idx, len(all_z) - 2) - try: - z1 = all_z[idx] - z2 = all_z[idx + 1] + idx = np.where(idx <= len(all_z) - 1, idx, len(all_z) - 1) + idx = np.where(idx > 0, idx, 1) - f = (z - z1) / (z2 - z1) + try: + z1 = all_z[idx - 1] + z2 = all_z[idx] - out = f * all_abs[:, idx] + (1 - f) * all_abs[:, idx + 1] + f = np.divide(z - z1, z2 - z1, + out=np.zeros_like(z), where=np.abs(z2-z1) > 1e-12) + # this is to avoid divide by zero errors (|f| gets very large) when z1 = z2 + out = (1-f) * all_abs[:, idx - 1] + f * all_abs[:, idx] except IndexError: out = all_abs[:, idx] diff --git a/tests/data/substrate_presence_profile.csv b/tests/data/substrate_presence_profile.csv index 17476fa1..6f1b62f2 100644 --- a/tests/data/substrate_presence_profile.csv +++ b/tests/data/substrate_presence_profile.csv @@ -1,20 +1,100 @@ -4.678228132479161024e+07,4.315724724342442118e+06,9.390180792711710092e+05 -6.938265521053141856e+04,2.534442823504240718e+06,8.420748909911215305e+05 -1.028983513619573813e+02,1.488370939467134420e+06,7.551400045618610457e+05 -1.562434525670335639e-01,8.757319686803899240e+05,6.774454794315428007e+05 -2.317193728623085421e-04,5.142803222402390675e+05,6.075067536990798544e+05 -3.518454905298698402e-07,3.025937161331444513e+05,5.450018546981358668e+05 -5.218137457066760257e-10,1.777004910771587165e+05,4.887364639479226316e+05 -7.923214328189151629e-13,1.045559135467049928e+05,4.384515514561979799e+05 -1.175081342546242938e-15,6.140126854159259528e+04,3.931862973123444244e+05 +4.733311893484660983e+07,4.319855093890735880e+06,9.392019805473309243e+05 +1.433569343737236224e+07,3.918113517695821822e+06,9.206209699882965069e+05 +4.341823749259334989e+06,3.553733442558917683e+06,9.024075628909341758e+05 +1.314999465534352930e+06,3.223240292678207159e+06,8.845544866538200295e+05 +3.982711650564905722e+05,2.923482623681675177e+06,8.670546125553576276e+05 +1.206235357108666212e+05,2.651602071791600436e+06,8.499009529072890291e+05 +3.653298392024997156e+04,2.405006097680339590e+06,8.330866582645147573e+05 +1.106466166322542267e+04,2.181343265114484355e+06,8.166050146901180269e+05 +3.351128125255747364e+03,1.978480818655225448e+06,8.004494410744949710e+05 +1.014948084320520138e+03,1.794484346605473198e+06,7.846134865075241541e+05 +3.073948058514957893e+02,1.627599335278353188e+06,7.690908277027227450e+05 +9.309988028029249563e+01,1.476234438696490135e+06,7.538752664723662892e+05 +2.819691669149860758e+01,1.338946304189113434e+06,7.389607272525557783e+05 +8.539924294910480995e+00,1.214425809190372005e+06,7.243412546772505157e+05 +2.586463303800504487e+00,1.101485577998876572e+06,7.100110112002976239e+05 +7.833548239111411959e-01,9.990486594636142254e+05,6.959642747645021882e+05 +2.372524051650460597e-01,9.061382576315053739e+05,6.821954365168175427e+05 +7.185593187545590010e-02,8.218684174324397463e+05,6.686989985687333392e+05 +2.176278785497093651e-02,7.454355765844633570e+05,6.554695718009753618e+05 +6.591227949429429890e-03,6.761109031616884749e+05,6.425018737116323318e+05 +1.996264274110634996e-03,6.132333457592495251e+05,6.297907263068602188e+05 +6.046020811613651543e-04,5.562033299846298760e+05,6.173310540333124809e+05 +1.831138277329821655e-04,5.044770411677911761e+05,6.051178817514748080e+05 +5.545906510389837243e-05,4.575612387724643340e+05,5.931463327490995871e+05 +1.679669542660556649e-05,4.150085530610304559e+05,5.814116267939351965e+05 +5.087170163048927641e-06,3.764132261936101131e+05,5.699090786709993845e+05 +1.540741563820887816e-06,3.414072270042708842e+05,5.586340953931099502e+05 +4.666413652445578704e-07,3.096567451996232849e+05,5.475821744118933566e+05 +1.413307199748011406e-07,2.808590214768872829e+05,5.367489027025871910e+05 +4.280453764833539580e-08,2.547394527896714862e+05,5.261299545469122240e+05 +1.296411727363517687e-08,2.310489738498722145e+05,5.157210898058146704e+05 +3.926413128118128925e-09,2.095616821467235277e+05,5.055181522263857769e+05 +1.189183669608921656e-09,1.900726838361993723e+05,4.955170677822721191e+05 +3.601652277964010406e-10,1.723961399600986624e+05,4.857138430469194427e+05 +1.090823605497679310e-10,1.563634943643946026e+05,4.761045635990040610e+05 +3.303749989089396049e-11,1.418218664190394629e+05,4.666853924594046548e+05 +1.000598212948384980e-11,1.286325932128917484e+05,4.574525685591022484e+05 +3.030485232478055257e-12,1.166699073227596091e+05,4.484024052373884479e+05 +9.178348040705698116e-13,1.058197375483344222e+05,4.395312887697826372e+05 +2.779820624115061361e-13,9.597862117733397463e+04,4.308356769250796642e+05 +8.419163121982884747e-14,8.705271740868478082e+04,4.223120975509350537e+05 +2.549887241960870891e-14,7.895691252616398560e+04,4.139571471874407725e+05 +7.722766933308788340e-15,7.161400828982375970e+04,4.057674897081266972e+05 +2.338970698067280225e-15,6.495398580604165181e+04,3.977398549878461054e+05 +7.083967077266117678e-16,5.891333785677807464e+04,3.898710375970192836e+05 +2.145498362086242311e-16,5.343446332141642051e+04,3.821578955217087059e+05 +6.498000502259175298e-17,4.846511791665053897e+04,3.745973489090132643e+05 +1.968027640998048892e-17,4.395791601690936659e+04,3.671863788372886484e+05 +5.960498004782220598e-18,3.986987880488957308e+04,3.599220261106921826e+05 0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00 -4.678228132479161769e+07,4.279038096903639846e+06,4.800241887743351981e+05 -6.938265521053146222e+04,2.543409179176834412e+06,6.713720965915143024e+05 -1.028983513619573813e+02,1.526295994198802859e+06,9.697707346578538418e+05 -1.562434525670335916e-01,9.069480006673237076e+05,4.948071984873304609e+05 -2.317193728623079458e-04,5.107409987542123999e+05,2.564359704550021270e+05 -3.518454905301201384e-07,2.669226631366466754e+05,7.103941144934552722e+05 -5.218137468521799330e-10,1.386796889530614717e+05,6.512383755456801737e+05 -7.923233578995703218e-13,9.546609559986203385e+04,1.298926754515937500e+05 -1.177398333859714643e-15,9.478370102525260882e+04,3.286310229598983424e+05 +4.733311893484663963e+07,4.282667181705768220e+06,4.793209129266348318e+05 +1.433569343737237528e+07,3.935690679749053903e+06,6.397875156054251129e+05 +4.341823749259336852e+06,3.579456892333367839e+06,8.857589856832476798e+05 +1.314999465534353629e+06,3.200067677556885872e+06,1.049587588847929379e+06 +3.982711650564908050e+05,2.870063670169110410e+06,1.016821597180192708e+06 +1.206235357108667085e+05,2.633910158631002065e+06,8.000482200601405930e+05 +3.653298392025000066e+04,2.436398019059109502e+06,5.307684266531388275e+05 +1.106466166322542995e+04,2.202226655436693691e+06,3.745463389090991113e+05 +3.351128125255748728e+03,1.946344973785133101e+06,4.241806074784563971e+05 +1.014948084320520593e+03,1.748102449458159506e+06,6.398257432524125325e+05 +3.073948058514960167e+02,1.628281805793287233e+06,8.746982502157207346e+05 +9.309988028029258089e+01,1.514999489191112109e+06,9.697813607246498577e+05 +2.819691669149862179e+01,1.350031086154998746e+06,8.564229267356696073e+05 +8.539924294910484548e+00,1.174698827616490191e+06,6.006282428604851011e+05 +2.586463303800505820e+00,1.065377320947901346e+06,3.604692009400026873e+05 +7.833548239111418621e-01,1.016159878644248471e+06,2.840902633087652503e+05 +2.372524051650461985e-01,9.461603418066937011e+05,4.141198646578903426e+05 +7.185593187545594174e-02,8.203893510514490772e+05,6.601629516958636232e+05 +2.176278785497095386e-02,7.013399701677301200e+05,8.571230679587048944e+05 +6.591227949429432492e-03,6.532354586315717315e+05,8.717234394591208547e+05 +1.996264274110634996e-03,6.437729508706863271e+05,6.884975104752674233e+05 +6.046020811613676479e-04,5.919862430269103497e+05,4.197660375614132499e+05 +1.831138277329805663e-04,4.898159013529409422e+05,2.331632098128859361e+05 +5.545906510389810137e-05,4.134773818447795347e+05,2.434733165245923738e+05 +1.679669542660764002e-05,4.074220648051212775e+05,4.385671329386560828e+05 +5.087170163046724508e-06,4.164172903233512188e+05,6.872252515255699400e+05 +1.540741563821539820e-06,3.684833158367379219e+05,8.238009631949120667e+05 +4.666413652459648980e-07,2.832223115049914923e+05,7.550988686288984027e+05 +1.413307199724239955e-07,2.416605518477708101e+05,5.205035274522408145e+05 +4.280453764984695256e-08,2.632216799270879128e+05,2.665653751195846125e+05 +1.296411727415697067e-08,2.759304174853296718e+05,1.524413162346418831e+05 +3.926413125966113510e-09,2.249153460444410739e+05,2.472016048313304200e+05 +1.189183671728786215e-09,1.550734522304506972e+05,4.852626470146932406e+05 +3.601652273462324458e-10,1.428458973650341504e+05,7.087012693216985790e+05 +1.090823589825075966e-10,1.803060424017931800e+05,7.691483292787405662e+05 +3.303750226490449633e-11,1.868984017404905171e+05,6.235178151732992847e+05 +1.000598078259219666e-11,1.310681833910598798e+05,3.618051140898090089e+05 +3.030484511810292531e-12,7.778873336641103378e+04,1.490184962617135025e+05 +9.178370356013072330e-13,9.008972001731641649e+04,1.187108307550786994e+05 +2.779800438305800518e-13,1.335438444675833161e+05,2.868070297880258877e+05 +8.419187756419852140e-14,1.281945786781990901e+05,5.413274972636202583e+05 +2.550058869804806459e-14,6.937182014481653459e+04,7.146382330844267271e+05 +7.720414077230066799e-15,3.463484467252287141e+04,6.916103841645757202e+05 +2.340143486525631316e-15,6.633852481738760252e+04,4.839025086069555255e+05 +7.093128769140944253e-16,1.074652032336971315e+05,2.224370835260323947e+05 +2.122620547953067387e-16,8.778946629494507215e+04,7.292548347354170983e+04 +6.690399351656450453e-17,3.003277296075194317e+04,1.289463619939968776e+05 +1.971148805515392964e-17,1.542559183352572836e+04,3.515654477500335197e+05 +4.345357146007531094e-18,6.077048092252996867e+04,5.944553968313524965e+05 0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00 diff --git a/tests/test_absorption.py b/tests/test_absorption.py index d2a1ce25..f796b895 100755 --- a/tests/test_absorption.py +++ b/tests/test_absorption.py @@ -243,7 +243,7 @@ def test_substrate_presence_profile(): }, ) - z_pos = np.linspace(0, my_structure.width, 10) + z_pos = np.linspace(0, my_structure.width, 50) profile_subs = my_structure[0].absorbed(z_pos) diff --git a/tests/test_depletion_approximation.py b/tests/test_depletion_approximation.py index ababd226..1e119e89 100644 --- a/tests/test_depletion_approximation.py +++ b/tests/test_depletion_approximation.py @@ -1,89 +1,101 @@ -from pytest import approx, mark, raises +from pytest import approx, raises import numpy as np from solcore.constants import kb, q, vacuum_permittivity - def test_get_j_dark(): from solcore.analytic_solar_cells.depletion_approximation import get_j_dark x = np.power(10, np.random.uniform(-8, -5)) xi = np.power(10, np.random.uniform(-8, -5)) - l = np.power(10, np.random.uniform(-9, -6)) + L = np.power(10, np.random.uniform(-9, -6)) s = np.power(10, np.random.uniform(1, 3)) d = np.power(10, np.random.uniform(-5, 0)) - Vbi = np.random.uniform(0.1,5) + Vbi = np.random.uniform(0.1, 5) minor = np.power(10, np.random.uniform(-7, -4)) - T = np.random.uniform(0.1,400) - es = np.random.uniform(1,20)*vacuum_permittivity + T = np.random.uniform(0.1, 400) + es = np.random.uniform(1, 20) * vacuum_permittivity Na = np.power(10, np.random.uniform(22, 25)) Nd = np.power(10, np.random.uniform(22, 25)) V = np.linspace(-6, 4, 20) V = np.where(V < Vbi - 0.001, V, Vbi - 0.001) - wn = (-xi + np.sqrt(xi ** 2 + 2. * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / (1 + Nd / Na) - wp = (-xi + np.sqrt(xi ** 2 + 2. * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / (1 + Na / Nd) + wn = (-xi + np.sqrt(xi**2 + 2.0 * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / ( + 1 + Nd / Na + ) + wp = (-xi + np.sqrt(xi**2 + 2.0 * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / ( + 1 + Na / Nd + ) w = wn + wp + xi - expected = (q*d*minor/l)*(np.exp(q*V/(kb*T))-1)*((((s*l)/d)*np.cosh((x-w)/l)+np.sinh((x-w)/l))/ - (((s*l)/d)*np.sinh((x-w)/l)+np.cosh((x-w)/l))) + expected = ( + (q * d * minor / L) + * (np.exp(q * V / (kb * T)) - 1) + * ( + (((s * L) / d) * np.cosh((x - w) / L) + np.sinh((x - w) / L)) + / (((s * L) / d) * np.sinh((x - w) / L) + np.cosh((x - w) / L)) + ) + ) - result = get_j_dark(x, w, l, s, d, V, minor, T) + result = get_j_dark(x, w, L, s, d, V, minor, T) assert result == approx(expected, nan_ok=True) def test_factor(): from solcore.analytic_solar_cells.depletion_approximation import factor - T = np.random.uniform(0.1,400) - Vbi = np.random.uniform(0.1,5) + + T = np.random.uniform(0.1, 400) + Vbi = np.random.uniform(0.1, 5) tp = np.power(10, np.random.uniform(-10, -5)) tn = np.power(10, np.random.uniform(-10, -5)) dEt = 0 - kT = kb*T + kT = kb * T V = np.linspace(-6, 4, 20) V = np.where(V < Vbi - 0.001, V, Vbi - 0.001) m = V >= -1120 * kT / q V = V[m] - b= np.exp(-V*q/(2*kT))*np.cosh((q*dEt/kT)+(1/2)*np.log(tp/tn)) + b = np.exp(-V * q / (2 * kT)) * np.cosh((q * dEt / kT) + (1 / 2) * np.log(tp / tn)) - z1 = np.sqrt(tp/tn)*np.exp(-q*(Vbi - V)/(kT * 2)) - z2 = np.sqrt(tp/tn)*np.exp(q*(Vbi - V)/(kT * 2)) + z1 = np.sqrt(tp / tn) * np.exp(-q * (Vbi - V) / (kT * 2)) + z2 = np.sqrt(tp / tn) * np.exp(q * (Vbi - V) / (kT * 2)) expected = np.zeros(b.shape) # For b values < 1 inds = b < 1 # use tan-1 formulation without issue - upper = np.arctan((b[inds]+z2[inds])/np.sqrt(1-b[inds]**2)) - lower = np.arctan((b[inds]+z1[inds])/np.sqrt(1-b[inds]**2)) + upper = np.arctan((b[inds] + z2[inds]) / np.sqrt(1 - b[inds] ** 2)) + lower = np.arctan((b[inds] + z1[inds]) / np.sqrt(1 - b[inds] ** 2)) - expected[inds] = (upper - lower)/np.sqrt(1-b[inds]**2) + expected[inds] = (upper - lower) / np.sqrt(1 - b[inds] ** 2) # for b values >=1, log formulation inds = (b >= 1) & (b <= 1e7) - upper = np.log(abs(z2[inds] + b[inds] - np.sqrt(b[inds]**2-1))) - \ - np.log(abs(z2[inds] + b[inds] + np.sqrt(b[inds]**2-1))) + upper = np.log(abs(z2[inds] + b[inds] - np.sqrt(b[inds] ** 2 - 1))) - np.log( + abs(z2[inds] + b[inds] + np.sqrt(b[inds] ** 2 - 1)) + ) - lower = np.log(abs(z1[inds] + b[inds] - np.sqrt(b[inds]**2-1))) - \ - np.log(abs(z1[inds] + b[inds] + np.sqrt(b[inds]**2-1))) - expected[inds] = (upper - lower)/(2*np.sqrt(b[inds]**2 - 1)) + lower = np.log(abs(z1[inds] + b[inds] - np.sqrt(b[inds] ** 2 - 1))) - np.log( + abs(z1[inds] + b[inds] + np.sqrt(b[inds] ** 2 - 1)) + ) + expected[inds] = (upper - lower) / (2 * np.sqrt(b[inds] ** 2 - 1)) # limit as b gets very large: sqrt(b**2-1) = b inds = b > 1e7 - upper = np.log(z2[inds]) - np.log(z2[inds] + 2*b[inds]) - lower = np.log(z1[inds]) - np.log(z1[inds] + 2*b[inds]) + upper = np.log(z2[inds]) - np.log(z2[inds] + 2 * b[inds]) + lower = np.log(z1[inds]) - np.log(z1[inds] + 2 * b[inds]) - expected[inds] = (upper - lower) / (2*b[inds]) + expected[inds] = (upper - lower) / (2 * b[inds]) result = factor(V, Vbi, tp, tn, kT, dEt) @@ -93,18 +105,18 @@ def test_factor(): def test_forward(): from solcore.analytic_solar_cells.depletion_approximation import factor, forward - T = np.random.uniform(0.1,400) - Vbi = np.random.uniform(0.1,5) + T = np.random.uniform(0.1, 400) + Vbi = np.random.uniform(0.1, 5) tp = np.power(10, np.random.uniform(-10, -5)) tn = np.power(10, np.random.uniform(-10, -5)) dEt = 0 - kT = kb*T + kT = kb * T V = np.linspace(-6, 4, 20) V = np.where(V < Vbi - 0.001, V, Vbi - 0.001) ni = np.power(10, np.random.uniform(2, 9)) - es = np.random.uniform(1,20)*vacuum_permittivity + es = np.random.uniform(1, 20) * vacuum_permittivity Na = np.power(10, np.random.uniform(22, 25)) Nd = np.power(10, np.random.uniform(22, 25)) xi = np.power(10, np.random.uniform(-8, -5)) @@ -112,14 +124,26 @@ def test_forward(): m = V >= -1120 * kT / q V = V[m] - wn = (-xi + np.sqrt(xi ** 2 + 2. * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / (1 + Nd / Na) - wp = (-xi + np.sqrt(xi ** 2 + 2. * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / (1 + Na / Nd) + wn = (-xi + np.sqrt(xi**2 + 2.0 * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / ( + 1 + Nd / Na + ) + wp = (-xi + np.sqrt(xi**2 + 2.0 * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / ( + 1 + Na / Nd + ) w = wn + wp + xi f_b = factor(V, Vbi, tp, tn, kT, dEt) - expected = 2 * q * ni * w / np.sqrt(tn * tp) * \ - np.sinh(q*V / (2 * kT)) / (q * (Vbi - V) / kT) * f_b + expected = ( + 2 + * q + * ni + * w + / np.sqrt(tn * tp) + * np.sinh(q * V / (2 * kT)) + / (q * (Vbi - V) / kT) + * f_b + ) result = forward(ni, V, Vbi, tp, tn, w, kT, dEt) @@ -129,24 +153,28 @@ def test_forward(): def test_get_J_srh(): from solcore.analytic_solar_cells.depletion_approximation import forward, get_Jsrh - T = np.random.uniform(0.1,400) - Vbi = np.random.uniform(0.1,5) + T = np.random.uniform(0.1, 400) + Vbi = np.random.uniform(0.1, 5) tp = np.power(10, np.random.uniform(-10, -5)) tn = np.power(10, np.random.uniform(-10, -5)) dEt = 0 - kT = kb*T + kT = kb * T V = np.linspace(-6, 4, 20) V = np.where(V < Vbi - 0.001, V, Vbi - 0.001) ni = np.power(10, np.random.uniform(2, 9)) - es = np.random.uniform(1,20)*vacuum_permittivity + es = np.random.uniform(1, 20) * vacuum_permittivity Na = np.power(10, np.random.uniform(22, 25)) Nd = np.power(10, np.random.uniform(22, 25)) xi = np.power(10, np.random.uniform(-8, -5)) - wn = (-xi + np.sqrt(xi ** 2 + 2. * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / (1 + Nd / Na) - wp = (-xi + np.sqrt(xi ** 2 + 2. * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / (1 + Na / Nd) + wn = (-xi + np.sqrt(xi**2 + 2.0 * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / ( + 1 + Nd / Na + ) + wp = (-xi + np.sqrt(xi**2 + 2.0 * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / ( + 1 + Na / Nd + ) w = wn + wp + xi @@ -166,40 +194,40 @@ def test_get_J_sc_diffusion_top(): from solcore.interpolate import interp1d from solcore.light_source import LightSource - D = np.power(10, np.random.uniform(-5, 0)) # Diffusion coefficient - L = np.power(10, np.random.uniform(-9, -6)) # Diffusion length - minority = np.power(10, np.random.uniform(-7, -4))# minority carrier density - s = np.power(10, np.random.uniform(0, 3)) # surface recombination velocity + D = np.power(10, np.random.uniform(-5, 0)) # Diffusion coefficient + L = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + minority = np.power(10, np.random.uniform(-7, -4)) # minority carrier density + s = np.power(10, np.random.uniform(0, 3)) # surface recombination velocity light_source = LightSource(source_type="standard", version="AM1.5g") - wl = np.linspace(300, 1800, 50)*1e-9 - wl_ls, phg = light_source.spectrum(output_units='photon_flux_per_m', x=wl) + wl = np.linspace(300, 1800, 50) * 1e-9 + wl_ls, phg = light_source.spectrum(output_units="photon_flux_per_m", x=wl) xa_nm = np.random.uniform(1, 1000) - xa = xa_nm*1e-9 - xb = np.random.uniform(xa_nm+1, 1100)*1e-9 + xa = xa_nm * 1e-9 + xb = np.random.uniform(xa_nm + 1, 1100) * 1e-9 - ## make a simple Beer-Lambert profile + # make a simple Beer-Lambert profile dist = np.linspace(0, xb, 1000) alphas = np.linspace(1e8, 10, len(wl)) - expn = np.exp(- alphas[:, None] * dist[None,:]) + expn = np.exp(-alphas[:, None] * dist[None, :]) - output = alphas[:, None]*expn + output = alphas[:, None] * expn output = output.T - gen_prof = interp1d(dist, output, axis = 0) - + gen_prof = interp1d(dist, output, axis=0) zz = np.linspace(xa, xb, 1002)[:-1] gg = gen_prof(zz) * phg g_vs_z = np.trapz(gg, wl, axis=1) - A = lambda x: np.interp(x, zz, g_vs_z) / D + minority / L ** 2 + def A(x): + return np.interp(x, zz, g_vs_z) / D + minority / L**2 def fun(x, y): out1 = y[1] - out2 = y[0] / L ** 2 - A(x) + out2 = y[0] / L**2 - A(x) return np.vstack((out1, out2)) def bc(ya, yb): @@ -207,15 +235,16 @@ def bc(ya, yb): right = yb[0] return np.array([left, right]) - guess = minority * np.ones((2, zz.size)) guess[1] = np.zeros_like(guess[0]) - solution = solve_bvp(fun, bc, zz, guess) + solution = solve_bvp(fun, bc, zz, guess, max_nodes=2 * zz.shape[0]) expected = solution.y[1][-1] - result = get_J_sc_diffusion(xa, xb, gen_prof, D, L, minority, s, wl, phg, side='top') + result = get_J_sc_diffusion( + xa, xb, gen_prof, D, L, minority, s, wl, phg, side="top" + ) assert result == approx(expected) @@ -226,56 +255,57 @@ def test_get_J_sc_diffusion_bottom(): from solcore.interpolate import interp1d from solcore.light_source import LightSource - D = np.power(10, np.random.uniform(-5, 0)) # Diffusion coefficient - L = np.power(10, np.random.uniform(-9, -6)) # Diffusion length - minority = np.power(10, np.random.uniform(-7, -4))# minority carrier density - s = np.power(10, np.random.uniform(0, 3)) # surface recombination velocity + D = np.power(10, np.random.uniform(-5, 0)) # Diffusion coefficient + L = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + minority = np.power(10, np.random.uniform(-7, -4)) # minority carrier density + s = np.power(10, np.random.uniform(0, 3)) # surface recombination velocity light_source = LightSource(source_type="standard", version="AM1.5g") - wl = np.linspace(300, 1800, 50)*1e-9 - wl_ls, phg = light_source.spectrum(output_units='photon_flux_per_m', x=wl) + wl = np.linspace(300, 1800, 50) * 1e-9 + wl_ls, phg = light_source.spectrum(output_units="photon_flux_per_m", x=wl) xa_nm = np.random.uniform(1, 1000) - xa = xa_nm*1e-9 - xb = np.random.uniform(xa_nm+1, 1100)*1e-9 + xa = xa_nm * 1e-9 + xb = np.random.uniform(xa_nm + 1, 1100) * 1e-9 - ## make a simple Beer-Lambert profile + # make a simple Beer-Lambert profile dist = np.linspace(0, xb, 1000) alphas = np.linspace(1e8, 10, len(wl)) - expn = np.exp(- alphas[:, None] * dist[None,:]) + expn = np.exp(-alphas[:, None] * dist[None, :]) - output = alphas[:, None]*expn + output = alphas[:, None] * expn output = output.T - gen_prof = interp1d(dist, output, axis = 0) - + gen_prof = interp1d(dist, output, axis=0) zz = np.linspace(xa, xb, 1002)[:-1] gg = gen_prof(zz) * phg g_vs_z = np.trapz(gg, wl, axis=1) - A = lambda x: np.interp(x, zz, g_vs_z) / D + minority / L ** 2 + def A(x): + return np.interp(x, zz, g_vs_z) / D + minority / L**2 def fun(x, y): out1 = y[1] - out2 = y[0] / L ** 2 - A(x) + out2 = y[0] / L**2 - A(x) return np.vstack((out1, out2)) def bc(ya, yb): left = ya[0] - right = yb[1] - s / D * (yb[0] - minority) + right = yb[1] + s / D * (yb[0] - minority) return np.array([left, right]) - guess = minority * np.ones((2, zz.size)) guess[1] = np.zeros_like(guess[0]) - solution = solve_bvp(fun, bc, zz, guess) + solution = solve_bvp(fun, bc, zz, guess, max_nodes=2 * zz.shape[0]) expected = solution.y[1][0] - result = get_J_sc_diffusion(xa, xb, gen_prof, D, L, minority, s, wl, phg, side='bottom') + result = get_J_sc_diffusion( + xa, xb, gen_prof, D, L, minority, s, wl, phg, side="bottom" + ) assert result == approx(expected) @@ -286,22 +316,22 @@ def test_get_J_sc_SCR(): from solcore.analytic_solar_cells.depletion_approximation import get_J_sc_SCR light_source = LightSource(source_type="standard", version="AM1.5g") - wl = np.linspace(300, 1800, 50)*1e-9 - wl_ls, phg = light_source.spectrum(output_units='photon_flux_per_m', x=wl) + wl = np.linspace(300, 1800, 50) * 1e-9 + wl_ls, phg = light_source.spectrum(output_units="photon_flux_per_m", x=wl) xa_nm = np.random.uniform(1, 1000) - xa = xa_nm*1e-9 - xb = np.random.uniform(xa_nm+1, 1100)*1e-9 + xa = xa_nm * 1e-9 + xb = np.random.uniform(xa_nm + 1, 1100) * 1e-9 - ## make a simple Beer-Lambert profile + # make a simple Beer-Lambert profile dist = np.linspace(0, xb, 1000) alphas = np.linspace(1e5, 10, len(wl)) - expn = np.exp(- alphas[:, None] * dist[None,:]) + expn = np.exp(-alphas[:, None] * dist[None, :]) - output = alphas[:, None]*expn + output = alphas[:, None] * expn output = output.T - gen_prof = interp1d(dist, output, axis = 0) + gen_prof = interp1d(dist, output, axis=0) zz = np.linspace(xa, xb, 1002)[:-1] gg = gen_prof(zz) * phg @@ -309,7 +339,7 @@ def test_get_J_sc_SCR(): result = get_J_sc_SCR(xa, xb, gen_prof, wl, phg) - assert expected == approx(result) + assert expected == approx(result) def test_get_J_sc_SCR_vs_WL(): @@ -318,22 +348,22 @@ def test_get_J_sc_SCR_vs_WL(): from solcore.analytic_solar_cells.depletion_approximation import get_J_sc_SCR_vs_WL light_source = LightSource(source_type="standard", version="AM1.5g") - wl = np.linspace(300, 1800, 50)*1e-9 - wl_ls, phg = light_source.spectrum(output_units='photon_flux_per_m', x=wl) + wl = np.linspace(300, 1800, 50) * 1e-9 + wl_ls, phg = light_source.spectrum(output_units="photon_flux_per_m", x=wl) xa_nm = np.random.uniform(1, 1000) - xa = xa_nm*1e-9 - xb = np.random.uniform(xa_nm+1, 1100)*1e-9 + xa = xa_nm * 1e-9 + xb = np.random.uniform(xa_nm + 1, 1100) * 1e-9 - ## make a simple Beer-Lambert profile + # make a simple Beer-Lambert profile dist = np.linspace(0, xb, 1000) alphas = np.linspace(1e5, 10, len(wl)) - expn = np.exp(- alphas[:, None] * dist[None,:]) + expn = np.exp(-alphas[:, None] * dist[None, :]) - output = alphas[:, None]*expn + output = alphas[:, None] * expn output = output.T - gen_prof = interp1d(dist, output, axis = 0) + gen_prof = interp1d(dist, output, axis=0) zz = np.linspace(xa, xb, 1002)[:-1] gg = gen_prof(zz) * phg @@ -341,97 +371,103 @@ def test_get_J_sc_SCR_vs_WL(): result = get_J_sc_SCR_vs_WL(xa, xb, gen_prof, wl, phg) - assert expected == approx(result) + assert expected == approx(result) def test_get_J_sc_diffusion_vs_WL_top(): - from solcore.analytic_solar_cells.depletion_approximation import get_J_sc_diffusion_vs_WL + from solcore.analytic_solar_cells.depletion_approximation import ( + get_J_sc_diffusion_vs_WL, + ) from scipy.integrate import solve_bvp from solcore.interpolate import interp1d from solcore.light_source import LightSource - D = np.power(10, np.random.uniform(-5, 0)) # Diffusion coefficient - L = np.power(10, np.random.uniform(-9, -6)) # Diffusion length - minority = np.power(10, np.random.uniform(-7, -4))# minority carrier density - s = np.power(10, np.random.uniform(0, 3)) # surface recombination velocity + D = np.power(10, np.random.uniform(-5, 0)) # Diffusion coefficient + L = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + minority = np.power(10, np.random.uniform(-7, -4)) # minority carrier density + s = np.power(10, np.random.uniform(0, 3)) # surface recombination velocity light_source = LightSource(source_type="standard", version="AM1.5g") - wl = np.linspace(300, 1800, 50)*1e-9 - wl_ls, phg = light_source.spectrum(output_units='photon_flux_per_m', x=wl) + wl = np.linspace(300, 1800, 50) * 1e-9 + wl_ls, phg = light_source.spectrum(output_units="photon_flux_per_m", x=wl) xa_nm = np.random.uniform(1, 1000) - xa = xa_nm*1e-9 - xb = np.random.uniform(xa_nm+1, 1100)*1e-9 + xa = xa_nm * 1e-9 + xb = np.random.uniform(xa_nm + 1, 1100) * 1e-9 - ## make a simple Beer-Lambert profile + # make a simple Beer-Lambert profile dist = np.linspace(0, xb, 1000) alphas = np.linspace(1e8, 10, len(wl)) - expn = np.exp(- alphas[:, None] * dist[None,:]) + expn = np.exp(-alphas[:, None] * dist[None, :]) - output = alphas[:, None]*expn + output = alphas[:, None] * expn output = output.T - gen_prof = interp1d(dist, output, axis = 0) - + gen_prof = interp1d(dist, output, axis=0) - zz = np.linspace(xa, xb, 1002)[:-1] + zz = np.linspace(xa, xb, 1002)[:-1] gg = gen_prof(zz) * phg expected = np.zeros_like(wl) for i in range(len(wl)): - - if np.all(gg[:,i] == 0): # no reason to solve anything if no generation at this wavelength + if np.all( + gg[:, i] == 0 + ): # no reason to solve anything if no generation at this wavelength expected[i] = 0 - A = lambda x: np.interp(x, zz, gg[:, i]) / D + minority / L ** 2 + def A(x): + return np.interp(x, zz, gg[:, i]) / D + minority / L**2 def fun(x, y): out1 = y[1] - out2 = y[0] / L ** 2 - A(x) + out2 = y[0] / L**2 - A(x) return np.vstack((out1, out2)) def bc(ya, yb): left = ya[1] - s / D * (ya[0] - minority) - right = yb[0] + right = yb[0] - minority return np.array([left, right]) guess = minority * np.ones((2, zz.size)) guess[1] = np.zeros_like(guess[0]) - solution = solve_bvp(fun, bc, zz, guess) + solution = solve_bvp(fun, bc, zz, guess, max_nodes=2 * zz.shape[0]) expected[i] = solution.y[1][-1] - result = get_J_sc_diffusion_vs_WL(xa, xb, gen_prof, D, L, minority, s, wl, phg, side='top') - + result = get_J_sc_diffusion_vs_WL( + xa, xb, gen_prof, D, L, minority, s, wl, phg, side="top" + ) assert result == approx(expected) def test_get_J_sc_diffusion_vs_WL_bottom(): - from solcore.analytic_solar_cells.depletion_approximation import get_J_sc_diffusion_vs_WL + from solcore.analytic_solar_cells.depletion_approximation import ( + get_J_sc_diffusion_vs_WL, + ) from scipy.integrate import solve_bvp from solcore.interpolate import interp1d from solcore.light_source import LightSource - D = np.power(10, np.random.uniform(-5, 0)) # Diffusion coefficient - L = np.power(10, np.random.uniform(-9, -6)) # Diffusion length - minority = np.power(10, np.random.uniform(-7, -4))# minority carrier density - s = np.power(10, np.random.uniform(0, 3)) # surface recombination velocity + D = np.power(10, np.random.uniform(-5, 0)) # Diffusion coefficient + L = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + minority = np.power(10, np.random.uniform(-7, -4)) # minority carrier density + s = np.power(10, np.random.uniform(0, 3)) # surface recombination velocity light_source = LightSource(source_type="standard", version="AM1.5g") wl = np.linspace(300, 1800, 50) * 1e-9 - wl_ls, phg = light_source.spectrum(output_units='photon_flux_per_m', x=wl) + wl_ls, phg = light_source.spectrum(output_units="photon_flux_per_m", x=wl) xa_nm = np.random.uniform(1, 1000) xa = xa_nm * 1e-9 xb = np.random.uniform(xa_nm + 1, 1100) * 1e-9 - ## make a simple Beer-Lambert profile + # make a simple Beer-Lambert profile dist = np.linspace(0, xb, 1000) alphas = np.linspace(1e8, 10, len(wl)) - expn = np.exp(- alphas[:, None] * dist[None, :]) + expn = np.exp(-alphas[:, None] * dist[None, :]) output = alphas[:, None] * expn output = output.T @@ -443,153 +479,192 @@ def test_get_J_sc_diffusion_vs_WL_bottom(): expected = np.zeros_like(wl) for i in range(len(wl)): - if np.all(gg[:, i] == 0): # no reason to solve anything if no generation at this wavelength + if np.all( + gg[:, i] == 0 + ): # no reason to solve anything if no generation at this wavelength expected[i] = 0 - A = lambda x: np.interp(x, zz, gg[:, i]) / D + minority / L ** 2 + def A(x): + return np.interp(x, zz, gg[:, i]) / D + minority / L**2 def fun(x, y): out1 = y[1] - out2 = y[0] / L ** 2 - A(x) + out2 = y[0] / L**2 - A(x) return np.vstack((out1, out2)) def bc(ya, yb): - left = ya[0] - right = yb[1] - s / D * (yb[0] - minority) + left = ya[0] - minority + right = yb[1] + s / D * (yb[0] - minority) return np.array([left, right]) guess = minority * np.ones((2, zz.size)) guess[1] = np.zeros_like(guess[0]) - solution = solve_bvp(fun, bc, zz, guess) + solution = solve_bvp(fun, bc, zz, guess, max_nodes=2 * zz.shape[0]) expected[i] = solution.y[1][0] - result = get_J_sc_diffusion_vs_WL(xa, xb, gen_prof, D, L, minority, s, wl, phg, side='bottom') + result = get_J_sc_diffusion_vs_WL( + xa, xb, gen_prof, D, L, minority, s, wl, phg, side="bottom" + ) assert result == approx(expected) def test_get_J_sc_diffusion_green_top(): - from solcore.analytic_solar_cells.depletion_approximation import get_J_sc_diffusion_green + from solcore.analytic_solar_cells.depletion_approximation import ( + get_J_sc_diffusion_green, + ) from solcore.light_source import LightSource - D = np.power(10, np.random.uniform(-5, 0)) # Diffusion coefficient - L = np.power(10, np.random.uniform(-9, -6)) # Diffusion length - minority = np.power(10, np.random.uniform(-7, -4))# minority carrier density - s = np.power(10, np.random.uniform(0, 3)) # surface recombination velocity + D = np.power(10, np.random.uniform(-5, 0)) # Diffusion coefficient + L = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + s = np.power(10, np.random.uniform(0, 3)) # surface recombination velocity light_source = LightSource(source_type="standard", version="AM1.5g") wl = np.linspace(300, 1800, 50) * 1e-9 - wl_ls, phg = light_source.spectrum(output_units='photon_flux_per_m', x=wl) + wl_ls, phg = light_source.spectrum(output_units="photon_flux_per_m", x=wl) xa_nm = np.random.uniform(1, 1000) xa = xa_nm * 1e-9 - xb = np.random.uniform(xa_nm+1, 1100) * 1e-9 + xb = np.random.uniform(xa_nm + 1, 1100) * 1e-9 beta1 = (xb - xa) / L - if beta1 > 1.e2: - xb = xa + 99. * L + if beta1 > 1.0e2: + xb = xa + 99.0 * L beta1 = (xb - xa) / L - ## make a simple Beer-Lambert profile + # make a simple Beer-Lambert profile alphas = np.linspace(1e8, 10, len(wl)) # theoreical value from Fonash, Solar cell device physics eq 4.10 # (the absorption profile has to be offset for xa != 0.) beta2 = (xb - xa) * alphas beta3 = L * s / D - c1 = beta2*beta2/(beta2*beta2 - beta1*beta1) - c2 = (beta3*beta1/beta2 + 1.) / (beta3 * np.sinh(beta1) + np.cosh(beta1)) - expected = (beta3*np.cosh(beta1) + np.sinh(beta1)) / (beta3*np.sinh(beta1) + np.cosh(beta1)) - expected *= beta1/beta2 - expected += 1. + c1 = beta2 * beta2 / (beta2 * beta2 - beta1 * beta1) + c2 = (beta3 * beta1 / beta2 + 1.0) / (beta3 * np.sinh(beta1) + np.cosh(beta1)) + expected = (beta3 * np.cosh(beta1) + np.sinh(beta1)) / ( + beta3 * np.sinh(beta1) + np.cosh(beta1) + ) + expected *= beta1 / beta2 + expected += 1.0 expected *= np.exp(-beta2) - expected = (expected - c2) * c1 * phg / D * np.exp(-alphas * xa) - expected -= minority / L * (beta3*np.cosh(beta1) + np.sinh(beta1)) / (beta3*np.sinh(beta1) + np.cosh(beta1)) - - result = get_J_sc_diffusion_green(xa, xb, lambda x: alphas * np.exp(-alphas * x), - D, L, minority, s, phg, side='top') - - assert beta1 < 1.e2 - assert result == approx(expected, rel=1.e-5) + expected = (expected - c2) * c1 * phg / D * np.exp(-alphas * xa) + + result = get_J_sc_diffusion_green( + xa, + xb, + lambda x: alphas * np.exp(-alphas * x), + D, + L, + s, + phg, + side="top", + ) + + assert beta1 < 1.0e2 + assert result == approx(expected, rel=1.0e-5) # test version for extreme L values - Lmax = np.log10((xb - xa) * 1.e-2) - L = np.power(10, np.random.uniform(Lmax - 2., Lmax)) + Lmax = np.log10((xb - xa) * 1.0e-2) + L = np.power(10, np.random.uniform(Lmax - 2.0, Lmax)) # theoretical value for Beer-Lambert profile - alphas_miL = alphas - 1. / L - idx_c = np.abs(alphas_miL) < 1.e-3 + alphas_miL = alphas - 1.0 / L + idx_c = np.abs(alphas_miL) < 1.0e-3 idx_i = np.logical_not(idx_c) - expected = np.zeros_like(alphas) + expected = np.exp(-alphas_miL * xa - xb / L) - np.exp(-alphas * xb) expected[idx_i] /= alphas_miL[idx_i] expected[idx_c] = np.exp(-alphas[idx_c] * xb) expected *= -phg * alphas / D - expected -= minority / L - result = get_J_sc_diffusion_green(xa, xb, lambda x: alphas * np.exp(-alphas * x), - D, L, minority, s, phg, side='top') + result = get_J_sc_diffusion_green( + xa, + xb, + lambda x: alphas * np.exp(-alphas * x), + D, + L, + s, + phg, + side="top", + ) - assert (xb - xa) / L > 1.e2 - assert result == approx(expected, rel=1.e-5) + assert (xb - xa) / L > 1.0e2 + assert result == approx(expected, rel=1.0e-5) def test_get_J_sc_diffusion_green_bottom(): - from solcore.analytic_solar_cells.depletion_approximation import get_J_sc_diffusion_green + from solcore.analytic_solar_cells.depletion_approximation import ( + get_J_sc_diffusion_green, + ) from solcore.light_source import LightSource - D = np.power(10, np.random.uniform(-5, 0)) # Diffusion coefficient - L = np.power(10, np.random.uniform(-9, -6)) # Diffusion length - minority = np.power(10, np.random.uniform(-7, -4))# minority carrier density - s = np.power(10, np.random.uniform(0, 3)) # surface recombination velocityx + D = np.power(10, np.random.uniform(-5, 0)) # Diffusion coefficient + L = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + s = np.power(10, np.random.uniform(0, 3)) # surface recombination velocityx light_source = LightSource(source_type="standard", version="AM1.5g") - wl = np.linspace(300, 1800, 50)*1e-9 - wl_ls, phg = light_source.spectrum(output_units='photon_flux_per_m', x=wl) + wl = np.linspace(300, 1800, 50) * 1e-9 + wl_ls, phg = light_source.spectrum(output_units="photon_flux_per_m", x=wl) xa_nm = np.random.uniform(1, 1000) - xa = xa_nm*1e-9 - xb = np.random.uniform(xa_nm+1, 1100)*1e-9 + xa = xa_nm * 1e-9 + xb = np.random.uniform(xa_nm + 1, 1100) * 1e-9 beta1 = (xb - xa) / L - if beta1 > 1.e2: - xb = xa + 99. * L + if beta1 > 1.0e2: + xb = xa + 99.0 * L beta1 = (xb - xa) / L - ## make a simple Beer-Lambert profile + # make a simple Beer-Lambert profile alphas = np.linspace(1e8, 10, len(wl)) # theoreical value from Fonash, Solar cell device physics eq 4.15 beta2 = (xb - xa) * alphas - beta3 = -L * s / D - c1 = beta2*beta2/(beta2*beta2 - beta1*beta1) * np.exp(-alphas * xa) - c2 = (beta3*beta1/beta2 - 1.) / (beta3 * np.sinh(beta1) + np.cosh(beta1)) - expected = (beta3*np.cosh(beta1) + np.sinh(beta1)) / (beta3*np.sinh(beta1) + np.cosh(beta1)) - expected *= -beta1/beta2 - expected += 1. - expected = (expected + c2 * np.exp(-beta2)) * c1 * phg / D - expected += minority / L * (beta3*np.cosh(beta1) + np.sinh(beta1)) / (beta3*np.sinh(beta1) + np.cosh(beta1)) - - result = get_J_sc_diffusion_green(xa, xb, lambda x: alphas * np.exp(-alphas * x), - D, L, minority, s, phg, side='bottom') - - assert beta1 < 1.e2 - assert result == approx(expected, rel=1.e-5) + beta3 = L * s / D + c1 = beta2 * beta2 / (beta2 * beta2 - beta1 * beta1) * np.exp(-alphas * xa) + c2 = (beta3 * beta1 / beta2 - 1.0) / (beta3 * np.sinh(beta1) + np.cosh(beta1)) + expected = (beta3 * np.cosh(beta1) + np.sinh(beta1)) / ( + beta3 * np.sinh(beta1) + np.cosh(beta1) + ) + expected *= -beta1 / beta2 + expected += 1.0 + expected = (expected + c2 * np.exp(-beta2)) * c1 * phg / D + + result = get_J_sc_diffusion_green( + xa, + xb, + lambda x: alphas * np.exp(-alphas * x), + D, + L, + s, + phg, + side="bottom", + ) + + assert beta1 < 1.0e2 + assert result == approx(expected, rel=1.0e-5) # test version for extreme L values - Lmax = np.log10((xb - xa) * 1.e-2) - L = np.power(10, np.random.uniform(Lmax - 2., Lmax)) + Lmax = np.log10((xb - xa) * 1.0e-2) + L = np.power(10, np.random.uniform(Lmax - 2.0, Lmax)) # theoretical value for Beer-Lambert profile beta1 = (xb - xa) / L expected = np.exp(-beta1 - alphas * xb) - np.exp(-alphas * xa) - expected *= -phg * alphas / D / (alphas + 1. / L) - expected += minority / L + expected *= -phg * alphas / D / (alphas + 1.0 / L) - result = get_J_sc_diffusion_green(xa, xb, lambda x: alphas * np.exp(-alphas * x), - D, L, minority, s, phg, side='bottom') + result = get_J_sc_diffusion_green( + xa, + xb, + lambda x: alphas * np.exp(-alphas * x), + D, + L, + s, + phg, + side="bottom", + ) - assert beta1 > 1.e2 - assert result == approx(expected, rel=1.e-5) + assert beta1 > 1.0e2 + assert result == approx(expected, rel=1.0e-5) def test_identify_layers_exceptions(): @@ -600,40 +675,49 @@ def test_identify_layers_exceptions(): Na = np.power(10, np.random.uniform(22, 25)) Nd = np.power(10, np.random.uniform(22, 25)) - Lp = np.power(10, np.random.uniform(-9, -6))# Diffusion length - Ln = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + Lp = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + Ln = np.power(10, np.random.uniform(-9, -6)) # Diffusion length - GaAs_n = material("GaAs")(Nd = Nd, hole_diffusion_length=Ln) - GaAs_p = material("GaAs")(Na = Na, electron_diffusion_length=Lp) + GaAs_n = material("GaAs")(Nd=Nd, hole_diffusion_length=Ln) + GaAs_p = material("GaAs")(Na=Na, electron_diffusion_length=Lp) GaAs_i = material("GaAs")() - Ge_n = material("Ge")(Nd = Nd, hole_diffusion_length=Ln) + Ge_n = material("Ge")(Nd=Nd, hole_diffusion_length=Ln) - n_width = np.random.uniform(500, 1000)*1e-9 - p_width = np.random.uniform(3000, 5000)*1e-9 + n_width = np.random.uniform(500, 1000) * 1e-9 + p_width = np.random.uniform(3000, 5000) * 1e-9 i_width = np.random.uniform(300, 500) * 1e-9 - test_junc = Junction([Layer(n_width, GaAs_n,role="emitter"), - Layer(p_width, GaAs_p, role="neither")]) + test_junc = Junction( + [Layer(n_width, GaAs_n, role="emitter"), Layer(p_width, GaAs_p, role="neither")] + ) with raises(RuntimeError): identify_layers(test_junc) - test_junc = Junction([Layer(n_width, GaAs_n,role="emitter"), - Layer(i_width, GaAs_i, role="intrinsic"), - Layer(p_width, GaAs_p, role="nothing")]) + test_junc = Junction( + [ + Layer(n_width, GaAs_n, role="emitter"), + Layer(i_width, GaAs_i, role="intrinsic"), + Layer(p_width, GaAs_p, role="nothing"), + ] + ) with raises(RuntimeError): identify_layers(test_junc) - test_junc = Junction([Layer(n_width, Ge_n,role="emitter"), - Layer(p_width, GaAs_p, role="base")]) + test_junc = Junction( + [Layer(n_width, Ge_n, role="emitter"), Layer(p_width, GaAs_p, role="base")] + ) with raises(AssertionError): identify_layers(test_junc) def test_process_junction_np(): - from solcore.analytic_solar_cells.depletion_approximation import identify_layers, identify_parameters + from solcore.analytic_solar_cells.depletion_approximation import ( + identify_layers, + identify_parameters, + ) from solcore import material from solcore.structure import Layer, Junction from solcore.state import State @@ -644,36 +728,58 @@ def test_process_junction_np(): options = State() options.T = np.random.uniform(250, 350) - Lp = np.power(10, np.random.uniform(-9, -6))# Diffusion length - Ln = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + Lp = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + Ln = np.power(10, np.random.uniform(-9, -6)) # Diffusion length GaAs_window = material("GaAs")() - GaAs_n = material("GaAs")(Nd = Nd, hole_diffusion_length=Ln) - GaAs_p = material("GaAs")(Na = Na, electron_diffusion_length=Lp) + GaAs_n = material("GaAs")(Nd=Nd, hole_diffusion_length=Ln) + GaAs_p = material("GaAs")(Na=Na, electron_diffusion_length=Lp) - n_width = np.random.uniform(500, 1000)*1e-9 - p_width = np.random.uniform(3000, 5000)*1e-9 - window_width = np.random.uniform(25, 200)*1e-9 + n_width = np.random.uniform(500, 1000) * 1e-9 + p_width = np.random.uniform(3000, 5000) * 1e-9 + window_width = np.random.uniform(25, 200) * 1e-9 - test_junc = Junction([Layer(window_width, GaAs_window, role="window"), - Layer(n_width, GaAs_n,role="emitter"), - Layer(p_width, GaAs_p, role="base")]) + test_junc = Junction( + [ + Layer(window_width, GaAs_window, role="window"), + Layer(n_width, GaAs_n, role="emitter"), + Layer(p_width, GaAs_p, role="base"), + ] + ) id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np = identify_layers(test_junc) - xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es = identify_parameters(test_junc, options.T, pRegion, nRegion, iRegion) + xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es = identify_parameters( + test_junc, options.T, pRegion, nRegion, iRegion + ) ni_expect = GaAs_n.ni assert [id_top, id_bottom] == approx([1, 2]) - assert pn_or_np == 'np' - assert [xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es] == approx([n_width, p_width, 0, 0, 0, Lp, Ln, GaAs_n.hole_mobility * kb * options.T / q, - GaAs_p.electron_mobility * kb * options.T / q, Nd, Na, ni_expect, GaAs_n.permittivity]) - - + assert pn_or_np == "np" + assert [xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es] == approx( + [ + n_width, + p_width, + 0, + 0, + 0, + Lp, + Ln, + GaAs_n.hole_mobility * kb * options.T / q, + GaAs_p.electron_mobility * kb * options.T / q, + Nd, + Na, + ni_expect, + GaAs_n.permittivity, + ] + ) def test_process_junction_pn(): - from solcore.analytic_solar_cells.depletion_approximation import identify_layers, identify_parameters + from solcore.analytic_solar_cells.depletion_approximation import ( + identify_layers, + identify_parameters, + ) from solcore import material from solcore.structure import Layer, Junction from solcore.state import State @@ -690,24 +796,46 @@ def test_process_junction_pn(): GaAs_n = material("GaAs")(Nd=Nd, hole_diffusion_length=Ln) GaAs_p = material("GaAs")(Na=Na, electron_diffusion_length=Lp) - p_width = np.random.uniform(500, 1000)*1e-9 - n_width = np.random.uniform(3000, 5000)*1e-9 + p_width = np.random.uniform(500, 1000) * 1e-9 + n_width = np.random.uniform(3000, 5000) * 1e-9 - test_junc = Junction([Layer(p_width, GaAs_p, role="emitter"), Layer(n_width, GaAs_n, role="base")]) + test_junc = Junction( + [Layer(p_width, GaAs_p, role="emitter"), Layer(n_width, GaAs_n, role="base")] + ) id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np = identify_layers(test_junc) - xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es = identify_parameters(test_junc, options.T, pRegion, nRegion, iRegion) + xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es = identify_parameters( + test_junc, options.T, pRegion, nRegion, iRegion + ) ni_expect = GaAs_n.ni assert [id_top, id_bottom] == approx([0, 1]) - assert pn_or_np == 'pn' - assert [xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es] == approx([n_width, p_width, 0, 0, 0, Lp, Ln, GaAs_n.hole_mobility * kb * options.T / q, - GaAs_p.electron_mobility * kb * options.T / q, Nd, Na, ni_expect, GaAs_n.permittivity]) + assert pn_or_np == "pn" + assert [xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es] == approx( + [ + n_width, + p_width, + 0, + 0, + 0, + Lp, + Ln, + GaAs_n.hole_mobility * kb * options.T / q, + GaAs_p.electron_mobility * kb * options.T / q, + Nd, + Na, + ni_expect, + GaAs_n.permittivity, + ] + ) def test_process_junction_nip(): - from solcore.analytic_solar_cells.depletion_approximation import identify_layers, identify_parameters + from solcore.analytic_solar_cells.depletion_approximation import ( + identify_layers, + identify_parameters, + ) from solcore import material from solcore.structure import Layer, Junction from solcore.state import State @@ -718,35 +846,58 @@ def test_process_junction_nip(): options = State() options.T = np.random.uniform(250, 350) - Lp = np.power(10, np.random.uniform(-9, -6)) # Diffusion length - Ln = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + Lp = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + Ln = np.power(10, np.random.uniform(-9, -6)) # Diffusion length - GaAs_n = material("GaAs")(Nd = Nd, hole_diffusion_length=Ln) - GaAs_p = material("GaAs")(Na = Na, electron_diffusion_length=Lp) + GaAs_n = material("GaAs")(Nd=Nd, hole_diffusion_length=Ln) + GaAs_p = material("GaAs")(Na=Na, electron_diffusion_length=Lp) GaAs_i = material("GaAs")() - n_width = np.random.uniform(500, 1000)*1e-9 - p_width = np.random.uniform(3000, 5000)*1e-9 - i_width = np.random.uniform(100, 300)*1e-9 + n_width = np.random.uniform(500, 1000) * 1e-9 + p_width = np.random.uniform(3000, 5000) * 1e-9 + i_width = np.random.uniform(100, 300) * 1e-9 - test_junc = Junction([Layer(n_width, GaAs_n,role="emitter"), - Layer(i_width, GaAs_i, role="intrinsic"), - Layer(p_width, GaAs_p, role="base")]) + test_junc = Junction( + [ + Layer(n_width, GaAs_n, role="emitter"), + Layer(i_width, GaAs_i, role="intrinsic"), + Layer(p_width, GaAs_p, role="base"), + ] + ) id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np = identify_layers(test_junc) - xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es = identify_parameters(test_junc, options.T, pRegion, nRegion, iRegion) + xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es = identify_parameters( + test_junc, options.T, pRegion, nRegion, iRegion + ) ni_expect = GaAs_n.ni assert [id_top, id_bottom] == approx([0, 2]) - assert pn_or_np == 'np' - assert [xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es] == approx([n_width, p_width, i_width, 0, 0, Lp, Ln, GaAs_n.hole_mobility * kb * options.T / q, - GaAs_p.electron_mobility * kb * options.T / q, Nd, Na, ni_expect, GaAs_n.permittivity]) - + assert pn_or_np == "np" + assert [xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es] == approx( + [ + n_width, + p_width, + i_width, + 0, + 0, + Lp, + Ln, + GaAs_n.hole_mobility * kb * options.T / q, + GaAs_p.electron_mobility * kb * options.T / q, + Nd, + Na, + ni_expect, + GaAs_n.permittivity, + ] + ) def test_process_junction_pin(): - from solcore.analytic_solar_cells.depletion_approximation import identify_layers, identify_parameters + from solcore.analytic_solar_cells.depletion_approximation import ( + identify_layers, + identify_parameters, + ) from solcore import material from solcore.structure import Layer, Junction from solcore.state import State @@ -758,7 +909,7 @@ def test_process_junction_pin(): options.T = np.random.uniform(250, 350) Lp = np.power(10, np.random.uniform(-9, -6)) # Diffusion length - Ln = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + Ln = np.power(10, np.random.uniform(-9, -6)) # Diffusion length GaAs_n = material("GaAs")(Nd=Nd, hole_diffusion_length=Ln) GaAs_p = material("GaAs")(Na=Na, electron_diffusion_length=Lp) @@ -766,25 +917,49 @@ def test_process_junction_pin(): p_width = np.random.uniform(500, 1000) n_width = np.random.uniform(3000, 5000) - i_width = np.random.uniform(100, 300)*1e-9 + i_width = np.random.uniform(100, 300) * 1e-9 - test_junc = Junction([Layer(p_width, GaAs_p, role="emitter"), - Layer(i_width, GaAs_i, role="intrinsic"), - Layer(n_width, GaAs_n, role="base")]) + test_junc = Junction( + [ + Layer(p_width, GaAs_p, role="emitter"), + Layer(i_width, GaAs_i, role="intrinsic"), + Layer(n_width, GaAs_n, role="base"), + ] + ) id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np = identify_layers(test_junc) - xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es = identify_parameters(test_junc, options.T, pRegion, nRegion, iRegion) + xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es = identify_parameters( + test_junc, options.T, pRegion, nRegion, iRegion + ) ni_expect = GaAs_n.ni assert [id_top, id_bottom] == approx([0, 2]) - assert pn_or_np == 'pn' - assert [xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es] == approx([n_width, p_width, i_width, 0, 0, Lp, Ln, GaAs_n.hole_mobility * kb * options.T / q, - GaAs_p.electron_mobility * kb * options.T / q, Nd, Na, ni_expect, GaAs_n.permittivity]) + assert pn_or_np == "pn" + assert [xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd_c, Na_c, ni, es] == approx( + [ + n_width, + p_width, + i_width, + 0, + 0, + Lp, + Ln, + GaAs_n.hole_mobility * kb * options.T / q, + GaAs_p.electron_mobility * kb * options.T / q, + Nd, + Na, + ni_expect, + GaAs_n.permittivity, + ] + ) def test_process_junction_set_in_junction(): - from solcore.analytic_solar_cells.depletion_approximation import identify_layers, identify_parameters + from solcore.analytic_solar_cells.depletion_approximation import ( + identify_layers, + identify_parameters, + ) from solcore import material from solcore.structure import Layer, Junction from solcore.state import State @@ -792,13 +967,13 @@ def test_process_junction_set_in_junction(): options = State() options.T = np.random.uniform(250, 350) - Lp = np.power(10, np.random.uniform(-9, -6)) # Diffusion length - Ln = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + Lp = np.power(10, np.random.uniform(-9, -6)) # Diffusion length + Ln = np.power(10, np.random.uniform(-9, -6)) # Diffusion length sn = np.power(10, np.random.uniform(0, 3)) sp = np.power(10, np.random.uniform(0, 3)) - se = np.random.uniform(1, 20)*vacuum_permittivity + se = np.random.uniform(1, 20) * vacuum_permittivity GaAs_n = material("GaAs")() GaAs_p = material("GaAs")() @@ -806,39 +981,67 @@ def test_process_junction_set_in_junction(): p_width = np.random.uniform(500, 1000) n_width = np.random.uniform(3000, 5000) - i_width = np.random.uniform(100, 300)*1e-9 + i_width = np.random.uniform(100, 300) * 1e-9 mun = np.power(10, np.random.uniform(-5, 0)) mup = np.power(10, np.random.uniform(-5, 0)) Vbi = np.random.uniform(0, 3) - test_junc = Junction([Layer(p_width, GaAs_p, role="emitter"), - Layer(i_width, GaAs_i, role="intrinsic"), - Layer(n_width, GaAs_n, role="base")], - sn = sn, sp = sp, permittivity = se, - ln = Ln, lp= Lp, - mup = mup, mun = mun, Vbi = Vbi) + test_junc = Junction( + [ + Layer(p_width, GaAs_p, role="emitter"), + Layer(i_width, GaAs_i, role="intrinsic"), + Layer(n_width, GaAs_n, role="base"), + ], + sn=sn, + sp=sp, + permittivity=se, + ln=Ln, + lp=Lp, + mup=mup, + mun=mun, + Vbi=Vbi, + ) id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np = identify_layers(test_junc) - xn, xp, xi, sn_c, sp_c, ln, lp, dn, dp, Nd_c, Na_c, ni, es = identify_parameters(test_junc, options.T, pRegion, nRegion, iRegion) + xn, xp, xi, sn_c, sp_c, ln, lp, dn, dp, Nd_c, Na_c, ni, es = identify_parameters( + test_junc, options.T, pRegion, nRegion, iRegion + ) ni_expect = GaAs_n.ni assert [id_top, id_bottom] == approx([0, 2]) - assert pn_or_np == 'pn' - assert [xn, xp, xi, sn_c, sp_c, ln, lp, dn, dp, Nd_c, Na_c, ni, es] == approx([n_width, p_width, i_width, sn, sp, Ln, Lp, mun * kb * options.T / q, - mup * kb * options.T / q, 1, 1, ni_expect, se]) + assert pn_or_np == "pn" + assert [xn, xp, xi, sn_c, sp_c, ln, lp, dn, dp, Nd_c, Na_c, ni, es] == approx( + [ + n_width, + p_width, + i_width, + sn, + sp, + Ln, + Lp, + mun * kb * options.T / q, + mup * kb * options.T / q, + 1, + 1, + ni_expect, + se, + ] + ) def test_get_depletion_widths(): - from solcore.analytic_solar_cells.depletion_approximation import get_depletion_widths + from solcore.analytic_solar_cells.depletion_approximation import ( + get_depletion_widths, + ) from solcore.structure import Junction - xi = np.power(10, np.random.uniform(-10, -6)) + xi = np.power(10, np.random.uniform(-10, -6)) Vbi = np.random.uniform(0, 3) - es = np.random.uniform(1, 20)*vacuum_permittivity + es = np.random.uniform(1, 20) * vacuum_permittivity Na = np.power(10, np.random.uniform(22, 25)) Nd = np.power(10, np.random.uniform(22, 25)) @@ -847,8 +1050,12 @@ def test_get_depletion_widths(): test_junc = Junction() - wn_e = (-xi + np.sqrt(xi ** 2 + 2. * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / (1 + Nd / Na) - wp_e = (-xi + np.sqrt(xi ** 2 + 2. * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / (1 + Na / Nd) + wn_e = (-xi + np.sqrt(xi**2 + 2.0 * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / ( + 1 + Nd / Na + ) + wp_e = (-xi + np.sqrt(xi**2 + 2.0 * es * (Vbi - V) / q * (1 / Na + 1 / Nd))) / ( + 1 + Na / Nd + ) wn_r, wp_r = get_depletion_widths(test_junc, es, Vbi, V, Na, Nd, xi) @@ -857,13 +1064,15 @@ def test_get_depletion_widths(): def test_get_depletion_widths_onesided(): - from solcore.analytic_solar_cells.depletion_approximation import get_depletion_widths + from solcore.analytic_solar_cells.depletion_approximation import ( + get_depletion_widths, + ) from solcore.structure import Junction - xi = np.power(10, np.random.uniform(-10, -6)) + xi = np.power(10, np.random.uniform(-10, -6)) Vbi = np.random.uniform(0, 3) - es = np.random.uniform(1, 20)*vacuum_permittivity + es = np.random.uniform(1, 20) * vacuum_permittivity Na = np.power(10, np.random.uniform(22, 25)) Nd = np.power(10, np.random.uniform(22, 25)) @@ -882,11 +1091,13 @@ def test_get_depletion_widths_onesided(): def test_get_depletion_widths_set_in_junction(): - from solcore.analytic_solar_cells.depletion_approximation import get_depletion_widths + from solcore.analytic_solar_cells.depletion_approximation import ( + get_depletion_widths, + ) from solcore.structure import Junction - wn = np.random.uniform(1,100) - wp = np.random.uniform(1,100) + wn = np.random.uniform(1, 100) + wp = np.random.uniform(1, 100) test_junc = Junction(wn=wn, wp=wp) wn_r, wp_r = get_depletion_widths(test_junc, 0, 0, 0, 0, 0, 0) @@ -896,15 +1107,26 @@ def test_get_depletion_widths_set_in_junction(): def test_dark_iv_depletion_pn(pn_junction): - from solcore.analytic_solar_cells.depletion_approximation import iv_depletion, get_depletion_widths, get_j_dark, get_Jsrh, identify_layers, identify_parameters + from solcore.analytic_solar_cells.depletion_approximation import ( + iv_depletion, + get_depletion_widths, + get_j_dark, + get_Jsrh, + identify_layers, + identify_parameters, + ) from scipy.interpolate import interp1d test_junc, options = pn_junction options.light_iv = False T = options.T - id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np = identify_layers(test_junc[0]) - xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd, Na, ni, es = identify_parameters(test_junc[0], T, pRegion, nRegion, iRegion) + id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np = identify_layers( + test_junc[0] + ) + xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd, Na, ni, es = identify_parameters( + test_junc[0], T, pRegion, nRegion, iRegion + ) niSquared = ni**2 @@ -928,12 +1150,14 @@ def test_dark_iv_depletion_pn(pn_junction): min_top, min_bot = niSquared / Na, niSquared / Nd JtopDark = get_j_dark(x_top, w_top, l_top, s_top, d_top, V, min_top, T) - JbotDark = get_j_dark(x_bottom, w_bottom, l_bottom, s_bottom, d_bottom, V, min_bot, T) + JbotDark = get_j_dark( + x_bottom, w_bottom, l_bottom, s_bottom, d_bottom, V, min_bot, T + ) JnDark, JpDark = JbotDark, JtopDark - lifetime_n = ln ** 2 / dn - lifetime_p = lp ** 2 / dp # Jenny p163 + lifetime_n = ln**2 / dn + lifetime_p = lp**2 / dp # Jenny p163 Jrec = get_Jsrh(ni, V, Vbi, lifetime_p, lifetime_n, w, kbT) @@ -941,17 +1165,33 @@ def test_dark_iv_depletion_pn(pn_junction): J_sc_bot = 0 J_sc_scr = 0 - current = Jrec + JnDark + JpDark + V / 1e14- J_sc_top - J_sc_bot - J_sc_scr - iv = interp1d(test_junc[0].voltage, current, kind='linear', bounds_error=False, assume_sorted=True, - fill_value=(current[0], current[-1]), copy=True) + current = Jrec + JnDark + JpDark + V / 1e14 - J_sc_top - J_sc_bot - J_sc_scr + iv = interp1d( + test_junc[0].voltage, + current, + kind="linear", + bounds_error=False, + assume_sorted=True, + fill_value=(current[0], current[-1]), + copy=True, + ) iv_depletion(test_junc[0], options) - assert test_junc[0].iv(options.internal_voltages) == approx(iv(options.internal_voltages), nan_ok=True) + assert test_junc[0].iv(options.internal_voltages) == approx( + iv(options.internal_voltages), nan_ok=True + ) def test_dark_iv_depletion_np(np_junction): - from solcore.analytic_solar_cells.depletion_approximation import iv_depletion, get_depletion_widths, get_j_dark, get_Jsrh, identify_layers, identify_parameters + from solcore.analytic_solar_cells.depletion_approximation import ( + iv_depletion, + get_depletion_widths, + get_j_dark, + get_Jsrh, + identify_layers, + identify_parameters, + ) from scipy.interpolate import interp1d test_junc, options = np_junction @@ -960,8 +1200,12 @@ def test_dark_iv_depletion_np(np_junction): test_junc[0].voltage = options.internal_voltages - id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np = identify_layers(test_junc[0]) - xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd, Na, ni, es = identify_parameters(test_junc[0], T, pRegion, nRegion, iRegion) + id_top, id_bottom, pRegion, nRegion, iRegion, pn_or_np = identify_layers( + test_junc[0] + ) + xn, xp, xi, sn, sp, ln, lp, dn, dp, Nd, Na, ni, es = identify_parameters( + test_junc[0], T, pRegion, nRegion, iRegion + ) niSquared = ni**2 @@ -983,12 +1227,14 @@ def test_dark_iv_depletion_np(np_junction): min_bot, min_top = niSquared / Na, niSquared / Nd JtopDark = get_j_dark(x_top, w_top, l_top, s_top, d_top, V, min_top, T) - JbotDark = get_j_dark(x_bottom, w_bottom, l_bottom, s_bottom, d_bottom, V, min_bot, T) + JbotDark = get_j_dark( + x_bottom, w_bottom, l_bottom, s_bottom, d_bottom, V, min_bot, T + ) JpDark, JnDark = JbotDark, JtopDark - lifetime_n = ln ** 2 / dn - lifetime_p = lp ** 2 / dp # Jenny p163 + lifetime_n = ln**2 / dn + lifetime_p = lp**2 / dp # Jenny p163 Jrec = get_Jsrh(ni, V, Vbi, lifetime_p, lifetime_n, w, kbT) @@ -996,17 +1242,25 @@ def test_dark_iv_depletion_np(np_junction): J_sc_bot = 0 J_sc_scr = 0 - current = Jrec + JnDark + JpDark + V / 1e14- J_sc_top - J_sc_bot - J_sc_scr - iv = interp1d(test_junc[0].voltage, current, kind='linear', bounds_error=False, assume_sorted=True, - fill_value=(current[0], current[-1]), copy=True) + current = Jrec + JnDark + JpDark + V / 1e14 - J_sc_top - J_sc_bot - J_sc_scr + iv = interp1d( + test_junc[0].voltage, + current, + kind="linear", + bounds_error=False, + assume_sorted=True, + fill_value=(current[0], current[-1]), + copy=True, + ) iv_depletion(test_junc[0], options) - assert test_junc[0].iv(options.internal_voltages) == approx(iv(options.internal_voltages), nan_ok=True) + assert test_junc[0].iv(options.internal_voltages) == approx( + iv(options.internal_voltages), nan_ok=True + ) def test_qe_depletion_np(np_junction): - from solcore.analytic_solar_cells import qe_depletion test_junc, options = np_junction @@ -1019,23 +1273,35 @@ def test_qe_depletion_np(np_junction): assert np.all(test_junc[0].eqe_emitter(wl) < 1) assert np.all(test_junc[0].eqe_base(wl) < 1) assert np.all(test_junc[0].eqe_scr(wl) < 1) - assert np.all(test_junc[0].eqe(wl)[test_junc[0].eqe(wl) > 1e-3] - <= test_junc.absorbed[test_junc[0].eqe(wl) > 1e-3]) - assert np.all(test_junc[0].eqe_emitter(wl) + test_junc[0].eqe_base(wl) + - test_junc[0].eqe_scr(wl) == approx(test_junc[0].eqe(wl))) + assert np.all( + test_junc[0].eqe(wl)[test_junc[0].eqe(wl) > 1e-3] + <= test_junc.absorbed[test_junc[0].eqe(wl) > 1e-3] + ) + assert np.all( + test_junc[0].eqe_emitter(wl) + + test_junc[0].eqe_base(wl) + + test_junc[0].eqe_scr(wl) + == approx(test_junc[0].eqe(wl)) + ) assert np.all(test_junc[0].iqe(wl) >= test_junc[0].eqe(wl)) - options['da_mode'] = 'green' + options["da_mode"] = "green" qe_depletion(test_junc[0], options) assert np.all(test_junc[0].eqe(wl) < 1) assert np.all(test_junc[0].eqe_emitter(wl) < 1) assert np.all(test_junc[0].eqe_base(wl) < 1) assert np.all(test_junc[0].eqe_scr(wl) < 1) - assert np.all(test_junc[0].eqe(wl)[test_junc[0].eqe(wl) > 1e-3] - <= test_junc.absorbed[test_junc[0].eqe(wl) > 1e-3]) - assert np.all(test_junc[0].eqe_emitter(wl) + test_junc[0].eqe_base(wl) + - test_junc[0].eqe_scr(wl) == approx(test_junc[0].eqe(wl))) + assert np.all( + test_junc[0].eqe(wl)[test_junc[0].eqe(wl) > 1e-3] + <= test_junc.absorbed[test_junc[0].eqe(wl) > 1e-3] + ) + assert np.all( + test_junc[0].eqe_emitter(wl) + + test_junc[0].eqe_base(wl) + + test_junc[0].eqe_scr(wl) + == approx(test_junc[0].eqe(wl)) + ) assert np.all(test_junc[0].iqe(wl) >= test_junc[0].eqe(wl)) @@ -1052,27 +1318,39 @@ def test_qe_depletion_pn(pn_junction): assert np.all(test_junc[0].eqe_emitter(wl) < 1) assert np.all(test_junc[0].eqe_base(wl) < 1) assert np.all(test_junc[0].eqe_scr(wl) < 1) - assert np.all(test_junc[0].eqe(wl)[test_junc[0].eqe(wl) > 1e-3] - <= test_junc.absorbed[test_junc[0].eqe(wl) > 1e-3]) - assert np.all(test_junc[0].eqe_emitter(wl) + test_junc[0].eqe_base(wl) + - test_junc[0].eqe_scr(wl) == approx(test_junc[0].eqe(wl))) + assert np.all( + test_junc[0].eqe(wl)[test_junc[0].eqe(wl) > 1e-3] + <= test_junc.absorbed[test_junc[0].eqe(wl) > 1e-3] + ) + assert np.all( + test_junc[0].eqe_emitter(wl) + + test_junc[0].eqe_base(wl) + + test_junc[0].eqe_scr(wl) + == approx(test_junc[0].eqe(wl)) + ) assert np.all(test_junc[0].iqe(wl) >= test_junc[0].eqe(wl)) - options['da_mode'] = 'green' + options["da_mode"] = "green" qe_depletion(test_junc[0], options) assert np.all(test_junc[0].eqe(wl) < 1) assert np.all(test_junc[0].eqe_emitter(wl) < 1) assert np.all(test_junc[0].eqe_base(wl) < 1) assert np.all(test_junc[0].eqe_scr(wl) < 1) - assert np.all(test_junc[0].eqe(wl)[test_junc[0].eqe(wl) > 1e-3] - <= test_junc.absorbed[test_junc[0].eqe(wl) > 1e-3]) - assert np.all(test_junc[0].eqe_emitter(wl) + test_junc[0].eqe_base(wl) + - test_junc[0].eqe_scr(wl) == approx(test_junc[0].eqe(wl))) + assert np.all( + test_junc[0].eqe(wl)[test_junc[0].eqe(wl) > 1e-3] + <= test_junc.absorbed[test_junc[0].eqe(wl) > 1e-3] + ) + assert np.all( + test_junc[0].eqe_emitter(wl) + + test_junc[0].eqe_base(wl) + + test_junc[0].eqe_scr(wl) + == approx(test_junc[0].eqe(wl)) + ) assert np.all(test_junc[0].iqe(wl) >= test_junc[0].eqe(wl)) -def test_iv_depletion_np(np_junction): +def test_iv_depletion_np(np_junction): from solcore.analytic_solar_cells import iv_depletion test_junc, options = np_junction @@ -1082,32 +1360,33 @@ def test_iv_depletion_np(np_junction): iv_depletion(test_junc[0], options) - wl_sp, ph = options.light_source.spectrum(output_units='photon_flux_per_m', x=wl) - Jph = q*np.trapz(test_junc.absorbed*ph, wl) + wl_sp, ph = options.light_source.spectrum(output_units="photon_flux_per_m", x=wl) + Jph = q * np.trapz(test_junc.absorbed * ph, wl) approx_Voc = V[np.argmin(abs(test_junc[0].iv(V)))] quadrant = (V > 0) * (V < approx_Voc) - power = abs(test_junc[0].iv(V[quadrant])*V[quadrant])[:-1] + power = abs(test_junc[0].iv(V[quadrant]) * V[quadrant])[:-1] assert abs(test_junc[0].iv(0)) <= Jph - assert approx_Voc < test_junc[0][1].material.band_gap/q + assert approx_Voc < test_junc[0][1].material.band_gap / q assert np.all(power < options.light_source.power_density) - options['da_mode'] = 'green' + options["da_mode"] = "green" iv_depletion(test_junc[0], options) approx_Voc = V[np.argmin(abs(test_junc[0].iv(V)))] quadrant = (V > 0) * (V < approx_Voc) - power = abs(test_junc[0].iv(V[quadrant])*V[quadrant])[:-1] + power = abs(test_junc[0].iv(V[quadrant]) * V[quadrant])[:-1] assert abs(test_junc[0].iv(0)) <= Jph - assert approx_Voc < test_junc[0][1].material.band_gap/q + assert approx_Voc < test_junc[0][1].material.band_gap / q assert np.all(power < options.light_source.power_density) + def test_iv_depletion_pn(pn_junction): from solcore.analytic_solar_cells import iv_depletion @@ -1118,29 +1397,28 @@ def test_iv_depletion_pn(pn_junction): iv_depletion(test_junc[0], options) - wl_sp, ph = options.light_source.spectrum(output_units='photon_flux_per_m', x=wl) - Jph = q*np.trapz(test_junc.absorbed*ph, wl) + wl_sp, ph = options.light_source.spectrum(output_units="photon_flux_per_m", x=wl) + Jph = q * np.trapz(test_junc.absorbed * ph, wl) approx_Voc = V[np.argmin(abs(test_junc[0].iv(V)))] quadrant = (V > 0) * (V < approx_Voc) - power = abs(test_junc[0].iv(V[quadrant])*V[quadrant])[:-1] + power = abs(test_junc[0].iv(V[quadrant]) * V[quadrant])[:-1] assert abs(test_junc[0].iv(0)) <= Jph - assert approx_Voc < test_junc[0][1].material.band_gap/q + assert approx_Voc < test_junc[0][1].material.band_gap / q assert np.all(power < options.light_source.power_density) - options['da_mode'] = 'green' + options["da_mode"] = "green" iv_depletion(test_junc[0], options) approx_Voc = V[np.argmin(abs(test_junc[0].iv(V)))] quadrant = (V > 0) * (V < approx_Voc) - power = abs(test_junc[0].iv(V[quadrant])*V[quadrant])[:-1] + power = abs(test_junc[0].iv(V[quadrant]) * V[quadrant])[:-1] assert abs(test_junc[0].iv(0)) <= Jph - assert approx_Voc < test_junc[0][1].material.band_gap/q + assert approx_Voc < test_junc[0][1].material.band_gap / q assert np.all(power < options.light_source.power_density) -