import numpy as np
import matplotlib.pyplot as plt

global muVacio
global muFeRel
global kDisp
global cteRoz
global g
global resCu

densidadHierro = 7800 # kg/m3
muVacio = 4*np.pi*1e-7
muFeRel = 5000
resCu = 1.7e-8 # ohm/m
kDisp = 2 # cte de dispersion
cteRoz = 0.1
g = 9.8

__author__ = 'Oscar Suescun' 

def calcular_diametro_ext(longitudBobina, diametroCuBobina, espirasBobina, diametroInteriorBobina):
    
    espirasCapa = int(longitudBobina/diametroCuBobina)
    numeroCapas = int(espirasBobina/espirasCapa + .5)

    return 2 * numeroCapas * diametroCuBobina + diametroInteriorBobina

def reluctancia_funcX(longFe, diamFe, longC, diamCint, diamCu, espiras):
    '''
    Esta funcion devuelve un vector que contenga la reluctancia total en
    funcion de la distancia que recorre el vastago dentro de la bobina.

    Para esto se utiliza un vector de 0 -> longC con un estep de longC/100

    Entrada:
        - longFe --> longitud del vastago
        - diamFe --> diametro del vastago
        - longC --> longitud de la bobina
        - diamCint --> diametro interior de la bobina
        - diamCu --> diametro del hilo de cobre
        - espiras --> numero de espiras en la bobina

    Salida:
        - relTotal --> devuelve un vector con la reluctancia total en 
                       funcion de la distancia recorrida en la bobina.

    '''

    diamCext = calcular_diametro_ext(longC, diamCu, espiras, diamCint)

    seccionFe = (np.pi * diamFe**2) / 4
    seccionSc = (np.pi * diamCext**2) / 4
    seccionDisp = (np.pi * (kDisp * diamCint)**2) / 4

    x = np.arange(0,longC + longC/100, longC/100)
    longFe = np.full(len(x), longFe)
    longC = np.full(len(x), longC)


    relFe = longFe / (muVacio * muFeRel * seccionFe)
    relDisp = longC / (muVacio * seccionDisp)
    relAire = (longC - x) / (muVacio * seccionSc)
    relFlujo = (longC + longFe - x) / (muVacio * seccionDisp)

    relTotal = relFe + relDisp + relAire + relFlujo

    return relTotal


def induccion(corriente, reluctancia, espiras, diamFe):
    flujo = espiras * corriente / reluctancia

    secFe = np.pi * diamFe**2 / 4

    return flujo / secFe

def fuerza_magnetica(corriente, reluctancia, espiras, diamFe):
    '''
    Solo se calcula la fuerza magetica

    '''

    b = induccion(corriente, reluctancia, espiras, diamFe)

    secFe = np.pi * diamFe**2 / 4

    return  np.sign(corriente) * 0.5 * b**2 * secFe / muVacio

def calcular_aceleracion(masa, anguloDisparo, fuerzaMagnetica):
    '''
    Retocar esto, tiene en cuenta el principio solo pero hay cosas raras
    
    '''

    anguloDisparo = anguloDisparo * np.pi / 180
    acel = fuerzaMagnetica/masa - g*(cteRoz * np.cos(anguloDisparo) + np.sin(anguloDisparo))
    return acel
    
    
def calcular_resistencia(diamCu, longC, diamC, espiras):

    espirasCapa = int(longC/diamCu)
    nCapasEnteras = int(espiras/espirasCapa) # solo capas completametne llenas
    espirasUltimaCapa = espiras - espirasCapa * nCapasEnteras # las espiras de la ultima capa

    longitud = 0

    for i in range(nCapasEnteras): longitud += espirasCapa*np.pi*(diamC + (1+i)*diamCu) # Calculo longitud de total de las capas llenas
    
    longitud += espirasUltimaCapa*np.pi*(diamC + (1 + nCapasEnteras)*diamCu)
    
    seccionCu = np.pi * diamCu**2 / 4

    return  resCu * longitud / seccionCu
        
def calcular_masa(diamFe, longFe): return np.pi * diamFe**2 / 4 * longFe * densidadHierro



if __name__ == '__main__':
    from funcionesSimELectrica import *
    from scipy.integrate import cumulative_trapezoid

    tiempoSimulacion = 10e-3
    step = tiempoSimulacion / 1e6

    numeroModulos = 2

    longitudBobina = 53.21e-3 # m
    diametroInteriorBobina = 6.035e-3 * 2 # m
    espirasBobina = 500
    diametroCuBobina = 0.8e-3 # m

    longitudVastago = 96e-3 # m
    diametroVastago = 3.045e-3 * 2 # m
        
    Tension = 30 # V
    Capacitancia = 500e-3 # F

    masa = calcular_masa(diametroVastago, longitudVastago)
    anguloDisparo = 40


    resistencia = calcular_resistencia(diametroCuBobina, longitudBobina,
                                       diametroInteriorBobina, espirasBobina)
    
    reluctancia = np.mean(reluctancia_funcX(longitudVastago, diametroVastago,
                                    longitudBobina, diametroInteriorBobina,
                                    diametroCuBobina, espirasBobina))
    
    inductancia = espirasBobina**2 / reluctancia

    tiempo = np.arange(0,tiempoSimulacion, step)

    corriente = solver_manual(Tension, resistencia, inductancia, Capacitancia, tiempo)
    fuerza = fuerza_magnetica(corriente, reluctancia, espirasBobina, diametroVastago)
    aceleracion = calcular_aceleracion(masa, anguloDisparo, fuerza)
    aceleracion[aceleracion<0]=0
    velocidad = cumulative_trapezoid(aceleracion, tiempo, initial=0)
    posicion = cumulative_trapezoid(velocidad, tiempo, initial=0)

    idx = posicion <= longitudBobina

    tiempo = tiempo[idx]
    posicion = posicion[idx]
    velocidad = velocidad[idx]
    aceleracion = aceleracion[idx]

    plt.figure()
    #plt.plot(t, corriente, label='corriente')
    #plt.plot(t, fuerza, label = 'Fuerza')
    plt.plot(tiempo, aceleracion, label='aceleracion')
    # plt.plot(tiempo, velocidad, label = 'velocidad')
    #plt.plot(posicion, velocidad, label = 'posicion')
    plt.legend()
    plt.show()