%matplotlib inline

import scipy as np
from matplotlib import pyplot as plt

"""
Explanation: Difusión en una dimensión
End of explanation
"""


data = np.loadtxt('initial_f.csv' , delimiter=',' )

x = data[ : , 0 ]
u0 = data[ : , 1 ]

L  = x[-1] - x[0]      # longitud del sistema 1D
nx = x.size            # nodos espaciales
dx = x[ 1] - x[0]      # equiespaciados

T= 0.03            # tiempo total
dt = 1e-5     # intervalo temporal
nt = int(T / dt)            # pasos temporales


D = 1              # difusividad

"""
Explanation: Parámetros
End of explanation
"""


Co = D * dt / dx**2
Co

"""
Explanation: ¡Número difusivo
End of explanation
"""


plt.plot( x , u0 , 'o-')

"""
Explanation: Condiciones iniciales
End of explanation
"""


u = np.zeros( nx + 2  )

u[ 1 : nx +1 ] = u0

nnx = u.size

un = u.copy()         # distribución actual

for i in range( 1 , nx +1 ):
   u[i] = un[i] + Co * (un[i+1] + un[i-1] - 2 * un[i] )

        

"""
Explanation: Un paso en el tiempo
Recordemos que queremos implementar $u_i^{n+1} = u_i^n + \mathrm{Co} (u_{i+1}^n +  u_{i-1}^n  - 2 u_i^n)$
Aumentamos nuestro sistema con dos extremos adicionales, que ponemos 0 (condiciones de Dirichlet homogéneas)
End of explanation
"""


uf = u[ 1 : nx +1 ]

plt.plot( x , uf)

"""
Explanation: El vector uf recibe la solución final, y le quitamos los dos extremos:
End of explanation
"""


u[ 0 ] = u[ nx + 1 ] = 0
u[ 1 : nx +1 ] = u0


for n in range(nt ):
    un = u.copy()
    for i in range( 1 , nx +1 ): 
           u[i] = un[i] + (Co / 2.0) * (un[i+1] + un[i-1] - 2 * un[i] )

        

uf = u[ 1 : nx +1 ]

plt.plot(x , uf , x , u0 , 'r')

"""
Explanation: Tiempo completo
End of explanation
"""


import scipy.integrate as integrate

I0 = integrate.trapz( x, u0)
print( I0 )

integrate.trapz( x, uf)

"""
Explanation: Es interesante comprobar si se pierde "masa" (energía... lo que sea )
End of explanation
"""


u[ 0 ] = u0[0]
u[-1]  = u0[-1]
u[ 1 : nx +1 ] = u0

I = [I0]

for n in range( nt ):
    un = u.copy()
    un [0]  = un[1]   # derivada cero a la izquierda
    un [-1] = un[-2]  # derivada cero a la derecha
    
    # una alternativa sería hacer un bucle de 2 a nx sólo, con expresiones
    # especiales para 1 y nx+1
    for i in range( 1 , nx +1 ): 
           u[i] = un[i] + (Co / 2.0) * (un[i+1] + un[i-1] - 2 * un[i] )
    
    I.append(  integrate.trapz( x, u[1 : nx+1 ] ) )

uf = u[ 1 : nx +1 ]

plt.plot(x , uf , x , u0 , 'r')

integrate.trapz( x, uf )

I0

plt.plot( (-1)*np.array(I) )

"""
Explanation: ... probablemente, se pierde por los bordes sobre todo
Condiciones "adiabáticas" (Neumann homogéneas)
End of explanation
"""


