c1 = [1]*100 + [0]*100
nc1=c1
DtildeList=[]
DCuList=[]
DNiList=[]
DCuS=2.03
DtildeS=1.45
DNiS=2
 
 
dt = 100
dx=0.01
xlist=np.linspace(-1,1,200)
 
 
def DtildeL(X,DA,DB):
    return X*DB+(1-X)*DA


    
 
Xlist=np.linspace(0,1,200)
for i in range(200):
    DCuList.append(10**(-12.92-Xlist[i]*DCuS))
    DNiList.append(10**(-13.5-Xlist[i]*DNiS))
    DtildeY=DtildeL(Xlist[i], DNiList[i],DCuList[i])
    DtildeList.append(DtildeY)
 
def DtildeFunc(X):
    return DtildeL(X,10**(-12.92-(1-X)*DCuS),10**(-13.5-(1-X)*DNiS)) #mm^2/sec
   
TlistN=np.linspace(973,1273,1080)
for k in range(10800):
   
    for i in range(1, len(c1)): #This loop calculates the numerical solution for the concentration
        if i== len(c1)-1:
            nc1[i] =c1[i] +(dt * DtildeFunc(c1[i])*10**6)/(dx**2) * (c1[i] - 2 * c1[i] + c1[i-1])+((DtildeFunc(c1[i])*10**6 - DtildeFunc(c1[i-1])*10**6)*(c1[i]-c1[i-1])*(1/(dx**2)))
        else:
            nc1[i] = c1[i]+(dt * DtildeFunc(c1[i])*10**6)/(dx**2) * (c1[i+1] - 2 * c1[i] + c1[i-1])+((DtildeFunc(c1[i])*10**6 - DtildeFunc(c1[i-1])*10**6)*(c1[i]-c1[i-1])*(1/(dx**2)))
    c1=nc1.copy()
    if k==1080:
        c30=c1
    if k==3600:
        c100=c1
    if k==7200:
        c200=c1

plt.figure(figsize=(12,8))
plt.plot(xlist,c30, label="Diffusion profile @30h")
plt.plot(xlist,c100, label="Diffusion profile @100h")
plt.plot(xlist,c200, label="Diffusion profile @200h")
plt.plot(xlist,c1, label="Diffusion profile @300h")
plt.legend()
plt.ylabel('Concentration of Cu',fontsize=12)
plt.xlabel('position [mm]', fontsize = 12)
plt.title(r'Diffusion profile with use of $\tilde{D}$')
plt.grid()

;