diff --git a/docs/Installation/installation.rst b/docs/Installation/installation.rst index b50316e..b4cf084 100644 --- a/docs/Installation/installation.rst +++ b/docs/Installation/installation.rst @@ -81,6 +81,9 @@ Downloading and installing S4 (start from the directory in which you want to dow Here we download my fork of S4 (modified to be compatible with Python3), move to the downloaded directory, and then we make the Python extension. If you activated the virtual environment before running this, S4 should automatically install into that virtual environment. +**Note for users of new Mac devices with Apple M1/ARM chips:** You can install packages using Homebrew as described above, +but instead of :literal:`make S4_pyext`, you should run :literal:`make S4_pyext --file="Makefile.m1"`. + .. _install: Installing RayFlare @@ -128,8 +131,8 @@ things like 'pip install' will point to the right place. As an alternative to the venv method above, PyCharm also has built-in support for making a new virtualenv environment based on one of your 'system interpreters' (these are the versions of Python installed -- if you followed the -instructions above and are on Ubuntu 18, these are probably Python2.7, Python3.6 and the newly-installed -Python3.8). You can create one of those environments by going to Project Interpreter in the same way, adding a new +instructions above and are on Ubuntu 20, these are probably Python2.7, Python3.8 and possibly the newly-installed +Python3.10). You can create one of those environments by going to Project Interpreter in the same way, adding a new environment and selecting 'Virtualenv environment' with the relevant base interpreter. Installing S4 on Windows @@ -166,4 +169,4 @@ Troubleshooting .. _official website: https://releases.ubuntu.com/18.04/ .. _here: https://itsfoss.com/install-linux-in-virtualbox/ .. _dual boot: https://linuxconfig.org/how-to-install-ubuntu-20-04-alongside-windows-10-dual-boot -.. _homebrew: https://brew.sh/ +.. _homebrew: https://brew.sh/ \ No newline at end of file diff --git a/docs/Installation/python_install.rst b/docs/Installation/python_install.rst index 4466c91..cf0ce99 100644 --- a/docs/Installation/python_install.rst +++ b/docs/Installation/python_install.rst @@ -1,6 +1,9 @@ -======================================== -Python versions and installing Python3.8 -======================================== +========================================== +Python versions and installing Python3.10 +========================================== + +Ubuntu +------- On Ubuntu, Python is an integral part of the operating system, so there will already be a version (likely more than one version, usually a version of Python2 and a version of Python3) installed @@ -10,39 +13,50 @@ working in a virtual environment. This also avoids having to use 'sudo' to insta There are various packages for creating virtual environments, all with confusingly similar names -- I use venv (not to be confused with virtualenv or pipenv...). -Currently, Ubuntu18 does not come with Python3.8 (the latest stable Python version), by default, but you can -install it pretty easily. From a terminal window (command line): +On Ubuntu20, Python3.8 is the default Python3. You can install the ability to make venvs with: .. code-block:: console - sudo apt install software-properties-common - sudo add-apt-repository ppa:deadsnakes/ppa - sudo apt install python3.8 python3.8-dev python3.8-venv - -This installs Python3.8, and also the ability to make virtual environments with python3.8-venv. -You can now use your new version of Python3.8 by just typing 'python3.8' on the command line. + sudo apt install python3-dev python3-venv -Now make the venv: +And make and activate the venv: .. code-block:: console - python3.8 -m venv mainenv + python3 -m venv mainenv source mainenv/bin/activate -This makes a virtual environment called 'mainenv' (first line) and then activates it (second line). If you close the terminal and open it again you would have to reactivate the virtual environment to make -sure you are using the virtual environment version of Python3.8. +This makes a virtual environment called 'mainenv' (first line) and then activates it (second line). If you close the +terminal and open it again you would have to reactivate the virtual environment to make +sure you are using the virtual environment version of Python3.10. -On Ubuntu20, Python3.8 is the default Python3. You can install the ability to make venvs with: +Ubuntu20 comes with Python 3.8 by default. Python 3.8 is sufficient for running RayFlare and necessary supporting packages, +but if you want to update to a more recent Python version (3.10), this is easy: .. code-block:: console - sudo apt install python3-dev python3-venv + sudo apt install software-properties-common + sudo add-apt-repository ppa:deadsnakes/ppa + sudo apt install python3.10 python3.10-dev python3.10-venv -And make and activate the venv: +This installs Python3.10, and also the ability to make virtual environments with python3.10-venv. +You can now use your new version of Python3.10 by just typing 'python3.10' on the command line. + +Now make the venv: .. code-block:: console - python3 -m venv mainenv + python3.10 -m venv mainenv source mainenv/bin/activate -On MacOS, Homebrew Python (default version Python3.8) also comes with the venv ability. \ No newline at end of file + +MacOS +------- + + +On MacOS, Homebrew Python (current default version Python3.10) also comes with the venv ability. This can be installed +with: + +.. code-block:: console + + brew install python \ No newline at end of file diff --git a/examples/anisotropic_material.py b/examples/anisotropic_material.py new file mode 100644 index 0000000..87cd1bf --- /dev/null +++ b/examples/anisotropic_material.py @@ -0,0 +1,84 @@ +import numpy as np +from solcore import material +from solcore.structure import Layer + +from rayflare.rigorous_coupled_wave_analysis import rcwa_structure +from rayflare.options import default_options + +import matplotlib.pyplot as plt + +wavelengths = np.linspace(400, 800, 20)*1e-9 + +options = default_options() + +options.wavelengths = wavelengths +options.orders = 1 +options.parallel = False +options.A_per_order = True +options.pol = "p" +options.theta_in = 0 +options.phi_in = 0 + +# [width of the layer in nm, wavelengths, n at these wavelengths, k at these wavelengths, geometry] + +Air = material('Air')() +GaAs = material("GaAs")() +Si = material("Si")() +Ag = material("Ag")() + +GaAs_n = GaAs.n(wavelengths) +GaAs_k = GaAs.k(wavelengths) + +Si_n = Si.n(wavelengths) +Si_k = Si.k(wavelengths) +# +# n_tensor = np.array([[[2, 0, 0], [0, 3, 0], [0, 0, 1]], [[1, 0, 0], [0, 2, 0], [0, 0, 2]]]) +# k_tensor = np.array([[[0.1, 0, 0], [0, 0.1, 0], [0, 0, 0.1]], [[0.2, 0, 0], [0, 0.2, 0], [0, 0, 0.2]]]) + +GaAs_n_tensor = np.array([np.diag([x,x,x]) for x in GaAs_n]) +GaAs_k_tensor = np.array([np.diag([x,x,x]) for x in GaAs_k]) + +Si_n_tensor = np.array([np.diag([x,x,x]) for x in Si_n]) +Si_k_tensor = np.array([np.diag([x,x,x]) for x in Si_k]) + +test_mat = [100, wavelengths*1e9, GaAs_n_tensor, GaAs_k_tensor, [{'type': 'rectangle', 'mat': Air, + 'center': (0, 0), 'angle': 0, 'halfwidths': (150, 150)}]] +test_mat2 = [1000, wavelengths*1e9, Si_n_tensor, Si_k_tensor, []] + +rcwa_setup = rcwa_structure([Layer(100e-9, GaAs, geometry=[{'type': 'rectangle', 'mat': Air, + 'center': (0, 0), 'angle': 0, 'halfwidths': (150, 150)}]), + Layer(1e-6, Si)], size=((400, 0), (0, 400)), options=options, incidence=Air, transmission=Ag) + +rcwa_setup_AS = rcwa_structure([test_mat, test_mat2], size=((400, 0), (0, 400)), options=options, incidence=Air, transmission=Ag) + +options.pol = "s" +RAT_s_AS = rcwa_setup_AS.calculate(options) +RAT_s = rcwa_setup.calculate(options) + +options.pol = "p" +RAT_p_AS = rcwa_setup_AS.calculate(options) +RAT_p = rcwa_setup.calculate(options) + +plt.figure() +plt.plot(wavelengths*1e9, RAT_s_AS["A_per_layer"], '--') +plt.plot(wavelengths*1e9, RAT_s_AS["R"] + RAT_s_AS["T"] + np.sum(RAT_s_AS["A_per_layer"], 1)) + +plt.plot(wavelengths*1e9, RAT_p_AS["A_per_layer"]) +plt.plot(wavelengths*1e9, RAT_p_AS["R"] + RAT_p_AS["T"] + np.sum(RAT_p_AS["A_per_layer"], 1)) +plt.show() + +plt.figure() +plt.plot(wavelengths*1e9, RAT_s["A_per_layer"]) +plt.plot(wavelengths*1e9, RAT_s["R"] + RAT_s["T"] + np.sum(RAT_s["A_per_layer"], 1)) + +plt.plot(wavelengths*1e9, RAT_p["A_per_layer"]) +plt.plot(wavelengths*1e9, RAT_p["R"] + RAT_p["T"] + np.sum(RAT_p["A_per_layer"], 1)) +plt.show() + +prof = rcwa_setup.calculate_profile(options) +prof_AS = rcwa_setup_AS.calculate_profile(options) + +plt.figure() +plt.plot(prof['profile'][[0,20,50,100,249],:].T) +plt.plot(prof_AS['profile'][[0,5, 10],:].T, '--') +plt.show() diff --git a/examples/full_RT_Sipyramids.py b/examples/full_RT_Sipyramids.py index 99a46cc..e2a8971 100644 --- a/examples/full_RT_Sipyramids.py +++ b/examples/full_RT_Sipyramids.py @@ -6,6 +6,7 @@ from solcore import material from solcore import si +from solcore.absorption_calculator.nk_db import search_db # imports for plotting import matplotlib.pyplot as plt @@ -17,7 +18,8 @@ # setting up Solcore materials Air = material('Air')() -Si = material('566', nk_db=True)() +pageid = search_db("Si/Green-1995")[0][0] +Si = material(str(pageid), nk_db=True)() # number of x and y points to scan across nxy = 25 @@ -31,7 +33,7 @@ options.ny = nxy options.n_rays = 2 * nxy ** 2 options.depth_spacing = si('1um') -options.parallel = True +options.parallel = False PVlighthouse = np.loadtxt('data/RAT_data_300um_2um_55.csv', delimiter=',', skiprows=1) diff --git a/examples/profile_pass_matrix.py b/examples/profile_pass_matrix.py index de2ede6..1f08a25 100644 --- a/examples/profile_pass_matrix.py +++ b/examples/profile_pass_matrix.py @@ -7,12 +7,12 @@ from rayflare.angles import make_angle_vector from rayflare.utilities import make_absorption_function - from solcore import material, si from solcore.solar_cell import SolarCell, Layer, Junction from solcore.solar_cell_solver import solar_cell_solver from solcore.light_source import LightSource from solcore.constants import q +from solcore.absorption_calculator import search_db # imports for plotting import matplotlib.pyplot as plt @@ -62,8 +62,12 @@ InAlP = material("AlInP")(Al=0.5) GaInP = material("GaInP")(In=0.5) Si = material("Si")() -MgF2 = material("203", nk_db=True)() -Ta2O5 = material("410", nk_db=True)() + +Ta2O5_pageid = str(search_db("Ta2O5/Rodriguez-de Marcos")[0][0]) +MgF2_pageid = str(search_db("MgF2/Rodriguez-de Marcos")[0][0]) + +Ta2O5 = material(Ta2O5_pageid, nk_db=True)() +MgF2 = material(MgF2_pageid, nk_db=True)() GaAs_1_th = 120e-9 GaAs_2_th = 1200e-9 diff --git a/examples/rcwa_examples.ipynb b/examples/rcwa_examples.ipynb index 0efad97..70db197 100644 --- a/examples/rcwa_examples.ipynb +++ b/examples/rcwa_examples.ipynb @@ -80,7 +80,7 @@ "metadata": {}, "outputs": [], "source": [ - "grating_circles = geometry=[{'type': 'circle', 'mat': Ag, 'center': (0, 0), 'radius': 115}]\n", + "grating_circles = [{'type': 'circle', 'mat': Ag, 'center': (0, 0), 'radius': 115}]\n", "\n", "grating_squares = [{'type': 'rectangle', 'mat': Ag, 'center': (0, 0),\n", " 'halfwidths': [115, 115], 'angle': 20}]\n", diff --git a/examples/tmm_rcwa_structures.py b/examples/tmm_rcwa_structures.py index e3f3322..6e05177 100644 --- a/examples/tmm_rcwa_structures.py +++ b/examples/tmm_rcwa_structures.py @@ -24,7 +24,7 @@ options = default_options() options.wavelengths = wavelengths -options.orders = 2 +options.orders = 1 size = ((100, 0), (0, 100)) diff --git a/rayflare/matrix_formalism/multiply_matrices.py b/rayflare/matrix_formalism/multiply_matrices.py index 7a4d4f4..041d74b 100644 --- a/rayflare/matrix_formalism/multiply_matrices.py +++ b/rayflare/matrix_formalism/multiply_matrices.py @@ -302,7 +302,6 @@ def matrix_multiplication(bulk_mats, bulk_thick, options, layer_names, calc_prof len_calcs = np.array([len(x) if x is not None else 0 for x in calc_prof_list]) - if np.any(len_calcs > 0) or options.bulk_profile: a = [[] for _ in range(n_interfaces)] @@ -420,6 +419,8 @@ def matrix_multiplication(bulk_mats, bulk_thick, options, layer_names, calc_prof else: + print("b") + a = [[] for _ in range(n_interfaces)] vr = [[] for _ in range(n_bulks)] vt = [[] for _ in range(n_bulks)] diff --git a/rayflare/rigorous_coupled_wave_analysis/rcwa.py b/rayflare/rigorous_coupled_wave_analysis/rcwa.py index 1a0e89f..e91c596 100644 --- a/rayflare/rigorous_coupled_wave_analysis/rcwa.py +++ b/rayflare/rigorous_coupled_wave_analysis/rcwa.py @@ -11,16 +11,14 @@ import matplotlib.pyplot as plt from joblib import Parallel, delayed -from solcore.absorption_calculator import OptiStack - from rayflare.angles import make_angle_vector, overall_bin -from rayflare.utilities import get_matrices_or_paths +from rayflare.utilities import get_matrices_or_paths, OptiStack from sparse import COO, save_npz, stack try: import S4 -except Exception as err: +except Exception as err: # pragma: no cover print('WARNING: The RCWA solver will not be available because an S4 installation has not been found.') @@ -49,7 +47,7 @@ def RCWA(structure, size, orders, options, structpath, incidence, transmission, :return: """ # TODO: when doing unpolarized, why not just set s=0.5 p=0.5 in S4? (Maybe needs to be normalised differently). Also don't know if this is faster, - # or if internally it will still do s & p separately + # or if internally it will still do s & p separately # TODO: if incidence angle is zero, s and p polarization are the same so no need to do both existing_mats, path_or_mats = get_matrices_or_paths(structpath, surf_name, front_or_rear, prof_layers) @@ -649,8 +647,6 @@ def __init__(self, structure, size, options, incidence, transmission): shapes_names = [str(x) for x in shape_mats] - # depth_spacing = options['depth_spacing'] - # RCWA options S4_options = dict(LatticeTruncation='Circular', DiscretizedEpsilon=False, @@ -676,22 +672,39 @@ def update_oc(self, wavelengths): shapes_oc = np.zeros((len(wavelengths), len(self.shape_mats)), dtype=complex) for i1, x in enumerate(self.shape_mats): - if isinstance(x, list): - if len(x) == 3: - shapes_oc[:, i1] = np.ones_like(wavelengths)*(x[1] + 1j*x[2])**2 - - if len(x) == 4: - shapes_oc[:, i1] = (x[2] + 1j*x[3])**2 - - else: - shapes_oc[:, i1] = (x.n(wavelengths) + 1j * x.k(wavelengths)) ** 2 + # if isinstance(x, list): # this doesn't work (because of necessary_materials) + # if len(x) == 3: + # shapes_oc[:, i1] = np.ones_like(wavelengths)*(x[1] + 1j*x[2])**2 + # + # if len(x) == 4: + # shapes_oc[:, i1] = (x[2] + 1j*x[3])**2 + # + # else: + shapes_oc[:, i1] = (x.n(wavelengths) + 1j * x.k(wavelengths)) ** 2 # prepare to pass to OptiStack. stack_OS = OptiStack(self.list_for_OS, no_back_reflection=False, substrate=self.transmission, incidence=self.incidence) - layers_oc = (np.array(stack_OS.get_indices(wavelengths*1e9))**2).T + layers_oc_list = stack_OS.get_indices(wavelengths*1e9) + is_anisotropic = np.array([x.size for x in layers_oc_list]) == 9*len(wavelengths) + + if np.any(is_anisotropic): + for i1 in np.where(~is_anisotropic)[0]: + + layers_oc_list[i1] = np.array([np.diag([x,x,x]) for x in layers_oc_list[i1]]) + + + layers_oc = np.transpose(layers_oc_list, [1,0,2,3])**2 + + else: + layers_oc = (np.array(layers_oc_list)**2).T + # at least one element in the structure has a dielectric tensor; convert all to this format + + # dielectric tensor: can't put in an array because not all same size if not all materials are anisotropic + # if ANY are anisotropic: expand all into dielectric tensor. + widths = stack_OS.get_widths() self.widths = widths @@ -752,7 +765,7 @@ def calculate(self, options): A_order = np.real(np.stack([item[3] for item in allres])) - S_for_orders = initialise_S(self.size, options['orders'], self.geom_list, self.layers_oc[0], + S_for_orders = initialise_S(self.size, options['orders'], self.geom_list, np.zeros(len(wl)), self.shapes_oc[0], self.shapes_names, self.widths, options['S4_options']) basis_set = S_for_orders.GetBasisSet() @@ -890,7 +903,7 @@ def get_fourier_epsilon(self, layer_index, wavelength, options, extent=None, n_p a_r = a_r.reshape((len(xs), len(ys))) a_i = a_i.reshape((len(xs), len(ys))) - if plot: + if plot: # pragma: no cover fig, axs = plt.subplots(1, 2, figsize=(7, 2.6)) im1 = axs[0].pcolormesh(xs, ys, a_r.T, cmap='magma') fig.colorbar(im1, ax=axs[0]) @@ -979,7 +992,7 @@ def vs_pol(s, p): E_mag = np.real(np.sqrt(np.sum(np.abs(E) ** 2, 2))) H_mag = np.real(np.sqrt(np.sum(np.abs(H) ** 2, 2))) - if plot: + if plot: # pragma: no cover fig, axs = plt.subplots(1, 2, figsize=(7, 2.6)) im1 = axs[0].pcolormesh(xs, ys, E_mag.T, cmap='magma', shading='auto') fig.colorbar(im1, ax=axs[0]) @@ -1108,7 +1121,7 @@ def vs_pol(s, p): E_mag = np.real(np.sqrt(np.sum(np.abs(E) ** 2, 2))) H_mag = np.real(np.sqrt(np.sum(np.abs(H) ** 2, 2))) - if plot: + if plot: # pragma: no cover fig, axs = plt.subplots(1, 2, figsize=(7, 2.6)) im1 = axs[0].pcolormesh(xs, ys, E_mag.T, cmap='magma') fig.colorbar(im1, ax=axs[0]) @@ -1155,15 +1168,22 @@ def edit_geom_list(self, layer_index, geom_index, geom_entry): def RCWA_structure_wl(wl, geom_list, layers_oc, shapes_oc, s_names, pol, theta, phi, widths, size, orders, A_per_order, S4_options): + if len(layers_oc.shape) > 1: + layers_oc = [tuple([tuple(y) for y in array]) for array in layers_oc] + n_inc = np.real(np.sqrt(layers_oc[0][0][0])) + + else: + n_inc = np.real(np.sqrt(layers_oc[0])) + def vs_pol(s, p): S.SetExcitationPlanewave((theta, phi), s, p, 0) S.SetFrequency(1 / wl) - out = rcwa_rat(S, len(widths), theta, np.sqrt(layers_oc[0])) + out = rcwa_rat(S, len(widths), theta, n_inc) R = out[0]['R'] T = out[0]['T'] - A_layer = rcwa_absorption_per_layer(S, len(widths), theta, np.sqrt(layers_oc[0])) + A_layer = rcwa_absorption_per_layer(S, len(widths), theta, n_inc) if A_per_order: - A_per_layer_order = rcwa_absorption_per_layer_order(S, len(widths), theta, np.sqrt(layers_oc[0])) + A_per_layer_order = rcwa_absorption_per_layer_order(S, len(widths), theta, n_inc) return R, T, A_layer, A_per_layer_order else: return R, T, A_layer @@ -1199,6 +1219,13 @@ def vs_pol(s, p): def RCWA_wl_prof(wl, rat_output_A, dist, geom_list, layers_oc, shapes_oc, s_names, pol, theta, phi, widths, size, orders, S4_options): + if len(layers_oc.shape) > 1: + layers_oc = [tuple([tuple(y) for y in array]) for array in layers_oc] + n_inc = np.real(np.sqrt(layers_oc[0][0][0])) + + else: + n_inc = np.real(np.sqrt(layers_oc[0])) + S = initialise_S(size, orders, geom_list, layers_oc, shapes_oc, s_names, widths, S4_options) profile_data = np.zeros(len(dist)) @@ -1212,7 +1239,7 @@ def RCWA_wl_prof(wl, rat_output_A, dist, geom_list, layers_oc, shapes_oc, s_name layer, d_in_layer = tmm.find_in_structure_with_inf(widths, d) # don't need to change this layer_name = 'layer_' + str(layer + 1) # layer_1 is air above so need to add 1 - data = rcwa_position_resolved(S, layer_name, d_in_layer, A, theta, np.sqrt(layers_oc[0])) + data = rcwa_position_resolved(S, layer_name, d_in_layer, A, theta, n_inc) profile_data[j] = data @@ -1234,7 +1261,7 @@ def RCWA_wl_prof(wl, rat_output_A, dist, geom_list, layers_oc, shapes_oc, s_name layer, d_in_layer = tmm.find_in_structure_with_inf(widths, d) # don't need to change this layer_name = 'layer_' + str(layer + 1) # layer_1 is air above so need to add 1 - data = rcwa_position_resolved(S, layer_name, d_in_layer, A, theta, np.sqrt(layers_oc[0])) + data = rcwa_position_resolved(S, layer_name, d_in_layer, A, theta, n_inc) profile_data[j] = data else: @@ -1248,9 +1275,9 @@ def RCWA_wl_prof(wl, rat_output_A, dist, geom_list, layers_oc, shapes_oc, s_name d) # don't need to change this layer_name = 'layer_' + str(layer + 1) # layer_1 is air above so need to add 1 S.SetExcitationPlanewave((theta, phi), 0, 1, 0) # p-polarization - data_p = rcwa_position_resolved(S, layer_name, d_in_layer, A, theta, np.sqrt(layers_oc[0])) + data_p = rcwa_position_resolved(S, layer_name, d_in_layer, A, theta, n_inc) S.SetExcitationPlanewave((theta, phi), 1, 0, 0) # p-polarization - data_s = rcwa_position_resolved(S, layer_name, d_in_layer, A, theta, np.sqrt(layers_oc[0])) + data_s = rcwa_position_resolved(S, layer_name, d_in_layer, A, theta, n_inc) profile_data[j] = np.real(0.5*(data_s + data_p)) return profile_data diff --git a/rayflare/transfer_matrix_method/tmm.py b/rayflare/transfer_matrix_method/tmm.py index 3f56173..728f291 100644 --- a/rayflare/transfer_matrix_method/tmm.py +++ b/rayflare/transfer_matrix_method/tmm.py @@ -477,6 +477,10 @@ def calculate_profile(layers): output['A'] = 1 - out['R'] - out['T'] output['T'] = out['T'] output['A_per_layer'] = A_per_layer[1:-1] + output['vw_list'] = out['vw_list'] + output['kz_list'] = out['kz_list'] + output['th_list'] = out['th_list'] + else: out = tmm.inc_tmm(pol, stack.get_indices(wavelength), stack.get_widths(), coherency_list, angle, wavelength) A_per_layer = np.array(tmm.inc_absorp_in_each_layer(out)) diff --git a/rayflare/utilities.py b/rayflare/utilities.py index 38a6982..21a2962 100644 --- a/rayflare/utilities.py +++ b/rayflare/utilities.py @@ -9,7 +9,121 @@ import os from sparse import load_npz import xarray as xr -from solcore.interpolate import interp1d +from solcore import si +from solcore.interpolate import interp1d as interp1d_SC +from solcore.absorption_calculator import OptiStack as OptiStack_SC +from solcore.absorption_calculator.transfer_matrix import np_cache + + +class interp1d(interp1d_SC): + def __call__(self, x_new: np.ndarray) -> np.ndarray: + """Find interpolated y_new = f(x_new). + + Parameters + ---------- + x_new : number or array + New independent variable(s). + + Returns + ------- + y_new : ndarray + Interpolated value(s) corresponding to x_new. + + """ + + # 1. Handle values in x_new that are outside of x. Throw error, + # or return a list of mask array indicating the outofbounds values. + # The behavior is set by the bounds_error variable. + x_new = np.asarray(x_new) + out_of_bounds = self._check_bounds(x_new) + + y_new = self._call(x_new) + + # Rotate the values of y_new back so that they correspond to the + # correct x_new values. For N-D x_new, take the last (for linear) + # or first (for other splines) N axes + # from y_new and insert them where self.axis was in the list of axes. + nx = x_new.ndim + ny = y_new.ndim + + # 6. Fill any values that were out of bounds with fill_value. + # and + # 7. Rotate the values back to their proper place. + + if nx == 0: + # special case: x is a scalar + if out_of_bounds: + if ny == 0: + return np.asarray(self.fill_value) + else: + y_new[...] = self.fill_value + return np.asarray(y_new) + elif self._kind in ('linear', 'nearest'): + if ny == 1: + y_new[..., out_of_bounds] = self.fill_value + axes = list(range(ny - nx)) + axes[self.axis:self.axis] = list(range(ny - nx, ny)) + return y_new.transpose(axes) + else: + y_new[..., out_of_bounds] = self.fill_value[..., None] + axes = list(range(ny - nx)) + axes[self.axis:self.axis] = list(range(ny - nx, ny)) + return y_new.transpose(axes) + else: # pragma: no cover + y_new[out_of_bounds] = self.fill_value + axes = list(range(nx, ny)) + axes[self.axis:self.axis] = list(range(nx)) + return y_new.transpose(axes) + + +class OptiStack(OptiStack_SC): # pragma: no cover + + 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. + + :param layer: The new layer to add as [thickness, wavelength, n, k] + :return: None + """ + # We make sure that the wavelengths are increasing, reversing the arrays otherwise. + + if len(layer[1]) > 1: + + if layer[1][0] > layer[1][-1]: + layer[1] = layer[1][::-1] + + layer[2] = layer[2][::-1] + layer[3] = layer[3][::-1] + + self.widths.append(layer[0]) + + if len(layer) >= 5: + self.models.append(layer[4]) + c = si(layer[5][0], 'nm') + w = 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=si(layer[1], 'nm'), y=layer[2], fill_value=layer[2][-1])(x)) + k_data = np_cache(lambda x: interp1d(x=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=si(layer[1], 'nm'), y=layer[2], fill_value=layer[2][-1], axis=0))) + self.k_data.append(np_cache(interp1d(x=si(layer[1], 'nm'), y=layer[3], fill_value=layer[3][-1], axis=0))) + + def get_matrices_or_paths(structpath, surf_name, front_or_rear, prof_layers=None): @@ -108,14 +222,14 @@ def make_absorption_function(profile_result, structure, options, matrix_method=F back_positions/1e9 + np.sum(widths_front)/1e9 + width_bulk/1e6)) all_absorption = np.hstack((prof_front, prof_bulk, prof_back)) - diff_absorb_fn = interp1d(all_positions, all_absorption) + diff_absorb_fn = interp1d_SC(all_positions, all_absorption) return all_positions, diff_absorb_fn else: positions = np.arange(0, structure.width, options.depth_spacing) - diff_absorb_fn = interp1d(positions, 1e9 * profile_result) + diff_absorb_fn = interp1d_SC(positions, 1e9 * profile_result) return positions, diff_absorb_fn diff --git a/setup.py b/setup.py index 8b61fcf..1ea1c70 100644 --- a/setup.py +++ b/setup.py @@ -24,26 +24,6 @@ def gen_data_files(*dirs): config = ConfigParser() config.read([default_config]) -# Option for updating the manifest -if "update_manifest" in sys.argv: - # Update the MANIFEST.in file with all the data from within solcore - include = "rayflare" - exclude = ["__pycache__", "egg", "darwin", "cpython"] - with open("MANIFEST.in", "w", encoding="utf-8") as f: - for root, dirs, files in os.walk("."): - if not any(sub in root for sub in exclude) and root[2:9] == include: - try: - files.remove(".DS_Store") - except ValueError: - pass - for file in files: - if any(sub in file for sub in exclude): - continue - include_line = "include " + os.path.join(root[2:], file) + "\n" - f.write(include_line) - - sys.exit() - # Get the long description from the README file with open(os.path.join(here, "README.md"), encoding="utf-8") as f: long_description = f.read() diff --git a/tests/test_rigorous_coupled_wave.py b/tests/test_rigorous_coupled_wave.py index d8bd039..754a96d 100644 --- a/tests/test_rigorous_coupled_wave.py +++ b/tests/test_rigorous_coupled_wave.py @@ -497,3 +497,174 @@ def test_matrix_generation(): assert full_mat.shape == (len(wavelengths), 6, options.n_theta_bins) assert A_mat.shape == (len(wavelengths), 4, options.n_theta_bins) + + +@mark.skipif(sys.platform == "win32", reason="S4 (RCWA) only installed for tests under Linux and macOS") +def test_rcwa_profile_circ(): + from rayflare.options import default_options + from solcore import material + from solcore import si + from rayflare.rigorous_coupled_wave_analysis import rcwa_structure + from solcore.structure import Layer + + Air = material('Air')() + Si = material('Si')() + GaAs = material('GaAs')() + Ge = material('Ge')() + + options = default_options() + + options.wavelengths = np.linspace(700, 1400, 4)*1e-9 + + options.depth_spacing = 10e-9 + options.theta_in = 0 + options.phi_in = 0 + options.pol = [0.707+0.707*1j, 0.707-0.707*1j] + + stack = [Layer(si('500nm'), GaAs), Layer(si('1.1um'), Si), Layer(si('0.834um'), Ge)] + + strt_rcwa = rcwa_structure(stack, ((100, 0), (0, 100)), options, Air, Air) + strt_rcwa.calculate(options) + output_rcwa_c = strt_rcwa.calculate_profile(options)['profile'] + # circular pol; not sure what to compare this against? + + assert np.all(output_rcwa_c >= 0) + + +@mark.skipif(sys.platform == "win32", reason="S4 (RCWA) only installed for tests under Linux and macOS") +def test_anisotropic_basic(): + from solcore import material + from solcore.structure import Layer + + from rayflare.rigorous_coupled_wave_analysis import rcwa_structure + from rayflare.options import default_options + + wavelengths = np.linspace(300, 1000, 25) * 1e-9 + + options = default_options() + + options.wavelengths = wavelengths + options.orders = 10 + options.A_per_order = True + options.theta_in = 1.3 + options.phi_in = 0.9 + + Air = material('Air')() + GaAs = material("GaAs")() + Si = material("Si")() + Ag = material("Ag")() + + GaAs_n = GaAs.n(wavelengths) + GaAs_k = GaAs.k(wavelengths) + + Si_n = Si.n(wavelengths) + Si_k = Si.k(wavelengths) + + GaAs_n_tensor = np.array([np.diag([x, x, x]) for x in GaAs_n]) + GaAs_k_tensor = np.array([np.diag([x, x, x]) for x in GaAs_k]) + + Si_n_tensor = np.array([np.diag([x, x, x]) for x in Si_n]) + Si_k_tensor = np.array([np.diag([x, x, x]) for x in Si_k]) + + test_mat = [100, wavelengths * 1e9, GaAs_n_tensor, GaAs_k_tensor, [{'type': 'rectangle', 'mat': Air, + 'center': (0, 0), 'angle': 0, + 'halfwidths': (150, 150)}]] + test_mat2 = [1000, wavelengths * 1e9, Si_n_tensor, Si_k_tensor, []] + + rcwa_setup = rcwa_structure([Layer(100e-9, GaAs, geometry=[{'type': 'rectangle', 'mat': Air, + 'center': (0, 0), 'angle': 0, + 'halfwidths': (150, 150)}]), + Layer(1e-6, Si)], size=((400, 0), (0, 400)), options=options, incidence=Air, + transmission=Ag) + + rcwa_setup_AS = rcwa_structure([test_mat, test_mat2], size=((400, 0), (0, 400)), options=options, incidence=Air, + transmission=Ag) + + options.pol = "s" + RAT_s_AS = rcwa_setup_AS.calculate(options) + RAT_s = rcwa_setup.calculate(options) + + options.pol = "p" + RAT_p_AS = rcwa_setup_AS.calculate(options) + RAT_p = rcwa_setup.calculate(options) + + prof = rcwa_setup.calculate_profile(options) + prof_AS = rcwa_setup_AS.calculate_profile(options) + + for key in RAT_s_AS.keys(): + assert RAT_s_AS[key] == approx(np.array(RAT_s[key])) + assert RAT_p_AS[key] == approx(np.array(RAT_p[key])) + + assert prof["profile"] == approx(prof_AS["profile"], abs=4e-5) + + + +@mark.skipif(sys.platform == "win32", reason="S4 (RCWA) only installed for tests under Linux and macOS") +def test_anisotropic_pol(): + from solcore import material + from rayflare.rigorous_coupled_wave_analysis import rcwa_structure + from rayflare.options import default_options + + wavelengths = np.linspace(400, 800, 20)*1e-9 + + options = default_options() + + options.wavelengths = wavelengths + options.orders = 1 + options.A_per_order = True + + Air = material("Air")() + + # [width of the layer in nm, wavelengths, n at these wavelengths, k at these wavelengths, geometry] + + n_tensor = np.array([[[2, 0, 0], [0, 3, 0], [0, 0, 1]], [[1, 0, 0], [0, 2, 0], [0, 0, 2]]]) + k_tensor = np.array([[[0.1, 0, 0], [0, 0.1, 0], [0, 0, 0.1]], [[0.2, 0, 0], [0, 0.2, 0], [0, 0, 0.2]]]) + + test_mat = [100, np.array([400, 900]), n_tensor, k_tensor, []] + test_mat2 = [1000, np.array([400, 900]), n_tensor, k_tensor, []] + + rcwa_setup_AS = rcwa_structure([test_mat, test_mat2], size=((400, 0), (0, 400)), options=options, incidence=Air, transmission=Air) + + options.pol = "s" + RAT_s_AS = rcwa_setup_AS.calculate(options) + + options.pol = "p" + RAT_p_AS = rcwa_setup_AS.calculate(options) + + options.pol = "u" + RAT_u_AS = rcwa_setup_AS.calculate(options) + + assert np.all(RAT_s_AS["A_per_layer"] > RAT_p_AS["A_per_layer"]) + assert np.all(RAT_s_AS["T"] > RAT_p_AS["T"]) + assert np.all(RAT_s_AS["R"] < RAT_p_AS["R"]) + + assert 0.5*(RAT_s_AS["A_per_layer"] + RAT_p_AS["A_per_layer"]) == approx(RAT_u_AS["A_per_layer"]) + assert 0.5*(RAT_s_AS["R"] + RAT_p_AS["R"]) == approx(RAT_u_AS["R"]) + assert 0.5*(RAT_s_AS["T"] + RAT_p_AS["T"]) == approx(RAT_u_AS["T"]) + + +@mark.skipif(sys.platform == "win32", reason="S4 (RCWA) only installed for tests under Linux and macOS") +def test_list_mats(): + from rayflare.rigorous_coupled_wave_analysis import rcwa_structure + from rayflare.options import default_options + from solcore import material + + opts = default_options() + opts.wavelengths = np.linspace(300, 400, 3)*1e-9 + + Air = material("Air")() + GaAs = material("GaAs")() + + const_n_mat = [100, 2, 0] + + test_mat = const_n_mat + [[]] + test_mat2 = const_n_mat + [[{'type': 'rectangle', 'mat': GaAs, + 'center': (0, 0), 'halfwidths': [300, 300], 'angle': 20}]] + + strt = rcwa_structure([test_mat], ((500, 0), (0, 500)), opts, Air, Air) + strt2 = rcwa_structure([test_mat2], ((500, 0), (0, 500)), opts, Air, Air) + + RAT = strt.calculate(opts) + RAT2 = strt2.calculate(opts) + + assert np.all(RAT2["A_per_layer"] > RAT["A_per_layer"]) \ No newline at end of file diff --git a/tests/test_utilities.py b/tests/test_utilities.py index 0e4555d..fc8c736 100644 --- a/tests/test_utilities.py +++ b/tests/test_utilities.py @@ -137,8 +137,8 @@ def test_matrix_method_profile(): InAlP = material("AlInP")(Al=0.5) GaInP = material("GaInP")(In=0.5) Si = material("Si")() - MgF2 = material("203", nk_db=True)() - Ta2O5 = material("410", nk_db=True)() + MgF2 = material("MgF2")() + Ta2O5 = material("TaOX2")() GaAs_1_th = 120e-9 GaAs_2_th = 1200e-9 @@ -193,4 +193,20 @@ def test_matrix_method_profile(): # Do one for a selection of layers, one for all layers profile calculated assert positions == approx(positions_expected) - assert generated == approx(generated_expected, abs=1e-4) \ No newline at end of file + assert generated == approx(generated_expected, abs=1e-4) + + +def test_interp1d(): + from rayflare.utilities import interp1d + + sc = interp1d(np.array([1,2]), np.array([2, 3]), fill_value=8) + + assert sc(1.5) == 2.5 + assert sc(2.5) == 8 + + assert np.all(sc(np.array([1.25, 1.75, 3])) == np.array([2.25, 2.75, 8])) + + sc = interp1d(np.array([1,2]), np.array([2, 3]), kind='zero', fill_value=8) + + assert sc(1.5) == 2 + assert sc(3) == 8