From 72fa19f97b2894c9ec815a0a8bca9306bbb1204a Mon Sep 17 00:00:00 2001
From: BrokenVoodooDoll <44895553+BrokenVoodooDoll@users.noreply.github.com>
Date: Thu, 21 Jan 2021 17:09:59 +0300
Subject: [PATCH] Add files

Added scripts for visualizing paths of rays reflected from a spherical mirror
---
 spherical_collimated.py   | 56 ++++++++++++++++++++++++++++
 spherical_point_source.py | 77 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 133 insertions(+)
 create mode 100644 spherical_collimated.py
 create mode 100644 spherical_point_source.py

diff --git a/spherical_collimated.py b/spherical_collimated.py
new file mode 100644
index 0000000..fdec5af
--- /dev/null
+++ b/spherical_collimated.py
@@ -0,0 +1,56 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+radius = 100 # curvature radius of the mirror in mm
+angle_deg = 0 # angle of incidence of the incident beam in degrees
+rays = 21 # number of rays
+
+focal_length = radius / 2 # focal length of the mirror
+a = 1.1 * focal_length  # mirror field
+inc_ang = -angle_deg * np.pi / 180 if angle_deg > 0.000001 else 0.000001 * np.pi / 180  # incident ray angle in radians
+var = np.arange(-a, a, 0.1)
+
+# mirror equation
+def surface(y):
+    return np.sqrt(radius ** 2 - y ** 2) - radius
+
+# reflection angle
+def refl_ang(y, inc_ang):
+    return inc_ang - 2 * np.arcsin((-y / np.tan(inc_ang) - surface(y)) / radius * np.sin(inc_ang))
+
+# incident ray vector (y_start, y_end)
+# x_vec is vector (x_start, x_end)
+def inc_vec(y, inc_ang, x_vec):
+    return np.tan(-inc_ang) * (x_vec - surface(y)) + y
+
+# reflected ray vector (y_start, y_end)
+# x_vec is vector (x_start, x_end)
+def refl_vec(y, inc_ang, x_vec):
+    r = refl_ang(y, inc_ang)
+    return np.tan(r) * (x_vec - surface(y) + y / np.tan(r))
+
+
+plt.figure(figsize=(13, 8))
+plt.plot(surface(var), var) # mirror surface visualization
+plt.plot([-2 * radius, 0], [0, 0]) # axis of the mirror
+plt.plot([-focal_length], [0], 'o') # focal point
+
+for y in np.linspace(-focal_length, focal_length, rays):
+    x_vec = np.array([-radius, surface(y)])
+    
+    plt.plot(x_vec, inc_vec(y, inc_ang, x_vec), 'k', lw=1)
+
+    r = refl_ang(y, inc_ang)
+    if (r < np.pi / 2 and r > -np.pi / 2):
+        plt.plot(x_vec, refl_vec(y, inc_ang, x_vec), 'r', lw=1)
+    else:
+        x_vec_out = np.array([surface(y), 0])
+        plt.plot(x_vec_out, refl_vec(y, inc_ang, x_vec_out), 'r', lw=1)
+
+plt.title("Radius = {:.1f} mm. Focal length = {:.1f} mm. Incident angle = {:.1f} deg. Number of rays = {}".format(radius, focal_length, angle_deg, rays))
+plt.xlabel("z, mm")
+plt.ylabel("r, mm")
+plt.ylim(-a, a)
+plt.xlim(-radius, 0)
+plt.grid()
+plt.show()
\ No newline at end of file
diff --git a/spherical_point_source.py b/spherical_point_source.py
new file mode 100644
index 0000000..aaf059e
--- /dev/null
+++ b/spherical_point_source.py
@@ -0,0 +1,77 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from numpy.lib.function_base import angle
+
+radius = 100 # curvature radius of the mirror in mm (must be positive)
+angle_d = 30 # maximum angle of incidence of the incident beam in degrees
+num_rays = 21 # number of rays
+source_pos = 80 # source position in mm (must be positive)
+
+focal_length = radius / 2 # focal length of the mirror
+y = np.linspace(-radius, radius, 1000)
+
+# mirror equation z = sqrt(R^2 - y^2) - R
+def surface(y):
+    return np.sqrt(radius ** 2 - y ** 2) - radius
+
+# angle between the incident ray and the line connecting the point of incidence
+# of the ray on the mirror and the center of curvature of the mirror
+def epsilon(inc_angle):
+    q = radius - source_pos
+    return np.arcsin(q / radius * np.sin(inc_angle))
+
+# angle of reflected ray
+def ref_angle(inc_angle):
+    return inc_angle - 2 * epsilon(inc_angle)
+
+# the z-coordinate of the intersection of the reflected ray with the axis
+def ref_z(inc_angle):
+    q = radius * np.sin(-epsilon(inc_angle)) / np.sin(ref_angle(inc_angle))
+    return radius - q
+
+# the y-coordinate of the intersection of the incident ray with the mirror 
+def height(inc_angle):
+    phi = ref_angle(inc_angle) + epsilon(inc_angle)
+    return radius * np.sin(phi)
+
+# line equation for extension of the reflected ray
+def line(inc_angle, z, z0):
+    return np.tan(inc_angle) * (z - z0)
+
+
+plt.figure(figsize=(13, 8))
+plt.plot(surface(y), y) # mirror surface visualization
+plt.plot([-2 * radius, 0], [0, 0]) # axis of the mirror
+plt.plot([-focal_length], [0], 'o') # focal point
+
+for ang in np.linspace(-angle_d, angle_d, num_rays):
+    inc_angle = ang * np.pi / 180
+    h = height(inc_angle)
+
+    z_inc = np.array([-source_pos, surface(h)])
+    y_inc = np.array([0, h])
+    plt.plot(z_inc, y_inc, 'k', lw=1) # draw incident beam
+
+    z_0 = ref_z(inc_angle)
+    if np.isnan(z_0):
+        z_0 = -2 * radius
+
+    if source_pos >= focal_length:
+        z_0 = -z_0 if z_0 > 0 else z_0
+    else:
+        z_0 = z_0 if z_0 > 0 else -z_0
+        
+    z_ref = np.array([surface(h), -2 * radius])
+    y_ref = np.array([h, line(ref_angle(inc_angle), -2 * radius, z_0)])   
+    if abs(source_pos) < abs(2 * focal_length) and abs(source_pos) > abs(focal_length) and abs(z_0) > abs(2 * radius):
+        z_ref = np.array([surface(h), z_0])
+        y_ref = np.array([h, 0])
+    plt.plot(z_ref, y_ref, 'r', lw=1)
+
+plt.title("Radius = {:.1f} mm. Focal length = {:.1f} mm. Source position = {:.1f} mm.\nMaximum incident angle = {:.1f} deg. Number of rays = {}".format(radius, focal_length, -source_pos, angle_d, num_rays))
+plt.xlabel("z, mm")
+plt.ylabel("r, mm")
+plt.ylim(-radius, radius)
+plt.xlim(-2 * radius, 0)
+plt.grid()
+plt.show()
\ No newline at end of file