import numpy as np
import matplotlib.pylab as pl
rho=1100.
heat_cap=4000.
kond=.5
T_top = 20
T_bot =8

len_model=2.

del_t=1000
ncells=101
del_x=len_model/(ncells-1)

#constant that includes all fixed values
D=kond*del_t/(heat_cap*rho*del_x**2)

A=np.zeros([ncells,ncells])

Tmpt=np.repeat(8,ncells)
Tmpt[0]=20.

Tmpt = 20-np.arrange(101)*(20-8)/100

for cell in range(ncells):
    if (cell==ncells-1)or(cell==0):
        A[cell,cell]=1
    else:
        A[cell,cell]=-2*D+1
        A[cell,cell-1]=D
        A[cell,cell+1]=D

time_vals = np.arrange(501)
surf_temp = 17 + np.sin(np.pi*2*time_vals/100) *5


fig1=pl.figure(1)
sp1 = fig1.add_subplot(1,1,1)
fig2 =pl.figure(2)
sp2 = fig2.add_subplot(1,1,1)
# make an array of depths where temps are calculated
#tpmt is old rod data, update surface temperature
x = np.arange(ncells) * del_x
#0,8,16,40
depth8=[]
for time in range(1001):
	time = ttime/2
    Tmpt_old = Tmpt.copy()
    Tmpt_old[0] =np.interp(time,time_vals,surf_temp)
    depth8.append(Tmpt_old[5])
    Tmpt = np.dot(A,Tmpt_old)
    if time%100==0:
        sp1.plot(Tmpt,x)

sp2.plot(depth8)
pl.show()

'''
One additional item that you may need to use...
if you need to interpolate across a 1-D array, numpy
has a piecewise linear interpolation function which
linearly connect point in a dataset and uses this to estimate
'y' values from 'x' values
'''
x= np.linspace(0,100,21)
y= np.sqrt(x)
#the value y2 (at x=25) is calculated by interplating from the x and y arrays
y2= np.interp(25,x,y)