numpy version < 2.0

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,17 +30,15 @@ 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
         self.max_angle = np.pi - target_diffraction_angle + angle_tolerance
         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")

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,31 +151,51 @@ 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(
         "sphere_raytrace_thetas_2e" + str(exp_points) + "_" + str(nxs) + "_points_" + str(options.n_rays) + "_rays.txt",
         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()

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):

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",