# compare results with regular & random ray positions, and also compare results between RT + TMM and # RT + TMM + angular redistribution matrix method import seaborn as sns from cycler import cycler from pytest import approx import numpy as np from sparse import load_npz import xarray as xr from solcore.structure import Layer from solcore import material # rayflare imports from rayflare.textures.standard_rt_textures import planar_surface, regular_pyramids from rayflare.structure import Interface, BulkLayer, Structure from rayflare.matrix_formalism import process_structure, calculate_RAT from rayflare.options import default_options from rayflare.ray_tracing import rt_structure # Thickness of bottom Ge layer bulkthick = 0.3e-6 wavelengths = np.linspace(300, 1500, 6) * 1e-9 # wavelengths = wavelengths # set options options = default_options() options.wavelengths = wavelengths options.project_name = "rt_tmm_comparisons_planar7" options.n_rays = 2000000 options.n_theta_bins = 120 options.lookuptable_angles = 400 options.parallel = True options.I_thresh = 1e-3 options.bulk_profile = True options.randomize_surface = True options.c_azimuth = 1 options.periodic = True options.bulk_profile = True options.depth_spacing_bulk = 0.1e-9 options.depth_spacing = 0.1e-9 options.only_incidence_angle = False options.parallel = True options.theta_in = 0.3 # set up Solcore materials Ge = material("Ge")() GaAs = material("GaAs")() GaInP = material("GaInP")(In=0.5) Ag = material("Ag")() SiN = material("Si3N4")() Air = material("Air")() Ta2O5 = material("TaOx1")() # Ta2O5 (SOPRA database) MgF2 = material("MgF2")() # MgF2 (SOPRA database) front_materials = [ # Layer(120e-9, MgF2), # Layer(74e-9, Ta2O5), Layer(464e-9, GaInP), Layer(1682e-9, GaAs), Layer(200e-9, Ge), ] back_materials = [Layer(200e-9, Ge), Layer(100e-9, SiN)] # RT/TMM, matrix framework bulk_Ge = BulkLayer(bulkthick, Ge, name="Ge_bulk") # bulk thickness in m ## RT with TMM lookup tables surf_pyr = regular_pyramids() # [texture, flipped texture] surf_pyr_inv = regular_pyramids(upright=False, elevation_angle=20) surf_planar = planar_surface() front_surf = Interface( "RT_TMM", layers=front_materials, texture=surf_planar, name="GaInP_GaAs_RT_planar", coherent=True, prof_layers=[2, 3] ) # front_surf = Interface( # "TMM", # layers=front_materials, # name="GaInP_GaAs_RT_planar_TMM", # coherent=True, # prof_layers=[2, 3] # ) back_surf = Interface( "RT_TMM", layers=back_materials, texture=surf_pyr_inv, name="SiN_Ag_RT", coherent=True, prof_layers=[1, 2], ) SC = Structure([front_surf, bulk_Ge, back_surf], incidence=Air, transmission=Ag) process_structure(SC, options) results_RT = calculate_RAT(SC, options) prof_front = results_RT[2][0] prof_back = results_RT[2][1] prof_Ge = results_RT[3][0] results_per_pass_RT = results_RT[1] # sum over passes results_per_layer_front_RT = np.sum(results_per_pass_RT["a"][0], 0) front_surf = planar_surface( interface_layers=front_materials, prof_layers=[2, 3], ) # pyramid size in microns back_surf = regular_pyramids(upright=False, elevation_angle=20, interface_layers=back_materials, prof_layers=[1, 2], ) # pyramid size in microns rtstr = rt_structure( [front_surf, back_surf], [Ge], [bulkthick], Air, Ag, options, use_TMM=True ) # RT + TMM options.n_rays = 10000 result_RT_only = rtstr.calculate(options) rt_front = result_RT_only["interface_profiles"][0] rt_back = result_RT_only["interface_profiles"][1] rt_Ge = result_RT_only["profile"] import matplotlib.pyplot as plt pal = sns.cubehelix_palette(n_colors=len(wavelengths)) cols = cycler("color", pal) params = { "axes.prop_cycle": cols, } plt.rcParams.update(params) plt.figure() plt.plot(wavelengths * 1e9, result_RT_only['R'], label='R RT') plt.plot(wavelengths * 1e9, results_RT[0]['R'][0], '--', label='R ARM') plt.plot(wavelengths * 1e9, result_RT_only['T'], label='T RT') plt.plot(wavelengths * 1e9, results_RT[0]['T'][0], '--', label='T ARM') plt.plot(wavelengths * 1e9, result_RT_only["A_per_interface"][0]) plt.plot(wavelengths * 1e9, results_per_layer_front_RT, '--') plt.plot(wavelengths * 1e9, result_RT_only["A_per_layer"][:, 0]) plt.plot(wavelengths * 1e9, results_RT[0]["A_bulk"][0], '--') plt.xlabel("Wavelength (nm)") plt.legend() plt.show() pos = np.arange(0, options.depth_spacing*1e9*len(rt_front.T), options.depth_spacing*1e9) plt.figure() plt.plot(pos, rt_front.T) plt.plot(pos, prof_front.T, '--') plt.xlabel("Position (nm)") plt.ylabel("a(z)") plt.legend([str(np.round(wl*1e9, 0)) for wl in wavelengths]) plt.tight_layout() plt.show() pos = np.arange(0, options.depth_spacing*1e9*len(rt_back.T), options.depth_spacing*1e9) plt.figure() plt.semilogy(pos, rt_back.T) plt.semilogy(pos, prof_back.T, '--') plt.tight_layout() plt.xlabel("Position (nm)") plt.ylabel("a(z)") plt.show() pos = np.arange(0, options.depth_spacing_bulk*1e9*len(rt_Ge.T), options.depth_spacing_bulk*1e9) plt.figure() plt.plot(pos,rt_Ge.T) plt.plot(pos,prof_Ge.T / 1e9, '--') plt.xlabel("Position (nm)") plt.ylabel("a(z)") plt.show() int_front_rt = np.trapz(rt_front, dx=0.1, axis=1) int_front_arm = np.trapz(prof_front, dx=0.1, axis=1) plt.figure() plt.plot(wavelengths*1e9, int_front_arm, label="Integrated ARM") plt.plot(wavelengths*1e9, int_front_rt, '-.', label="Integrated RT") plt.plot(wavelengths*1e9, np.sum(results_per_layer_front_RT, axis=1), '--', label="ARM per layer total") plt.ylabel("A") plt.xlabel("Wavelength (nm)") plt.show() # total A is correct, I think, but integrated is notß= # assert int_front_rt == approx(np.sum(result_RT_only["A_per_interface"][0], axis=1)) # assert int_front_arm == approx(np.sum(results_per_layer_front_RT, axis=1)) # # assert result_RT_only["R"] == approx(results_RT[0]["R"][0], abs=0.05) # assert result_RT_only["T"] == approx(results_RT[0]["T"][0], abs=0.05) # assert result_RT_only["A_per_interface"][0] == approx(results_per_layer_front_RT, abs=0.05) # assert result_RT_only["A_per_layer"][:, 0] == approx(results_RT[0]["A_bulk"][0], abs=0.05) # # RT_matrix = load_npz("/Users/phoebe/RayFlare_results/rt_tmm_comparisons_matrix/GaInP_GaAs_RTrearA.npz").todense() # TMM_matrix = load_npz("/Users/phoebe/RayFlare_results/rt_tmm_comparisons_matrix/GaInP_GaAs_TMMrearA.npz").todense() # # plt.figure() # plt.plot(RT_matrix[1,2,:], 'o') # plt.plot(TMM_matrix[1,2,:], 'x') # plt.show() # RT_matrix = load_npz("/Users/phoebe/RayFlare_results/rt_tmm_comparisons_planar3/GaInP_GaAs_RTrearA.npz").todense() # TMM_matrix = load_npz("/Users/phoebe/RayFlare_results/rt_tmm_comparisons_planar2/GaInP_GaAs_RTrearA.npz").todense() # # plt.figure() # plt.plot(RT_matrix[1,2,:], 'o') # plt.plot(TMM_matrix[1,2,:], 'x') # plt.show() # RT_matrix = load_npz("/Users/phoebe/RayFlare_results/rt_tmm_comparisons_planar3/GaInP_GaAs_RTfrontA.npz").todense() # TMM_matrix = load_npz("/Users/phoebe/RayFlare_results/rt_tmm_comparisons_planar2/GaInP_GaAs_RTfrontA.npz").todense() # # plt.figure() # plt.plot(RT_matrix[1,2,:], 'o') # plt.plot(TMM_matrix[1,2,:], 'x') # plt.ylim(0, 0.7) # plt.show() lookup_RT = xr.open_dataset("/Users/phoebe/RayFlare_results/rt_tmm_comparisons_planar7/GaInP_GaAs_RT_planarrearprofmat.nc") lookup_TMM = xr.open_dataset("/Users/phoebe/RayFlare_results/rt_tmm_comparisons_planar7/GaInP_GaAs_RT_planar_TMMrearprofmat.nc") wl_select = wavelengths[-1] fig, (ax1, ax2) = plt.subplots(1,2, figsize=(10,4)) lookup_RT.profile.sel(wl=wl_select*1e9).plot.imshow(ax=ax1) lookup_TMM.profile.sel(wl=wl_select).plot.imshow(ax=ax2) plt.show() wl_select = wavelengths[-4] fig, (ax1, ax2) = plt.subplots(1,2, figsize=(10,4)) lookup_RT.profile.sel(wl=wl_select*1e9).plot.imshow(ax=ax1) lookup_TMM.profile.sel(wl=wl_select).plot.imshow(ax=ax2) plt.show() # from rayflare.transfer_matrix_method import tmm_structure # # tmm_str = tmm_structure(front_materials[::-1], incidence=Ge, transmission=Air) # # options.theta_in = 0 # # res = tmm_str.calculate(options)