############  functions for performing cross correlations in R (splus)
###  source("/home/lees/Progs/R_stuff/HEAT.EQN.R")
###
erf<-function(x)
  {
    ## copied from phillpots
    
    t = 1/(1+0.47047*x)
    a1 = 0.34802
    a2 = -0.09587
    a3 = 0.74785

    y = 1 - (a1*t + a2*t^2 + a3*t^3)* exp(-(x^2))

    return(y)

  }
### j = seq(from=0, to=2.5, by=0.1)
###  y = erf(j)
###   plot(c(rev(-j), j),c(rev(-y),y))

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

    for(tim in t)
      {
        t1 = erf(x/(2*sqrt(k*tim)))
        t2 = rev(-t1)
        T = T0*(0.5+0.5*c(t2,t1))
        
        plot(c(rev(-x),x),T, xlab="Distance, m", ylab="Temp, K")
        abline(v=0, lty=2)
        atim = tim/(24*3600)
        
        title(paste(sep=' ', "time=",atim, " days") )
        locator()
      }
    return(T)
  }
###   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)

###  DO.DYKE()

###  source("/home/lees/Progs/R_stuff/HEAT.EQN.R")
###
DO.HALFSPACE<-function()
  {
    k = 1e-6
    T0 = 1000
    x = seq(0,20, length=1000)
    y =  heat.sol(x, T0, k, seq(from=0,to=1000*24*3600, by=100*24*3600))
  }
###

DO.DYKE<-function(a=a, x=x, t=t, k=k, T0=T0, NDIM=TRUE)
  {
    if(missing(a)) { a = 6 }
    if(missing(x)) { x = seq(0,20, length=1000) }
    if(missing(T0)) { T0 = 1000 }
    if(missing(k)) { k = 1e-6  }
    if(missing(t)) { t = seq(from=0,to=1000*24*3600, by=100*24*3600)  }
       if(missing(NDIM)) { NDIM = TRUE }

    TMAT = matrix(ncol=length(t), nrow=length(x))
    
    AX = a-x
    flag1 = AX<0
    
   
    AX = a+x
    flag2 = AX<0

    i = 0
    for(tim in t)
      {
        i = i+1
        tim = t[i]
        e1 = erf(abs(a-x)/(2*sqrt(k*tim)))
        e1[flag1] = -e1[flag1]
        
        
        
        e2 = erf(abs(a+x)/(2*sqrt(k*tim)))
        e2[flag2] = -e2[flag2]
    
        
        Temp  =  0.5*(e1 + e2)
        TMAT[,i] = Temp 
      }

    
   ###  matplot(x/a,TMAT, type='l')

    
    BMAT =   rbind( TMAT[ rev(1:length(x)), ], TMAT)

    EX = c( rev(-x/a), x/a)
   
    if(NDIM==TRUE)
      {
        matplot(EX , BMAT , type='l', xlab="Distance, x/a", ylab="Temperature, T/T0" )
      }
    else
      {
        matplot(a*EX , T0*BMAT , type='l', xlab="Distance, m", ylab="Temperature, K" )
      }
    grid(col='darkgray')

    px = rep(0,length(t))
    py = rep(0,length(t))
    
    tx = a*seq(from=0.01, to=.99, length=length(t))
    i =0
    for(tim in t)
      {
        i = i+1
        tim = t[i]
        e1 = erf(abs(a-tx[i])/(2*sqrt(k*tim)))
        e2 = erf(abs(a+tx[i])/(2*sqrt(k*tim)))
         py[i] =  0.5*(e1 + e2)
      }
    
   #   labs = k*t/a^2
    if(NDIM==TRUE)
      {
        labs = format.default(k*t/a^2, digits=3)
        text(tx/a, py, labels=labs)
      }
    else
      {
        tlabs = t/(24*3600)
        labs = format.default(tlabs, digits=3)
        text(tx, T0*py, labels=labs)

      }
  }

###  source("/home/lees/Progs/R_stuff/HEAT.EQN.R")

###    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)



###