D:\Lees\CLASSES
source("D:/Lees/CLASSES/HEAT.EQN.R")


###  gaussian
### postscript(file="gauss_err.ps", paper="letter", horizontal=TRUE, print.it=FALSE,  onefile=FALSE)

j = seq(from=-2.5, to=2.5, by=0.01)
gaus = exp(-j^2)
plot(j,gaus, type='n', xlab='x', ylab="gauss")
lines(j,gaus)
lines(j[j<=1.2],gaus[j<=1.2], type='h', col=gray(.75) )

lines(j[j>=0&j<=1.2],gaus[j>=0&j<=1.2], type='h')
title(main="Illustration of Gaussian and Error Function")

dev.off()

################################
################################
################################

j = seq(from=0, to=2.5, by=0.1)
y = erf(j)
plot(j,y, type='n', xlab='x', ylab="erf")
abline(v=0, h=0, col=gray(.8) )
lines(c(rev(-j), j),c(rev(-y),y), type='b')
title(main="Error Function")


################################
################################

###   erf is anti-symmetirc about the origin
plot(c(rev(-j), j),c(rev(-y),y), type='n', xlab='x', ylab="erf")
abline(v=0, h=0, col=gray(.8) )
lines(c(rev(-j), j),c(rev(-y),y), type='b')
title(main="Error Function")

################################

###   erfc is the complement of erf:   erfc=1-erf


plot(c(rev(-j), j),1-c(rev(-y),y), type='n', xlab='x', ylab="erfc")
abline(v=0, h=0, col=gray(.8) )
lines(c(rev(-j), j),1-c(rev(-y),y), type='b')
title(main="Error Function Complement")



################################
################################
  k = 1e-6
  t  = 1*365*24*60*60
   T0 = 1000
   x = seq(0,20, length=1000)
   T =  heat.sol(x, T0, k, t)


   T =  heat.sol(x, 1000, 10e-6, 3600)

  y =  erf(x/(2*sqrt(k*t)))

   y =  heat.sol(x, T0, k, seq(from=0,to=100*24*3600, by=10*24*3600))

   y =  heat.sol(x, T0, k, seq(from=0,to=1000*24*3600, by=100*24*3600))
    

   y =  heat.sol(x, T0, k, t)
   plot(c(rev(-x),x),y)

##########################
    a = 10; k=10^6

   DO.DYKE(NDIM=TRUE)
   DO.DYKE(NDIM=FALSE)
   DO.DYKE(t = seq(from=0,to=500*24*3600, by=50*24*3600) , NDIM=FALSE)
   DO.DYKE(t = seq(from=0,to=500*24*3600, by=50*24*3600) , NDIM=TRUE)
   DO.DYKE(t = a^2*c(0, 0.05, 0.2, 0.5, 2.0)/k  , NDIM=TRUE)
   DO.DYKE(t = a^2*c(0, 0.05, 0.2, 0.5, 2.0)/k  , NDIM=FALSE)