Link here: https://projecteuler.net/problem=144.
Algorithm:
Modelling:
from collections import namedtuple
Ray = namedtuple('Ray', 'point, direction')
r = Ray((0.1, 0.2), (0.2, 0.4))
r.point
(0.1, 0.2)
r.direction
(0.2, 0.4)
def iterate(ray, on_ellipse=True):
"Computes the next state of the ray."
new_point = compute_intersection(ray, on_ellipse=on_ellipse)
normal = compute_normal(new_point)
new_direction = compute_reflexion(ray, normal)
return Ray(new_point, new_direction)
#TODO: compute_intersection, compute_normal, compute_reflexion
import sympy
x, y = sympy.symbols('x, y')
e = 4 * x**2 + y ** 2 - 100
e
4*x**2 + y**2 - 100
Let's compute the first direction:
import numpy as np
d = np.array([1.4, -9.6]) - np.array([0, 10.1])
d
array([ 1.4, -19.7])
p = np.array([0, 10.1])
l = (x - p[0])*(-d[1]) + (y - p[1])*d[0]
Let's now solve the system of equations:
solutions = sympy.solve([e, l])
solutions
[{y: 9.99998989720405, x: 0.00710731694996590}, {y: -9.60000000000000, x: 1.40000000000000}]
How can we choose the right solution? Easy, just use the one that has a positive scalar product.
np.sum(d * np.array([0.00710731694996590, 9.99998989720405]))
-196.98985073118985
np.sum(d * np.array([1.40000000000000, -9.60000000000000]))
191.07999999999998
solutions[0][x]
0.00710731694996590
solutions[0][y]
9.99998989720405
sympy.lambdify()
def compute_intersection(ray, on_ellipse=True):
"Computes intersection point between ray and ellipse, assuming ray starts on ellipse."
x, y = sympy.symbols('x, y')
e = 4 * x**2 + y ** 2 - 100
p = ray.point
d = ray.direction
l = (x - p[0])*(-d[1]) + (y - p[1])*d[0]
solutions = sympy.solve([e, l])
if on_ellipse:
for solution in solutions:
dist = np.linalg.norm(np.array([solution[x], solution[y]]) - np.array(p))
if dist > 1e-1:
return (solution[x], solution[y])
else:
for solution in solutions:
scalarprod = np.sum(np.array(d) * np.array([solution[x], solution[y]]))
if scalarprod > 0:
return (solution[x], solution[y])
raise ValueError('no intersection found!')
r = Ray(p, d)
compute_intersection(r)
(0.00710731694996590, 9.99998989720405)
compute_intersection(r, on_ellipse=False)
(1.40000000000000, -9.60000000000000)
We already know the slope of the tangent to the ellipse:
$$ m = {-4x \over y} $$def compute_normal(point):
"Returns normal vector to ellipse at given point."
x, y = point
n = -np.array((x / 25., y / 100.))
n /= np.linalg.norm(n)
return n.astype(np.float64)
compute_normal((0, 10))
array([-0., -1.])
Let's plot some normals here:
import matplotlib.pyplot as plt
%matplotlib notebook
from matplotlib.patches import Ellipse
plt.figure()
plt.gca().add_patch(Ellipse((0, 0), 10, 20, alpha=0.1))
plt.xlim(-12, 12)
plt.ylim(-12, 12)
def plot_normal(pt):
x, y = pt
n = compute_normal(pt)
plt.arrow(x, y, n[0], n[1], head_width=0.05, head_length=0.1)
plot_normal((0, 10))
plot_normal((0, -10))
plot_normal((5, 0))
plot_normal((-5, 0))
compute_intersection(r)
(1.40000000000000, -9.60000000000000)
compute_normal(compute_intersection(r))
array([-0.00284292, -0.99999596])
def compute_reflexion(ray, normal):
"Returns the new direction vector."
i = np.array(ray.direction)
n = np.array(normal)
i -= 2 * np.sum(i * n) * n
return (i[0], i[1])
Let's test this:
compute_reflexion(r, compute_normal(compute_intersection(r)))
(1.5119878928830672, 19.691721423272657)
r
Ray(point=array([ 0. , 10.1]), direction=array([ 1.4, -19.7]))
iterate(r)
Ray(point=(0.00710731694996590, 9.99998989720405), direction=(1.5119878928830672, 19.691721423272657))
def plot_ray(ray):
"Plots a ray."
x, y = ray.point[0], ray.point[1]
u, v = ray.direction[0], ray.direction[1]
plt.arrow(x, y , u, v, head_width=0.05, head_length=0.1)
plt.figure()
plt.gca().add_patch(Ellipse((0, 0), 10, 20, alpha=0.1))
plt.xlim(-12, 12)
plt.ylim(-12, 12)
rr = iterate(r, on_ellipse=False)
for i in range(10):
plot_ray(rr)
rr = iterate(rr)
Let's just plot the points of intersection as they get generated:
points = []
rr = iterate(r, on_ellipse=False)
for i in range(10):
rr = iterate(rr)
points.append(rr.point)
plt.figure()
plt.gca().add_patch(Ellipse((0, 0), 10, 20, alpha=0.1))
plt.xlim(-12, 12)
plt.ylim(-12, 12)
plt.plot(np.array(points)[:, 0], np.array(points)[:, 1])
[<matplotlib.lines.Line2D at 0xb7ee160>]
from ipywidgets import interact
%matplotlib inline
@interact
def plot_incremental_rays(index=(1, len(points) - 1)):
plt.figure()
plt.gca().add_patch(Ellipse((0, 0), 10, 20, alpha=0.1))
plt.xlim(-12, 12)
plt.ylim(-12, 12)
plt.plot(np.array(points)[:index, 0], np.array(points)[:index, 1], '-o')
def is_out(ray):
"Checks whether the ray has escaped the ellipse."
x, y = ray.point
return (-0.01 <= x <= 0.01) and y > 0
is_out(Ray((0.011, 0.1), (0, 0)))
False
r
Ray(point=array([ 0. , 10.1]), direction=array([ 1.4, -19.7]))
points = []
rr = iterate(r, on_ellipse=False)
points.append(rr.point)
while True:
rr = iterate(rr)
points.append(rr.point)
if is_out(rr):
break
print('.', end='')
.................................................................................................................................................................................................................................................................................................................................................................
print("solution is {}".format(len(points) - 1))
solution is 354
points[-1]
(-0.00984876093143555, 9.99998060036281)
%matplotlib notebook
plt.figure()
plt.plot(np.array(points)[:, 0], 'o')
plt.hlines([0.01, -0.01], 0, len(points))
<matplotlib.collections.LineCollection at 0x11b2c828>
plt.figure()
plt.plot(np.array(points)[:, 1])
[<matplotlib.lines.Line2D at 0x1144c048>]
interact(plot_incremental_rays,
index=(1, len(points) - 1))
points
[(1.40000000000000, -9.60000000000000), (-3.99059761936164, -6.02499149886379), (0.324582878083967, 9.97890694520291), (0.416490435780230, -9.96524675397526), (-4.47344349721357, 4.46690195873498), (1.13233786470162, 9.74018705367879), (0.0536506662265838, -9.99942430463143), (-1.83444237340102, 9.30264933847790), (3.28824950976850, 7.53323706291186), (-0.226491477766366, -9.98973505364366), (-0.566604327511133, 9.93558443898408), (4.90431042704293, -1.94703799161500), (-0.852425522672158, -9.85360253476810), (-0.118726109219474, 9.99718042469768), (2.37474089458685, -8.80013765428174), (-2.61551238787611, -8.52269791764712), (0.146551577819212, 9.99570360405684), (0.761081905005412, -9.88347192718699), (-4.98153079674417, -0.858722122899896), (0.636896896911528, 9.91854068756175), (0.193346746067536, -9.99252061009335), (-3.01616727994218, 7.97564666704993), (2.03349686334379, 9.13562050585968), (-0.0783865432841545, -9.99877103444858), (-1.01380874000134, 9.79228100877357), (4.67723035500788, 3.53469444574370), (-0.470867662547838, -9.95555797419044), (-0.283481011270524, 9.98391476651299), (3.71796643012843, -6.68632204562060), (-1.55836389076242, -9.50189496552510), (0.0165277996782151, 9.99994536621832), (1.34102998909082, -9.63361584626647), (-4.09588844700140, -5.73539809593735), (0.342025820357107, 9.97657623400113), (0.396277651636882, -9.96854332845339), (-4.38245889950231, 4.81416825388269), (1.18296241992587, 9.71608973055377), (0.0440005983245893, -9.99961278197252), (-1.75978004931405, 9.36016542119556), (3.39838127580374, 7.33512227689943), (-0.240440910613546, -9.98843094154499), (-0.540546115451592, 9.94139022412262), (4.85767643187498, -2.36894886665485), (-0.891498903134685, -9.83976213243179), (-0.108070164764551, 9.99766389502823), (2.28324356785164, -8.89647094298950), (-2.71553654711948, -8.39663295881353), (0.158114718975740, 9.99499869647680), (0.727244533834306, -9.89365764275440), (-4.99562966766892, -0.418015423189514), (0.666956461311776, 9.91063450616851), (0.180830009788430, -9.99345796159866), (-2.91036949781890, 8.13135889901936), (2.11741654361105, 9.05903906180832), (-0.0884886341801505, -9.99843382967971), (-0.969852325025090, 9.81007369343236), (4.74687176861991, 3.14146998220848), (-0.494071151321741, -9.95105897830610), (-0.268104970890267, 9.98561359648648), (3.60755272469153, -6.92403447090515), (-1.62568009484008, -9.45667261339648), (0.0259805442601294, 9.99986500135276), (1.28429990014651, -9.66448628049803), (-4.19842294972815, -5.43074386624741), (0.360130134227017, 9.97402752882137), (0.376828658409305, -9.97155958959308), (-4.28686157996732, 5.14696718240948), (1.23563680193946, 9.68983006944762), (0.0344364069660533, -9.99976282396243), (-1.68769922019544, 9.41311241665619), (3.50900316348113, 7.12373407664251), (-0.254857510375878, -9.98700108128662), (-0.515518313569446, 9.94670616201655), (4.80259877058855, -2.78211793333165), (-0.932201502026572, -9.82466291729531), (-0.0976249773269686, 9.99809369105969), (2.19433903802941, -8.98551638720451), (-2.81772592991015, -8.26085239764323), (0.169985883419199, 9.99421928905669), (0.694766497365913, -9.90298934951218), (-4.99998521047336, 0.0243224215598438), (0.698271455063897, 9.90200322662903), (0.168665330301781, -9.99430878177259), (-2.80639310294019, 8.27626914781527), (2.20398261805475, 8.97607054769792), (-0.0987633478281681, -9.99804896989923), (-0.927646161653755, 9.82638745394645), (4.80905074810866, 2.73717438400661), (-0.518220179848705, -9.94614454855701), (-0.253249256594594, 9.98716472559340), (3.49683261536947, -7.14763224014451), (-1.69549665746593, -9.40750574478099), (0.0354840005576638, 9.99974817397007), (1.22974493715681, -9.69282773797971), (-4.29757421533640, -5.11110785003630), (0.378930110539142, 9.97124104037743), (0.358106857584838, -9.97431891981617), (-4.18730715652254, 5.46496432813064), (1.29042762667322, 9.66122073866828), (0.0249394296559529, -9.99987560419597), (-1.61815982456346, 9.46183043225120), (3.61971182648870, 6.89860458155933), (-0.269769002777829, -9.98543432908960), (-0.491475347622385, 9.95157313849011), (4.73957207068452, -3.18525137895726), (-0.974596863552674, -9.80819268847289), (-0.0873702097139703, 9.99847317273083), (2.10806309748528, -9.06776046816870), (-2.92191269324341, -8.11478310568140), (0.182188093620622, 9.99335929476027), (0.663592115892426, -9.91153782290627), (-4.99455495818489, 0.466565192318248), (0.730896146376426, 9.89258122498109), (0.156829136626174, -9.99507971392026), (-2.70443404506700, 8.41095392827295), (2.29317419847386, 8.88624827369881), (-0.109230698260387, -9.99761344612954), (-0.887125720158343, 9.84134298896925), (4.86322345361817, 2.32297881193796), (-0.543358979139638, -9.94077683479281), (-0.238885344391391, 9.98858023789868), (3.38624560431714, -7.35753782382201), (-1.76785939086567, -9.35407358836310), (0.0450567152213289, 9.99959397023968), (1.17729912307893, -9.71884083104515), (-4.39269674693032, -4.77669981870654), (0.398461240287576, 9.96819515057535), (0.340076892365315, -9.97684272849465), (-4.08443943785959, 5.76796479824807), (1.34740076652591, 9.63005943374502), (0.0154911320973545, -9.99995200485009), (-1.55111540299715, 9.50663789287990), (3.73006331830606, 6.65931757507556), (-0.285204012930471, -9.98371847980668), (-0.468373150605604, 9.95602864435248), (4.66914565312826, -3.57719379954369), (-1.01874963053747, -9.79022965824189), (-0.0772858857194548, 9.99880530700914), (2.02443490527876, -9.14366738552731), (-3.02789728925990, -7.95784845437322), (0.194744990607333, 9.99241199883859), (0.633667296291042, -9.91936807616517), (-4.97939184263934, 0.906988153079881), (0.764886388221083, 9.88229706356050), (0.145298471717413, -9.99577677904355), (-2.60465682582633, 8.53598566532918), (2.38495353824062, 8.78908336982499), (-0.119911064678056, -9.99712385370268), (-0.848227585935610, 9.85505148894877), (4.90890191730165, 1.90019153383249), (-0.569533291134686, -9.93491456033504), (-0.224985604564191, 9.98987116588375), (3.27619056804565, -7.55422407977119), (-1.84280758110263, -9.29603360988560), (0.0547173670342382, 9.99940118401994), (1.12689578732265, -9.74271130322857), (-4.48313413815260, -4.42787005199146), (0.418760261366236, 9.96486624967952), (0.322704600896244, -9.97915061326572), (-3.97888287523214, 6.05590325721421), (1.40662085502927, 9.59612792123922), (0.00607307289600574, -9.99999262355440), (-1.48651419949427, 9.54783232670158), (3.83957376445349, 6.40552052758259), (-0.301192111491396, -9.98184017343002), (-0.446169128071070, 9.96010704945611), (4.59191252815097, -3.95693787356343), (-1.06472537429184, -9.77064171430702), (-0.0673523527424391, 9.99909269095553), (1.94345896446204, -9.21367836500757), (-3.13544619686423, -7.78947421757582), (0.207680875934373, 9.99137000691519), (0.604939543137603, -9.92653980985287), (-4.95464293979329, 1.34389484582172), (0.800299593073777, 9.87107300374705), (0.134050951276196, -9.99640542244300), (-2.50719688237237, 8.65192782980124), (2.47926391087429, 8.68406597400925), (-0.130825229378654, -9.99657636580805), (-0.810889594777716, 9.86761532794651), (4.94566432847514, 1.47024399342228), (-0.596790398685127, -9.92851272246498), (-0.211523256163942, 9.99104757512480), (3.16702592938540, -7.73820308924508), (-1.92037274658044, -9.23302085217858), (0.0644848023956064, 9.99916830746638), (1.07846798674574, -9.76461096031267), (-4.56822785051538, -4.06511712292554), (0.439865204478405, 9.96122855914654), (0.305956969121586, -9.98126050818151), (-3.87123619901044, 6.32883252787470), (1.46815071645876, 9.55919106907308), (-0.00333313159223542, -9.99999777804651), (-1.42430014466149, 9.58569123181365), (3.94772092170207, 6.13693719190796), (-0.317763861309131, -9.97978479295940), (-0.424822120011745, 9.96383985546711), (4.50849814659926, -4.32363017017461), (-1.11259038855275, -9.74928564096878), (-0.0575502448285888, 9.99933757192349), (1.86512696044980, -9.27821134085731), (-3.24428983693418, -7.60909546633911), (0.221020754406279, 9.99022518787673), (0.577357958838524, -9.93310782934844), (-4.92054600833228, 1.77564295947720), (0.837194692899198, 9.85882448290493), (0.123064723148217, -9.99697055590672), (-2.41216300761646, 8.75933093898997), (2.57602737444526, 8.57066694396852), (-0.141994416315457, -9.99596670377312), (-0.775050940372287, 9.87912871458370), (4.97316426560854, 1.03467325716718), (-0.625179150637130, -9.92152226820232), (-0.198472321964268, 9.99211864169238), (3.05907045985408, -7.91001590937671), (-2.00057726420937, -9.16464742582682), (0.0743780711636534, 9.99889351929102), (1.03194886614963, -9.78469857229184), (-4.64732695971651, -3.68909318369275), (0.461815438207198, 9.95725393891981), (0.289802083298979, -9.98318881971398), (-3.76206725748557, 6.58691124887986), (1.53205071336727, 9.51899587386629), (-0.0127458424066102, -9.99996750864748), (-1.36441372853888, 9.62047299822095), (4.05394673595617, 5.85338051455064), (-0.334950864556515, -9.97753635289451), (-0.404292359945486, 9.96725592882810), (4.41954926567458, -4.67657322759959), (-1.16241144287865, -9.72600630011410), (-0.0478604460918855, 9.99954186504559), (1.78941945693989, -9.33766095060748), (-3.35412075874800, -7.41616449001389), (0.234790377529844, 9.98896861114688), (0.550873234607410, -9.93912243196414), (-4.87742397191615, 2.20066844224902), (0.875632081808088, 9.84545955398902), (0.112318427414757, -9.99747659579418), (-2.31963975691627, 8.85872934412908), (2.67514244190120, 8.44834017201932), (-0.153440330199433, -9.99529010385757), (-0.740652258114627, 9.88967830266479), (4.99113851043795, 0.595101233819135), (-0.654749980579312, -9.91388974377492), (-0.185807584183565, 9.99309272280806), (2.95260459571357, -8.07022331819166), (-2.08343285655114, -9.09050219344194), (0.0844164627174124, 9.99857467058626), (0.987271965750077, -9.80312074099753), (-4.71979867408089, -3.30060629348261), (0.484651712407504, 9.95291167802879), (0.274209082732614, -9.98495055149437), (-3.65190903446066, 6.83039102951643), (1.59837800272488, 9.47527049965440), (-0.0221834329985889, -9.99990157857574), (-1.30679277250081, 9.65241786284445), (4.15766099181677, 5.55476554937298), (-0.352785809955450, -9.97507737760346), (-0.384541432073931, 9.97038171526417), (4.32572325056847, -5.01522416627265), (-1.21425548960755, -9.70063577420737), (-0.0382640545712798, 9.99970716813803), (1.71630744998415, -9.39239878564127), (-3.46459219581351, -7.21015975321161), (0.249016287576240, 9.98759047789240), (0.525437633797612, -9.94462976545415), (-4.82567758072196, 2.61750712466459), (0.915673538559779, 9.83087828645669), (0.101791157137347, -9.99792749730236), (-2.22968980020821, 8.95063872466037), (2.77648168635953, 8.31652558351387), (-0.165185196199823, -9.99454128031025), (-0.707635688092924, 9.89934376268190), (4.99941312362711, 0.153211217670698), (-0.685554918188670, -9.90555691602392), (-0.173504540875650, 9.99397742128639), (2.84787215485984, -8.21939759096103), (-2.16893893721101, -9.01015069499950), (0.0946195423887688, 9.99820926810359), (0.944371480803901, -9.82001273038773), (-4.78503937275883, -2.90061938292344), (0.508416199714271, 9.94816826714709), (0.259148112608650, -9.98655941868497), (-3.54125678234443, 7.05960350203880), (1.66718569332479, 9.42772334425934), (-0.0316643243286346, -9.99979947210239), (-1.25137310753798, 9.68174888038948), (4.25824612596106, 5.24112199084799), (-0.371302520202977, -9.97238876869347), (-0.365532226675534, 9.97324143721822), (4.22767808059384, -5.33919021832577), (-1.26818931938236, -9.67299247393577), (-0.0287423463408198, 9.99983477414039), (1.64575378193389, -9.44277384866339), (-3.57531710172166, -6.99059586076513), (0.263725862386027, 9.98608004564529), (0.501004968582454, -9.94967215971575), (-4.76577655642840, 3.02481325982118), (0.957382126244436, 9.81497212718359), (0.0914624196591103, -9.99832678517563), (-2.14235618764001, 9.03555420885306), (2.87988931985245, 8.17465289915109), (-0.177251800321505, -9.99371438440839), (-0.675944919945838, 9.90819831557674), (4.99790752463220, -0.289277549940218), (-0.717647591865551, -9.89646036396642), (-0.161539363019156, 9.99477964423336), (2.74508235242527, -8.35811530870527), (-2.25708081355513, -8.92313536848597), (0.105007188279201, 9.99779445486047), (0.903182479097403, -9.83549925717072), (-4.84248585648797, -2.49024555392757), (0.533152534456932, 9.94298715175719), (0.244590277259554, -9.98802795275826), (-3.43056618346342, 7.27494760417620), (1.73852189528069, 9.37604215426310), (-0.0412070199730003, -9.99966039053426), (-1.19808916791956, 9.70867289503850), (4.35506323990654, 4.91260589765340), (-0.390535999172610, -9.96944966050790), (-0.347228894440679, 9.97585727541558), (4.12606327329984, -5.64822161913861), (-1.32427915834316, -9.64288021511788), (-0.0192767399391997, 9.99992568118330), (1.57771442089055, -9.48911317376160), (-3.68586778195413, -6.75703446608126), (0.278947360735102, 9.98442554580661), (0.477530571302519, -9.95428843332792), (-4.69824957463656, 3.42137453923269), (1.00082206546679, 9.79762322061321), (0.0813120984906048, -9.99867758108822), (-2.05766450591623, 9.11394904113305), (2.98517879316551, 8.02214623971289), (-0.189663530436945, -9.99280295917461), (-0.645525222104020, 9.91630923027868), (4.98663639916843, -0.730636496456409), (-0.751083221078908, -9.88653103874432), (-0.149888852275934, 9.99550565643648), (2.64441202356069, -8.48695117215780), (-2.34782774964667, -8.82897612591382), (0.115599628568141, 9.99732698792530), (0.863641080781321, -9.84969524068372), (-4.89162646219149, -2.07073953397139), (0.558905848586973, 9.93732846441442), (0.230507593719100, -9.98936759744836), (-3.32025243056125, 7.47687736888926), (1.81242865455203, 9.31989321229782), (-0.0508301414595587, -9.99948324599210), (-1.14687450711368, 9.73338150180763), (4.44745930892862, 4.56951012491456), (-0.410522478990221, -9.96623726272734), (-0.329596799794617, 9.97824953577834), (4.02151187711018, -5.94220236015629), (-1.38259020114319, -9.61008727030152), (-0.00984876093143555, 9.99998060036281)]
len(points)
355
%matplotlib inline
from moviepy.video.io.bindings import mplfig_to_npimage
import moviepy.editor as mpy
duration = 5
fig, ax = plt.subplots(figsize=(6, 4), facecolor='white', dpi=150)
plt.tight_layout()
def make_frame_mpl(t):
index = int(t / duration * (len(points) - 1))
ax.clear()
ax.add_patch(Ellipse((0, 0), 10, 20, alpha=0.1))
ax.set_xlim(-12, 12)
ax.set_ylim(-12, 12)
ax.plot(np.array(points)[:index, 0], np.array(points)[:index, 1], '-o')
return mplfig_to_npimage(fig) # RGB image of the figure
animation = mpy.VideoClip(make_frame_mpl, duration=duration)
plt.close(fig)
animation.ipython_display(fps=25, width=800)
99%|███████████████████████████████████████▋| 125/126 [00:10<00:00, 11.57it/s]