# -*- coding: utf-8 -*-
"""
Created on Wed Apr  1 03:11:31 2026

@author: pjaub
"""

import numpy as np
import matplotlib.pyplot as plt

plt.close('all')

# ==================================================
# 1) Fonction et gradient
# ==================================================
def f(x, y):
    return 2*np.sin(x)*np.cos(y) + 0.5*np.cos(2*x) - 0.7*np.sin(2*y)

def grad_f(x, y):
    # dérivées partielles
    gx = 2*np.cos(x)*np.cos(y) - np.sin(2*x)
    gy = -2*np.sin(x)*np.sin(y) - 1.4*np.cos(2*y)
    return gx, gy

# ==================================================
# 2) Intégration d'une ligne de plus grande pente
# ==================================================
def ligne_pente(x0, y0, h=0.06, n=250, sens=1,
                xmin=-4, xmax=4, ymin=-4, ymax=4):
    """
    sens = 1  : ligne de plus grande montée
    sens = -1 : ligne de plus grande descente
    """
    xs = [x0]
    ys = [y0]

    x, y = x0, y0

    for _ in range(n):
        gx, gy = grad_f(x, y)
        norme = np.sqrt(gx**2 + gy**2)

        # arrêt près d'un point critique
        if norme < 1e-8:
            break

        # direction normalisée
        x = x + sens * h * gx / norme
        y = y + sens * h * gy / norme

        # arrêt si on sort de la fenêtre
        if not (xmin <= x <= xmax and ymin <= y <= ymax):
            break

        xs.append(x)
        ys.append(y)

    xs = np.array(xs)
    ys = np.array(ys)
    zs = f(xs, ys)
    return xs, ys, zs

# ==================================================
# 3) Grille
# ==================================================
xmin, xmax = -4, 4
ymin, ymax = -4, 4

x = np.linspace(xmin, xmax, 250)
y = np.linspace(ymin, ymax, 250)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)

# gradient sur la grille
Ux, Uy = grad_f(X, Y)

# plan de projection sous la surface
z0 = Z.min() - 1.0

# ==================================================
# 4) Points de départ des lignes de pente
# ==================================================
points_depart = [
    (-3.2, -2.8),
    (-2.4,  0.8),
    (-1.4, -1.5),
    (-0.3,  2.6),
    ( 1.2, -2.4),
    ( 2.0,  1.3),
    ( 3.0, -0.8),
]

# ==================================================
# 5) Sous-échantillonnage pour le champ de gradient
# ==================================================
step = 14
Xq = X[::step, ::step]
Yq = Y[::step, ::step]
Uq = Ux[::step, ::step]
Vq = Uy[::step, ::step]

# normalisation pour afficher surtout la direction
Nq = np.sqrt(Uq**2 + Vq**2)
Nq[Nq == 0] = 1
Uq = Uq / Nq
Vq = Vq / Nq

# ==================================================
# 6) Figure avec deux subplots
# ==================================================
fig = plt.figure(figsize=(14, 6))

# --------------------------------------------------
# SUBPLOT 1 : surface 3D
# --------------------------------------------------
ax1 = fig.add_subplot(1, 2, 1, projection='3d')

# surface
ax1.plot_surface(X, Y, Z, cmap='viridis', alpha=0.55, edgecolor='none')

# lignes de niveau projetées sous la surface
ax1.contour(X, Y, Z, levels=18, zdir='z', offset=z0, cmap='plasma')

# champ de gradient projeté sous la surface
ax1.quiver(
    Xq, Yq, z0*np.ones_like(Xq),
    Uq, Vq, np.zeros_like(Uq),
    length=0.45, color='black', normalize=False
)

# lignes de plus grande pente
for (x0, y0) in points_depart:
    # montée
    xs, ys, zs = ligne_pente(x0, y0, sens=1,
                             xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
    ax1.plot(xs, ys, zs, color='red', linewidth=2.5)
    ax1.plot(xs, ys, z0*np.ones_like(xs), color='red', linewidth=1.5)

    # descente
    xs, ys, zs = ligne_pente(x0, y0, sens=-1,
                             xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
    ax1.plot(xs, ys, zs, color='blue', linewidth=2.5, linestyle='--')
    ax1.plot(xs, ys, z0*np.ones_like(xs), color='blue', linewidth=1.5, linestyle='--')

ax1.set_title("Surface 3D")
ax1.set_xlabel("x")
ax1.set_ylabel("y")
ax1.set_zlabel("z")
ax1.set_zlim(z0, Z.max())
ax1.view_init(elev=30, azim=-60)

# --------------------------------------------------
# SUBPLOT 2 : plan (x,y)
# --------------------------------------------------
ax2 = fig.add_subplot(1, 2, 2)

# lignes de niveau
ax2.contour(X, Y, Z, levels=18, cmap='plasma')

# gradient
ax2.quiver(Xq, Yq, Uq, Vq, color='black')

# lignes de plus grande pente
for (x0, y0) in points_depart:
    # montée
    xs, ys, _ = ligne_pente(x0, y0, sens=1,
                            xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
    ax2.plot(xs, ys, color='red', linewidth=2.5)

    # descente
    xs, ys, _ = ligne_pente(x0, y0, sens=-1,
                            xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
    ax2.plot(xs, ys, color='blue', linewidth=2.5, linestyle='--')

# points de départ
for (x0, y0) in points_depart:
    ax2.plot(x0, y0, 'ko', markersize=4)

ax2.set_title("Plan (x,y)")
ax2.set_xlabel("x")
ax2.set_ylabel("y")
ax2.set_aspect('equal')
ax2.set_xlim(xmin, xmax)
ax2.set_ylim(ymin, ymax)
ax2.grid(True, alpha=0.25)

plt.tight_layout()
plt.show()