In [1]:
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
# %matplotlib widget
In [2]:
def get_refraction_points(p1, radius, center1, center2, pre_refractive_index, lens_refractive_index, post_refractive_index, ray_length = 8):
length = p1[1]
refractive_index_ratio1 = pre_refractive_index/lens_refractive_index
refractive_index_ratio2 = lens_refractive_index/post_refractive_index
theta1 = np.arcsin(length/r)
x1 = center1[0] + radius*np.cos(np.pi - theta1)
y1 = center1[1] + radius*np.sin(np.pi - theta1)
A = np.array([x1, y1])
# ray_length = np.linalg.norm((c2[0] - A[0], c2[1] - A[1]))
theta2 = np.arcsin(refractive_index_ratio1 * np.sin(theta1))
x2 = A[0] + ray_length*np.cos(-theta1 + theta2)
y2 = A[1] + ray_length*np.sin(-theta1 + theta2)
B = np.array([x2, y2])
from sympy import symbols, Eq, solve, sin, cos, pi
x, y, a = symbols('x, y, a')
m = (A[1] - B[1])/(A[0] - B[0])
b = A[1] - m*A[0]
eq1 = Eq(y, m*x + b)
eq2 = Eq(x, center2[0]+radius*cos(a))
eq3 = Eq(y, center2[1]+radius*sin(a))
soln = solve((eq1, eq2, eq3))
alpha = 0
# soln = np.array(soln).astype(np.float64)
for items in soln:
if items[a] <= theta or items[a] >= 2*pi-theta:
# print(items[a], items[x], items[y])
B = np.array([items[x], items[y]]).astype(np.float64)
alpha = np.float64(items[a])
theta3 = np.arcsin(refractive_index_ratio2* np.sin(theta1 - theta2 + alpha) ) - alpha
x3 = B[0] + ray_length*np.cos(-theta3)
y3 = B[1] + ray_length*np.sin(-theta3)
C = np.array([x3, y3])
# print(theta3 , (theta1 - theta2 + alpha))
return A, B, C
In [18]:
# plotter = Plotter()
c1 = np.array([10, 0])
c2 = np.array([2, 0])
r = 5
length = 2
ray_length = 15
theta = np.arccos(((c1[0]-c2[0])/2)/r)
obj_x = 0
obj_y = 0
pu0 = np.array([obj_x, obj_y])
pu1 = np.array([obj_x, obj_y + length])
pre_refractive_index = 1
lens_refractive_index = 1.8
post_refractive_index = 1
list_of_values = []
no_of_rays = 4
for i in np.linspace(0, length, no_of_rays):
A, B, C = get_refraction_points(p1 = np.array([obj_x, obj_y + i]),
radius = r,
center1 = c1,
center2 = c2,
pre_refractive_index = pre_refractive_index,
lens_refractive_index = lens_refractive_index,
post_refractive_index = post_refractive_index,
ray_length = ray_length,
)
list_of_values.append([A, B, C])
In [19]:
fig = plt.figure(figsize = (9, 9))
plt.plot(*c1, color = 'b', marker = 'o')
plt.plot(*c2, color = 'b', marker = 'o')
# plt.plot(*pu0, color = 'b', marker = 'o')
# plt.plot(pu0[0] + pu0[1]+length, color = 'r', marker = 'o')
angle_circle1 = np.linspace(np.pi - theta, np.pi + theta, 50)
circle1_x = r*np.cos(angle_circle1) + c1[0]
circle1_y = r*np.sin(angle_circle1) + c1[1]
plt.plot(circle1_x, circle1_y, color = 'black')
angle_circle2 = np.linspace(- theta, theta, 50)
circle2_x = r*np.cos(angle_circle2) + c2[0]
circle2_y = r*np.sin(angle_circle2) + c2[1]
plt.plot(circle2_x, circle2_y, color = 'black')
plt.arrow(*pu0, 0, length, width = 0.1, color = 'blue', length_includes_head = True)
for n, i in enumerate(np.linspace(0, length, no_of_rays)):
plt.plot([obj_x, list_of_values[n][0][0]], [obj_y + i, list_of_values[n][0][1]], color = 'r', alpha = 0.5)
plt.plot([list_of_values[n][0][0], list_of_values[n][1][0]], [list_of_values[n][0][1], list_of_values[n][1][1]], color = 'r', alpha = 0.5)
plt.plot([list_of_values[n][1][0], list_of_values[n][2][0]], [list_of_values[n][1][1], list_of_values[n][2][1]], color = 'r', alpha = 0.5)
plt.gca().set_aspect('equal')
In [ ]: