k = 1*10^(-6)
dt = 3600
 dz = 20*10^(-2)

T0=25

T1 = 1200

##  initialize

###  at t=0

Ntsteps  =5
Nlays = 6
t = seq(from=1, to=Ntsteps)
lay =  seq(from=1, to=Nlays-1)
T = rep(1200,length=Nlays)
dTi = rep(0,length=Nlays)

celltemp = matrix(NA, nrow=Ntsteps+1, ncol=Nlays+2)

celltemp[1, ] = c(0, T0, T)


for(j in t)
{
for(i in lay)
{

if(i==1)
{
Tip1 = T[i+1]
Ti = T[i]
Tim1 = T0
}
else
{
Tip1 = T[i+1]
Ti = T[i]
Tim1 = T[i-1]
}

dTi[i] = k*dt*(Tip1-2*Ti+Tim1)/dz^2


}

T = T + dTi
print(paste(sep=" ", "##### TIME ######", j))
## print(dTi)
## print(T)


celltemp[j+1, ] = c(j, T0, T)

dTi = rep(0,length=Nlays)

}

colnames(celltemp)<-c("time", "surf", paste(sep="", "Cell", 1:Nlays))

print(celltemp)

#####################################
#####  D:\LEES\R_Programs\R_stuff
source("D:/LEES/R_Programs/R_stuff/HEAT.EQN.R")
source("c:/Documents and Settings/leesj/Desktop/GEOL_154/HEAT.EQN.R")


###  heat.sol<-function(x, T0, k, t)

get.heat<-function(x, T0, k, t)
  {
	tim = t
        t1 = erf(x/(2*sqrt(k*tim)))
        T = T0*(0.5+0.5*c(t1))
    return(T)
  }


###     we use get.heat2 because the 
##   temperature at the boundary remains constant
###    for all time
get.heat2<-function(x, T0, k, t)
  {
	tim = t
        t1 = erf(x/(2*sqrt(k*tim)))
        T = T0*(c(t1))
    return(T)
  }




##  distance in meters
x =  seq(from=0, to=80, by=0.5)/100

i = 1
Tx =  get.heat2(x, T1-T0, k, i*dt)
Tx = Tx+T0

plot(Tx, max(x)-x, type='n', xlim=c(700, 1200) , axes=FALSE, xlab="Temp", ylab="Depth, cm")
axis(3)
axis(2, at=pretty(x), labels=100*(max(x)-pretty(x)))
## points(Tx, max(x)-x)

for(j in 1:5)
{
Tx =  get.heat2(x, T1-T0, k, j*dt)
Tx = Tx+T0
lines(Tx, max(x)-x, lty=2, col=j)
}


z = c(0, (dz)*lay-dz/2)

for(i in 1:Ntsteps)
{

tem = celltemp[i+1,2:(Nlays+1)]
lines(tem, max(x)-z, type='b', col=i)

}

grid()