From b6012f727afdc119df7a068028924ffdf3201735 Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Tue, 19 Apr 2022 11:07:06 +1000 Subject: [PATCH 01/13] start setting up anisotropic materials --- examples/anisotropic_material.py | 29 +++++++++++++++++++++++++++++ examples/tmm_rcwa_structures.py | 2 +- 2 files changed, 30 insertions(+), 1 deletion(-) create mode 100644 examples/anisotropic_material.py diff --git a/examples/anisotropic_material.py b/examples/anisotropic_material.py new file mode 100644 index 0000000..4556be6 --- /dev/null +++ b/examples/anisotropic_material.py @@ -0,0 +1,29 @@ +import numpy as np +from solcore import si, material +from solcore.structure import Layer +from solcore.constants import q, h, c +from solcore.interpolate import interp1d +from solcore.solar_cell import SolarCell + +from rayflare.rigorous_coupled_wave_analysis import rcwa_structure +from rayflare.transfer_matrix_method import tmm_structure +from rayflare.options import default_options + +wavelengths = np.linspace(300, 1000, 6)*1e-9 + +options = default_options() + +options.wavelengths = wavelengths +options.orders = 1 + + +# [width of the layer in nm, wavelengths, n at these wavelengths, k at these wavelengths, geometry] + +Air = material('Air')() + +n_tensor = np.array([[[2, 0, 0], [0, 2, 0], [0, 0, 2]], [[1, 0, 0], [0, 1, 0], [0, 0, 1]]]) +k_tensor = np.array([[[0, 0, 0], [0, 0, 0], [0, 0, 0]], [[0, 0, 0], [0, 0, 0], [0, 0, 0]]]) + +test_mat = [100, np.array([400, 900]), n_tensor, k_tensor, []] + +rcwa_setup = rcwa_structure([test_mat], size=((100, 0), (0, 100)), options=options, incidence=Air, transmission=Air) \ No newline at end of file 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)) From 9bfd1ccfb4fb4f9d6200486cb58f6755ab7a83d4 Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Wed, 20 Apr 2022 13:18:05 +1000 Subject: [PATCH 02/13] anisotropic material seems to be working for basic rcwa_structure calculation --- examples/anisotropic_material.py | 73 ++++++++++- .../rigorous_coupled_wave_analysis/rcwa.py | 41 +++++- rayflare/utilities.py | 122 +++++++++++++++++- 3 files changed, 224 insertions(+), 12 deletions(-) diff --git a/examples/anisotropic_material.py b/examples/anisotropic_material.py index 4556be6..77f3b4d 100644 --- a/examples/anisotropic_material.py +++ b/examples/anisotropic_material.py @@ -9,21 +9,88 @@ from rayflare.transfer_matrix_method import tmm_structure from rayflare.options import default_options +import matplotlib.pyplot as plt + wavelengths = np.linspace(300, 1000, 6)*1e-9 options = default_options() options.wavelengths = wavelengths options.orders = 1 - +options.parallel = False # [width of the layer in nm, wavelengths, n at these wavelengths, k at these wavelengths, geometry] Air = material('Air')() +Si = material("Si")() n_tensor = np.array([[[2, 0, 0], [0, 2, 0], [0, 0, 2]], [[1, 0, 0], [0, 1, 0], [0, 0, 1]]]) -k_tensor = np.array([[[0, 0, 0], [0, 0, 0], [0, 0, 0]], [[0, 0, 0], [0, 0, 0], [0, 0, 0]]]) +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 = [200, np.array([400, 900]), n_tensor, k_tensor, []] + +rcwa_setup = rcwa_structure([test_mat, test_mat2], size=((100, 0), (0, 100)), options=options, incidence=Air, transmission=Si) + +RAT = rcwa_setup.calculate(options) + +plt.figure() +plt.plot(wavelengths*1e9, RAT["A_per_layer"]) +plt.plot(wavelengths*1e9, RAT["R"] + RAT["T"] + np.sum(RAT["A_per_layer"], 1)) +plt.show() -rcwa_setup = rcwa_structure([test_mat], size=((100, 0), (0, 100)), options=options, incidence=Air, transmission=Air) \ No newline at end of file +# import S4 +# +# S = S4.New(Lattice=((1,0),(0,1)), NumBasis=100) +# +# S.SetMaterial(Name = 'Silicon', Epsilon = ( +# (12+0.01j, 0, 0), +# (0, 12+0.01j, 0), +# (0, 0, 12+0.01j) +# )) +# +# +# S.SetMaterial(Name = 'Vacuum', Epsilon = ( +# (1, 0, 0), +# (0, 1, 0), +# (0, 0, 1) +# )) +# +# S.SetMaterial(Name = 'mat2', Epsilon = ( +# (10+0.02j, 0, 0), +# (0, 12+0.01j, 0), +# (0, 0, 13+0.03j) +# )) +# +# +# S.AddLayer('inc', 0, 'Vacuum') +# S.AddLayer(Name = 'slab', Thickness = 0.6, Material = 'Silicon') +# S.AddLayer(Name = 'slab2', Thickness = 0.6, Material = 'mat2') +# S.AddLayer('trn', 0, 'Vacuum') +# +# S.SetRegionCircle( +# Layer = 'slab', +# Material = 'Vacuum', +# Center = (0,0), +# Radius = 0.2 +# ) +# +# S.SetExcitationPlanewave( +# IncidenceAngles=( +# 10, # polar angle in [0,180) +# 30 # azimuthal angle in [0,360) +# ), +# sAmplitude = 0.707+0.707j, +# pAmplitude = 0.707-0.707j, +# Order = 0 +# ) +# +# S.SetFrequency(1.2) +# +# from time import time +# +# start = time() +# +# (forw,back) = S.GetAmplitudes(Layer = 'trn', zOffset = 0) +# +# print(time()-start) \ No newline at end of file diff --git a/rayflare/rigorous_coupled_wave_analysis/rcwa.py b/rayflare/rigorous_coupled_wave_analysis/rcwa.py index 1a0e89f..7130cd4 100644 --- a/rayflare/rigorous_coupled_wave_analysis/rcwa.py +++ b/rayflare/rigorous_coupled_wave_analysis/rcwa.py @@ -11,11 +11,12 @@ 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 solcore.absorption_calculator import OptiStack +from rayflare.utilities import OptiStack + from sparse import COO, save_npz, stack try: @@ -691,7 +692,26 @@ def update_oc(self, wavelengths): 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 + + print(layers_oc.shape) + + # 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? Does this slow down the calculation? Test. (very slightly slower but negligible) + widths = stack_OS.get_widths() self.widths = widths @@ -1155,15 +1175,24 @@ 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): + print('a', len(layers_oc.shape)) + + if len(layers_oc.shape) > 1: + layers_oc = [tuple([tuple(y) for y in array]) for array in layers_oc] + n_inc = np.sqrt(layers_oc[0][0][0]) + + else: + n_inc = 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 diff --git a/rayflare/utilities.py b/rayflare/utilities.py index 38a6982..46dcbb8 100644 --- a/rayflare/utilities.py +++ b/rayflare/utilities.py @@ -9,7 +9,123 @@ 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 + + print('n', nx, ny) + + # 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: + 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): + + 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 +224,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 From 7050b6cf83e02d5bd4cd0468559fa01b2bbd6ba1 Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Wed, 20 Apr 2022 14:26:02 +1000 Subject: [PATCH 03/13] anisotropic material working in rcwa_structure (not yet implemented in matrix calcs) --- examples/anisotropic_material.py | 64 ++++++++++++++++--- .../rigorous_coupled_wave_analysis/rcwa.py | 34 +++++----- rayflare/utilities.py | 2 - 3 files changed, 71 insertions(+), 29 deletions(-) diff --git a/examples/anisotropic_material.py b/examples/anisotropic_material.py index 77f3b4d..c8a7afb 100644 --- a/examples/anisotropic_material.py +++ b/examples/anisotropic_material.py @@ -11,34 +11,80 @@ import matplotlib.pyplot as plt -wavelengths = np.linspace(300, 1000, 6)*1e-9 +wavelengths = np.linspace(300, 1000, 250)*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.4 # [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")() -n_tensor = np.array([[[2, 0, 0], [0, 2, 0], [0, 0, 2]], [[1, 0, 0], [0, 1, 0], [0, 0, 1]]]) -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 = GaAs.n(wavelengths) +GaAs_k = GaAs.k(wavelengths) -test_mat = [100, np.array([400, 900]), n_tensor, k_tensor, []] -test_mat2 = [200, np.array([400, 900]), n_tensor, k_tensor, []] +Si_n = Si.n(wavelengths) +Si_k = Si.k(wavelengths) -rcwa_setup = rcwa_structure([test_mat, test_mat2], size=((100, 0), (0, 100)), options=options, incidence=Air, transmission=Si) +# 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]]]) -RAT = rcwa_setup.calculate(options) +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, []] +test_mat2 = [1000, wavelengths*1e9, Si_n_tensor, Si_k_tensor, []] + +rcwa_setup = rcwa_structure([Layer(100e-9, GaAs), Layer(1e-6, Si)], size=((100, 0), (0, 100)), options=options, incidence=Air, transmission=Ag) + +rcwa_setup_AS = rcwa_structure([test_mat, test_mat2], size=((100, 0), (0, 100)), 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(wavelengths*1e9, RAT["A_per_layer"]) -plt.plot(wavelengths*1e9, RAT["R"] + RAT["T"] + np.sum(RAT["A_per_layer"], 1)) +plt.plot(prof['profile'][[0,20,50,100,249],:].T) +plt.plot(prof_AS['profile'][[0,20,50,100,249],:].T, '--') plt.show() + # import S4 # # S = S4.New(Lattice=((1,0),(0,1)), NumBasis=100) diff --git a/rayflare/rigorous_coupled_wave_analysis/rcwa.py b/rayflare/rigorous_coupled_wave_analysis/rcwa.py index 7130cd4..70f4486 100644 --- a/rayflare/rigorous_coupled_wave_analysis/rcwa.py +++ b/rayflare/rigorous_coupled_wave_analysis/rcwa.py @@ -12,10 +12,7 @@ from joblib import Parallel, delayed from rayflare.angles import make_angle_vector, overall_bin -from rayflare.utilities import get_matrices_or_paths - -# from solcore.absorption_calculator import OptiStack -from rayflare.utilities import OptiStack +from rayflare.utilities import get_matrices_or_paths, OptiStack from sparse import COO, save_npz, stack @@ -50,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) @@ -650,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, @@ -707,10 +702,8 @@ def update_oc(self, wavelengths): 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 - print(layers_oc.shape) - - # 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? Does this slow down the calculation? Test. (very slightly slower but negligible) + # 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() @@ -772,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() @@ -1175,8 +1168,6 @@ 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): - print('a', len(layers_oc.shape)) - if len(layers_oc.shape) > 1: layers_oc = [tuple([tuple(y) for y in array]) for array in layers_oc] n_inc = np.sqrt(layers_oc[0][0][0]) @@ -1228,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.sqrt(layers_oc[0][0][0]) + + else: + n_inc = 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)) @@ -1241,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 @@ -1263,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: @@ -1277,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/utilities.py b/rayflare/utilities.py index 46dcbb8..5ed3905 100644 --- a/rayflare/utilities.py +++ b/rayflare/utilities.py @@ -46,8 +46,6 @@ def __call__(self, x_new: np.ndarray) -> np.ndarray: nx = x_new.ndim ny = y_new.ndim - print('n', nx, ny) - # 6. Fill any values that were out of bounds with fill_value. # and # 7. Rotate the values back to their proper place. From b7f69160423fc662a9e0d5e2d0a86350bd495da7 Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Fri, 22 Apr 2022 13:43:24 +1000 Subject: [PATCH 04/13] :white_check_mark: add tests for anisotropic materials in rcwa_structure calculatioins --- examples/anisotropic_material.py | 46 ++++--- examples/rcwa_examples.ipynb | 2 +- .../rigorous_coupled_wave_analysis/rcwa.py | 8 +- tests/test_rigorous_coupled_wave.py | 116 ++++++++++++++++++ 4 files changed, 147 insertions(+), 25 deletions(-) diff --git a/examples/anisotropic_material.py b/examples/anisotropic_material.py index c8a7afb..072fe4b 100644 --- a/examples/anisotropic_material.py +++ b/examples/anisotropic_material.py @@ -11,7 +11,7 @@ import matplotlib.pyplot as plt -wavelengths = np.linspace(300, 1000, 250)*1e-9 +wavelengths = np.linspace(400, 800, 20)*1e-9 options = default_options() @@ -21,7 +21,7 @@ options.A_per_order = True options.pol = "p" options.theta_in = 0 -options.phi_in = 0.4 +options.phi_in = 0 # [width of the layer in nm, wavelengths, n at these wavelengths, k at these wavelengths, geometry] @@ -36,8 +36,8 @@ 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]]]) +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]) @@ -45,43 +45,49 @@ 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, []] +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), Layer(1e-6, Si)], size=((100, 0), (0, 100)), options=options, incidence=Air, transmission=Ag) +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=((100, 0), (0, 100)), options=options, incidence=Air, transmission=Ag) +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) +# RAT_s = rcwa_setup.calculate(options) options.pol = "p" RAT_p_AS = rcwa_setup_AS.calculate(options) -RAT_p = rcwa_setup.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["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() +# 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 = 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,20,50,100,249],:].T, '--') +# 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/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/rayflare/rigorous_coupled_wave_analysis/rcwa.py b/rayflare/rigorous_coupled_wave_analysis/rcwa.py index 70f4486..6891165 100644 --- a/rayflare/rigorous_coupled_wave_analysis/rcwa.py +++ b/rayflare/rigorous_coupled_wave_analysis/rcwa.py @@ -1170,10 +1170,10 @@ def RCWA_structure_wl(wl, geom_list, layers_oc, shapes_oc, s_names, pol, theta, if len(layers_oc.shape) > 1: layers_oc = [tuple([tuple(y) for y in array]) for array in layers_oc] - n_inc = np.sqrt(layers_oc[0][0][0]) + n_inc = np.real(np.sqrt(layers_oc[0][0][0])) else: - n_inc = np.sqrt(layers_oc[0]) + n_inc = np.real(np.sqrt(layers_oc[0])) def vs_pol(s, p): S.SetExcitationPlanewave((theta, phi), s, p, 0) @@ -1221,10 +1221,10 @@ def RCWA_wl_prof(wl, rat_output_A, dist, geom_list, layers_oc, shapes_oc, s_name if len(layers_oc.shape) > 1: layers_oc = [tuple([tuple(y) for y in array]) for array in layers_oc] - n_inc = np.sqrt(layers_oc[0][0][0]) + n_inc = np.real(np.sqrt(layers_oc[0][0][0])) else: - n_inc = np.sqrt(layers_oc[0]) + 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)) diff --git a/tests/test_rigorous_coupled_wave.py b/tests/test_rigorous_coupled_wave.py index d8bd039..3272701 100644 --- a/tests/test_rigorous_coupled_wave.py +++ b/tests/test_rigorous_coupled_wave.py @@ -497,3 +497,119 @@ 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_anisotropic_basic(): + 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 + + 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=2e-5) + + + +@mark.skipif(sys.platform == "win32", reason="S4 (RCWA) only installed for tests under Linux and macOS") +def test_anisotropic_pol(): + import numpy as np + 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"]) + + From 9a81a41b7d956f55989a933029d2794315111a13 Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Tue, 26 Apr 2022 15:27:59 +1000 Subject: [PATCH 05/13] make Codacy happy, fix failing test --- examples/anisotropic_material.py | 94 +++++------------------------ rayflare/utilities.py | 2 +- tests/test_rigorous_coupled_wave.py | 4 +- tests/test_utilities.py | 18 +++++- 4 files changed, 34 insertions(+), 84 deletions(-) diff --git a/examples/anisotropic_material.py b/examples/anisotropic_material.py index 072fe4b..87cd1bf 100644 --- a/examples/anisotropic_material.py +++ b/examples/anisotropic_material.py @@ -1,12 +1,8 @@ import numpy as np -from solcore import si, material +from solcore import material from solcore.structure import Layer -from solcore.constants import q, h, c -from solcore.interpolate import interp1d -from solcore.solar_cell import SolarCell from rayflare.rigorous_coupled_wave_analysis import rcwa_structure -from rayflare.transfer_matrix_method import tmm_structure from rayflare.options import default_options import matplotlib.pyplot as plt @@ -35,9 +31,9 @@ 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]]]) +# +# 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]) @@ -49,9 +45,6 @@ 'center': (0, 0), 'angle': 0, 'halfwidths': (150, 150)}]] test_mat2 = [1000, wavelengths*1e9, Si_n_tensor, Si_k_tensor, []] -test_mat = [100, np.array([400, 900]), n_tensor, k_tensor, []] -test_mat2 = [1000, np.array([400, 900]), n_tensor, 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) @@ -60,11 +53,11 @@ options.pol = "s" RAT_s_AS = rcwa_setup_AS.calculate(options) -# RAT_s = rcwa_setup.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) +RAT_p = rcwa_setup.calculate(options) plt.figure() plt.plot(wavelengths*1e9, RAT_s_AS["A_per_layer"], '--') @@ -74,75 +67,18 @@ 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() +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 = 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['profile'][[0,20,50,100,249],:].T) plt.plot(prof_AS['profile'][[0,5, 10],:].T, '--') plt.show() - - -# import S4 -# -# S = S4.New(Lattice=((1,0),(0,1)), NumBasis=100) -# -# S.SetMaterial(Name = 'Silicon', Epsilon = ( -# (12+0.01j, 0, 0), -# (0, 12+0.01j, 0), -# (0, 0, 12+0.01j) -# )) -# -# -# S.SetMaterial(Name = 'Vacuum', Epsilon = ( -# (1, 0, 0), -# (0, 1, 0), -# (0, 0, 1) -# )) -# -# S.SetMaterial(Name = 'mat2', Epsilon = ( -# (10+0.02j, 0, 0), -# (0, 12+0.01j, 0), -# (0, 0, 13+0.03j) -# )) -# -# -# S.AddLayer('inc', 0, 'Vacuum') -# S.AddLayer(Name = 'slab', Thickness = 0.6, Material = 'Silicon') -# S.AddLayer(Name = 'slab2', Thickness = 0.6, Material = 'mat2') -# S.AddLayer('trn', 0, 'Vacuum') -# -# S.SetRegionCircle( -# Layer = 'slab', -# Material = 'Vacuum', -# Center = (0,0), -# Radius = 0.2 -# ) -# -# S.SetExcitationPlanewave( -# IncidenceAngles=( -# 10, # polar angle in [0,180) -# 30 # azimuthal angle in [0,360) -# ), -# sAmplitude = 0.707+0.707j, -# pAmplitude = 0.707-0.707j, -# Order = 0 -# ) -# -# S.SetFrequency(1.2) -# -# from time import time -# -# start = time() -# -# (forw,back) = S.GetAmplitudes(Layer = 'trn', zOffset = 0) -# -# print(time()-start) \ No newline at end of file diff --git a/rayflare/utilities.py b/rayflare/utilities.py index 5ed3905..035e7f8 100644 --- a/rayflare/utilities.py +++ b/rayflare/utilities.py @@ -76,7 +76,7 @@ def __call__(self, x_new: np.ndarray) -> np.ndarray: return y_new.transpose(axes) -class OptiStack(OptiStack_SC): +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 diff --git a/tests/test_rigorous_coupled_wave.py b/tests/test_rigorous_coupled_wave.py index 3272701..7e2f1f1 100644 --- a/tests/test_rigorous_coupled_wave.py +++ b/tests/test_rigorous_coupled_wave.py @@ -501,7 +501,6 @@ def test_matrix_generation(): @mark.skipif(sys.platform == "win32", reason="S4 (RCWA) only installed for tests under Linux and macOS") def test_anisotropic_basic(): - import numpy as np from solcore import material from solcore.structure import Layer @@ -564,13 +563,12 @@ def test_anisotropic_basic(): 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=2e-5) + 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(): - import numpy as np from solcore import material from rayflare.rigorous_coupled_wave_analysis import rcwa_structure from rayflare.options import default_options diff --git a/tests/test_utilities.py b/tests/test_utilities.py index 0e4555d..1731bc9 100644 --- a/tests/test_utilities.py +++ b/tests/test_utilities.py @@ -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 From 40253931687a75be02df1bb221d6e5a5beb65552 Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Tue, 26 Apr 2022 16:15:20 +1000 Subject: [PATCH 06/13] :white_check_mark: add test for vector polarization --- tests/test_rigorous_coupled_wave.py | 32 +++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/tests/test_rigorous_coupled_wave.py b/tests/test_rigorous_coupled_wave.py index 7e2f1f1..f22c0d2 100644 --- a/tests/test_rigorous_coupled_wave.py +++ b/tests/test_rigorous_coupled_wave.py @@ -499,6 +499,38 @@ def test_matrix_generation(): 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_s >= 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 cddb4ca96608071248d572e50501e3b39bcb7cd7 Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Tue, 26 Apr 2022 16:24:33 +1000 Subject: [PATCH 07/13] fix typo in test --- tests/test_rigorous_coupled_wave.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_rigorous_coupled_wave.py b/tests/test_rigorous_coupled_wave.py index f22c0d2..6f12072 100644 --- a/tests/test_rigorous_coupled_wave.py +++ b/tests/test_rigorous_coupled_wave.py @@ -528,7 +528,7 @@ def test_rcwa_profile_circ(): output_rcwa_c = strt_rcwa.calculate_profile(options)['profile'] # circular pol; not sure what to compare this against? - assert np.all(output_rcwa_s >= 0) + assert np.all(output_rcwa_c >= 0) @mark.skipif(sys.platform == "win32", reason="S4 (RCWA) only installed for tests under Linux and macOS") From 48e019119fdc0089fd58b2dbfe0743d94db14caf Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Wed, 27 Apr 2022 14:43:36 +1000 Subject: [PATCH 08/13] :white_check_mark: add test for constant n/k in RCWA. Remove functionality which doesn't work for shape constant n/k --- .../rigorous_coupled_wave_analysis/rcwa.py | 20 +++++++-------- tests/test_rigorous_coupled_wave.py | 25 +++++++++++++++++++ 2 files changed, 35 insertions(+), 10 deletions(-) diff --git a/rayflare/rigorous_coupled_wave_analysis/rcwa.py b/rayflare/rigorous_coupled_wave_analysis/rcwa.py index 6891165..1e44ff7 100644 --- a/rayflare/rigorous_coupled_wave_analysis/rcwa.py +++ b/rayflare/rigorous_coupled_wave_analysis/rcwa.py @@ -18,7 +18,7 @@ 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.') @@ -672,15 +672,15 @@ 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. diff --git a/tests/test_rigorous_coupled_wave.py b/tests/test_rigorous_coupled_wave.py index 6f12072..05f54ec 100644 --- a/tests/test_rigorous_coupled_wave.py +++ b/tests/test_rigorous_coupled_wave.py @@ -643,3 +643,28 @@ def test_anisotropic_pol(): 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}]] + + str = rcwa_structure([test_mat], ((500, 0), (0, 500)), opts, Air, Air) + str2 = rcwa_structure([test_mat2], ((500, 0), (0, 500)), opts, Air, Air) + + RAT = str.calculate(opts) + RAT2 = str2.calculate(opts) + + assert np.all(RAT2["A_per_layer"] > RAT["A_per_layer"]) \ No newline at end of file From 5cf221262ce6484336a7214781f9a62222099abe Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Wed, 27 Apr 2022 14:52:34 +1000 Subject: [PATCH 09/13] make Codacy happy --- tests/test_rigorous_coupled_wave.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_rigorous_coupled_wave.py b/tests/test_rigorous_coupled_wave.py index 05f54ec..754a96d 100644 --- a/tests/test_rigorous_coupled_wave.py +++ b/tests/test_rigorous_coupled_wave.py @@ -661,10 +661,10 @@ def test_list_mats(): test_mat2 = const_n_mat + [[{'type': 'rectangle', 'mat': GaAs, 'center': (0, 0), 'halfwidths': [300, 300], 'angle': 20}]] - str = rcwa_structure([test_mat], ((500, 0), (0, 500)), opts, Air, Air) - str2 = rcwa_structure([test_mat2], ((500, 0), (0, 500)), opts, Air, Air) + 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 = str.calculate(opts) - RAT2 = str2.calculate(opts) + 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 From fb2e1d8c29f7fc52bb26b9512373403adfdf5fd3 Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Wed, 27 Apr 2022 15:57:52 +1000 Subject: [PATCH 10/13] ignore plotting stuff in tests --- rayflare/rigorous_coupled_wave_analysis/rcwa.py | 6 +++--- rayflare/utilities.py | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/rayflare/rigorous_coupled_wave_analysis/rcwa.py b/rayflare/rigorous_coupled_wave_analysis/rcwa.py index 1e44ff7..e91c596 100644 --- a/rayflare/rigorous_coupled_wave_analysis/rcwa.py +++ b/rayflare/rigorous_coupled_wave_analysis/rcwa.py @@ -903,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]) @@ -992,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]) @@ -1121,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]) diff --git a/rayflare/utilities.py b/rayflare/utilities.py index 035e7f8..21a2962 100644 --- a/rayflare/utilities.py +++ b/rayflare/utilities.py @@ -69,7 +69,7 @@ def __call__(self, x_new: np.ndarray) -> np.ndarray: axes = list(range(ny - nx)) axes[self.axis:self.axis] = list(range(ny - nx, ny)) return y_new.transpose(axes) - else: + 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)) From f87451e724bf31d7de73b49165e49a1a0b3c9584 Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Wed, 1 Jun 2022 10:35:27 +1000 Subject: [PATCH 11/13] updates for using refractiveindex.info database more consistently --- examples/full_RT_Sipyramids.py | 6 ++++-- examples/profile_pass_matrix.py | 10 +++++++--- tests/test_utilities.py | 4 ++-- 3 files changed, 13 insertions(+), 7 deletions(-) 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/tests/test_utilities.py b/tests/test_utilities.py index 1731bc9..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 From 1c6fb8eb566aa5e0769822444e39f8fa338e67ca Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Fri, 30 Sep 2022 12:09:57 +0200 Subject: [PATCH 12/13] update docs --- docs/Installation/installation.rst | 7 ++- docs/Installation/python_install.rst | 54 ++++++++++++------- .../matrix_formalism/multiply_matrices.py | 3 +- rayflare/transfer_matrix_method/tmm.py | 4 ++ setup.py | 20 ------- 5 files changed, 45 insertions(+), 43 deletions(-) diff --git a/docs/Installation/installation.rst b/docs/Installation/installation.rst index b50316e..26fa5aa 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 diff --git a/docs/Installation/python_install.rst b/docs/Installation/python_install.rst index 4466c91..a2a46a9 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 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/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/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() From 90023b18cff5d4c86711a7f4791170b4cdeb7a23 Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Fri, 30 Sep 2022 12:15:29 +0200 Subject: [PATCH 13/13] :memo: Update documentation for more recent Python versions --- docs/Installation/installation.rst | 2 +- docs/Installation/python_install.rst | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/Installation/installation.rst b/docs/Installation/installation.rst index 26fa5aa..b4cf084 100644 --- a/docs/Installation/installation.rst +++ b/docs/Installation/installation.rst @@ -169,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 a2a46a9..cf0ce99 100644 --- a/docs/Installation/python_install.rst +++ b/docs/Installation/python_install.rst @@ -59,4 +59,4 @@ with: .. code-block:: console - brew install python + brew install python \ No newline at end of file