This repository has been archived on 2023-03-25. You can view files and clone it, but cannot push or open issues or pull requests.
mightyscape-1.1-deprecated/extensions/fablabchemnitz/raytracing/lens.py
2021-10-24 19:19:52 +02:00

229 lines
7.8 KiB
Python

"""
Module to add a lens object in the document
"""
from math import cos, pi, sin, sqrt, acos, tan
import inkex
class Lens(inkex.GenerateExtension):
"""
Produces a PathElement corresponding to the shape of the lens calculated
from user parameters.
"""
@property
def style(self):
return {
"stroke": "#000000",
"fill": "#b7c2dd",
"stroke-linejoin": "round",
"stroke-width": str(self.svg.unittouu("1px")),
}
def add_arguments(self, pars):
pars.add_argument("--focal_length", type=float, default=100.0)
pars.add_argument("--focal_length_unit", default="mm")
pars.add_argument("--diameter", type=float, default=1.0)
pars.add_argument("--diameter_unit", default="in")
pars.add_argument("--edge_thickness", type=float, default=2.0)
pars.add_argument("--edge_thickness_unit", default="mm")
pars.add_argument("--optical_index", type=float, default=1.5168)
pars.add_argument("--lens_type", default="plano_con")
def to_document_units(self, value: float, unit: str) -> float:
c1x, c1y, c2x, c2y = self.svg.get_viewbox()
document_width = self.svg.unittouu(self.document.getroot().get("width"))
scale_factor = (c2x - c1x) / document_width
return self.svg.unittouu(str(value) + unit) * scale_factor
def generate(self):
opts = self.options
d = self.to_document_units(opts.diameter, opts.diameter_unit)
f = self.to_document_units(opts.focal_length, opts.focal_length_unit)
e = self.to_document_units(opts.edge_thickness, opts.edge_thickness_unit)
optical_index = opts.optical_index
lens_path = []
if opts.lens_type == "plano_con":
# Radius of curvature from Lensmaker's equation
roc = (optical_index - 1) * abs(f)
if 2 * roc < d:
inkex.utils.errormsg(
"Focal length is too short or diameter is too large."
)
return None
elif (roc ** 2 - (d / 2) ** 2) ** 0.5 - roc < -e and f < 0:
inkex.utils.errormsg("Edge thickness is too small.")
return None
else:
sweep = 1 if f < 0 else 0
lens_path = arc_to_path(
[-d / 2, 0], [roc, roc, 0.0, 0, sweep, +d / 2, 0]
)
lens_path += [
[[+d / 2, 0], [+d / 2, 0], [+d / 2, -e]],
[[+d / 2, -e], [+d / 2, -e], [-d / 2, -e]],
[[+d / 2, -e], [-d / 2, -e], [-d / 2, +e]],
]
# no need to close the path correctly as it's done after
elif opts.lens_type == "bi_con":
roc = (
(optical_index - 1) * abs(f) * (1 + (1 - e / f / optical_index) ** 0.5)
)
if 2 * roc < d:
inkex.utils.errormsg(
"Focal length is too short or diameter is too large."
)
return None
elif (roc ** 2 - (d / 2) ** 2) ** 0.5 - roc < -e / 2 and f < 0:
inkex.utils.errormsg("Edge thickness is too small.")
return None
else:
sweep = 1 if f < 0 else 0
lens_path = arc_to_path(
[-d / 2, 0], [roc, roc, 0.0, 0, sweep, +d / 2, 0]
)
lens_path += [
[[+d / 2, 0], [+d / 2, 0], [+d / 2, -e]],
[[+d / 2, -e], [+d / 2, -e], [+d / 2, -e]],
]
lens_path += arc_to_path(
[+d / 2, -e], [roc, roc, 0.0, 0, sweep, -d / 2, -e]
)
lens_path += [
[[-d / 2, -e], [-d / 2, -e], [-d / 2, 0]],
[[-d / 2, -e], [-d / 2, 0], [-d / 2, 0]],
]
lens = inkex.PathElement()
lens.style = self.style
closed_path = inkex.Path(inkex.CubicSuperPath([lens_path]))
closed_path.close()
lens.path = closed_path.transform(inkex.Transform("rotate(90)"))
lens.desc = (
f"L{opts.focal_length}{opts.focal_length_unit}\n"
f"optics:glass:{optical_index:.4f}"
)
yield lens
def arc_to_path(point, params):
"""Approximates an arc with cubic bezier segments.
Arguments:
point: Starting point (absolute coords)
params: Arcs parameters as per
https://www.w3.org/TR/SVG/paths.html#PathDataEllipticalArcCommands
Returns a list of triplets of points : [control_point_before, node, control_point_after]
(first and last returned triplets are [p1, p1, *] and [*, p2, p2])
"""
A = point[:]
rx, ry, theta, long_flag, sweep_flag, x2, y2 = params[:]
theta = theta * pi / 180.0
B = [x2, y2]
# Degenerate ellipse
if rx == 0 or ry == 0 or A == B:
return [[A[:], A[:], A[:]], [B[:], B[:], B[:]]]
# turn coordinates so that the ellipse morph into a *unit circle* (not 0-centered)
mat = mat_prod(
(rot_mat(theta), [[1.0 / rx, 0.0], [0.0, 1.0 / ry]], rot_mat(-theta))
)
apply_mat(mat, A)
apply_mat(mat, B)
k = [-(B[1] - A[1]), B[0] - A[0]]
d = k[0] * k[0] + k[1] * k[1]
k[0] /= sqrt(d)
k[1] /= sqrt(d)
d = sqrt(max(0, 1 - d / 4.0))
# k is the unit normal to AB vector, pointing to center O
# d is distance from center to AB segment (distance from O to the midpoint of AB)
# for the last line, remember this is a unit circle, and kd vector is orthogonal to AB (Pythagorean thm)
if (
long_flag == sweep_flag
): # top-right ellipse in SVG example https://www.w3.org/TR/SVG/images/paths/arcs02.svg
d *= -1
O = [(B[0] + A[0]) / 2.0 + d * k[0], (B[1] + A[1]) / 2.0 + d * k[1]]
OA = [A[0] - O[0], A[1] - O[1]]
OB = [B[0] - O[0], B[1] - O[1]]
start = acos(OA[0] / norm(OA))
if OA[1] < 0:
start *= -1
end = acos(OB[0] / norm(OB))
if OB[1] < 0:
end *= -1
# start and end are the angles from center of the circle to A and to B respectively
if sweep_flag and start > end:
end += 2 * pi
if (not sweep_flag) and start < end:
end -= 2 * pi
nb_sectors = int(abs(start - end) * 2 / pi) + 1
d_theta = (end - start) / nb_sectors
v = 4 * tan(d_theta / 4.0) / 3.0
# I would use v = tan(d_theta/2)*4*(sqrt(2)-1)/3 ?
p = []
for i in range(0, nb_sectors + 1, 1):
angle = start + i * d_theta
v1 = [
O[0] + cos(angle) - (-v) * sin(angle),
O[1] + sin(angle) + (-v) * cos(angle),
]
pt = [O[0] + cos(angle), O[1] + sin(angle)]
v2 = [O[0] + cos(angle) - v * sin(angle), O[1] + sin(angle) + v * cos(angle)]
p.append([v1, pt, v2])
p[0][0] = p[0][1][:]
p[-1][2] = p[-1][1][:]
# go back to the original coordinate system
mat = mat_prod((rot_mat(theta), [[rx, 0], [0, ry]], rot_mat(-theta)))
for pts in p:
apply_mat(mat, pts[0])
apply_mat(mat, pts[1])
apply_mat(mat, pts[2])
return p
def mat_prod(m_list):
"""Get the product of the mat"""
prod = m_list[0]
for mat in m_list[1:]:
a00 = prod[0][0] * mat[0][0] + prod[0][1] * mat[1][0]
a01 = prod[0][0] * mat[0][1] + prod[0][1] * mat[1][1]
a10 = prod[1][0] * mat[0][0] + prod[1][1] * mat[1][0]
a11 = prod[1][0] * mat[0][1] + prod[1][1] * mat[1][1]
prod = [[a00, a01], [a10, a11]]
return prod
def rot_mat(theta):
"""Rotate the mat"""
return [[cos(theta), -sin(theta)], [sin(theta), cos(theta)]]
def apply_mat(mat, point):
"""Apply the given mat"""
x = mat[0][0] * point[0] + mat[0][1] * point[1]
y = mat[1][0] * point[0] + mat[1][1] * point[1]
point[0] = x
point[1] = y
def norm(point):
"""Normalise"""
return sqrt(point[0] * point[0] + point[1] * point[1])
if __name__ == "__main__":
Lens().run()