diff --git a/src/doptics/functions.py b/src/doptics/functions.py index 0f2b5ca..c32822d 100644 --- a/src/doptics/functions.py +++ b/src/doptics/functions.py @@ -3,6 +3,7 @@ import scipy as sp from typing import Callable, List from numpy.typing import ArrayLike +from icecream import ic YL = 0 YR = 2 @@ -10,14 +11,15 @@ XL = 1 XR = 9 -uniform = lambda x: 1 -triangle = lambda x: x -normal_10 = lambda x: np.exp(-0.5*((x-0)/SIGMA)**2) + 0.1 + +def uniform(x): return 1 +def triangle(x): return x +def normal_source(x, xl=XL, xr=XR, sigma=SIGMA): return np.exp(-0.5*((x-((xl+xr)*0.5))/sigma)**2) + 0.01 +def normal_target(x, yl=YL, yr=YR, sigma=SIGMA): return np.exp(-0.5*((x-((yl+yr)*0.5))/sigma)**2) + 0.01 +normal_10 = lambda x: np.exp(-0.5 * ((x - 0) / SIGMA) ** 2) + 0.1 G1 = lambda y1: 1 G2 = lambda y2: 1 - np.abs(y2 - 12) if 11 < y2 < 13 else 0 E = lambda x: 1 / (np.exp(10 * (x - 0.5)) + np.exp(-10 * (x - 0.5))) -normal_target = lambda x: np.exp(-0.5*((x-((YR+YL)*0.5))/SIGMA)**2) + 0.01 -normal_source = lambda x: np.exp(-0.5*((x-((XR+XL)*0.5))/SIGMA)**2) + 0.01 def cdf_sampling_source(source_density: Callable, linear_samples: ArrayLike, total_samples: int = 1000): @@ -48,14 +50,13 @@ def cdf_sampling_source(source_density: Callable, linear_samples: ArrayLike, tot Phi_arr[i] = Phi_arr[i - 1] + elt # abuses np.zeros() Phi_arr = np.roll(Phi_arr, 1) Phi_arr[0] = 0 - # return Phi_arr print(Phi_arr) distr_samples = np.zeros(len(linear_samples)) - inverted_Phi = reversed(Phi_arr) + reversed_Phi = reversed(Phi_arr) quantiles = np.linspace(0, 1, len(linear_samples)) for i in range(len(distr_samples)): - # ladder_index = next(x[0] for x in enumerate(inverted_Phi) if x[1] < delta*i) + # ladder_index = next(x[0] for x in enumerate(reversed_Phi) if x[1] < delta*i) # ladder_index = len([x for x in Phi_arr if x <= + i*delta]) - 1 ladder_index = len([x for x in Phi_arr if x <= + quantiles[i]]) - 1 distr_samples[i] = fine_ladder[ladder_index] @@ -95,25 +96,31 @@ def construct_target_density_intervals_from_angular(angle_density: Callable, np.array([np.cos(large_angle), np.cos(small_angle) * np.sin(abs_large_angle)/np.sin(abs_small_angle)]) ) - print(y1_span) - print(y2_span) + ic(y1_span) + ic(y2_span) if l1 == l2: l2 = 1.4 * l2 y2_span[0] = 1.4 * y2_span[0] y2_span[1] = 1.4 * y2_span[1] - print(y2_span) + ic(y2_span) + + def y1_density(y1): + return (angle_density(np.arccos(y1 / np.linalg.norm(np.array( + [y1 - center[0], l1 - center[1]], dtype=object)))) / + np.linalg.norm(np.array([y1 - center[0], l1 - center[1]], dtype=object)) + ) + + def y2_density(y2): + return (angle_density(np.arccos(y2 / np.linalg.norm(np.array( + [y2 - center[0], l2 - center[1]], dtype=object)))) / + np.linalg.norm(np.array([y2 - center[0], l2 - center[1]], dtype=object)) + ) - y1_density = lambda y1: (angle_density(np.arccos(y1 / np.linalg.norm(np.array( - [y1-center[0], l1-center[1]], dtype=object)))) / - np.linalg.norm(np.array([y1-center[0], l1-center[1]], dtype=object))) - y2_density = lambda y2: (angle_density(np.arccos(y2 / np.linalg.norm(np.array( - [y2-center[0], l2-center[1]], dtype=object)))) / - np.linalg.norm(np.array([y2-center[0], l2-center[1]], dtype=object))) + # for i in np.linspace(small_angle, large_angle, 100): + # print(f'y2({i}) = {y2_density(i)}') - for i in np.linspace(small_angle, large_angle, 100): - print(f'y2({i}) = {y2_density(i)}') y1_span = y1_span + center[0] y2_span = y2_span + center[0] for j in np.linspace(y2_span[0], y2_span[1], 100): @@ -123,15 +130,13 @@ def construct_target_density_intervals_from_angular(angle_density: Callable, return y1_density, y2_density, y1_span, y2_span, l1 + center[1], l2 + center[1] -def f(x): - return 1 +def f(x): return 1 -def g(x): - return x**2 +def g(x, mu=0): return (x-mu)**2 -def rescaling_target_distribution(): +def rescaling_target_distribution() -> None: xl = -10 xr = 4 # f is the density on x diff --git a/src/doptics/solver.py b/src/doptics/solver.py index 4b4cddd..0f60287 100644 --- a/src/doptics/solver.py +++ b/src/doptics/solver.py @@ -47,7 +47,7 @@ def solve_single_mirror_parallel_source(starting_distribution: Callable, target_ # plt.savefig(f'solution-mirror-{FILENAME[result_type]}.png') -def f(x:float, u:float) -> float: +def f(x: float, u: float) -> float: return u**2 diff --git a/src/doptics/symbolic.py b/src/doptics/symbolic.py index d806032..42dd36d 100644 --- a/src/doptics/symbolic.py +++ b/src/doptics/symbolic.py @@ -19,4 +19,3 @@ def u_prime_mv(): dudx = solve(Eq(diff(sol, x), 0), diff(u))[0].simplify() return lambdify([x, y, u, V, t1, t2, L1], dudx), lambdify([x, y, u, V, t1, t2, L1], sol) - diff --git a/src/doptics/two_mirrors.py b/src/doptics/two_mirrors.py index 99f3d9b..390984a 100644 --- a/src/doptics/two_mirrors.py +++ b/src/doptics/two_mirrors.py @@ -37,7 +37,7 @@ def solve_two_mirrors_parallel_source_two_targets(starting_density: Callable, ta """ x_discrete = np.linspace(x_span[0], x_span[1], number_rays) ax = sp.integrate.quad(starting_density, x_span[0], x_span[1])[0] - starting_density_rescaled = lambda x: 1 / ax * starting_density(x) + def starting_density_rescaled(x): return 1 / ax * starting_density(x) res = [] @@ -156,7 +156,7 @@ def solve_two_mirrors_parallel_source_point_target(starting_density: Callable, t small_angle=-0.7 * np.pi, large_angle=-0.3 * np.pi) # new stuff ends - print(f'back in main program') + print('back in main program') print(G1(y1_span[1])) x_discrete = np.linspace(x_span[0], x_span[1], number_rays) @@ -283,4 +283,3 @@ def solve_two_mirrors_parallel_source_point_target(starting_density: Callable, t plt.show() return None - diff --git a/src/ray/ray.py b/src/ray/ray.py index 0fd9cfa..5ec0c93 100644 --- a/src/ray/ray.py +++ b/src/ray/ray.py @@ -4,4 +4,4 @@ class Ray: def __init__(self, x: NDArray[float], y: NDArray[float]): self.x = x - self.y = y \ No newline at end of file + self.y = y diff --git a/src/ray/raytracing/raytracing.py b/src/ray/raytracing/raytracing.py index 4451c4b..5472497 100644 --- a/src/ray/raytracing/raytracing.py +++ b/src/ray/raytracing/raytracing.py @@ -4,6 +4,7 @@ import matplotlib.pyplot as plt import scipy as sp import matplotlib as mpl +from functools import lru_cache mpl.rcParams['figure.dpi'] = 100 COLORS = {'szegedblue': '#3f3fa6'} @@ -27,8 +28,6 @@ def intersection(l1, l2): return np.nan - - """ Example of ray tracing Variable names are same as in the report @@ -36,35 +35,36 @@ def intersection(l1, l2): L1 = 8 L2 = 14 -y1_span = (10.5, 13.5) -y2_span = (11, 13) -x_span = (0, 3) +y1_SPAN = (10.5, 13.5) +y2_SPAN = (11, 13) +x_SPAN = (0, 3) u0 = 4 w0 = 6 # E = lambda x: np.exp(-((x - 1.5) / 1)**2 / 2) / (2 * np.pi)**.5 -G1 = lambda y1: 1 -G2 = lambda y2: 1 +def G1(y1): return 1 +def G2(y2): return 1 + + E = G2 -pm1 = -1 -pm2 = 1 +PM1 = -1 +PM2 = 1 + + +@lru_cache(maxsize=100) +def integrate(f_, left: float, right: float) -> float: + return sp.integrate.quad(f_, left, right)[0] if __name__ == '__main__': # calculate theoretical target from the mapping between emitting density E and target density G1 - G1_int = sp.integrate.quad(G1, y1_span[0], y1_span[1])[0] - G1_fixed = lambda y1: G1(y1) / G1_int - E_int = sp.integrate.quad(E, x_span[0], x_span[1])[0] - E_fixed = lambda x: E(x) / E_int - # a2 = G1_int / sp.integrate.quad(G2, y2_span[0], y2_span[1])[0] - # print(f'calculated a2 as {a2}') - # G2_fixed = lambda y2: a2 * G2(y2) - # a1 = G1_int / sp.integrate.quad(E, x_span[0], x_span[1])[0] - # E_fixed = lambda x: a1 * E(x) - - m1_prime = lambda x, m: pm1 * E_fixed(x) / G1_fixed(m) - y1_0 = y1_span[1] - m1_solved = sp.integrate.solve_ivp(m1_prime, x_span, [y1_0], dense_output=1) - m1 = lambda x: m1_solved.sol(x)[0] + def G1_fixed(y1, y1_span=y1_SPAN): return G1(y1) / integrate(G1, y1_span[0], y1_span[1]) + def E_fixed(x, x_span=x_SPAN): return E(x) / integrate(x, x_span[0], x_span[1]) + def m1_prime(x, m, pm1=PM1): return pm1 * E_fixed(x) / G1_fixed(m) + + y1_0 = y1_SPAN[1] + m1_solved = sp.integrate.solve_ivp(m1_prime, x_SPAN, [y1_0], dense_output=1) + def m1(x): return m1_solved.sol(x)[0] + # m1 = lambda x: m1_solved.sol(x)[0] # absue that solution was calculated already by importing mirror points reflectors = np.load("reflectors.npz") @@ -79,7 +79,7 @@ def intersection(l1, l2): # spline_B = sp.interpolate.CubicSpline(np.flipud(B[:][0]), np.flipud(B[:][1])) # spline_B = sp.interpolate.CubicSpline(np.flipud(B[:][0]), B[:][1]) - # x_arr_u = np.linspace(x_span[0], x_span[1], 1000) + # x_arr_u = np.linspace(x_SPAN[0], x_SPAN[1], 1000) # x_arr_w = np.linspace(B[0][0], B[-1][0], 1000) # plt.plot(x_arr_u, spline_A(x_arr_u), "r") # plt.plot(A[:, 0], A[:, 1], "k") @@ -161,7 +161,6 @@ def intersection(l1, l2): # plt.plot([xis, xis, xis + 10], [np.zeros(len(xis)), hits_U, hits_U - 10 * U_reflection]) # plt.show() - # Plotting # for i in range(n_rays): # plt.plot([Ps[i][0], As[i][0], Bs[i][0], Bs[i][0] + 4 * u2s[i][0]], @@ -181,7 +180,7 @@ def intersection(l1, l2): ray_low = (Bs[i][0], Bs[i][1]) ray_upp = (Bs[i][0] + 4 * u2s[i][0], Bs[i][1] + 4 * u2s[i][1]) ray = line(ray_low, ray_upp) - target_line = line((y1_span[0], L1), (y1_span[1], L1)) + target_line = line((y1_SPAN[0], L1), (y1_SPAN[1], L1)) intersec = intersection(target_line, ray) # print(intersec) try: diff --git a/test2.py b/test2.py index 873bcb2..7b3277c 100644 --- a/test2.py +++ b/test2.py @@ -8,10 +8,6 @@ import doptics.two_mirrors as tms from doptics.functions import cdf_sampling_source -# Enable ctrl+c -# import signal -# signal.signal(signal.SIGINT, signal.SIG_DFL) - if __name__ == '__main__': l1 = 5 l2 = 10 @@ -35,20 +31,18 @@ number_rays=15 ) - G1 = lambda y1: 1 - G2 = lambda y2: 1 - np.abs(y2 - 12) if 11 < y2 < 13 else 0 - E = lambda x: 1 / (np.exp(10 * (x - 0.5)) + np.exp(-10 * (x - 0.5))) + def G1(y1): return 1 + def G2(y2): return 1 - np.abs(y2 - 12) if 11 < y2 < 13 else 0 + def E(x): return 1 / (np.exp(10 * (x - 0.5)) + np.exp(-10 * (x - 0.5))) # new stuff # E = func.uniform angle_result = tms.solve_two_mirrors_parallel_source_point_target(starting_density=func.uniform, - target_distribution_1=lambda y1: 1, - target_distribution_2=lambda y2: 1 - - np.abs(y2 - 12) - if 11 < y2 < 13 else 0, + target_distribution_1=G1, + target_distribution_2=G2, x_span=x_span, y1_span=y1_span, y2_span=y2_span, u0=u0, w0=w0, l1=l1, l2=l2, # color='#a69f3f', - number_rays=16) \ No newline at end of file + number_rays=16)