# -*- coding: utf-8 -*-
"""
Created on Fri Sep  5 16:55:22 2025

@author: AKourgli
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.constants import speed_of_light, pi
from scipy import special
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D

def reflection_coefficient(epsilon_r, theta_i, polarization='vertical'):
    """
    Coefficient de réflexion selon les équations de Fresnel
    epsilon_r: constante diélectrique relative
    theta_i: angle d'incidence (radians)
    polarization: 'vertical' ou 'horizontal'
    """
    sin_theta_t = np.sin(theta_i) / np.sqrt(epsilon_r)
    cos_theta_t = np.sqrt(1 - sin_theta_t**2)
    
    if polarization == 'vertical':
        num = epsilon_r * np.cos(theta_i) - np.sqrt(epsilon_r - np.sin(theta_i)**2)
        den = epsilon_r * np.cos(theta_i) + np.sqrt(epsilon_r - np.sin(theta_i)**2)
    else:  # horizontal
        num = np.cos(theta_i) - np.sqrt(epsilon_r - np.sin(theta_i)**2)
        den = np.cos(theta_i) + np.sqrt(epsilon_r - np.sin(theta_i)**2)
    
    return num / den

# Simulation de la réflexion
theta_i = np.linspace(0, pi/2, 100)
epsilon_r_values = [4, 9, 16]  # béton, verre, métal

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
for eps in epsilon_r_values:
    Gamma_v = reflection_coefficient(eps, theta_i, 'vertical')
    plt.plot(np.degrees(theta_i), np.abs(Gamma_v), label=f'εr={eps} (vertical)')

plt.xlabel('Angle d\'incidence (degrés)')
plt.ylabel('|Γ|')
plt.title('Coefficient de réflexion (polarisation verticale)')
plt.legend()
plt.grid(True)

plt.subplot(1, 2, 2)
for eps in epsilon_r_values:
    Gamma_h = reflection_coefficient(eps, theta_i, 'horizontal')
    plt.plot(np.degrees(theta_i), np.abs(Gamma_h), label=f'εr={eps} (horizontal)')

plt.xlabel('Angle d\'incidence (degrés)')
plt.ylabel('|Γ|')
plt.title('Coefficient de réflexion (polarisation horizontale)')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

def free_space_path_loss(d, f, Gt=1, Gr=1):
    """
    Calcul de l'affaiblissement en espace libre
    d: distance (m)
    f: fréquence (Hz)
    Gt, Gr: gains d'antenne (linéaires)
    """
    lambda_ = speed_of_light / f
    Pr_Pt = Gt * Gr * (lambda_ / (4 * pi * d))**2
    return 10 * np.log10(Pr_Pt)  # en dB

# Simulation
distances = np.linspace(1, 10000, 1000)  # de 1m à 10km
frequencies = [900e6, 2.4e9, 5e9]  # 900 MHz, 2.4 GHz, 5 GHz

plt.figure(figsize=(10, 6))
for f in frequencies:
    loss = free_space_path_loss(distances, f)
    plt.plot(distances, loss, label=f'{f/1e9} GHz')

plt.xlabel('Distance (m)')
plt.ylabel('Affaiblissement (dB)')
plt.title('Affaiblissement en espace libre (Formule de Friis)')
plt.legend()
plt.grid(True)
plt.show()

def reflection_coefficient(epsilon_r, theta_i, polarization='vertical'):
    """
    Coefficient de réflexion selon les équations de Fresnel
    epsilon_r: constante diélectrique relative
    theta_i: angle d'incidence (radians)
    polarization: 'vertical' ou 'horizontal'
    """
    sin_theta_t = np.sin(theta_i) / np.sqrt(epsilon_r)
    cos_theta_t = np.sqrt(1 - sin_theta_t**2)
    
    if polarization == 'vertical':
        num = epsilon_r * np.cos(theta_i) - np.sqrt(epsilon_r - np.sin(theta_i)**2)
        den = epsilon_r * np.cos(theta_i) + np.sqrt(epsilon_r - np.sin(theta_i)**2)
    else:  # horizontal
        num = np.cos(theta_i) - np.sqrt(epsilon_r - np.sin(theta_i)**2)
        den = np.cos(theta_i) + np.sqrt(epsilon_r - np.sin(theta_i)**2)
    
    return num / den

# Simulation de la réflexion
theta_i = np.linspace(0, pi/2, 100)
epsilon_r_values = [4, 9, 16]  # béton, verre, métal

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
for eps in epsilon_r_values:
    Gamma_v = reflection_coefficient(eps, theta_i, 'vertical')
    plt.plot(np.degrees(theta_i), np.abs(Gamma_v), label=f'εr={eps} (vertical)')

plt.xlabel('Angle d\'incidence (degrés)')
plt.ylabel('|Γ|')
plt.title('Coefficient de réflexion (polarisation verticale)')
plt.legend()
plt.grid(True)

plt.subplot(1, 2, 2)
for eps in epsilon_r_values:
    Gamma_h = reflection_coefficient(eps, theta_i, 'horizontal')
    plt.plot(np.degrees(theta_i), np.abs(Gamma_h), label=f'εr={eps} (horizontal)')

plt.xlabel('Angle d\'incidence (degrés)')
plt.ylabel('|Γ|')
plt.title('Coefficient de réflexion (polarisation horizontale)')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

def knife_edge_diffraction(v):
    """
    Modèle de diffraction par un obstacle coupant (couteau de Kirchhoff)
    v: paramètre de Fresnel-Kirchhoff
    """
    # Fonctions de Fresnel
    C, _ = special.fresnel(v)
    S, _ = special.fresnel(v)
    
    # Formule de diffraction
    F = ((1 + 1j) / 2) * (0.5 - C - 1j*(0.5 - S))
    return 20 * np.log10(np.abs(F))

# Simulation de la diffraction
v = np.linspace(-3, 3, 200)  # paramètre de Fresnel-Kirchhoff
diffraction_loss = knife_edge_diffraction(v)

plt.figure(figsize=(10, 6))
plt.plot(v, diffraction_loss)
plt.xlabel('Paramètre de Fresnel-Kirchhoff (v)')
plt.ylabel('Perte par diffraction (dB)')
plt.title('Diffraction par un obstacle coupant')
plt.grid(True)
plt.show()

def okumura_hata(f, d, hb, hm, environment='urban'):
    """
    Modèle d'Okumura-Hata pour la propagation en milieu urbain
    f: fréquence (MHz)
    d: distance (km)
    hb: hauteur de l'antenne BS (m)
    hm: hauteur de l'antenne MS (m)
    environment: 'urban', 'suburban', 'open'
    """
    # Calcul de la perte de base
    L = 69.55 + 26.16 * np.log10(f) - 13.82 * np.log10(hb)
    L += (44.9 - 6.55 * np.log10(hb)) * np.log10(d)
    
    # Correction pour la hauteur de l'antenne mobile
    a_hm = (1.1 * np.log10(f) - 0.7) * hm - (1.56 * np.log10(f) - 0.8)
    
    # Correction selon l'environnement
    if environment == 'suburban':
        L -= 2 * (np.log10(f/28))**2 - 5.4
    elif environment == 'open':
        L -= 4.78 * (np.log10(f))**2 + 18.33 * np.log10(f) - 40.94
    
    return L - a_hm

# Simulation pour différents environnements
f = 900  # MHz
d = np.linspace(1, 20, 100)  # km
hb, hm = 30, 1.5  # m

environments = ['urban', 'suburban', 'open']
colors = ['red', 'blue', 'green']

plt.figure(figsize=(10, 6))
for env, color in zip(environments, colors):
    loss = okumura_hata(f, d, hb, hm, env)
    plt.plot(d, loss, color=color, label=env.capitalize())

plt.xlabel('Distance (km)')
plt.ylabel('Perte de propagation (dB)')
plt.title('Modèle d\'Okumura-Hata (f=900 MHz)')
plt.legend()
plt.grid(True)
plt.show()

def multipath_channel(direct_path, reflected_paths, frequency, velocity=0):
    """
    Simulation d'un canal multi-trajets
    direct_path: longueur du trajet direct (m)
    reflected_paths: liste de tuples (longueur, amplitude, déphasage)
    frequency: fréquence (Hz)
    velocity: vitesse du mobile (m/s) pour l'effet Doppler
    """
    lambda_ = speed_of_light / frequency
    t = np.linspace(0, 1e-5, 1000)  # temps
    
    # Signal du trajet direct
    signal = np.exp(1j * 2 * pi * frequency * t)
    
    # Ajout des trajets réfléchis
    for path_length, amplitude, phase_shift in reflected_paths:
        delay = (path_length - direct_path) / speed_of_light
        doppler_shift = (velocity / speed_of_light) * frequency if velocity else 0
        reflected_signal = amplitude * np.exp(1j * (2 * pi * (frequency + doppler_shift) * (t - delay) + phase_shift))
        signal += reflected_signal
    
    return t, signal

# Simulation
direct_path = 1000  # m
reflected_paths = [
    (1100, 0.8, pi/4),   # trajet 1: 100m plus long, amplitude 0.8, déphasage π/4
    (1200, 0.6, pi/2),   # trajet 2: 200m plus long, amplitude 0.6, déphasage π/2
]

t, signal = multipath_channel(direct_path, reflected_paths, 900e6)

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.plot(t, np.real(signal))
plt.xlabel('Temps (s)')
plt.ylabel('Amplitude')
plt.title('Signal résultant (partie réelle)')
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(t, 20 * np.log10(np.abs(signal)))
plt.xlabel('Temps (s)')
plt.ylabel('Puissance (dB)')
plt.title('Évanouissement (fading) du signal')
plt.grid(True)

plt.tight_layout()
plt.show()

def visualize_propagation_3d():
    """Visualisation 3D de la propagation d'ondes"""
    fig = plt.figure(figsize=(10, 8))
    ax = fig.add_subplot(111, projection='3d')
    
    # Source
    x0, y0, z0 = 0, 0, 10
    ax.scatter(x0, y0, z0, color='red', s=100, label='Source')
    
    # Ondes sphériques
    r = np.linspace(1, 50, 100)
    theta = np.linspace(0, 2*pi, 100)
    R, Theta = np.meshgrid(r, theta)
    
    X = R * np.cos(Theta)
    Y = R * np.sin(Theta)
    Z = np.sqrt(np.maximum(0, R**2 - X**2 - Y**2))  # Demi-sphère
    
    # Affichage des ondes
    ax.plot_surface(X, Y, Z, alpha=0.3, color='blue')
    ax.plot_surface(X, Y, -Z, alpha=0.3, color='blue')
    
    # Obstacle
    x_obs = np.linspace(-20, 20, 10)
    y_obs = np.linspace(-20, 20, 10)
    X_obs, Y_obs = np.meshgrid(x_obs, y_obs)
    Z_obs = np.zeros_like(X_obs) + 15
    
    ax.plot_wireframe(X_obs, Y_obs, Z_obs, color='gray', alpha=0.5, label='Obstacle')
    
    ax.set_xlabel('X (m)')
    ax.set_ylabel('Y (m)')
    ax.set_zlabel('Z (m)')
    ax.set_title('Propagation d\'ondes électromagnétiques en 3D')
    ax.legend()
    
    plt.show()

visualize_propagation_3d()

def animate_wave_propagation():
    """Animation de la propagation d'une onde"""
    fig, ax = plt.subplots(figsize=(10, 6))
    
    x = np.linspace(-10, 10, 1000)
    t_values = np.linspace(0, 2*np.pi, 100)
    
    line, = ax.plot(x, np.sin(x))
    ax.set_ylim(-1.5, 1.5)
    ax.set_xlabel('Distance')
    ax.set_ylabel('Amplitude')
    ax.set_title('Propagation d\'une onde électromagnétique')
    ax.grid(True)
    
    def animate(i):
        line.set_ydata(np.sin(x - t_values[i]))
        return line,
    
    ani = animation.FuncAnimation(fig, animate, frames=len(t_values), 
                                  interval=50, blit=True)
    plt.close()
    return ani

# Pour afficher l'animation dans un notebook Jupyter
# from IPython.display import HTML
# HTML(animate_wave_propagation().to_jshtml())

# Pour sauvegarder l'animation
# animate_wave_propagation().save('wave_propagation.gif', writer='pillow', fps=20)