import numpy as np rng = np.random.default_rng() x_dim = 7 y_dim = 7.3 d = 20.6 - 6.8 - 7 - 2.5 def generate_ray(): x_start = rng.uniform(0, x_dim) y_start = rng.uniform(0, y_dim) phi = rng.uniform(0, 1) * 2 * np.pi theta = np.arccos(np.sqrt(rng.uniform(0, 1))) x_final = x_start + np.cos(phi) * d * np.tan(theta) y_final = y_start + np.sin(phi) * d * np.tan(theta) return (0 <= x_final and x_dim >= x_final and 0 <= y_final and y_dim >= y_final) num = 1000000 count_in = 0 for i in range(num): if generate_ray(): count_in += 1 print(count_in)