From 97edb30331f3ec33ba661c69f9bef92363c358ea Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Tue, 30 Jul 2024 11:14:28 +1000 Subject: [PATCH] Analytical and Lambertian methods working (top incidence): --- examples/pvk_Si_analytical_RT.py | 36 ++- examples/pvk_Si_analytical_RT_lambertian.py | 279 ++++++++++++++++++++ rayflare/ray_tracing/analytical_rt.py | 67 +++-- rayflare/ray_tracing/rt.py | 173 +++++------- 4 files changed, 414 insertions(+), 141 deletions(-) create mode 100644 examples/pvk_Si_analytical_RT_lambertian.py diff --git a/examples/pvk_Si_analytical_RT.py b/examples/pvk_Si_analytical_RT.py index 540bc54..5019e47 100644 --- a/examples/pvk_Si_analytical_RT.py +++ b/examples/pvk_Si_analytical_RT.py @@ -33,11 +33,12 @@ d = 100e-6 -wavelengths = np.linspace(300, 1150, 10) * 1e-9 +wavelengths = np.linspace(300, 1200, 50) * 1e-9 AM15G = LightSource(source_type='standard', version='AM1.5g', x=wavelengths, output_units='photon_flux_per_m') +lambert_approx = 30 options = default_options() options.analytical_ray_tracing = 0 @@ -50,17 +51,19 @@ options.randomize_surface = True options.I_thresh = 1e-3 options.depth_spacing_bulk = 1e-8 -options.n_rays = 500 +options.n_rays = n_rays +options.lambertian_approximation = lambert_approx # n_bounces = np.arange(1, 11, dtype=int) # n_bounces = [1, 2, 3, 5, 10, 15] # 10, 15, 20, 30, 70] -front_text = regular_pyramids(10, True, - interface_layers=[Layer(100e-9, coverglass_ARC)],) - -# front_text = planar_surface( +# front_text = regular_pyramids(10, True, # interface_layers=[Layer(100e-9, coverglass_ARC)],) +front_text = planar_surface( + # interface_layers=[Layer(100e-9, coverglass_ARC)], +) + front_text_2 = regular_pyramids(52, True, 1, interface_layers=[Layer(100e-9, MgF2), Layer(1000e-9, Pvk)] ) @@ -74,8 +77,9 @@ # options.n_rays = 10 # options.analytical_ray_tracing = 0 # result = rt_str.calculate(options) -# -# options.n_rays = n_rays + +options.n_rays = n_rays +options.lambertian_approximation = 0 # # start = time() result_1 = rt_str.calculate(options) @@ -83,7 +87,7 @@ A_layer = result_1['A_per_layer'] A_per_interface = result_1['A_per_interface'] -total = result_1['R'] + result_1['T'] + np.sum(A_layer, axis=1) + A_per_interface[1][:,1] +total_1 = result_1['R'] + result_1['T'] + np.sum(A_layer, axis=1) + A_per_interface[1][:,1] # # plt.figure() # plt.plot(wavelengths*1e9, result_1['R'], label='R') @@ -96,8 +100,9 @@ # plt.legend() # plt.show() -# run again without analytical RT + options.analytical_ray_tracing = 2 +options.lambertian_approximation = lambert_approx start = time() result = rt_str.calculate(options) @@ -132,11 +137,20 @@ plt.plot(wavelengths*1e9, A_layer[:,0], 'b--', label='glass') plt.plot(wavelengths*1e9, A_per_interface[1][:,1], 'y--') # plt.plot(wavelengths*1e9, A_per_interface[2][:,0], label='Ge') -# plt.plot(wavelengths*1e9, total, 'k--', label='total', alpha=0.5) +plt.plot(wavelengths*1e9, total_1, 'k--', label='total', alpha=0.5) # plt.axhline(1) plt.legend() plt.show() +# number of passes: +plt.figure() +plt.plot(wavelengths*1e9, np.mean(result_1['n_passes'],axis=1), 'k--') +plt.plot(wavelengths*1e9, np.mean(result['n_passes'], axis=1), 'k') + +plt.plot(wavelengths*1e9, np.max(result_1['n_passes'],axis=1), 'r--') +plt.plot(wavelengths*1e9, np.max(result['n_passes'], axis=1), 'r') +plt.show() + # result_1 is just RT, result is analytical RT fig, axes = plt.subplots(5, 2, figsize=(8, 15)) diff --git a/examples/pvk_Si_analytical_RT_lambertian.py b/examples/pvk_Si_analytical_RT_lambertian.py new file mode 100644 index 0000000..171878f --- /dev/null +++ b/examples/pvk_Si_analytical_RT_lambertian.py @@ -0,0 +1,279 @@ +from rayflare.ray_tracing.analytical_rt import lambertian_scattering + +from solcore.light_source import LightSource +from solcore.constants import q +import numpy as np +from rayflare.textures import regular_pyramids, planar_surface +import matplotlib.pyplot as plt + +from solcore import material +from solcore.structure import Layer +from rayflare.options import default_options +from rayflare.ray_tracing import rt_structure +from rayflare.ray_tracing.analytical_rt import analytical_start +from time import time +import seaborn as sns + +SiN = material("Si3N4")() +Si = material("Si")() +Air = material("Air")() +Ag = material("Ag")() +MgF2 = material("MgF2")() +Ge = material("Ge")() + +coverglass_ARC = material("coverglass_AR_JJ")() +coverglass = material("coverglass_JJ")() + +GaAs = material("GaAs")() + +Pvk = material("Pvk_Ox_165")() +epoxy = material("BK7")() + +n_rays = 800 + +d = 100e-6 + +wavelengths = np.linspace(300, 1180, 50) * 1e-9 + +AM15G = LightSource(source_type='standard', version='AM1.5g', x=wavelengths, + output_units='photon_flux_per_m') + +lambert_approx = 30 +options = default_options() + +options.analytical_ray_tracing = 0 +options.wavelength = wavelengths +options.project_name = 'integration_testing' +options.nx = 10 +options.ny = 10 +options.theta_in = 0.1 +options.parallel = True +# options.n_jobs = -3 +options.randomize_surface = True +options.I_thresh = 5e-3 +options.depth_spacing_bulk = 1e-8 +options.n_rays = n_rays +options.lambertian_approximation = lambert_approx + +# n_bounces = np.arange(1, 11, dtype=int) +# n_bounces = [1, 2, 3, 5, 10, 15] # 10, 15, 20, 30, 70] + +# front_text = regular_pyramids(10, True, +# interface_layers=[Layer(100e-9, coverglass_ARC)],) + +front_text = planar_surface( + interface_layers=[Layer(100e-9, coverglass_ARC)], +) + +front_text_2 = regular_pyramids(52, True, 1, + interface_layers=[Layer(100e-9, MgF2), Layer(400e-9, Pvk)] + ) +rear_text = regular_pyramids(52, False, 1) + +rt_str = rt_structure(textures=[front_text, front_text_2, rear_text], materials=[coverglass, Si], + widths=[10e-6, d], incidence=Air, transmission=Ag, + options=options, use_TMM=True, save_location='current', + overwrite=True) + +# options.n_rays = 10 +# options.analytical_ray_tracing = 0 +# result = rt_str.calculate(options) + +options.n_rays = n_rays +options.lambertian_approximation = 0 +# # +start = time() +result_1 = rt_str.calculate_profile(options) +normal_time = time() - start +print('Elapsed time: ', normal_time) + +A_layer = result_1['A_per_layer'] +A_per_interface = result_1['A_per_interface'] +total_1 = result_1['R'] + result_1['T'] + np.sum(A_layer, axis=1) + A_per_interface[1][:,1] +# +# plt.figure() +# plt.plot(wavelengths*1e9, result_1['R'], label='R') +# plt.plot(wavelengths*1e9, result_1['T'], label='T') +# plt.plot(wavelengths*1e9, A_layer, label='A') +# plt.plot(wavelengths*1e9, A_per_interface[1][:,1], label='GaAs') +# # plt.plot(wavelengths*1e9, A_per_interface[2][:,0], label='Ge_back') +# plt.plot(wavelengths*1e9, total, 'k-', label='total') +# plt.axhline(1) +# plt.legend() +# plt.show() + +n_bounces = [3, 5, 10, 20, 30, 40, 60, 80, 100, 120, 140, 200, 0] +options.analytical_ray_tracing = 2 + + +timing = np.zeros(len(n_bounces)) +result_list = [] + +for i1, n_b in enumerate(n_bounces): + + options.lambertian_approximation = n_b + start = time() + result = rt_str.calculate_profile(options) + result_list.append(result) + timing[i1] = time() - start + print(f'{n_b} passes elapsed time: ', timing[i1]) + +n_bounces_plot = n_bounces[:] +n_bounces_plot[-1] = max(n_bounces) + 10 + +plt.figure() +alphas = np.linspace(0.2, 1, len(n_bounces)) + +for i1, result in enumerate(result_list): + A_layer_2 = result['A_per_layer'] + A_per_interface_2 = result['A_per_interface'] + # import matplotlib.pyplot as plt + plt.plot(wavelengths*1e9, result['R'], 'k-', alpha=alphas[i1]) + plt.plot(wavelengths*1e9, A_layer_2[:,0], 'r-', alpha=alphas[i1]) + plt.plot(wavelengths*1e9, A_layer_2[:,1], 'g-', alpha=alphas[i1]) + plt.plot(wavelengths*1e9, A_per_interface_2[1][:,1], 'y-', alpha=alphas[i1]) + +plt.show() + +for i1, result in enumerate(result_list[:-1]): + ratio = result['A_per_layer'][:,1]/result_list[-1]['A_per_layer'][:,1] + plt.plot(wavelengths*1e9, ratio, 'k-', alpha=alphas[i1], label=f'{n_bounces[i1]} passes') + +plt.legend() +plt.show() + +max_diff = np.zeros(len(result_list)-1) +for i1, result in enumerate(result_list[:-1]): + difference = result['A_per_layer'][:,1] - result_list[-1]['A_per_layer'][:,1] + plt.plot(wavelengths*1e9, difference, 'k-', alpha=alphas[i1], label=f'{n_bounces[i1]} passes') + max_diff[i1] = np.max(np.abs(difference)) +plt.legend() +plt.show() + +plt.figure() +plt.plot(n_bounces_plot, timing, 'o', mfc='none') +plt.plot(max(n_bounces) + 20, normal_time, 'o', mfc='none') +plt.xlabel("Allowed passes") +plt.ylabel("Time (s)") +ax2 = plt.gca().twinx() +ax2.plot(n_bounces_plot[:-1], max_diff, 'o', mfc='none', color='r') +ax2.set_ylabel("Max difference in A") +plt.show() + +result_i = 3 + +int_A = 1e9*np.trapz(result_1['profile'], dx=options.depth_spacing_bulk) +int_A_lamb = 1e9*np.trapz(result_list[result_i]['profile'], dx=options.depth_spacing_bulk) + +total_2 = (result_list[result_i]['R'] + result_list[result_i]['T'] + + np.sum(result_list[result_i]['A_per_layer'], axis=1) + + result_list[result_i]['A_per_interface'][1][:,1]) + +plt.figure() +plt.plot(wavelengths*1e9, int_A, label='full RT') +plt.plot(wavelengths*1e9, result_1['A_per_layer'][:,1]) +plt.plot(wavelengths*1e9, total_1) +plt.plot(wavelengths*1e9, result_1['R']) + +plt.plot(wavelengths*1e9, int_A_lamb, '-r', label='lambertian') +plt.plot(wavelengths*1e9, result_list[result_i]['A_per_layer'][:,1], '--k') +plt.plot(wavelengths*1e9, total_2, '--') +plt.plot(wavelengths*1e9, result_list[result_i]['R'], '--') +plt.legend() +plt.show() + +cols = sns.color_palette('viridis', len(result_list)) + +# check if integrated A makes sense (bulk) +for i1, result in enumerate(result_list): + int_A = 1e9*np.trapz(result['profile'], dx=options.depth_spacing_bulk) + plt.plot(wavelengths*1e9, int_A, label=f'{n_bounces[i1]} passes', color=cols[i1]) + plt.plot(wavelengths*1e9, np.sum(result['A_per_layer'], axis=1), '--', color=cols[i1]) + +plt.xlim(500, 1180) +plt.legend() +plt.show() + + + + + +plt.figure() +plt.semilogy(result_1['profile'][-2]) +plt.semilogy(result_list[0]['profile'][-2]) +plt.show() + +# +# A_layer_2 = result['A_per_layer'] +# A_per_interface_2 = result['A_per_interface'] +# total = result['R'] + result['T'] + np.sum(A_layer_2, axis=1) + A_per_interface_2[1][:,1] +# # import matplotlib.pyplot as plt +# plt.figure() +# plt.plot(wavelengths*1e9, result['R'], 'k-', label='R') +# plt.plot(wavelengths*1e9, result['T'], 'r-', label='T') +# plt.plot(wavelengths*1e9, A_layer_2, 'g-', label='A') +# plt.plot(wavelengths*1e9, A_per_interface_2[1][:,1], 'y-', label='GaAs') +# plt.plot(wavelengths*1e9, total, 'k-', label='total', alpha=0.5) +# plt.axhline(1) +# plt.legend() +# plt.show() +# +# plt.figure() +# plt.plot(wavelengths*1e9, result['R'], 'k-', label='R') +# plt.plot(wavelengths*1e9, result['T'], 'r-', label='T') +# plt.plot(wavelengths*1e9, A_layer_2[:,0], 'b-', label='glass') +# plt.plot(wavelengths*1e9, A_per_interface_2[1][:,1], 'y-', label='Pvk') +# plt.plot(wavelengths*1e9, A_layer_2[:,1], 'g-', label='Si') +# +# plt.plot(wavelengths*1e9, total, 'k-', label='total', alpha=0.5) +# +# plt.plot(wavelengths*1e9, result_1['R'], 'k--') +# plt.plot(wavelengths*1e9, result_1['T'], 'r--') +# plt.plot(wavelengths*1e9, A_layer[:,1], 'g--') +# plt.plot(wavelengths*1e9, A_layer[:,0], 'b--', label='glass') +# plt.plot(wavelengths*1e9, A_per_interface[1][:,1], 'y--') +# # plt.plot(wavelengths*1e9, A_per_interface[2][:,0], label='Ge') +# plt.plot(wavelengths*1e9, total_1, 'k--', label='total', alpha=0.5) +# # plt.axhline(1) +# plt.legend() +# plt.show() +# +# # number of passes: +# plt.figure() +# plt.plot(wavelengths*1e9, np.mean(result_1['n_passes'],axis=1), 'k--') +# plt.plot(wavelengths*1e9, np.mean(result['n_passes'], axis=1), 'k') +# +# plt.plot(wavelengths*1e9, np.max(result_1['n_passes'],axis=1), 'r--') +# plt.plot(wavelengths*1e9, np.max(result['n_passes'], axis=1), 'r') +# plt.show() +# +# # result_1 is just RT, result is analytical RT +# +# fig, axes = plt.subplots(5, 2, figsize=(8, 15)) +# for i, ax in enumerate(axes.flatten()): +# ax.hist(result['n_interactions'][i], color='k', bins=51, alpha=0.5, label='Analytical RT', histtype='step', range=[0,50]) +# ax.hist(result_1['n_interactions'][i], color='r', bins=51, alpha=0.5, label='RT', histtype='step', range=[0, 50]) +# plt.title('n_interactions') +# plt.show() +# +# fig, axes = plt.subplots(5, 2, figsize=(8, 15)) +# for i, ax in enumerate(axes.flatten()): +# ax.hist(result['n_passes'][i], color='k', bins=11, alpha=0.5, label='Analytical RT', histtype='step', range=[0,10]) +# ax.hist(result_1['n_passes'][i], color='r', bins=11, alpha=0.5, label='RT', histtype='step', range=[0,10]) +# plt.title('n_passes') +# plt.show() +# +# +# fig, axes = plt.subplots(5, 2, figsize=(8, 15)) +# for i, ax in enumerate(axes.flatten()): +# ax.hist(result['thetas'][i], color='k', bins=30, alpha=0.5, label='Analytical RT', histtype='step', range=[0, np.pi]) +# ax.hist(result_1['thetas'][i], color='r', bins=30, alpha=0.5, label='RT', histtype='step', range=[0, np.pi]) +# plt.title('thetas') +# plt.show() +# +# fig, axes = plt.subplots(5, 2, figsize=(8, 15)) +# for i, ax in enumerate(axes.flatten()): +# ax.hist(result['phis'][i], color='k', bins=30, alpha=0.5, label='Analytical RT', histtype='step', range=[0, 2*np.pi]) +# ax.hist(result_1['phis'][i], color='r', bins=30, alpha=0.5, label='RT', histtype='step', range=[0, 2*np.pi]) +# plt.title('phis') +# plt.show() \ No newline at end of file diff --git a/rayflare/ray_tracing/analytical_rt.py b/rayflare/ray_tracing/analytical_rt.py index 096a180..183ae7d 100644 --- a/rayflare/ray_tracing/analytical_rt.py +++ b/rayflare/ray_tracing/analytical_rt.py @@ -7,7 +7,7 @@ from copy import deepcopy theta_lamb = np.linspace(0, 0.999 * np.pi / 2, 100) -def traverse_vectorised(width, theta, alpha, I_i, positions, I_thresh, direction): +def traverse_vectorised(width, theta, alpha, I_i, positions, direction): ratio = alpha[None, :] / np.real(np.abs(np.cos(theta))) DA_u = I_i[:, :, None] * ratio[:, :, None] * np.exp((-ratio[:, :, None] * positions[None, None, :])) @@ -220,7 +220,7 @@ def analytical_start(nks, I_rem_data = I_remaining.data[0] - if np.all(np.abs(normals[:,2]) > 0.99): + if np.all(np.abs(normals[:,2]) > 0.999): # if the surface is planar, can just use TMM or Fresnel equations directly # should already have a lookuptable, if necessary @@ -259,7 +259,7 @@ def analytical_start(nks, R_total = xr.DataArray(I_remaining*R[None, :], dims=["unique_direction", "wl"]) # INTERFACE (not bulk!) absorption - A_per_interface[surf_index] = xr.DataArray((I_rem_data[:,None]*A_per_int_layer)[None, :, :], dims=["unique_direction", "wl", "layer"]) + A_per_interface[surf_index] = xr.DataArray((I_rem_data[None, :]*A_per_int_layer.T)[None, :, :], dims=["unique_direction", "layer", "wl"]) else: @@ -271,7 +271,7 @@ def analytical_start(nks, T = 1 - R # no interface absorption: - A_per_interface[surf_index] = xr.DataArray(np.zeros((1, n_wl, 1)), dims=["unique_direction", "wl", "layer"]) + A_per_interface[surf_index] = xr.DataArray(np.zeros((1, 1, n_wl)), dims=["unique_direction", "layer", "wl"]) n_interactions += 1 # can only have one interaction with planar surface regardless of # angle of incidence @@ -280,7 +280,7 @@ def analytical_start(nks, { "I": R_total, "direction": final_R_directions, - "n_interactions": n_interactions, + "n_interactions": np.array([n_interactions]), } ) @@ -288,7 +288,7 @@ def analytical_start(nks, { "I": I_remaining * T, "direction": final_T_directions, - "n_interactions": n_interactions, + "n_interactions": np.array([n_interactions]), "theta_t": theta_t, } ) @@ -330,7 +330,6 @@ def analytical_start(nks, alphas[surf_index + 1], # units? np.ones_like(T_data.theta_t), depths[surf_index + 1], - I_thresh, initial_dir, ) @@ -686,8 +685,6 @@ def analytical_per_face(current_surf, "I": R_total, "direction": final_R_directions, "n_interactions": n_interactions, - # "theta_out_R": theta_out_R, - # "phi_out_R": phi_out_R, } ) @@ -697,7 +694,6 @@ def analytical_per_face(current_surf, final_T_directions = xr.DataArray(final_T_directions, dims=["unique_direction", "xyz", "wl"]) final_T_n_interactions = xr.DataArray(final_T_n_interactions, dims=["unique_direction"]) theta_out_T = xr.DataArray(theta_out_T, dims=["unique_direction", "wl"]) - # phi_out_T = xr.DataArray(phi_out_T, dims=["outgoing", "wl"]) T_data = xr.Dataset( { @@ -705,11 +701,9 @@ def analytical_per_face(current_surf, "direction": final_T_directions, "n_interactions": final_T_n_interactions, "theta_t": theta_out_T, - # "phi_out_T": phi_out_T, } ) - return R_data, A_data, T_data @@ -860,7 +854,7 @@ def lambertian_scattering(strt, save_location, options): A_polar_rear = np.sum(np.mean(A_surf_weighted_rear, 2), 1) # calculate travel distance for each ray - I_rear = I_theta[:, None] * np.exp(-strt.widths[0] * strt.mats[1].alpha(options.wavelength[None, :]) / np.cos(theta_lamb)[:, None]) + I_rear = I_theta[:, None] * np.exp(-strt.widths[mat_index - 1] * strt.mats[mat_index].alpha(options.wavelength[None, :]) / np.cos(theta_lamb)[:, None]) R_1 = np.sum(I_theta[:, None]*R_polar_front, axis=0) R_2 = np.sum(I_theta[:, None]*R_polar_rear, axis=0) @@ -921,18 +915,42 @@ def lambertian_scattering(strt, save_location, options): # Concatenate the two xarrays along the new dimension merged = xr.concat([initial_down, initial_up], dim="direction") - return merged, front_surf_P, rear_surf_P, [R_1, R_2] + result_list.append([merged, front_surf_P, rear_surf_P, [R_1, R_2]]) + + # stack merged, front_surf_P and rear_surf_P arrays along a nr + + return result_list def calculate_lambertian_profile(strt, I_wl, options, initial_direction, - lambertian_results, alphas, position): + lambertian_R_results, alphas, position, + total_A): + def traverse_lambertian(width, theta, alpha, I_i, positions, direction): + + ratio = alpha / np.real(np.abs(np.cos(theta))) + DA_u = I_i[:, None] * ratio[:, None] * np.exp((-ratio[:, None] * positions[None, :])) + I_back = I_i * np.exp(-ratio * width) + + if direction == -1: + DA_u = np.flip(DA_u) + + intgr = np.trapz(DA_u, positions, axis=1) + + DA = np.divide( + ((I_i[:, None] - I_back[:, None]) * DA_u).T, intgr, + ).T + + DA[intgr == 0] = 0 + + return DA, I_back + I_theta = np.cos(theta_lamb) I_theta = I_theta / np.sum(I_theta) profile_wl = np.zeros((len(I_wl), len(position))) - [R_top, R_bot] = lambertian_results + [R_top, R_bot] = lambertian_R_results if initial_direction == 1: R1 = R_bot # CHECK @@ -942,9 +960,11 @@ def calculate_lambertian_profile(strt, I_wl, options, initial_direction, R1 = R_top R2 = R_bot - for i1, I0 in enumerate(I_wl): + cont_wl_ind = np.where(I_wl > options.I_thresh)[0] + + for i1 in cont_wl_ind: - I = I0 + I = I_wl[i1] I_angular = I * I_theta direction = -initial_direction # first thing that happens is reflection! DA = np.zeros((len(theta_lamb), len(position))) @@ -956,13 +976,12 @@ def calculate_lambertian_profile(strt, I_wl, options, initial_direction, # absorption - DA_pass, _, I_angular = traverse_vectorised( + DA_pass, I_angular = traverse_lambertian( strt.widths[0]*1e6, theta_lamb, alphas[i1], I_angular, position, - options.I_thresh, direction, ) @@ -972,21 +991,23 @@ def calculate_lambertian_profile(strt, I_wl, options, initial_direction, direction = -direction - DA_pass, _, I_angular = traverse_vectorised( + DA_pass, I_angular = traverse_lambertian( strt.widths[0], theta_lamb, alphas[i1], I_angular, position, - options.I_thresh, direction, ) + DA += DA_pass direction = -direction I = np.sum(I_angular) - profile_wl[i1] = np.sum(DA, axis=0) + sum_over_angles = np.sum(DA, axis=0) + int_I = np.trapz(sum_over_angles, position) + profile_wl[i1] = (total_A[i1]/int_I)*np.sum(DA, axis=0) return profile_wl diff --git a/rayflare/ray_tracing/rt.py b/rayflare/ray_tracing/rt.py index ec86600..e5ec136 100644 --- a/rayflare/ray_tracing/rt.py +++ b/rayflare/ray_tracing/rt.py @@ -784,7 +784,7 @@ def __init__(self): pass def isel(self, wl): - return None + return State(I=0) class rt_structure: """Set up structure for RT calculations. @@ -910,14 +910,14 @@ def calculate(self, options): lambertian_approximation = options.lambertian_approximation analytical_rt = options.analytical_ray_tracing - if options.lambertian_approximation: - if self.lambertian_results is None: - raise ValueError("Lambertian approximation results not found - must provide" - "options argument when initialising rt_structure.") - - if len(self.mats) > 3: - raise ValueError("Lambertian approximation currently only implemented for structures" - "with a single bulk material and two interface textures.") + # if options.lambertian_approximation: + # if self.lambertian_results is None: + # raise ValueError("Lambertian approximation results not found - must provide" + # "options argument when initialising rt_structure.") + # + # if len(self.mats) > 3: + # raise ValueError("Lambertian approximation currently only implemented for structures" + # "with a single bulk material and two interface textures.") if not options["parallel"]: n_jobs = 1 @@ -1077,9 +1077,9 @@ def calculate(self, options): # plt.legend() # plt.show() - remaining_I_per_wl = np.sum(prop_rays.I, axis=0) + # remaining_I_per_wl = np.sum(prop_rays.I, axis=0) # pr.enable() - wl_continuing = np.where(remaining_I_per_wl > 1e-9)[0] + # wl_continuing = np.where(remaining_I_per_wl > 1e-9)[0] R_to_add = np.sum(overall_R.I, axis=0) @@ -1096,7 +1096,7 @@ def calculate(self, options): A_interface_to_add.append(0) else: - wl_continuing = np.arange(len(wavelengths)) + # wl_continuing = np.arange(len(wavelengths)) prop_rays = dummy_prop_rays() # R_to_add = 0 # T_to_add = 0 @@ -1131,12 +1131,10 @@ def calculate(self, options): initial_dir, periodic, lambertian_approximation, - analytical_rt, tmm_args + [wavelengths[i1]], - # self.lambertian_results, prop_rays.isel(wl=i1), ) - for i1 in wl_continuing) + for i1 in range(len(wavelengths))) # pr.disable() @@ -1183,7 +1181,7 @@ def calculate(self, options): A_per_interface = 0 interface_profiles = 0 - non_abs = np.logical_and(~np.isnan(thetas), np.abs(thetas) != 10) + non_abs = np.logical_and(~np.isnan(thetas), np.abs(thetas) < 10) refl = np.logical_and( non_abs, @@ -1206,43 +1204,59 @@ def calculate(self, options): absorption_profiles[absorption_profiles < 0] = 0 - A_layer[:, 1] = A_layer[:, 1] - # process A_interfaces - if np.any(np.abs(thetas) == 10): # Lambertian scattering + if np.any(np.abs(thetas) >= 10): # Lambertian scattering # +ve if travelling down, about to hit rear surface. # stop normal RT right before next interaction with surface, but AFTER taking into account bulk # absorption on this pass - direction = np.sign(thetas[np.abs(thetas)==10])[0] - lambertian_RAT = self.lambertian_results[0].sel(direction=direction) - lambertian_A1 = self.lambertian_results[1].sel(direction=direction) - lambertian_A2 = self.lambertian_results[2].sel(direction=direction) + # per wavelength, find how many rays need to be accounted for with the Lambertian model: + # which direction it's travelling in, which material it is in, + # and the total intensity in those rays + + # indexing: wavelengths, material, direction + I_transmitted = np.zeros(len(wavelengths)) + + for i1 in range(len(mats)-2): - # where np.abs(thetas) == 10, want to divide remaining intensity in Is at the same location correctly between - # reflection, interface absorption, and transmission - n_lambertian = np.sum(np.abs(thetas) == 10, axis=1) - I_lambertian = np.array([np.sum(Is[i1][np.abs(thetas[i1]) == 10]) for i1 in - range(len(wavelengths))]) / (nx * ny * n_reps) - I_RAT = lambertian_RAT * I_lambertian + direction_in_mat = np.sign(thetas[np.abs(thetas/10) == (i1+1)])[0] + I_per_mat = [np.sum(Is[l][np.all((np.abs(thetas[l]/10) == (i1+1), np.sign(thetas[l]) == direction_in_mat), axis=0)]) for l in range(len(wavelengths))] + I_lambertian = np.array(I_per_mat)/(n_reps*nx*ny) + I_transmitted - R += I_RAT.sel(event='R') - T += I_RAT.sel(event='T') - A_layer[:, 1] += I_RAT.sel(event='A_bulk') + lambertian_RAT = self.lambertian_results[i1][0].sel(direction=direction_in_mat) + lambertian_A1 = self.lambertian_results[i1][1].sel(direction=direction_in_mat) + lambertian_A2 = self.lambertian_results[i1][2].sel(direction=direction_in_mat) + + # where np.abs(thetas) == 10, want to divide remaining intensity in Is at the same location correctly between + # reflection, interface absorption, and transmission + + I_RAT = lambertian_RAT * I_lambertian + + R += I_RAT.sel(event='R').data + + if i1 == len(mats) - 3: # final non-transmission bulk + T += I_RAT.sel(event='T').data + else: + I_transmitted = I_RAT.sel(event='T').data - if sum(self.tmm_or_fresnel) > 0: - add_frontsurf = (I_lambertian * lambertian_A1).data.T - add_backsurf = (I_lambertian * lambertian_A2).data.T + A_layer[:, i1+1] += I_RAT.sel(event='A_bulk').data - A_per_interface[0] = A_per_interface[0] + add_frontsurf - A_per_interface[1] = A_per_interface[-1] + add_backsurf + add_profile = calculate_lambertian_profile(self, I_lambertian, options, + direction_in_mat, + self.lambertian_results[i1][3], alphas[i1+1], + depths[i1+1], + I_RAT.sel(event='A_bulk').data) - add_profile = calculate_lambertian_profile(self, I_lambertian, options, direction, - self.lambertian_results[3], alphas[1], depths[1]) + absorption_profiles[:, depth_indices[i1+1]] += add_profile - absorption_profiles += add_profile + if sum(self.tmm_or_fresnel) > 0: + add_frontsurf = (I_lambertian * lambertian_A1).data.T + add_backsurf = (I_lambertian * lambertian_A2).data.T + + A_per_interface[i1] = A_per_interface[i1] + add_frontsurf + A_per_interface[i1+1] = A_per_interface[i1+1] + add_backsurf # add results from analytical ray tracing: if analytical_rt > 0: @@ -1266,7 +1280,7 @@ def calculate(self, options): for wl_ind in wl_ind_anlt: # TODO: transmission not included here - A_interf = np.sum([np.sum(A_interf[wl_ind]) if not isinstance(A_interf, int) else 0 for A_interf in A_interface_to_add]) + A_interf = np.sum([np.sum(A_intf[wl_ind]) if not isinstance(A_intf, int) else 0 for A_intf in A_interface_to_add]) total_A = np.sum(A_bulk_to_add[wl_ind]) + A_interf n_absorbed = int(np.round(total_A*max_rays, 0)) @@ -1315,8 +1329,13 @@ def calculate(self, options): scale_I_R = R_weights*max_rays/n_rays_per_R_direction - theta_R = np.arccos(overall_R.direction.isel(wl=wl_ind, xyz=2)) - phi_R = np.arctan(overall_R.direction.isel(wl=wl_ind, xyz=1)/overall_R.direction.isel(wl=wl_ind, xyz=0)) + if 'wl' in overall_R.direction.dims: + theta_R = np.arccos(overall_R.direction.isel(wl=wl_ind, xyz=2)) + phi_R = np.arctan(overall_R.direction.isel(wl=wl_ind, xyz=1)/overall_R.direction.isel(wl=wl_ind, xyz=0)) + + else: + theta_R = np.arccos(overall_R.direction.isel(xyz=2)) + phi_R = np.arctan(overall_R.direction.isel(xyz=1)/overall_R.direction.isel(xyz=0)) current_start = n_remaining[wl_ind] + n_absorbed @@ -1421,8 +1440,7 @@ def parallel_inner( initial_mat, initial_dir, periodic, - lambertian_approximation, - analytical_rt, + lambertian_approximation=0, tmm_args=None, existing_rays=None, ): @@ -1526,7 +1544,7 @@ def parallel_inner( end_ind = nx*ny*np.ones(n_reps, dtype=int) - if existing_rays is not None: + if np.sum(existing_rays.I) > 1e-9: ds, i_mats, i_dirs, surf_inds, n_remaining, I_in, n_inter_in, n_passes_in = make_rt_args(existing_rays, xs, ys, n_reps) stop_before = int(np.ceil(n_remaining/(nx*ny))) @@ -2029,73 +2047,13 @@ def single_ray_stack( # should end when either material = final material (len(materials)-1) & direction == 1 or # material = 0 & direction == -1 - # print("NEW RAY") - profile = np.zeros(len(z_pos)) # do everything in microns A_per_layer = np.zeros(len(widths)) - # direction: start travelling downwards; 1 = down, -1 = up - - # surf_below = mat_i - # surf_above = mat_i - 1 - A_interface_array = 0 A_interface_index = 0 - # max_below = np.max(surfaces[surf_below].Points[:, 2]) - # min_above = np.min(surfaces[surf_above].Points[:, 2]) - - # all of this can be done outside - will be the same for every ray! - # print(max_below, min_above) - - ######## - - # same for all wavelengths & rays - - # generally same across wavelengths, but can be changed by analytical - # ray tracing happening first - # if direction == 1 and mat_i > 0: - # surf_index = mat_i - # z_offset = -cum_width[mat_i - 1] - 1e-8 - # # print('z_offset', z_offset, r_a_0) - # - # elif direction == 1 and mat_i == 0: - # surf_index = 0 - # z_offset = r_a_0[2] - # - # else: - # surf_index = mat_i - 1 - # z_offset = -cum_width[mat_i] + 1e-8 - - # different for each ray, but same across wavelengths - should pregenerate as an array - # r_a = r_a_0 + np.array([x, y, 0]) - # r_b = np.array([x, y, 0]) - # - # d = (r_b - r_a) / np.linalg.norm( - # r_b - r_a - # ) # direction (unit vector) of ray. Always downwards! - - # print(d, -r_a_0/np.linalg.norm(r_a_0)) - # d = -r_a_0/np.linalg.norm(r_a_0) - # want to translate along d so r_a is in between the correct surfaces - # first translate to z = 0: - - # nd = -r_a[2] / d[2] - # - # # print(r_a[2], d[2], nd) - # - # r_a = r_a + nd * d - # - # r_a[2] = z_offset - # - # if direction != 1: - # d[2] = -d[2] - - ##### above here outside - - # don't need direction_i, mat_i after this. - stop = False I = I_in @@ -2221,7 +2179,7 @@ def single_ray_stack( if lambertian_approximation and n_passes >= lambertian_approximation: # choose a direction randomly, with probability determined by Lambertian scattering stop = True - theta = 10*direction # +ve if travelling down, about to hit rear surface. + theta = 10*direction*mat_i # +ve if travelling down, about to hit rear surface. # stop right before next interaction with surface, but AFTER taking into account bulk # absorption on this pass @@ -2340,6 +2298,7 @@ def traverse(width, theta, alpha, x, y, I_i, positions, I_thresh, direction): def decide_RT_Fresnel(n0, n1, theta, d, N, side, pol, rnd, wl=None, lookuptable=None): + ratio = np.clip(np.real(n1) / np.real(n0), -1, 1) if abs(theta) > np.arcsin(ratio):