import numpy as np
import matplotlib.pyplot as plt
import math
import cmath
im = complex(0.0,1.0)
# SS Stainless Steel
kappa_SS = 22
rho_SS = 7200
c_SS = 700
D = kappa_SS/(rho_SS*c_SS)
f = 100
w = 2.0*math.pi*f
gama2 = im*(w/D)
gama = cmath.sqrt(gama2)
mu = math.sqrt((2*D)/w)
Nx = 100
xi = -10*mu
xf = 10*mu
dx = (xf-xi)/Nx
x = np.zeros(Nx+1)
TTn = np.zeros(Nx+1,dtype=complex)
TTap = np.zeros(Nx+1,dtype=complex)
TTam = np.zeros(Nx+1,dtype=complex)
for ix in range(0,Nx+1):
x[ix] = xi+ix*dx
TTap[ix]= cmath.exp(+im*gama*x[ix])
TTam[ix]= cmath.exp(-im*gama*x[ix])
plt.xlim(0,10)
plt.ylim(-1,1)
plt.plot(x/mu,TTap.real,'-b',x/mu,TTap.imag,'--b')
Nt = 10000;
dt = 0.2*(dx*dx)/(2.0*D);
t = np.zeros(Nt+1)
fuente = np.zeros(Nt+1)
for it in range(0,Nt+1):
t[it] = it*dt
fuente[it]= math.sin(w*t[it])
plt.plot(t,fuente)
ic1=0
ic2=0
for it in range(0,Nt+1):
for ix in range(1,Nx):
TTn[ix]= TTn[ix] + ((D*dt)/(dx*dx))*(TTn[ix+1]-2.0*TTn[ix]+TTn[ix-1])
TTn[int(Nx/2)]=fuente[it]
if ic1==0:
ic2=ic2+1
print(ic2)
plt.ylim(-1,1)
plt.xlim(0,10)
plt.plot(x/mu,TTn.real,'-b',x/mu,TTap.real,'--b')
if ic2 < 10: plt.savefig("T_0000{}.png".format(ic2))
if ic2 >= 10 and ic2<100: plt.savefig( "T_000{}.png".format(ic2))
if ic2 >= 100 and ic2<1000: plt.savefig( "T_00{}.png".format(ic2))
if ic2 >= 1000 and ic2<10000: plt.savefig( "T_0{}.png".format(ic2))
if ic2 >= 10000 and ic2<100000: plt.savefig( "T_{}.png".format(ic2))
plt.close()
ic1=ic1+1
if ic1==50: ic1=0
