Lancer de rayons (sujet Centrale 2021)

Sujet : https://www.concours-centrale-supelec.fr/CentraleSupelec/2021/Multi/sujets/I007.pdf

import numpy as np
import plotly
import matplotlib.pyplot as plt

plotly.offline.init_notebook_mode(connected=True)  # pour export html

def point(x, y, z):
    return np.array([x, y, z])

I Géométrie

# Q1
def vec(A, B):
    return np.array(B) - np.array(A)

# Q2 
def ps(v1, v2):
    return np.sum(v1 * v2)

# Q3
def norme(v):
    return np.sqrt(np.sum(v**2))

# Q4
def unitaire(v):
    return v / norme(v)

# Q5
def pt(r, t):
    S, u = r
    return S + t * u

def dir(A, B):
    return unitaire(vec(A, B))

def ra(A, B):
    return A, dir(A, B)

# Q6
def sp(A, B):
    return (A, norme(B - A))

# Q8
def intersection(r, s):
    S, u = r
    C, r_s = s
    CS = vec(S, C)
    b = 2*ps(u, CS)
    delta = b**2 - 4*(norme(CS)**2 - r_s**2)
    if delta < 0:
        return None
    t1, t2 = -b - delta**0.5, -b + delta**0.5
    t = min(abs(t1), abs(t2))/2
    return S + t*u, t

II Optique

Q 9. P est le premier point d'intersection d'un rayon partant de S et allant en P. (ou : Produit scalaire positif avec la normale au plan tangent)

# Q10
def au_dessus(s, P, src):
    return ps(vec(P, src), vec(s[0], P)) > 0

# Q11
def visible(obj, j, P, src):
    for i, o in enumerate(obj):
        inter = intersection(ra(src, P), o)
        if i != j and inter != None and inter[1] < norme(vec(src, P)):
            return False
    return au_dessus(obj[j], P, src)

# Q12
def couleur_diffusée(r, Cs, N, kd):
    cos = ps(-r[1], N)
    return cos*(Cs*kd)

# Q13
def rayon_réfléchi(s, P, src):
    C, r = s
    u = dir(src, P)
    N = dir(C, P)
    return P, 2*N + u

Digression : affichage d'un rayon avec Plotly

import plotly.graph_objects as go

def line(P, Q):
    return go.Scatter3d(x=[P[0], Q[0]], y=[P[1], Q[1]], z=[P[2], Q[2]])

theta, phi = np.mgrid[0:2*np.pi:100j, 0:2*np.pi:100j]
s = (point(0, 0, 0), 1)
src = point(3, 3, 2)
P = point(1, 0, 0)
_, u = rayon_réfléchi(s, P, src)

fig = go.Figure(data=[ 
    line(src, P),
    go.Surface(x=np.cos(theta)*np.cos(phi), y=np.sin(theta)*np.cos(phi), z=np.sin(phi)),
    line(P, P + 3*u),
])
fig.show(renderer="notebook")