# -*- coding: utf-8 -*-
"""
Created on Thu Oct 12 13:42:08 2017

@author: Louis
"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

#==============================================================================
# Définition de plusieurs variables et d'exemples
#==============================================================================
p = 0.01   #Précision pour les dérivées
preci = 5   #Nombre de décimales de précision pour la comparaison dans verification_surface_reglee()
Lu = np.linspace(-10,10,40) #Discrétisation d'un intervalle des u
Lv = np.linspace(-10,10,40) #Discrétisation d'un intervalle des v

#Définition de plusieurs surfaces connues :
def f1(u,v): #Hélice
    x=1*v*np.cos(u)
    y=1*v*np.sin(u)
    z=1*u
    return np.array([x,y,z])
def f2(u,v): #Sphère
    x=np.cos(u)*np.cos(v)
    y=np.sin(u)*np.cos(v)
    z=np.sin(v)
    return np.array([x,y,z])
def f3(u,v): #Cône
    x=v*np.cos(u)
    y=v*np.sin(u)
    z=v
    return np.array([x,y,z])
def f4(u,v): #Hyperboloïde
    a,b,c,d = 1,1,1,1
    x=a*np.sqrt(u**2 + d)*np.cos(v)
    y=b*np.sqrt(u**2 + d)*np.sin(v)
    z=c*u
    return np.array([x,y,z])

#==============================================================================
# Mise en place de la méthode
#==============================================================================

#Définition des fonctions de dérivation partielles et de la fonction d'intrégration de Gp pour obtenir G et K
def derive_u(f,u,v,pas=p):
    return (f(u+pas,v)-f(u,v))/pas
def derive_v(f,u,v,pas=p):
    return (f(u,v+pas)-f(u,v))/pas
def derive_u_u(f,u,v,pas=p):
    return (derive_u(f,u+pas,v,pas)-derive_u(f,u,v,pas))/pas
def derive_v_v(f,u,v,pas=p):
    return (derive_v(f,u,v+pas,pas)-derive_v(f,u,v,pas))/pas
def derive_u_v(f,u,v,pas=p):
    return (derive_v(f,u+pas,v,pas)-derive_v(f,u,v,pas))/pas
def derive_v_u(f,u,v,pas=p):
    return (derive_u(f,u,v+pas,pas)-derive_u(f,u,v,pas))/pas
def primitive_u(Gp,f,u,v,pas=p):
    S=0
    P = np.arange(0,u+pas,pas)
    for p in P:
        S += pas * Gp(f,p,v)[0]
    return np.array([S,Gp(f,u,v)[1],Gp(f,u,v)[2]])

#Obtention de K et de G
def K(f,u,v):
    return derive_v(f,u,v)

def Gp(f,u,v):
    return derive_u(f,u,v) - v*derive_u_v(f,u,v)

def G(f,u,v):
    return primitive_u(Gp,f,u,v)

def fonct_surf_reglee(f,u,v):   #Recomposition de l'hypothétique surface reglée sous forme f(u,v)=G(u)+v.K(u)
    return G(f,u,v) + v*K(f,u,v)

def verification_surface_reglee(f,Lu,Lv):   #Vérification de l'égalité de la fonction réelle et de la fonction recomposée pour tout u et tout v
    b=True
    for u in Lu:
        for v in Lv:
            if not(np.array_equal(np.around(f(u,v),preci),np.around(fonct_surf_reglee(f,u,v),preci))):
                b=False
#                print("Pas égalité pour u={} et v={}".format(u,v))
#                print(f(u,v))
#                print(fonct_surf_reglee(f,u,v))
            else:
#                print("OK pour {} ; {}".format(u,v))
                pass
    return b


#==============================================================================
# Essais
#==============================================================================
f = lambda u,v : f4(v,u) #Fonction à utiliser, renommée f()

#Affichage de la surface réelle et de la surface recomposée
#ax1 = Axes3D(plt.figure())
#ax2 = Axes3D(plt.figure())
#Lx = [[f(u,v)[0] for u in Lu] for v in Lv]
#Ly = [[f(u,v)[1] for u in Lu] for v in Lv]
#Lz = [[f(u,v)[2] for u in Lu] for v in Lv]
#ax1.plot_surface(Lx, Ly, Lz, rstride=1, cstride=1,cmap=cm.jet)
#
#fsurf = lambda u,v : fonct_surf_reglee(f,u,v)
#Lx = [[fsurf(u,v)[0] for u in Lu] for v in Lv]
#Ly = [[fsurf(u,v)[1] for u in Lu] for v in Lv]
#Lz = [[fsurf(u,v)[2] for u in Lu] for v in Lv]
#ax2.plot_surface(Lx, Ly, Lz, rstride=1, cstride=1,cmap=cm.jet)
#
#ax1.set_xlim(-5,5);ax1.set_ylim(-5,5);ax1.set_zlim(-5,5)
#ax2.set_xlim(-5,5);ax2.set_ylim(-5,5);ax2.set_zlim(-5,5)
#
#plt.show()

print(verification_surface_reglee(f,Lu,Lv)) #Affiche si la surface est réglée ou non

#==============================================================================
# Recherche du rayon de courbure principal
#==============================================================================

def deriveeUnSurRn(f,u,v,mu):   #Renvoi la dérivée de 1/Rn
    d1 = derive_u(f,u,v)
    d2 = derive_v(f,u,v)
    d3 = derive_u_u(f,u,v)
    d4 = derive_v_v(f,u,v)
    d5 = derive_u_v(f,u,v)
    H = np.array([d1[1]*d2[2]-d1[2]*d2[1],d1[2]*d2[0]-d1[0]*d2[2],d1[0]*d2[1]-d1[1]*d2[0]]) #produit vectoriel
    h = H/(np.sqrt(H[0]**2 + H[1]**2 + H[2]**2))
    L = h[0]*d3[0] + h[1]*d3[1] + h[2]*d3[2]
    M = h[0]*d5[0] + h[1]*d5[1] + h[2]*d5[2]
    N = h[0]*d4[0] + h[1]*d4[1] + h[2]*d4[2]
    E = d1[0]**2 + d1[1]**2 + d1[2]**2
    F = d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]
    G = d2[0]**2 + d2[1]**2 + d2[2]**2
    return ((2*M+2*N*mu)*(E+2*F+G*mu*mu)-(L+2*M*mu+N*mu*mu)*(2*F+2*G*mu))/((E+2*F*mu+G*mu*mu)**2)

def recherche_zero(Lx,Ly):  #Renvoie l'abscisse auquel s'annule la fonction
    for p in range(len(Lx)-1):
        if Ly[p]*Ly[p+1]<0:
            return Lx[p]

Lmu = np.linspace(-100,100,10000)   #Intervalle de recherche du zéro
Ly=[deriveeUnSurRn(f,3,-8,mu) for mu in Lmu]

def rayon_courbure_mini(f,Lu,Lv):   #Renvoie le rayon de courbure minimal de la surface pour tout u et tout v
    mini=100
    for u in Lu:
        for v in Lv:
            Ly=[deriveeUnSurRn(f,u,v,mu) for mu in Lmu]
            mu=recherche_zero(Lmu,Ly)
            if(mu<mini):
                mini=mu
    return mu

def usinable(f,Lu,Lv,d):
    mu=rayon_courbure_mini(f,Lu,Lv)
    if mu>=d:
        return("La surface est usinable par une fraise de {}".format(d))
    else:
        return("La surface n'est pas usinable par une fraise de {}".format(d))

print(usinable(f,Lu,Lv,4))