#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 14 15:03:44 2026

@author: adomps
"""

# Étude du mouvement d'une charge ponctuelle autour d'un fil électrisé.
# J'ai renouvelé ce vieil exo à partir  X/ENS PSI 2024, mais les valeur numériques
# sont plutôt celles de la vidéo de l'expérience
# https://www.youtube.com/watch?v=qHrBhgwq__Q.

# On peut dans un premier temps observer les mouvemens sans frottement 
# (prendre eta = 0) et retrouver les résultats de l'étude par l'énergie
# potentielle effective. Dans un seconde temps, donner à eta la valeur de la
# viscosité de l'air, et reprendre le calcul sur une durée plus longue pour
# voir la goutte descendre lentement en spirale et s'écraser sur le fil.
# Pour ce mouvement, prendre alpha = 1 (vitesse de départ correspondant à 
# l'orbite circulaire).


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import fsolve



###############################################################################
###  Paramètres à modifier pour voir différentes situations  ##################
alpha = 1.2 # facteur v_0 / v_c
eta = 1.8e-5 
eta = 0  # Pendre eta = 0 pour supprimer les frottements
###############################################################################


R = 3e-3 # rayon de la goutte
r0 = 20e-3 # distance initiale
T = 3 # période de rotation très grossièrement déduite de la vidéo
v_c = 2 * np.pi * r0 / T # vitesse de la goutte en orbitre circulaire
rho_eau = 1000
m = 4 / 3 * np.pi * R**3 * rho_eau
rho_air = 1.2
C_Stokes = 6 * np.pi * eta * R
rf = 3e-3 # rayon du fil

# facteur beta déduit de la période observée, en supposant le mouvement
# circulaire à la vitesse v0
beta = m * v_c**2  # beta = - q * lambda / (2 * pi * epsilon0)

v_0 = alpha * v_c


def Epeff(r) :
    ''' Epeff au facteur m v_c**2 près. '''
    return 0.5 * alpha**2 * r0**2 / r**2 + np.log(r/r0)


def ecart_energie(r) :
    ''' Renvoie l'écart entre Em et Epeff, au facteur  m v_c**2 près. 
    L'annulation de cet écart définit les extremums radiaux du mouvement. '''
    return 0.5 * alpha**2 - Epeff(r)

r_extre = alpha * r0 # extremum de Epeff
# r_stop : extremum radial distinct de r0
r_stop = fsolve(ecart_energie, r_extre + 0.5 * (r_extre - r0)) #  on recherche du bon côté
print('r_stop = ', r_stop)
# r_gauche et r_droit : bornes pour le graphique
r_gauche = 0.6 * min(r0, r_stop)
r_droit = 1.5 * max(r0, r_stop)
r = np.linspace(r_gauche, r_droit, 1000)
e0 = 0.5 * alpha**2

plt.figure(0)
plt.clf()
plt.plot(r, Epeff(r))
plt.vlines(r0, 0, Epeff(r0), color = 'black', linestyle ='dashed')
plt.vlines(r_stop, 0, Epeff(r_stop), color = 'black', linestyle ='dashed')
plt.hlines(e0, r_gauche, r_droit, color = 'red')
plt.ylim(0, np.max(Epeff(r)))
plt.xlim(r_gauche, r_droit)
plt.grid()
plt.text(r0 + 0.0005, 0.02*e0, '$r_0$')
plt.xlabel('$r$ (m)')
plt.ylabel('$E_{peff}\,\, / \,\,m v_c^2$')



def systeme_diff(X, t) :
    r, theta, rp, thetap = X
    #### À compléter
    ####
    return [rp, thetap, rpp, thetapp]
    
X0 = [r0, 0, 0, v_0/r0]  
t_duree = 140
t = np.linspace(0, t_duree, 10000)    

# Utiliser odeint pour résoudre le système
soluce = ### À compléter
r = soluce[:, 0]
theta = soluce[:, 1]
    
# On trace r(t) et theta(t)
# plt.figure(1)
# plt.subplot(2, 1, 1)
# plt.plot(t, r)
# plt.subplot(2, 1, 2)
# plt.plot(t, theta)

# J'ai multiplié les rayons par 1000 pour afficher en mm
r0_cercle = [1000 * r0 for i in range(len(theta))]
r_stop_cercle = [1000 * r_stop for i in range(len(theta))]

plt.figure(2)
plt.clf()
ax = plt.subplot(111, polar=True)
ax.set_title('')
ax.plot(theta, 1000 * r, color='black', linewidth=2)
ax.plot(theta, r0_cercle, color = 'red', linewidth=2)
ax.plot(theta, r_stop_cercle, color = 'red', linewidth=2)

# On recherche le contact avec le fil.
# Le contact se produit quand r = R + rf
if eta > 0 :
    tau = m / (6 * np.pi * eta * R)
    tf = tau * np.log(r0 / (rf+R)) # instant de contact avec le fil par raisonnement énergétique
    print('instant de contact avec le fil dans le modèle : ', tf)
    i_contact = np.argmin(abs(r-(rf+R)))
    t_contact = t[i_contact]