From 6f332fcfdef0ed6de30cebe739c4cbb669dbd94e Mon Sep 17 00:00:00 2001
From: BrokenVoodooDoll <44895553+BrokenVoodooDoll@users.noreply.github.com>
Date: Wed, 27 Jan 2021 18:13:16 +0300
Subject: [PATCH] Added MATLAB/Octave versions of the Python scripts

---
 parabolic_collimated.m   |  68 ++++++++++++++++++++++++++
 parabolic_point_source.m | 102 +++++++++++++++++++++++++++++++++++++++
 spherical_collimated.m   |  64 ++++++++++++++++++++++++
 spherical_point_source.m |  90 ++++++++++++++++++++++++++++++++++
 4 files changed, 324 insertions(+)
 create mode 100644 parabolic_collimated.m
 create mode 100644 parabolic_point_source.m
 create mode 100644 spherical_collimated.m
 create mode 100644 spherical_point_source.m

diff --git a/parabolic_collimated.m b/parabolic_collimated.m
new file mode 100644
index 0000000..61dc476
--- /dev/null
+++ b/parabolic_collimated.m
@@ -0,0 +1,68 @@
+clc
+clear
+
+focal_length = 100; % focal length in mm
+angle_deg = 0; % angle of incidence of the incident beam in degrees
+rays = 21; % number of rays
+
+p = 2 * focal_length; % parameter of the parabola equation y**2 = 2*p*z
+a = 1.1 * focal_length;  % mirror field
+inc_ang = -angle_deg * pi / 180;
+if angle_deg < 0.000001
+  inc_ang = 0.000001 * pi / 180;  % incident ray angle in radians
+end
+var = -a:0.1:a;
+
+% mirror equation z = -y^2 / (2 * p)
+function s = surface(y, p)
+    s = -y .^ 2 / (2 * p);
+end
+
+% reflection angle
+function angle = refl_ang(y, inc_ang, p)
+    angle = 2 * atan(y / p) - inc_ang;
+end
+
+
+% incident ray vector (y_start, y_end)
+% x_vec is vector (x_start, x_end)
+function v = inc_vec(y, inc_ang, x_vec, p)
+    v = tan(-inc_ang) * (x_vec - surface(y, p)) + y;
+end
+
+
+% reflected ray vector (y_start, y_end)
+% x_vec is vector (x_start, x_end)
+function v = refl_vec(y, inc_ang, x_vec, p)
+    sigma = refl_ang(y, inc_ang, p);
+    v = tan(sigma) * (x_vec - surface(y, p) + y / tan(sigma));
+end
+
+
+figure
+hold on
+plot(surface(var, p), var) % mirror surface visualization
+plot([-p 0], [0 0]) % axis of the mirror
+plot([-focal_length], [0], 'o') % focal point
+
+y = linspace(-focal_length, focal_length, rays);
+for i = 1:length(y)
+    x_vec = [-p surface(y(i), p)];
+    
+    plot(x_vec, inc_vec(y(i), inc_ang, x_vec, p), 'k')
+
+    r = refl_ang(y, inc_ang, p);
+    if (r < pi / 2 & r > -pi / 2)
+        plot(x_vec, refl_vec(y(i), inc_ang, x_vec, p), 'r')
+    else
+        x_vec_out = [surface(y(i), p) 0];
+        plot(x_vec_out, refl_vec(y(i), inc_ang, x_vec_out, p), 'r')
+    end
+end
+
+title(sprintf("Focal length = %.1f mm. Incident angle = %.1f{\\deg}. Number of rays = %d", focal_length, angle_deg, rays))
+xlabel("z, mm")
+ylabel("r, mm")
+ylim([-a a])
+xlim([-p 0])
+grid on
\ No newline at end of file
diff --git a/parabolic_point_source.m b/parabolic_point_source.m
new file mode 100644
index 0000000..5b1a40d
--- /dev/null
+++ b/parabolic_point_source.m
@@ -0,0 +1,102 @@
+clc
+clear
+
+focal_length = 100; % focal length in mm
+angle_d = 25; % maximum angle of incidence of the incident beam in degrees
+num_rays = 21; % number of rays
+source_pos = 400; % source position in mm (must be positive)
+% there is a bug when source_pos = focal_length because of very small
+% angles of reflectance.
+% Instead, better to choose source_pos = focal_length + 0.01 or something like
+% this
+
+p = 2 * focal_length; % parameter of the parabola equation y**2 = 2*p*z
+y = linspace(-p, p, 1000);
+
+% mirror equation z = -y^2 / (2 * p)
+function s = surface(y, p)
+    s = -y .^ 2 / (2 * p);
+end
+
+function angle = phi(y, p)
+    angle = atan(y / p);
+end
+
+% 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
+function angle = epsilon(y, inc_angle, p)
+    angle = inc_angle - phi(y, p);
+end
+
+% angle of reflected ray
+function angle = ref_angle(y, inc_angle, p)
+    angle = phi(y, p) - epsilon(y, inc_angle, p);
+end
+
+% the z-coordinate of the intersection of the reflected ray with the axis
+function z = ref_z(y, inc_angle, p, s)
+    if inc_angle == 0
+        z = 0;
+        return
+    end
+    tg_sigma = tan(inc_angle);
+    z_0 = -(y - s * tg_sigma) / tg_sigma;
+    z = y ./ tan(ref_angle(y, inc_angle, p)) + z_0;
+end
+
+% the y-coordinate of the intersection of the incident ray with the mirror 
+function h = height(inc_angle, p, s)
+    if inc_angle == 0
+        h = 0;
+        return
+    end
+    tg_sigma = tan(inc_angle);
+    h = -p / tg_sigma * (1 - sqrt(1 + 2 * s / p * tg_sigma ^ 2));
+end
+
+% line equation for extension of the reflected ray
+function l = line(ref_angle, z, z0)
+    l = tan(ref_angle) * (z - z0);
+end
+
+figure
+hold on
+plot(surface(y, p), y) % mirror surface visualization
+plot([-2 * p 0], [0 0]) % axis of the mirror
+plot([-focal_length], [0], 'o') % focal point
+
+angles = linspace(-angle_d, angle_d, num_rays);
+for i  = 1:length(angles)
+    inc_angle = angles(i) * pi / 180;
+    h = height(inc_angle, p, source_pos);
+
+    z_inc = [-source_pos surface(h, p)];
+    y_inc = [0 h];
+    plot(z_inc, y_inc, 'k') % draw incident beam
+
+    z_0 = ref_z(h, inc_angle, p, source_pos);
+    if isnan(z_0)
+        z_0 = -2 * p;
+    end
+
+    if source_pos >= focal_length
+        if z_0 > 0
+          z_0 = -z_0;
+        end
+    else
+        if z_0 < 0
+          z_0 = -z_0;
+        end
+    end
+    
+    z_ref = [surface(h, p) -2 * p];
+    y_ref = [h line(ref_angle(h, inc_angle, p), -2 * p, z_0)];
+    plot(z_ref, y_ref, 'r')
+end
+
+title(sprintf("Focal length = %.1f mm. Source position = %.1f mm.\nMaximum incident angle = %.1f{\\deg}. Number of rays = %d", focal_length, -source_pos, angle_d, num_rays))
+xlabel("z, mm")
+ylabel("y, mm")
+ylim([-p p])
+xlim([-2 * p 0])
+grid on
\ No newline at end of file
diff --git a/spherical_collimated.m b/spherical_collimated.m
new file mode 100644
index 0000000..fe5450d
--- /dev/null
+++ b/spherical_collimated.m
@@ -0,0 +1,64 @@
+clc
+clear
+
+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 * pi / 180;
+if angle_deg < 0.000001
+  inc_ang = 0.000001 * pi / 180;  % incident ray angle in radians
+end
+var = -a:0.1:a;
+
+% mirror equation
+function s = surface(y, r)
+    s = sqrt(r ^ 2 - y .^ 2) - r;
+end
+
+% reflection angle
+function angle = refl_ang(y, inc_ang, r)
+    angle = inc_ang - 2 * asin((-y / tan(inc_ang) - surface(y, r)) / r * sin(inc_ang));
+end
+
+% incident ray vector (y_start, y_end)
+% x_vec is vector (x_start, x_end)
+function v = inc_vec(y, inc_ang, x_vec, r)
+    v = tan(-inc_ang) * (x_vec - surface(y, r)) + y;
+end
+
+% reflected ray vector (y_start, y_end)
+% x_vec is vector (x_start, x_end)
+function v = refl_vec(y, inc_ang, x_vec, r)
+    sigma =  refl_ang(y, inc_ang, r);
+    v = tan(sigma) .* (x_vec - surface(y, r) + y ./ tan(sigma));
+end
+
+figure
+hold on
+plot(surface(var, radius), var) % mirror surface visualization
+plot([-2 * radius 0], [0 0]) % axis of the mirror
+plot([-focal_length], [0], 'o') % focal point
+
+y = linspace(-focal_length, focal_length, rays);
+for i = 1:length(y)
+    x_vec = [-radius surface(y(i), radius)];
+    plot(x_vec, inc_vec(y(i), inc_ang, x_vec, radius), 'k')
+    
+    r = refl_ang(y(i), inc_ang, radius);
+    if (r < pi / 2 & r > -pi / 2)
+        plot(x_vec, refl_vec(y(i), inc_ang, x_vec, radius), 'r')
+    else
+        x_vec_out = [surface(y(i)) 0];
+        plot(x_vec_out, refl_vec(y(i), inc_ang, x_vec_out, radius), 'r')
+    end
+end
+
+title(sprintf("Radius = %.1f mm. Focal length = %.1f mm.\nIncident angle = %.1f{\\deg}. Number of rays = %d", radius, focal_length, angle_deg, rays))
+xlabel("z, mm")
+ylabel("r, mm")
+ylim([-a a])
+xlim([-radius 0])
+grid
\ No newline at end of file
diff --git a/spherical_point_source.m b/spherical_point_source.m
new file mode 100644
index 0000000..d6a2896
--- /dev/null
+++ b/spherical_point_source.m
@@ -0,0 +1,90 @@
+clc
+clear
+
+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 = linspace(-radius, radius, 1000);
+
+% mirror equation z = sqrt(R^2 - y^2) - R
+function s = surface(y, r)
+    s = sqrt(r ^ 2 - y .^ 2) - r;
+end
+
+% 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
+function e = epsilon(inc_angle, r, s)
+    q = r - s;
+    e = asin(q / r * sin(inc_angle));
+end
+
+% angle of reflected ray
+function r = ref_angle(inc_angle, r, s)
+    r = inc_angle - 2 * epsilon(inc_angle, r, s);
+end
+
+% the z-coordinate of the intersection of the reflected ray with the axis
+function z = ref_z(inc_angle, r, s)
+    q = r * sin(-epsilon(inc_angle, r, s)) / sin(ref_angle(inc_angle, r, s));
+    z = r - q;
+end
+
+% the y-coordinate of the intersection of the incident ray with the mirror 
+function h = height(inc_angle, r, s)
+    phi = ref_angle(inc_angle, r, s) + epsilon(inc_angle, r, s);
+    h = r * sin(phi);
+end    
+
+% line equation for extension of the reflected ray
+function l = line(inc_angle, z, z0)
+    l = tan(inc_angle) * (z - z0);
+end    
+
+figure
+hold on
+plot(surface(y, radius), y) % mirror surface visualization
+plot([-2 * radius 0], [0 0]) % axis of the mirror
+plot([-focal_length], [0], 'o') % focal point
+
+angles = linspace(-angle_d, angle_d, num_rays);
+for i = 1:length(angles)
+    inc_angle = angles(i) * pi / 180;
+    h = height(inc_angle, radius, source_pos);
+
+    z_inc = [-source_pos surface(h, radius)];
+    y_inc = [0 h];
+    plot(z_inc, y_inc, 'k') % draw incident beam
+
+    z_0 = ref_z(inc_angle, radius, source_pos);
+    if isnan(z_0)
+        z_0 = -2 * radius;
+    end 
+
+    if source_pos >= focal_length
+        if z_0 > 0
+            z_0 = -z_0; 
+        end
+    else
+        if z_0 < 0
+            z_0 = -z_0; 
+        end
+    end    
+        
+    z_ref = [surface(h, radius) -2 * radius];
+    y_ref = [h line(ref_angle(inc_angle, radius, source_pos), -2 * radius, z_0)];
+    if abs(source_pos) < abs(2 * focal_length) & abs(source_pos) > abs(focal_length) & abs(z_0) > abs(2 * radius)
+        z_ref = [surface(h, radius) z_0]
+        y_ref = [h 0]
+    end
+    plot(z_ref, y_ref, 'r')
+end
+
+title(sprintf("Radius = %.1f mm. Focal length = %.1f mm. Source position = %.1f mm.\nMaximum incident angle = %.1f{\\deg}. Number of rays = %d", radius, focal_length, -source_pos, angle_d, num_rays))
+xlabel("z, mm")
+ylabel("r, mm")
+ylim([-radius radius])
+xlim([-2 * radius 0])
+grid on
\ No newline at end of file