diff --git a/examples/Si_optimize_diffraction_angle.py b/examples/Si_optimize_diffraction_angle.py index 8e8bb88..8b70b9c 100644 --- a/examples/Si_optimize_diffraction_angle.py +++ b/examples/Si_optimize_diffraction_angle.py @@ -8,11 +8,20 @@ from rayflare.analytic.diffraction import get_order_directions import seaborn as sns -# Paper: https://doi.org/10.1002/adom.201700585 -MgF2 = material("MgF2")() # MgF2 (SOPRA database) + +# Optimize an Si diffraction grating of hexagonally-arranged pillars on bulk Si +# to optimally diffract as much incident light as possible into some target angle + Air = material("Air")() Si = material("Si")() + +# Class which defines the optimization problem. This must have a function called 'fitness' +# which returns the fitness of the solution, and a function called 'get_bounds' which returns +# the bounds of the solution space. The fitness function must take a single argument, which is +# the solution to be evaluated. The bounds function must return two lists, the first of which +# contains the lower bounds for each parameter, and the second of which contains the upper bounds +# for each parameter. class optimize_surface(): def __init__(self, @@ -21,9 +30,8 @@ def __init__(self, angle_tolerance=10*np.pi/180, pol='s', ): + self.wavelengths = wavelengths - # self.angle = target_diffraction_angle - # self.tolerance = angle_tolerance self.target = target_diffraction_angle self.tolerance = angle_tolerance self.min_angle = np.pi - target_diffraction_angle - angle_tolerance @@ -31,7 +39,6 @@ def __init__(self, self.pol = pol def calculate(self, x): - # wavelengths = np.array([1e-6]) # lattice vectors for the grating. Units are in nm! # Note: annoyingly have to set options here internally every time, since the options object @@ -67,28 +74,29 @@ def fitness(self, x): RAT, options, size = self.calculate(x) - # order_power = np.flip(np.argsort(RAT['T_per_order'], axis=1), axis=1) # dimensions: (wl, orders) + # diffraction orders retained in the Fourier transform: orders = np.array(RAT['basis_set']) # dimensions: (orders, 2) - # orders = orders[order_power] # dimensions: (wl, orders, 2) - # power_sorted = np.sort(RAT['T_per_order'], axis=0) # dimensions: (wl, orders - - # orders = orders[power_sorted > 1e-7] - # power_sorted = power_sorted[power_sorted > 1e-7] + # Calculate which directions in air and Si each of these diffraction orders correspond + # to diffracted_directions = get_order_directions(self.wavelengths * 1e9, size, 4, Air, Si, options.theta_in, 0, np.pi / 3) orders_analytical = diffracted_directions['order_index'] + # Match the diffraction orders from the RCWA calculation to the analytical calculation orders_analytical_ind = [np.where( np.all((orders_analytical[0].flatten() == x[0], orders_analytical[1].flatten() == x[1]), axis=0))[0][0] for x in orders] + # pick up angles in silicon (transmission) theta_analytical = diffracted_directions['theta_t'][:, orders_analytical_ind] + # mask to select only the orders which are within the target angle +/- the tolerance mask = np.all((theta_analytical > self.min_angle, theta_analytical < self.max_angle), axis=0) + # Sum the transmitted power at all wavelengths which falls within the allowed angles T_in_condition = np.sum(RAT['T_per_order'][mask]) return [-T_in_condition/len(self.wavelengths)] @@ -97,45 +105,52 @@ def get_bounds(self): # [lower bounds for all parameters], [upper bounds for all parameters] return [(100, 0, 200), (10000, 0.5, 3000)] - - -crit_angle = np.arcsin(Air.n(5e-6)/Si.n(5e-6)) +# wavelengths: 4.5 - 5.5 um, 10 points test_wl = np.linspace(4.5, 5.5,10)*1e-6 +# create an instance of the optimize_surface class which will be passed to Pygmo. We can choose +# the wavelengths, target angle, and angle tolerance here. Can also set the polarization ('s' by +# default) obj_class = optimize_surface(wavelengths=test_wl, target_diffraction_angle=45*np.pi/180, angle_tolerance=10*np.pi/180) + +# Create the pygmo problem to solve: prob = pg.problem(obj_class) -n_generations = 40 -pop_size = 5*len(prob.get_bounds()[0]) +# Number of parameters to optimize (length of candidate vectors) +n_params = len(prob.get_bounds()[0]) -# algo = pg.algorithm(pg.de(gen=n_generations)) -# pop = pg.population(prob, pop_size) -# pop = algo.evolve(pop) +# Number of generations and population size: +n_generations = 40 +pop_size = 5*n_params algo = pg.algorithm(pg.de(gen=1)) # note: I am running this one generation at a time so I can get the # best population at each step and plot how it evolves. If you just care about the final result, # you can just set gen=n_generations here and not have the for loop below. -n_params = len(prob.get_bounds()[0]) +# algo = pg.algorithm(pg.de(gen=n_generations)) +# pop = pg.population(prob, pop_size) +# pop = algo.evolve(pop) + +# arrays of zeros to store results best_f = np.zeros(n_generations) mean_f = np.zeros(n_generations) best_x = np.zeros((n_generations, n_params)) +# create initial populate, will be generated randomly within upper and lower bounds of each parameter pop = pg.population(prob, pop_size) for i1 in range(n_generations): pop = algo.evolve(pop) - print(i1, pop.champion_f[0]) best_f[i1] = pop.champion_f[0] mean_f[i1] = np.mean(pop.get_f()) best_x[i1] = pop.champion_x print(pop.champion_x, pop.champion_f) - +# Plot best and mean fitness in population at each generation fig, ax1 = plt.subplots(1, 1, figsize=(5, 4)) ax1.plot(np.arange(n_generations) + 1, best_f, 'ro-', label="Champion fitness") ax1.plot(np.arange(n_generations) + 1, mean_f, 'kx-', label="Mean fitness") diff --git a/examples/lens_hyperhemisphere_rt.py b/examples/lens_hyperhemisphere_rt.py index e18d4e4..21d5d9b 100644 --- a/examples/lens_hyperhemisphere_rt.py +++ b/examples/lens_hyperhemisphere_rt.py @@ -25,8 +25,8 @@ # (32 GB RAM, M1 chip, 10 cores) n_thetas = 50 # 100 used for paper data -exp_points = 15 # 2**exp_points on surface of WHOLE sphere. exp_points = 15 used for paper data -nxs = 70 # 70 points used for paper data +exp_points = 13 # 2**exp_points on surface of WHOLE sphere. exp_points = 15 used for paper data +nxs = 30 # 70 points used for paper data thetas = np.linspace(0, np.pi / 2 - 0.05, n_thetas) # 100 angles used for paper data # structure: @@ -57,7 +57,7 @@ options.periodic = 0 options.wavelength = np.array([6e-6]) options.parallel = False -options.n_rays = 4*nxs ** 2 # every emission point on the surface 4 times +options.n_rays = 5*nxs ** 2 # every emission point on the surface 5 times options.theta = 0.1 options.nx = nxs @@ -120,7 +120,7 @@ def parallel_theta(options, th): options.theta_in = th result = rtstr.calculate(options) - return result['T'][0], np.mean(result['n_interactions']), result['thetas'][0] + return result['T'][0], result['n_interactions'], result['thetas'][0] if os.path.isfile( "sphere_raytrace_totalT_2e" + str(exp_points) + "_" + str(nxs) + "_points_" + str(options.n_rays) + "_rays_2.txt" @@ -151,8 +151,8 @@ def parallel_theta(options, th): "sphere_raytrace_totalT_2e" + str(exp_points) + "_" + str(nxs) + "_points_" + str(options.n_rays) + "_rays.txt", T_total, ) - np.savetxt( - "sphere_raytrace_ninter_2e" + str(exp_points) + "_" + str(nxs) + "_points_" + str(options.n_rays) + "_rays.txt", + np.save( + "sphere_raytrace_ninter_2e" + str(exp_points) + "_" + str(nxs) + "_points_" + str(options.n_rays) + "_rays.npy", n_interactions, ) np.savetxt( @@ -160,22 +160,42 @@ def parallel_theta(options, th): theta_distribution, ) -min_angle_1 = np.pi - 17.5 * np.pi / 180 -min_angle_2 = np.pi - np.pi * 45 / 180 +min_angle_1 = np.pi - 20 * np.pi / 180 +min_angle_2 = np.pi - 45 * np.pi / 180 -T_175 = np.array([np.sum(x > min_angle_1) / options.n_rays for x in theta_distribution]) +T_20 = np.array([np.sum(x > min_angle_1) / options.n_rays for x in theta_distribution]) T_45 = np.array([np.sum(x > min_angle_2) / options.n_rays for x in theta_distribution]) plt.figure() -plt.plot(thetas * 180 / np.pi, T_total, color=pal[0]) +plt.plot(thetas * 180 / np.pi, T_total, color=pal[0], label=r"$\delta <$ 90°") plt.scatter(thetas * 180 / np.pi, T_total, edgecolors=pal[0], facecolors="none") -plt.plot(thetas * 180 / np.pi, T_45, color=pal[1]) +plt.plot(thetas * 180 / np.pi, T_45, color=pal[1], label=r"$\delta <$ 45°") plt.scatter(thetas * 180 / np.pi, T_45, edgecolors=pal[1], facecolors="none") -plt.plot(thetas * 180 / np.pi, T_175, color=pal[2]) -plt.scatter(thetas * 180 / np.pi, T_175, edgecolors=pal[2], facecolors="none") +plt.plot(thetas * 180 / np.pi, T_20, color=pal[2], label=r"$\delta <$ 20°") +plt.scatter(thetas * 180 / np.pi, T_20, edgecolors=pal[2], facecolors="none") plt.xlim(0, 90) plt.xlabel(r"$\beta$ (rads)") plt.ylabel("Transmission") +plt.legend() +plt.show() + +# selection n_interactions for rays escaping within a cone of 45 degrees: +options.n_rays = 1 +sample_result = result = rtstr.calculate(options) +n_interactions_escape = [n_interactions[i, 0, theta_distribution[i] > min_angle_2] for i in range(len(thetas))] +mean_interactions = [np.mean(x) for x in n_interactions_escape] +min_interactions = [np.min(x) for x in n_interactions_escape] +max_interactions = [np.max(x) for x in n_interactions_escape] + + +# plot the mean number of interactions for rays escaping in the forward direction, and shade the range: + +plt.figure() +plt.plot(thetas * 180 / np.pi, mean_interactions, color='k', label="Mean") +# plt.fill_between(thetas * 180 / np.pi, min_interactions, max_interactions, color='k', alpha=0.3) +plt.xlabel(r"$\beta$ (rads)") +plt.ylabel("Number of interactions before escape") +plt.legend() +plt.show() -plt.show() \ No newline at end of file diff --git a/rayflare/ray_tracing/rt.py b/rayflare/ray_tracing/rt.py index 5ce7f3c..a715012 100644 --- a/rayflare/ray_tracing/rt.py +++ b/rayflare/ray_tracing/rt.py @@ -1066,6 +1066,7 @@ def calculate(self, options): "A_per_interface": A_per_interface, "A_interfaces": A_interfaces, "interface_profiles": interface_profiles, + "xy": [xs, ys], } def calculate_profile(self, options): diff --git a/setup.py b/setup.py index 099c0dd..8f80b98 100644 --- a/setup.py +++ b/setup.py @@ -50,7 +50,7 @@ def gen_data_files(*dirs): install_requires = [ "matplotlib", "scipy", - "numpy", + "numpy<2.0.0", "solcore", "xarray", "sparse<=0.14.0",