Atelier 3 - Déviation des forages

import numpy as np
import plotly.graph_objs as go
import ipywidgets as widgets
from IPython.display import display, clear_output, Markdown

def d2r(deg):
    return np.deg2rad(deg)

def compute_cosine_direction(measurements, show_debug=False):
    """
    Retourne uniquement les cosinus directeurs (l, m, n)
    pour chaque station donnée (MD, azimut, inclinaison).
    
    measurements : liste de tuples (MD, azimut, inclinaison)
    """
    cos_dirs = []
    for i, (md, azi, plo) in enumerate(measurements):
        azi_rad, inc_rad = d2r(azi), d2r(plo)
        l = np.cos(inc_rad) * np.sin(azi_rad)  # X
        m = np.cos(inc_rad) * np.cos(azi_rad)  # Y
        n = -np.sin(inc_rad)                   # Z
        cos_dirs.append((md, l, m, n))

        if show_debug:
            print(f"Station {i+1} (M{i+1}={md}): azimut={azi}°, inclinaison={plo}°")
            print(f"  → cos_dir = (lx={l:.4f}, ly={m:.4f}, lz={n:.4f})\n")

    return cos_dirs

def compute_positions(measurements, start_pos=(0,0,0)):
    """
    Calcule la trajectoire cumulée (positions des stations)
    et les points milieux pour la méthode des points milieux.
    """
    depths = [m[0] for m in measurements]
    cos_dirs = compute_cosine_direction(measurements)

    # Segments effectifs
    seg_lengths = []
    n = len(measurements)
    for i in range(n):
        if i == 0:
            seg_len = (depths[1] - depths[0]) / 2 if n > 1 else 0
        elif i == n-1:
            seg_len = (depths[-1] - depths[-2]) / 2
        else:
            seg_len = (depths[i+1] - depths[i-1]) / 2
        seg_lengths.append(seg_len)

    # Positions cumulées
    positions = [start_pos]
    mids = []
    x, y, z = start_pos
    for i, ((md, _, _), (_, l, m, n_)) in enumerate(zip(measurements, cos_dirs)):
        dx, dy, dz = seg_lengths[i] * l, seg_lengths[i] * m, seg_lengths[i] * n_
        x, y, z = x + dx, y + dy, z + dz
        positions.append((x, y, z))
        mids.append((md, x, y, z))
    return positions, mids

def interpolate_composites(measurements, start_pos=(0,0,0), targets=[], show_debug=False):
    """
    Calcule la position (x,y,z) des composites pour des longueurs données (targets)
    à partir de mesures de déviation (MD, azimut, inclinaison).
    
    measurements : liste de tuples (MD, azimut, inclinaison)
    start_pos : coordonnées initiales (x,y,z)
    targets : longueurs où placer les composites
    """
    n_stations = len(measurements)

    # Cas particulier : une seule mesure
    if n_stations == 1:
        cos_dirs = compute_cosine_direction(measurements, show_debug=show_debug)
        l, m, n_ = cos_dirs[0][1:]  # ignorer MD
        interp_XYZ = []
        for t in targets:
            dx, dy, dz = t * l, t * m, -t * n_
            xt, yt, zt = start_pos[0] + dx, start_pos[1] + dy, start_pos[2] + dz
            interp_XYZ.append((t, xt, yt, zt))

        return interp_XYZ

    # Cas général : plusieurs mesures
    depths = [m[0] for m in measurements]

    # 1) longueurs "effectives" par méthode des points milieux
    seg_lengths = []
    for i in range(n_stations):
        if i == 0:
            seg_len = (depths[1] - depths[0]) / 2
        elif i == n_stations - 1:
            seg_len = (depths[-1] - depths[-2]) / 2
        else:
            seg_len = (depths[i+1] - depths[i-1]) / 2
        seg_lengths.append(seg_len)

    # 2) calculer tous les cosinus directeurs en une seule fois
    cos_dirs = compute_cosine_direction(measurements, show_debug=show_debug)

    # 3) interpolation des composites

    interp_XYZ = []
    seg_lengths[-1] = 8000000
    for t in targets:
        x, y, z = start_pos
        depth_accum = 0
        for i, seg_len in enumerate(seg_lengths):
            l, m, n_ = cos_dirs[i][1:]  # ignorer MD

            if t <= depth_accum + seg_len:
                remaining = t - depth_accum
                dx, dy, dz = remaining * l, remaining * m, remaining * n_
                x, y, z = x + dx, y + dy, z + dz
                break
            else:
                dx, dy, dz = seg_len * l, seg_len * m, seg_len * n_
                x, y, z = x + dx, y + dy, z + dz
                depth_accum += seg_len

        interp_XYZ.append((t, x, y, z))

    return interp_XYZ


def plotly_3d_path(positions, mids, interp, measurements):
    pos_arr = np.array(positions)
    fig = go.Figure()

    fig.add_trace(go.Scatter3d(
        x=pos_arr[:,0], y=pos_arr[:,1], z=pos_arr[:,2],
        mode='lines+markers',
        line=dict(color='blue', width=5),
        marker=dict(size=5, color='blue'),
        name='Trajectoire'
    ))

    md_vals = [m[0] for m in measurements]
    stations_coords = pos_arr[1:]
    fig.add_trace(go.Scatter3d(
        x=stations_coords[:,0], y=stations_coords[:,1], z=stations_coords[:,2],
        mode='markers+text',
        marker=dict(size=7, color='orange'),
        text=[f"MD {md:.1f} m" for md in md_vals],
        textposition='top center',
        name='Stations (mesures)'
    ))

    if mids:
        mids_np = np.array(mids)
        fig.add_trace(go.Scatter3d(
            x=mids_np[:,1], y=mids_np[:,2], z=mids_np[:,3],
            mode='markers+text',
            marker=dict(size=6, color='green'),
            text=[f"MD milieu {md:.1f} m" for md in mids_np[:,0]],
            textposition='bottom center',
            name='Points milieux'
        ))

    if interp:
        interp_np = np.array(interp)
        fig.add_trace(go.Scatter3d(
            x=interp_np[:,1], y=interp_np[:,2], z=interp_np[:,3],
            mode='markers+text',
            marker=dict(size=6, color='red'),
            text=[f"MD {md:.1f} m" for md in interp_np[:,0]],
            textposition='top center',
            name='Composites interpolés'
        ))

    fig.update_layout(
        scene=dict(
            xaxis_title='X (m)',
            yaxis_title='Y (m)',
            zaxis_title='Z (m, ↑)',
            aspectmode='data',
        ),
        title='Trajectoire 3D par méthode des points milieux',
        margin=dict(l=0, r=0, b=0, t=40),
        legend=dict(x=0.05, y=0.95)
    )
    fig.show()

# Widgets
num_points_w = widgets.BoundedIntText(value=3, min=1, max=50, description="Nombre stations:")
md_list_w = widgets.Text(value="0, 40, 75", description="Mesure (m):")
azi_list_w = widgets.Text(value="90,100,110", description="Azimut (\u00b0):")
inc_list_w = widgets.Text(value="40,35,25", description="Plongée (\u00b0):")
targets_w = widgets.Text(value="76.5", description="Composites (MD):")
collet_w = widgets.Text(value="0, 50, 32", description="Collet (X,Y,Z):")
show_debug_w = widgets.Checkbox(value=False, description="Afficher calculs intermédiaires")
button = widgets.Button(description="Calculer et Tracer", button_style='success', icon='play')
out = widgets.Output()

def on_button_clicked(b):
    with out:
        clear_output()
        try:
            md_list = list(map(float, md_list_w.value.strip().split(',')))
            azi_list = list(map(float, azi_list_w.value.strip().split(',')))
            inc_list = list(map(float, inc_list_w.value.strip().split(',')))
            collet = tuple(map(float, collet_w.value.strip().split(',')))

            targets_str = targets_w.value.strip()
            targets = list(map(float, targets_str.split(','))) if targets_str else []

            n = num_points_w.value
            if not (len(md_list) == len(azi_list) == len(inc_list) == n):
                print("⚠️ Erreur : Le nombre de valeurs pour MD, azimut et inclinaison doit correspondre au nombre de stations.")
                return

            measurements = list(zip(md_list, azi_list, inc_list))
            show_debug = show_debug_w.value

            # 🔹 Calcul trajectoire
            positions, mids = compute_positions(measurements, start_pos=collet)

            # 🔹 Interpolation des composites
            interp = interpolate_composites(measurements, collet, targets, show_debug)

            if interp:
                interp.sort(key=lambda x: x[0])
                md_output = "### Coordonnées des Composites\n\n"
                md_output += "| Profondeur (MD) | X | Y | Z |\n"
                md_output += "|:---:|:---:|:---:|:---:|\n"
                for t, x, y, z in interp:
                    md_output += f"| {t:.2f} m | {x:.2f} | {y:.2f} | {z:.2f} |\n"
                display(Markdown(md_output))
            elif targets:
                print("ℹ️ Aucun composite n'a pu être calculé.")

            # 🔹 Tracé 3D
            plotly_3d_path(positions, mids, interp, measurements)

        except Exception as e:
            print(f"❌ Erreur lors de l'exécution : {str(e)}")


button.on_click(on_button_clicked)

info_text = widgets.Output()
with info_text:
    display(Markdown("""
    ### ℹ️ Conventions et paramètres
    - **Azimut (°)** : 0-360°, mesuré depuis le Nord dans le sens horaire.
    - **Inclinaison (°)** : angle entre le plan XY et l'axe Z descendant. 90° = forage vertical descendant.
    - **MD (m)** : profondeur mesurée le long du forage.
    - **Composites (MD)** : profondeurs où calculer les coordonnées interpolées.
    - **Collet (X,Y,Z)** : coordonnées du départ du forage.
    ---
    """))

inputs = widgets.VBox([num_points_w, md_list_w, azi_list_w, inc_list_w, targets_w, collet_w, show_debug_w, button])
main_ui = widgets.VBox([info_text, inputs, out])

display(main_ui)