.packageName <- "spatialCovariance"
oneByoneHack <- function(info)
  {
    info$rowsep <- numeric(0)
    info$colsep <- numeric(0)
    info$nrows <- 1
    info$ncols <- 1
    info$lengths <- info$lengths[1]
    info$rowReps <- info$rowReps[1]
    info$locations <- matrix(info$locations[1,],1,)
    info$indices <- info$indices[1:3,]
    info$indices.preLimits$ax <- info$indices.preLimits$ax[1]
    info$indices.preLimits$bx <- info$indices.preLimits$bx[1]
    info$indices.preLimits$getV.i <- info$indices.preLimits$getV.i[1]
    info$indices.preLimits$getV.j <- info$indices.preLimits$getV.j[1]
    info$indices.preLimits$indices <- matrix(info$indices.preLimits$indices[1,],1,)
    info$indices.preLimits$evalFactor <- info$indices.preLimits$evalFactor[1]
    info
  }
defineK <- function(class)
  {
    ## current set of covariance functions
    ## analytic results exist for the power and ldt models
    ## the powerNI and ldtNI are for numerical integration
    
    if(class=="matern") K <- function(d,params) (params[1]*d)^params[2]*besselK(params[1]*d,params[2])
    if(class=="bess1") K <- function(d,params) (params[1]*d)*besselK(params[1]*d,1)
    if(class=="exp") K <- function(d,params) exp(-params[1]*d)
    ##if(class=="exp") K <- function(d,params) sqrt(params[1]*d)*besselK(params[1]*d,0.5)
    if(class=="bess0") K <- function(d,params) besselK(params[1]*d,0)

    if(class=="powerNI") K <- function(d,params) d^(-params[1])/params[1]
    ## params < 2 above
    if(class=="powerNI") K <- function(d,params) -d^(2*params[1])/params[1]
    ## range for params is params > -1 in power model

    if(class=="spline") K <- function(d,params) d^2*log(d)
    if(class=="ldtNI") K <- function(d,params) -log(d)
    if(class=="cauchy") K <- function(d,params) 1/(params[2])*(1+abs(params[1]*d)^params[3])^(-params[2]/params[3])
    ## params are (lambda, beta, alpha)

    cov.f <- function(d,rw,cw,ax,bx,i,j,params,K)
      {
        K(d,params)*f(d,rw,cw,ax,bx,i,j)
      }


    ## derivatives of covariance functions
    if(class=="dmatern1") K <- function(d,params) d * (d*params[1])^(-1+params[2])*params[2]*besselK(params[1]*d,params[2])+d/2*(d*params[1])^(params[2])*(-besselK(params[1]*d,-1+params[2])-besselK(params[1]*d,1+params[2]))

    if(class=="dbess1") K <- function(d,params) d*besselK(params[1]*d,1)+d/2*(d*params[1])*(-besselK(params[1]*d,0)-besselK(params[1]*d,2))
    if(class=="dexp") K <- function(d,params) -d*exp(-params[1]*d)
    if(class=="dexp2") K <- function(d,params) d*d*exp(-params[1]*d)

    if(class=="dbess0") K <- function(d,params) -d*besselK(params[1]*d,1)
    if(class=="dpowerNI") K <- function(d,params) d^(2*params[1])/params[1]^2-2*d^(2*params[1])*log(d)/params[1]

    list(K=K,cov.f=cov.f)
  }


Hypergeometric2F1 <- function(a,b,c,x)
  {

    ## Method 3
    ## numerical integration of results on pages 231 and 233 of Mathai 1999
    ## only used when b = 1/2, c=3/2
    ## or when b = 1, c = 5/2
    ret <- NA
    if(is.nan(x)) {
      ret <- 0
    } else if(x==0) ret <- 1
    if(is.na(ret))
      {
        ## this seems to be the best, in terms of speed and accuracy
        if(b==1 && c==5/2) {
          ## use Mathai 1999, page 231, lemma 2.4.1
          integrand1 <- function(u) u^(-2*a+1)*sqrt(u^2-(1-x))
          result <- integrate(integrand1, lower=sqrt(1-x),upper=1,rel.tol=1e-12)$value
          result <- result*3/(x*sqrt(x))
          ret <- result
        } else if(b==1/2 && c==3/2)
          {
            if(x==1) {
              ret <- sqrt(pi)*gamma(1-a)/(2*gamma(3/2-a))
            } else {
              ## use Mathai 1999, page 233, lemma 2.4.2
              integrand2 <- function(u) u^(2*a-2)*asin(sqrt(1-x)/u)
              result <- integrate(integrand2, lower=sqrt(1-x),upper=1,rel.tol=1e-12)$value
              
              result <- result*(2*a-1)/sqrt(1-x)^(2*a-1)+pi/2-1/sqrt(1-x)^(2*a-1)*asin(sqrt(1-x))
              result <- result/sqrt(x)
              ret <- result
            }
          }
      }
    ##cat("hypgeo =",ret,"\n")
    ##if(is.nan(ret)) cat(a,b,c,x,"error\n")
    if(is.nan(ret)) ret <- 0  ## if this happens 2F1 will be multiplied by 0 anyway
    ret
}
f.anal.ldt <- function(coords,info)
  {
    ## Copied over from ldt-analytic.r on Dec 10th
    ## trying to incorporate apply into my code
    
    ## created May 12th
    ## analytic results for the ldt model
    ## see write up and notes in binder for more details.

    a <- info$rowwidth
    b <- info$colwidth
    ax <- coords[1]
    bx <- coords[2]
    i <- coords[3]
    j <- coords[4]
    ## 5
    ## 6
    ## 7 - which rows to evaluate - evalFactor
    ## 8 - lengths
    ## 9 - indicator for which rows to evaluate analytic results
    evaluate <- coords[9]
    ## 10 - lower limit for numerical integration
    ## 11 - upper limit for numerical integration

    ret <- 0
    if(evaluate) {
      if(i !=1 && j !=1)
        {
          ## most common situation:  different row, different column
          ## simplified different row, column result
          ## this is bad due to 0/0 and log(0) see below
          ## also its old, note the double for loops
          if(FALSE) {
            as <- c(-ax,ax-a,ax,ax+a)
            bs <- c(-bx,bx-b,bx,bx+b)
            ret <- 100*a^2*b^2
            for(ii in 1:4)
              for(jj in 1:4)
                {
                  ret <- ret + (-1)^(ii+jj+1)*(8*as[ii]*bs[jj]*(bs[jj]^2*atan(as[ii]/bs[jj])+as[ii]^2*atan(bs[jj]/as[ii]))-(as[ii]^4-6*as[ii]^2*bs[jj]^2+bs[jj]^4)*log(as[ii]^2+bs[jj]^2))
                } 
            ret <- ret/(48*a^2*b^2)
          }
          
          ## different version, takes care of 0/0, or log(0)
          ## doesn't have for loops
          ret <- 100*a^2*b^2
          as <- c(-ax,ax-a,ax,ax+a)
          bs <- c(-bx,bx-b,bx,bx+b)
          asbs <- cbind(as.vector(matrix(as,4,4)),as.vector(matrix(1:4,4,4)),as.vector(t(matrix(bs,4,4))),as.vector(t(matrix(1:4,4,4))))
                  
          ret <- 100*a^2*b^2
          ret <- ret + sum(apply(asbs,1,sepRowsepCol.ldt.res))
          ret <- ret/(48*a^2*b^2)        
        } else if(i == 1 && j !=1)
          {
            ## same column, different row
            ax <- bx
            bx <- 0
            a <- b
            b <- info$rowwidth
            ret <- (50*a^2*b^2-8*b^3*(ax-a)*atan((ax-a)/b)+
                    16*ax*b^3*atan(ax/b)-8*b^3*(ax+a)*atan((ax+a)/b)+
                    16*ax^3*b*atan(b/ax)-8*b*(ax-a)^3*atan(b/(ax-a))-
                    8*b*(ax+a)^3*atan(b/(ax+a))+4*ax^4*log(ax)-
                    2*(ax+a)^4*log(ax+a)+
                    (((ax-a)^2-b^2)^2-4*b^2*(ax-a)^2)*log((ax-a)^2+b^2)-
                    2*((ax^2-b^2)^2-4*b^2*ax^2)*log(ax^2+b^2)+
                    (((ax+a)^2-b^2)^2-4*b^2*(ax+a)^2)*log((ax+a)^2+b^2))
            if(ax>a) ret <- ret - 2*(ax-a)^4*log(ax-a)
            ret <- ret/(24*a^2*b^2)
          } else if(i!=1 && j==1)
            {
              ## simplified same row result using mathematica
              ## problems when ax=a with log(ax-a) dealt with below
              ret <- (50*a^2*b^2-8*b^3*(ax-a)*atan((ax-a)/b)+
                      16*ax*b^3*atan(ax/b)-8*b^3*(ax+a)*atan((ax+a)/b)+
                      16*ax^3*b*atan(b/ax)-8*b*(ax-a)^3*atan(b/(ax-a))-
                      8*b*(ax+a)^3*atan(b/(ax+a))+4*ax^4*log(ax)-
                      2*(ax+a)^4*log(ax+a)+
                      (((ax-a)^2-b^2)^2-4*b^2*(ax-a)^2)*log((ax-a)^2+b^2)-
                      2*((ax^2-b^2)^2-4*b^2*ax^2)*log(ax^2+b^2)+
                      (((ax+a)^2-b^2)^2-4*b^2*(ax+a)^2)*log((ax+a)^2+b^2))
              if(ax>a) ret <- ret - 2*(ax-a)^4*log(ax-a)
              ret <- ret/(24*a^2*b^2)
            } else {
              ## diagonal result, the most rare
              ret <- (25*a^2*b^2-8*a*b^3*atan(a/b)-8*a^3*b*atan(b/a)-2*a^4*log(a)-2*b^4*log(b)+(a^4-6*a^2*b^2+b^4)*log(a^2+b^2))/(12*a^2*b^2)
            }
    }
    ret
  }
   
sepRowsepCol.ldt.res <- function(abVal)
  {
    aVal <- abVal[1]
    ii <- abVal[2]
    bVal <- abVal[3]
    jj <- abVal[4]
    
    numSteps <- 2-((aVal==0) + (bVal==0))
    ret <- 0
    if(numSteps>=1) ret <- ret + (-1)^(ii+jj+2)*(aVal^4-6*aVal^2*bVal^2+bVal^4)*log(aVal^2+bVal^2)
    if(numSteps>1) ret <- ret + (-1)^(ii+jj+1)*(8*aVal*bVal*(bVal^2*atan(aVal/bVal)+aVal^2*atan(bVal/aVal)))
    ret
  }
f.NI <- function(coords,params,eps,K,cov.f,info)
  {
    a <- info$rowwidth
    b <- info$colwidth
    ax <- coords[1]
    bx <- coords[2]
    i <- coords[3]
    j <- coords[4]
    ## 5 
    ## 6
    ## 7 - which rows to evaluate - evalFactor
    ## 8 - lengths
    ## 9 - indicator for which rows to evaluate analytic results
    minD <- coords[10]
    maxD <- coords[11]

    if(minD < maxD) {
      int.res <- try(integrate(cov.f,lower=minD, upper=maxD,rw=a,cw=b,ax=ax,bx=bx,i=i,j=j,rel.tol=eps,params=params,K=K,stop.on.error=F))
    } else { int.res <- list(value=0) }
    
    if(class(int.res)=="try-error") {

      cat("integration fail for plot",i,j,"with limits of integration",minD,maxD,"corresponding to",integrate(f,lower=minD, upper=maxD,rw=a,cw=b,ax=ax,bx=bx,i=i,j=j,rel.tol=eps)$value,"of the density and over this region the covariance function values are between",cov.f(minD,rw=a,cw=b,ax=ax,bx=bx,i=i,j=j,params=params,K=K),"and",cov.f(maxD,rw=a,cw=b,ax=ax,bx=bx,i=i,j=j,params=params,K=K),"\n")

      if(FALSE) {
        par(mfrow=c(3,2))
        gVals <- seq(0,1,length=101)
        gVals <- c(gVals^2,1-gVals^2)
        gVals <- gVals*(maxD-minD)+minD
        gVals <- sort(gVals)
        plot(gVals,K(gVals,params=params),xlab="",ylab="",main="Cov Function K")
        plot(gVals,f(gVals,rw=a,cw=b,ax=ax,bx=bx,i=i,j=j),xlab="",ylab="",main="Density f")
        plot(gVals,cov.f(gVals,rw=a,cw=b,ax=ax,bx=bx,i=i,j=j,params=params,K=K),xlab="",ylab="",main="K times f")
        gVals <- seq(0,1,length=101)
        gVals <- c(gVals^2,1-gVals^2)
        gVals <- gVals*(maxD-minD+4)+minD-4
        gVals <- sort(gVals)
        plot(gVals,K(gVals,params=params),xlab="",ylab="",main="Cov Function K, larger range")
        plot(gVals,f(gVals,rw=a,cw=b,ax=ax,bx=bx,i=i,j=j),xlab="",ylab="",main="Density f, larger range")
        plot(gVals,cov.f(gVals,rw=a,cw=b,ax=ax,bx=bx,i=i,j=j,params=params,K=K),xlab="",ylab="",main="K times f, larger range")
        locator(1)
      }
      
      int.res <- list(NULL)
      int.res$value <- 0
    }
    int.res$value
  }
## taken from numerics/Density/explicit.densities.r
## updated on January 9th 2003

f <- function(d, rowwidth, colwidth, ax, bx, i, j)
  {

    ## calculates the density of the distance (not its sqaure)
    ## one point is in rectangle of sides (rw by cw)
    ## rectangle is at location (i=1,j=1)
    ## rectangles have rs and cs seperations between rows and columns
    ## second point is in rectangle (i,j)

    ## depending on the values of i and j there are three different algorithms
    ## to compute the density of uu, the first is a formula calculated by Ghosh
    ## the other two take his technique and apply it more generally

    ## all three densities fcase1, fcase2 and fcase3 are for the square of the distance
    ## that is why we pass uu^2 and multiply the result by 2*uu

    ## only pass vectors uu with values on the postive real axis

    ## pass the dimensions of the plot and the coordinates of the lower left hand corner of the second plot, see report - July 20th 2003

    ## TP for temporary parameter
    uu <- d
    TPa <- rowwidth
    rw <- rowwidth
    TPb <- colwidth
    cw <- colwidth
    TPr <- ax
    TPs <- bx

    if(i != 1 && j != 1) {
      ret <- fcase3(uu^2,TPa=TPa,TPb=TPb,TPr=TPr,TPs=TPs)*2*uu
    } else if(j==1 && i!=1) {
      ret <- fcase2(uu^2,TPa=TPa,TPb=TPb,TPr=TPr)*2*uu
    } else if(j!=1 && i==1) {
      ret <- fcase2(uu^2,TPa=cw,TPb=rw,TPr=bx)*2*uu
    } else if(i==1 && j==1) {
      ret <- fcase1(uu^2,TPa=TPa,TPb=TPb)*2*uu
    }
    
    ret
  }

h <- function(uVect,alpha)
  {
    uTEMP <- uVect[uVect<=alpha^2]
    ret <- rep(alpha*pi/2,length(uTEMP))
    uTEMP <- uVect[uVect>alpha^2]
    ret <- c(ret,2*sqrt(uTEMP-alpha^2)-alpha*asin(1-2*alpha^2/uTEMP))
    ret
  }

## case 1

fcase1 <- function(uVect,TPa,TPb)
  {
    ## density for distance squared
    ## Ghosh, Bull Calcutta
    ## assumes a <= b

    lim1 <- 0
    lim2 <- min(TPa^2,TPb^2)
    lim3 <- max(TPa^2,TPb^2)
    lim4 <- TPa^2+TPb^2

    ret <- NULL
    uTEMP <- uVect[uVect < lim1]
    ret <- c(ret,rep(0,length(uTEMP)))
    uTEMP <- uVect[uVect >= lim1 & uVect < lim2]
    ret <- c(ret,TPa*TPb*pi-2*sqrt(uTEMP)*(TPa+TPb)+uTEMP)
    uTEMP <- uVect[uVect >= lim2 & uVect < lim3]
    ret <- c(ret,TPa*TPb*pi/2-lim2-2*max(TPa,TPb)*sqrt(uTEMP)+max(TPa,TPb)*h(uTEMP,min(TPa,TPb)))
    uTEMP <- uVect[uVect >= lim3 & uVect < lim4]
    ret <- c(ret,TPa*h(uTEMP,TPb)+TPb*h(uTEMP,TPa)-lim2-lim3-uTEMP)
    uTEMP <- uVect[uVect >=lim4]
    ret <- c(ret,rep(0,length(uTEMP)))
    
    ret <- ret/(lim2*lim3)
    ret
  }

## case 2
f2.1 <- function(u,TPa,TPb,TPr)
  {
    -(TPr-TPa)*TPb*pi/2-(TPr-TPa)^2-u+2*(TPr-TPa)*sqrt(u)+TPb*h(u,TPr-TPa)
  }

f2.2 <- function(u,TPa,TPb,TPr)
  {
    ((TPr+TPa)*TPb*pi/2 + TPr^2 + 2*TPr*TPa-TPa^2) + u - 2*(TPr+TPa)*sqrt(u) - 2*TPb*h(u,TPr)+TPb*h(u,TPr-TPa)
  }

f2.3 <- function(u,TPa,TPb,TPr)
  {
    -2*TPa^2-2*TPb*h(u,TPr)+TPb*h(u,TPr-TPa)+TPb*h(u,TPr+TPa)
  }

f2.4 <- function(u,TPa,TPb,TPr)
  {
    (TPr^2-2*TPr*TPa-TPa^2+TPb^2)+u-(TPr-TPa)*h(u,TPb)-2*TPb*h(u,TPr)+TPb*h(u,TPr+TPa)
  }

f2.5 <- function(u,TPa,TPb,TPr)
  {
    -(TPr+TPa)^2-TPb^2-u+(TPr+TPa)*h(u,TPb)+TPb*h(u,TPr+TPa)
  }

f2.6 <- function(u,TPa,TPb,TPr)
  {
    (TPb^2+(TPr+TPa)*TPb*pi/2+2*TPr^2)+2*u-2*(TPr+TPa)*sqrt(u)-2*TPb*h(u,TPr)-(TPr-TPa)*h(u,TPb)
  }

f2.7 <- function(u,TPa,TPb,TPr)
  {
    ((TPr+TPa)*TPb*pi/2-TPb^2)-2*(TPr+TPa)*sqrt(u)+(TPr+TPa)*h(u,TPb)
  }

f2.8 <- function(u,TPa,TPb,TPr)
  {
    (-(TPr-TPa)*TPb*pi/2+TPb^2)+2*(TPr-TPa)*sqrt(u)-(TPr-TPa)*h(u,TPb)
  }

fcase2 <- function(uVect,TPa,TPb,TPr)
  {
    ret <- NULL

    ## in the notes these are u1, u2, u3, u4, u5, u6
    ## but this will be confusing with uVect
    lim1 <- (TPr-TPa)^2
    lim2 <- TPr^2
    lim3 <- (TPr+TPa)^2
    lim4 <- (TPr-TPa)^2+TPb^2
    lim5 <- TPr^2+TPb^2
    lim6 <- (TPr+TPa)^2+TPb^2
    
    if(lim3<=lim4)
      {
        uTEMP <- uVect[uVect < lim1]
        ret <- c(ret,rep(0,length(uTEMP)))
        uTEMP <- uVect[uVect >= lim1 & uVect < lim2]
        ret <- c(ret,f2.1(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim2 & uVect < lim3]
        ret <- c(ret,f2.2(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim3 & uVect < lim4]
        ret <- c(ret,f2.3(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim4 & uVect < lim5]
        ret <- c(ret,f2.4(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim5 & uVect < lim6]
        ret <- c(ret,f2.5(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim6]
        ret <- c(ret,rep(0,length(uTEMP)))
      }

    if(lim4<lim3 && lim3<=lim5)
      {
        uTEMP <- uVect[uVect < lim1]
        ret <- c(ret,rep(0,length(uTEMP)))
        uTEMP <- uVect[uVect >= lim1 & uVect < lim2]
        ret <- c(ret,f2.1(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim2 & uVect < lim4]
        ret <- c(ret,f2.2(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim4 & uVect < lim3]
        ret <- c(ret,f2.6(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim3 & uVect < lim5]
        ret <- c(ret,f2.4(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim5 & uVect < lim6]
        ret <- c(ret,f2.5(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim6]
        ret <- c(ret,rep(0,length(uTEMP)))
      }

    if(lim2<=lim4 && lim5<lim3)
      {
        uTEMP <- uVect[uVect < lim1]
        ret <- c(ret,rep(0,length(uTEMP)))
        uTEMP <- uVect[uVect >= lim1 & uVect < lim2]
        ret <- c(ret,f2.1(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim2 & uVect < lim4]
        ret <- c(ret,f2.2(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim4 & uVect < lim5]
        ret <- c(ret,f2.6(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim5 & uVect < lim3]
        ret <- c(ret,f2.7(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim3 & uVect < lim6]
        ret <- c(ret,f2.5(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim6]
        ret <- c(ret,rep(0,length(uTEMP)))
      }

    if(lim4 < lim2)
      {
        uTEMP <- uVect[uVect < lim1]
        ret <- c(ret,rep(0,length(uTEMP)))
        uTEMP <- uVect[uVect >= lim1 & uVect < lim4]
        ret <- c(ret,f2.1(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim4 & uVect < lim2]
        ret <- c(ret,f2.8(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim2 & uVect < lim5]
        ret <- c(ret,f2.6(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim5 & uVect < lim3]
        ret <- c(ret,f2.7(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim3 & uVect < lim6]
        ret <- c(ret,f2.5(uTEMP,TPa,TPb,TPr))
        uTEMP <- uVect[uVect >= lim6]
        ret <- c(ret,rep(0,length(uTEMP)))
      }
    
    ret/(2*TPa^2*TPb^2)
  }

## CASE 3

## PARALLELOGRAM 1

f11 <- function(u,TPa,TPb,TPr,TPs)
  {
    ((TPr-TPa)^2+(TPs-TPb)^2)+u-2*(TPr-TPa)*sqrt(u-(TPs-TPb)^2)-2*(TPs-TPb)*sqrt(u-(TPr-TPa)^2)+(TPr-TPa)*(TPs-TPb)*(asin(1-2*(TPr-TPa)^2/u) + asin(1-2*(TPs-TPb)^2/u))
  }

f12 <- function(u,TPa,TPb,TPr,TPs)
  {
    TPa^2+2*(TPs-TPb)*(sqrt(u-TPr^2)-sqrt(u-(TPr-TPa)^2))-(TPr-TPa)*(TPs-TPb)*(asin(1-2*TPr^2/u)-asin(1-2*(TPr-TPa)^2/u))
  }

f13 <- function(u,TPa,TPb,TPr,TPs)
  {
    -(TPr*(TPr-2*TPa)+TPs*(TPs-2*TPb))-u+2*(TPr-TPa)*sqrt(u-TPs^2)+2*(TPs-TPb)*sqrt(u-TPr^2)-(TPr-TPa)*(TPs-TPb)*(asin(1-2*TPr^2/u)+asin(1-2*TPs^2/u))
  }

f14 <- function(u,TPa,TPb,TPr,TPs)
  {
    TPb^2-2*(TPr-TPa)*(sqrt(u-(TPs-TPb)^2)-sqrt(u-TPs^2))+(TPr-TPa)*(TPs-TPb)*(asin(1-2*(TPs-TPb)^2/u)-asin(1-2*TPs^2/u))
  }

f1 <- function(uVect,TPa,TPb,TPr,TPs)
  {
    lim <- c((TPr-TPa)^2+(TPs-TPb)^2,TPr^2+(TPs-TPb)^2,(TPr-TPa)^2+TPs^2,TPr^2+TPs^2)
    ret <- NULL
    if(lim[2]<=lim[3])
      {
        uTEMP <- uVect[uVect < lim[1]]
        ret <- c(ret,rep(0,length(uTEMP)))
        uTEMP <- uVect[uVect >= lim[1] & uVect < lim[2]]
        ret <- c(ret,f11(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[2] & uVect < lim[3]]
        ret <- c(ret,f12(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[3] & uVect < lim[4]]
        ret <- c(ret,f13(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >=lim[4]]
        ret <- c(ret,rep(0,length(uTEMP)))
      }

    if(lim[3] < lim[2])
      {
        uTEMP <- uVect[uVect < lim[1]]
        ret <- c(ret,rep(0,length(uTEMP)))
        uTEMP <- uVect[uVect >= lim[1] & uVect < lim[3]]
        ret <- c(ret,f11(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[3] & uVect < lim[2]]
        ret <- c(ret,f14(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[2] & uVect < lim[4]]
        ret <- c(ret,f13(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >=lim[4]]
        ret <- c(ret,rep(0,length(uTEMP)))
      }
    
    ret/(4*TPa^2*TPb^2)
  }
    
## PARALLELOGRAM 2

f21 <- function(u,TPa,TPb,TPr,TPs)
  {
    (-TPr^2-2*TPr*TPa-(TPs-TPb)^2)-u+2*(TPr+TPa)*sqrt(u-(TPs-TPb)^2)+2*(TPs-TPb)*sqrt(u-TPr^2)-(TPr+TPa)*(TPs-TPb)*(asin(1-2*(TPs-TPb)^2/u)+asin(1-2*TPr^2/u))
  }

f22 <- function(u,TPa,TPb,TPr,TPs)
  {
    TPa^2 - 2*(TPs-TPb)*(sqrt(u-(TPr+TPa)^2)-sqrt(u-TPr^2))+(TPr+TPa)*(TPs-TPb)*(asin(1-2*(TPr+TPa)^2/u)-asin(1-2*TPr^2/u))
  }

f23 <- function(u,TPa,TPb,TPr,TPs)
  {
    ((TPr+TPa)^2+TPs^2-2*TPs*TPb)+u-2*(TPr+TPa)*sqrt(u-TPs^2)-2*(TPs-TPb)*sqrt(u-(TPr+TPa)^2)+(TPr+TPa)*(TPs-TPb)*(asin(1-2*(TPr+TPa)^2/u)+asin(1-2*TPs^2/u))
  }

f24 <- function(u,TPa,TPb,TPr,TPs)
  {
    -TPb^2+2*(TPr+TPa)*(sqrt(u-(TPs-TPb)^2)-sqrt(u-TPs^2))-(TPr+TPa)*(TPs-TPb)*(asin(1-2*(TPs-TPb)^2/u)-asin(1-2*TPs^2/u))
  }

f2 <- function(uVect,TPa,TPb,TPr,TPs)
  {
    lim <- c((TPr)^2+(TPs-TPb)^2,(TPr+TPa)^2+(TPs-TPb)^2,(TPr)^2+(TPs)^2,(TPr+TPa)^2+(TPs)^2)
    ret <- NULL
    if(lim[2]<=lim[3])
      {
        uTEMP <- uVect[uVect < lim[1]]
        ret <- c(ret,rep(0,length(uTEMP)))
        uTEMP <- uVect[uVect >= lim[1] & uVect < lim[2]]
        ret <- c(ret,f21(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[2] & uVect < lim[3]]
        ret <- c(ret,f22(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[3] & uVect < lim[4]]
        ret <- c(ret,f23(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >=lim[4]]
        ret <- c(ret,rep(0,length(uTEMP)))
      }

    if(lim[3] < lim[2])
      {
        uTEMP <- uVect[uVect < lim[1]]
        ret <- c(ret,rep(0,length(uTEMP)))
        uTEMP <- uVect[uVect >= lim[1] & uVect < lim[3]]
        ret <- c(ret,f21(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[3] & uVect < lim[2]]
        ret <- c(ret,f24(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[2] & uVect < lim[4]]
        ret <- c(ret,f23(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >=lim[4]]
        ret <- c(ret,rep(0,length(uTEMP)))
      }
    
    ret/(4*TPa^2*TPb^2)
  }

## PARALLELOGRAM 3

f31 <- function(u,TPa,TPb,TPr,TPs)
  {
    (TPr*(TPr+2*TPa)+TPs*(TPs+2*TPb))+u-2*(TPr+TPa)*sqrt(u-TPs^2)-2*(TPs+TPb)*sqrt(u-TPr^2)+(TPr+TPa)*(TPs+TPb)*(asin(1-2*TPs^2/u)+asin(1-2*TPr^2/u))
  }

f32 <- function(u,TPa,TPb,TPr,TPs)
  {
    -TPa^2+2*(TPs+TPb)*(sqrt(u-(TPr+TPa)^2)-sqrt(u-TPr^2))-(TPr+TPa)*(TPs+TPb)*(asin(1-2*(TPr+TPa)^2/u)-asin(1-2*TPr^2/u))
  }

f33 <- function(u,TPa,TPb,TPr,TPs)
  {
    (-(TPr+TPa)^2-(TPs+TPb)^2)-u+2*(TPr+TPa)*sqrt(u-(TPs+TPb)^2)+2*(TPs+TPb)*sqrt(u-(TPr+TPa)^2)-(TPr+TPa)*(TPs+TPb)*(asin(1-2*(TPr+TPa)^2/u)+asin(1-2*(TPs+TPb)^2/u))
  }

f34 <- function(u,TPa,TPb,TPr,TPs)
  {
    -TPb^2+2*(TPr+TPa)*(sqrt(u-(TPs+TPb)^2)-sqrt(u-TPs^2))-(TPs+TPb)*(TPr+TPa)*(asin(1-2*(TPs+TPb)^2/u)-asin(1-2*TPs^2/u))
  }

f3 <- function(uVect,TPa,TPb,TPr,TPs)
  {
    lim <- c((TPr)^2+(TPs)^2,(TPr+TPa)^2+(TPs)^2,(TPr)^2+(TPs+TPb)^2,(TPr+TPa)^2+(TPs+TPb)^2)
    ret <- NULL
    if(lim[2]<=lim[3])
      {
        uTEMP <- uVect[uVect < lim[1]]
        ret <- c(ret,rep(0,length(uTEMP)))
        uTEMP <- uVect[uVect >= lim[1] & uVect < lim[2]]
        ret <- c(ret,f31(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[2] & uVect < lim[3]]
        ret <- c(ret,f32(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[3] & uVect < lim[4]]
        ret <- c(ret,f33(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >=lim[4]]
        ret <- c(ret,rep(0,length(uTEMP)))
      }

    if(lim[3] < lim[2])
      {
        uTEMP <- uVect[uVect < lim[1]]
        ret <- c(ret,rep(0,length(uTEMP)))
        uTEMP <- uVect[uVect >= lim[1] & uVect < lim[3]]
        ret <- c(ret,f31(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[3] & uVect < lim[2]]
        ret <- c(ret,f34(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[2] & uVect < lim[4]]
        ret <- c(ret,f33(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >=lim[4]]
        ret <- c(ret,rep(0,length(uTEMP)))
      }
    
    ret/(4*TPa^2*TPb^2)
  }

## PARALLELOGRAM 4

f41 <- function(u,TPa,TPb,TPr,TPs)
  {
    -(TPr-TPa)^2-TPs^2-2*TPb*TPs-u+2*(TPr-TPa)*sqrt(u-TPs^2)+2*(TPs+TPb)*sqrt(u-(TPr-TPa)^2)-(TPr-TPa)*(TPs+TPb)*(asin(1-2*TPs^2/u)+asin(1-2*(TPr-TPa)^2/u))
  }

f42 <- function(u,TPa,TPb,TPr,TPs)
  {
    -TPa^2-2*(TPs+TPb)*(sqrt(u-TPr^2)-sqrt(u-(TPr-TPa)^2))+(TPr-TPa)*(TPs+TPb)*(asin(1-2*TPr^2/u)-asin(1-2*(TPr-TPa)^2/u))
  }

f43 <- function(u,TPa,TPb,TPr,TPs)
  {
    (TPs+TPb)^2+TPr^2-2*TPa*TPr+u-2*(TPr-TPa)*sqrt(u-(TPs+TPb)^2)-2*(TPs+TPb)*sqrt(u-TPr^2)+(TPr-TPa)*(TPs+TPb)*(asin(1-2*TPr^2/u)+asin(1-2*(TPs+TPb)^2/u))
  }

f44 <- function(u,TPa,TPb,TPr,TPs)
  {
    TPb^2+2*(TPr-TPa)*(sqrt(u-TPs^2)-sqrt(u-(TPs+TPb)^2))-(TPs+TPb)*(TPr-TPa)*(asin(1-2*TPs^2/u)-asin(1-2*(TPs+TPb)^2/u))
  }

f4 <- function(uVect,TPa,TPb,TPr,TPs)
  {
    lim <- c((TPr-TPa)^2+(TPs)^2,(TPr)^2+(TPs)^2,(TPr-TPa)^2+(TPs+TPb)^2,(TPr)^2+(TPs+TPb)^2)
    ret <- NULL
    if(lim[2]<=lim[3])
      {
        uTEMP <- uVect[uVect < lim[1]]
        ret <- c(ret,rep(0,length(uTEMP)))
        uTEMP <- uVect[uVect >= lim[1] & uVect < lim[2]]
        ret <- c(ret,f41(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[2] & uVect < lim[3]]
        ret <- c(ret,f42(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[3] & uVect < lim[4]]
        ret <- c(ret,f43(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >=lim[4]]
        ret <- c(ret,rep(0,length(uTEMP)))
      }

    if(lim[3] < lim[2])
      {
        uTEMP <- uVect[uVect < lim[1]]
        ret <- c(ret,rep(0,length(uTEMP)))
        uTEMP <- uVect[uVect >= lim[1] & uVect < lim[3]]
        ret <- c(ret,f41(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[3] & uVect < lim[2]]
        ret <- c(ret,f44(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >= lim[2] & uVect < lim[4]]
        ret <- c(ret,f43(uTEMP,TPa,TPb,TPr,TPs))
        uTEMP <- uVect[uVect >=lim[4]]
        ret <- c(ret,rep(0,length(uTEMP)))
      }
    
    ret/(4*TPa^2*TPb^2)
  }

fcase3 <- function(u,TPa,TPb,TPr,TPs) f1(u,TPa,TPb,TPr,TPs)+f2(u,TPa,TPb,TPr,TPs)+f3(u,TPa,TPb,TPr,TPs)+f4(u,TPa,TPb,TPr,TPs)

if(!exists("aniso")) aniso <- 1

collapse.sum <- function(vect,positions)
      {
        res <- NULL
        pos <- 1
        for(count.i in 1:length(positions))
          {
            res <- c(res,sum(vect[pos:(pos+positions[count.i]-2)]))
            pos <- pos+positions[count.i]-1
          }
        res
      }

computeV <- function(info,class="ldt",params,eps=1e-8,cat.level=0)
  {
    getVTime <- proc.time()
    
    n <- info$nrows*info$ncols
    area <- info$rowwidth*info$colwidth

    ## density for distance between rectangles
    ptm <- proc.time()
    
    if(class != "ldt" && class != "power")
      {
        KandCov <- defineK(class)
        K <- KandCov$K  ## covariance function
        cov.f <- KandCov$cov.f  ## product of covariance function and density
      }

    if(cat.level>=1) cat("loading in functions",(proc.time()-ptm)[3],"seconds\n")
    
    ## this is what needs to be computed
    ptm <- proc.time()
    if(class=="ldt") {
      results <- apply(info$indices,1,f.anal.ldt,info=info)
      message <- "evaluating ldt analytic results"
    } else if(class=="power") {
      results <- apply(info$indices,1,f.anal.power,h=params[1],info=info)
      message <- "evaluating power analytic results"
    } else {
      results <- apply(info$indices,1,f.NI,params=params,eps=eps,K=K,cov.f=cov.f,info=info)
      message <- "evaluating numerical integrals"
    }
    
    ## collapse results back down
    results <- collapse.sum(results,info$lengths)
    if(cat.level>=1) cat(message,(proc.time()-ptm)[3],"seconds\n")

    ptm <- proc.time()
    results <- rep(results,info$rowReps)
    V <- matrix(0,n,n)
    
    ## now stick the results into V
    V[info$locations] <- results
    V[info$locations[,c(2,1)]] <- results
    diag(V) <- V[1,1]
    if(cat.level>=1) cat("Inserting values in V",(proc.time()-ptm)[3],"seconds\n")
    if(cat.level) cat("computing V takes",(proc.time()-getVTime)[3],"seconds\n")
    V
  }

if(FALSE) {
  ## compute info when we have at least one row or one column
  ## code fails when we only have one plot
  if(!exists("info") && nrows+ncols>2) {
    getV.prec <- getV.precompute(nrows,ncols,rowwidth,colwidth,rowsep,colsep,cat.level)
    info <- getV.prec$info
    rm(getV.prec)
  }
  
  ## this is a hack to get V when we have one row and one column
  ## this is hardly ever used, but just for a complete setup we have included it here
  ## it is used when checking the accuracy and timing of evaluating the power jobs
  
  if(!exists("info")) {
    getV.prec <- getV.precompute(nrows+1,ncols+1,rowwidth,colwidth,rowsep,colsep,cat.level)
    info <- getV.prec$info
    source("11case.R",local=T)
    rm(getV.prec)
  }
  
  if(info$aniso != aniso) {
    getV.prec.update <- getV.precompute.update(nrows,ncols,rowwidth,colwidth,rowsep,colsep,cat.level)
    info <- getV.prec.update$info
    rm(getV.prec.update)
  }
  
  if(exists("getV.precompute")) rm(getV.precompute)
  
  getVresult <- getV(info,class,eps,unitarea,cat.level,params)
  V <- getVresult$V
  area <- getVresult$area
  rm(getV,getVresult)
}
## April 24 2004

precompute <- function(nrows,ncols,rowwidth,colwidth,rowsep,colsep,cat.level=0)
  {
    if(nrows==1 && ncols==1) {
      ## this is for when there is one row and one column
      ## rarely used but there for completion
      
      info <- precompute(2,2,rowwidth,colwidth,rowsep,colsep,cat.level=0)
      info <- oneByoneHack(info)
    } else {
      preCompTime <- proc.time()
      
      ## April 15th 2004
      ## this code doesn't work when nrows=1 and ncols=1
      
      ## Jan 21st changes
      ## somethings I compute are integers, but R stores them as real numbers
      ## accounting for this in the initial steps at least cut back on memory allocation
      ##
      
      n <- nrows*ncols
      if(length(rowsep)==1) rowsep <- rep(rowsep,nrows-1)
      if(length(colsep)==1) colsep <- rep(colsep,ncols-1)
      nrow <- rep(seq(1:nrows), ncols) - 1
      ncol <- rep(seq(1:ncols), rep(nrows, ncols)) - 1
      t1 <- proc.time()
      num.ifail <- 0

      ## vector of x-coordinates:
      ax.vals <- 0
      for(gen.i in 1:(nrows-1))
        ax.vals <- c(ax.vals,ax.vals[gen.i]+rowwidth+rowsep[gen.i])
      if(cat.level>2) cat(ax.vals,"\n")

      bx.vals <- 0
      for(gen.i in 1:(ncols-1))
        bx.vals <- c(bx.vals,bx.vals[gen.i]+colwidth+colsep[gen.i])
      if(cat.level>2) cat(bx.vals,"\n")

      indices <- cbind(gl(n,n,n^2),gl(n,1,n^2))
      ## only use pairs where second is to the right of the first
      indices <- indices[indices[,2]>indices[,1],]
      indices <- rbind(as.integer(c(1,1)),indices)
      if(cat.level>1) {
        z <- sapply(ls(), function(x) object.size(get(x)))
        cat("Memory Allocation is approx:",sum(z)/1000000,"MB\n")
      }

      ## following what I used to do
      getV.i.bl <- as.integer(indices[,1]%%nrows)
      getV.i.bl[getV.i.bl==0] <- as.integer(nrows)
      getV.j.bl <- as.integer((indices[,1]-getV.i.bl)/nrows+1)
      getV.i.nt <- as.integer(indices[,2]%%nrows)
      getV.i.nt[getV.i.nt==0] <- as.integer(nrows)
      getV.j.nt <- as.integer((indices[,2]-getV.i.nt)/nrows+1)
      getV.i <- as.integer(abs(getV.i.nt-getV.i.bl)+1)
      getV.j <- as.integer(abs(getV.j.nt-getV.j.bl)+1)

      ax <- abs(ax.vals[getV.i.bl]-ax.vals[getV.i.nt])
      bx <- abs(bx.vals[getV.j.bl]-bx.vals[getV.j.nt])
      ax <- round(ax,12)
      bx <- round(bx,12)
      rm(getV.i.bl, getV.j.bl)
      rm(getV.i.nt, getV.j.nt)
      if(cat.level>1) {
        z <- sapply(ls(), function(x) object.size(get(x)))
        cat("Memory Allocation is approx:",sum(z)/1000000,"MB\n")
      }

      ## perhaps I can improve this step
      ## firstly getV.i, getV.j and indices are all integers
      ## but ax and bx are not, is it possible to store
      ## these in two groups to save memory
      ## if so, how do I use the apply command later
      ## later I will tack on another 5 columns,
      ## 3 of type integer and 2 of real
      ## now use lists
      indices <- list(ax=ax,bx=bx,getV.i=getV.i,getV.j=getV.j,indices=indices)
      rm(ax,bx,getV.i,getV.j)
      if(cat.level>1) {
        z <- sapply(ls(), function(x) object.size(get(x)))
        cat("Memory Allocation is approx:",sum(z)/1000000,"MB\n")
      }

      ## which of the rows need evaluation
      ## does not take row column symmetry into account
      evalFactor <- as.integer(as.factor(indices[[1]]):as.factor(indices[[2]]))

      ## now sort according to the numerical value of evalFactor
      indices$evalFactor <- evalFactor
      ord <- order(evalFactor)

      for(count.i in 1:length(indices))
        {
          if(is.matrix(indices[[count.i]])) {
            indices[[count.i]] <- indices[[count.i]][ord,]
          } else {
            indices[[count.i]] <- indices[[count.i]][ord]
          }
        }
      
      ## now make an old copy, and in the current one only keep the distinct rows
      ## all that is needed from the current version is the list of locations
      locations <- indices$indices
      if(cat.level>1) {
        z <- sapply(ls(), function(x) object.size(get(x)))
        cat("Memory Allocation is approx:",sum(z)/1000000,"MB\n")
      }

      keepRows <- (c(1,diff(indices$evalFactor))==1)
      keepRows <- (1:length(indices[[1]]))[keepRows]
      for(count.i in 1:length(indices))
        {
          if(is.matrix(indices[[count.i]])) {
            indices[[count.i]] <- indices[[count.i]][keepRows,]
          } else {
            indices[[count.i]] <- indices[[count.i]][keepRows]
          }
        }
      if(cat.level>1) {
        z <- sapply(ls(), function(x) object.size(get(x)))
        cat("Memory Allocation is approx:",sum(z)/1000000,"MB\n")
      }

      indices.preLimits <- indices

      ## now convert indices into a matrix again:
      indices <- cbind(indices[[1]],indices[[2]],indices[[3]],indices[[4]],indices[[5]],indices[[6]])

      ## limits of integration, may or may not be used
      ## how to incorporate aniso here...
      c.pts <- cbind((indices[,1]-rowwidth)^2+(indices[,2]-colwidth)^2, indices[,1]^2+(indices[,2]-colwidth)^2, indices[,1]^2+indices[,2]^2, (indices[,1]-rowwidth)^2+indices[,2]^2, indices[,1]^2+(indices[,2]-colwidth)^2, (indices[,1]+rowwidth)^2+(indices[,2]-colwidth)^2, (indices[,1]+rowwidth)^2+indices[,2]^2, indices[,1]^2+indices[,2]^2,  indices[,1]^2+indices[,2]^2, (indices[,1]+rowwidth)^2+indices[,2]^2, (indices[,1]+rowwidth)^2+(indices[,2]+colwidth)^2, indices[,1]^2+(indices[,2]+colwidth)^2, (indices[,1]-rowwidth)^2+indices[,2]^2, indices[,1]^2+indices[,2]^2, indices[,1]^2+(indices[,2]+colwidth)^2, (indices[,1]-rowwidth)^2+(indices[,2]+colwidth)^2)
      c.pts <- sqrt(c.pts)

      ## expand out the indices so each one corresponds to one integral,
      ## then stick them back together again
      extract.Limits <- function(vals) as.numeric(levels(as.factor(vals)))
      c.pts <- apply(c.pts,1,extract.Limits)
      lengths <- unlist(lapply(c.pts,length))  ## number of limits for each integral
      indices <- cbind(indices,lengths)

      ## expand the current indices to make numerical integration better
      indicesNEW <- NULL
      ptm <- proc.time()
      for(count.i in 1:length(c.pts))
        {
          c.pts[[count.i]] <- matrix(sort(c(c.pts[[count.i]],c.pts[[count.i]]))[2:(lengths[count.i]*2-1)],lengths[count.i]-1,2,byrow=T)
          ## indicator to say which rows are sufficient for analytic results
          ## 1 means evaluate analytic result
          ## 0 means do not, instead give value of 0
          c.pts[[count.i]] <- cbind(c(1,rep(0,dim(c.pts[[count.i]])[1]-1)),c.pts[[count.i]])
          indicesNEW <- rbind(indicesNEW,cbind(matrix(rep(indices[count.i,],lengths[count.i]-1),lengths[count.i]-1,dim(indices)[2],byrow=T),c.pts[[count.i]]))
        }
      if(cat.level>=1) cat("expanding the limits of integration",(proc.time()-ptm)[3],"seconds\n")

      ## function used to collapse the partial integrals into the corrent number of results
      collapse.sum <- function(vect,positions)
        {
          res <- NULL
          pos <- 1
          for(count.i in 1:length(positions))
            {
              res <- c(res,sum(vect[pos:(pos+positions[count.i]-2)]))
              pos <- pos+positions[count.i]-1
            }
          res
        }

      rowReps <- as.integer(summary.factor(as.factor(evalFactor),maxsum=length(evalFactor)))
      ## this indicates how many times a particular entry will have to be recorded in the matrix
      ## it has nothing to do with the rows in the lattice, more to do with the number of repetitions in the rows of the matrix that contains the results, unfortunate variable name

      ## what is needed from here:
      ## geometric info associated with the results below
      ## indicesNEW - information on what needs to be computed
      ## collapse.sum - converts the computed into the correct results
      ##              - by adding the parts from the different integrals
      ##              - or in the case analytic results, the result and some zeros
      ## rowReps      - needed to expand out the results
      ## locations
      ## lengths - know how much to expand the results by

      ## relabel these:
      ## info
      ## collapse.sum
      ## locations
      ## lengths as part of info, to be used in collapse.sum

      info <- list(rowwidth=rowwidth,colwidth=colwidth,rowsep=rowsep,colsep=colsep,nrows=nrows,ncols=ncols,lengths=lengths,rowReps=rowReps,locations=locations,indices=indicesNEW,aniso=1,indices.preLimits=indices.preLimits)
      
      if(cat.level) cat("precompute takes",(proc.time()-preCompTime)[3],"seconds\n")
    }
    info
  }
precompute.update <- function(info,cat.level=0,aniso=1)
  {
    preCompTime <- proc.time()
    
    ## aniso value has changed.
    ## we need to change the limits of integration

    ## divide out by the old aniso
    info$rowwidth <- info$rowwidth/info$aniso
    info$rowsep <- info$rowsep/info$aniso
    info$indices[,1] <- info$indices[,1]/info$aniso

    ## multiply by the new one
    info$aniso <- aniso
    info$rowwidth <- info$rowwidth*info$aniso
    info$rowsep <- info$rowsep*info$aniso
    info$indices[,1] <- info$indices[,1]*info$aniso

    ## compute new limits of integration
    ## the preLimits are always computed with aniso = 1
    indices <- info$indices.preLimits
    indices <- cbind(indices[[1]],indices[[2]],indices[[3]],indices[[4]],indices[[5]],indices[[6]])
    indices[,1] <- indices[,1]*info$aniso  ## rowsep information via ax

    c.pts <- cbind((indices[,1]-info$rowwidth)^2+(indices[,2]-info$colwidth)^2, indices[,1]^2+(indices[,2]-info$colwidth)^2, indices[,1]^2+indices[,2]^2, (indices[,1]-info$rowwidth)^2+indices[,2]^2, indices[,1]^2+(indices[,2]-info$colwidth)^2, (indices[,1]+info$rowwidth)^2+(indices[,2]-info$colwidth)^2, (indices[,1]+info$rowwidth)^2+indices[,2]^2, indices[,1]^2+indices[,2]^2,  indices[,1]^2+indices[,2]^2, (indices[,1]+info$rowwidth)^2+indices[,2]^2, (indices[,1]+info$rowwidth)^2+(indices[,2]+info$colwidth)^2, indices[,1]^2+(indices[,2]+info$colwidth)^2, (indices[,1]-info$rowwidth)^2+indices[,2]^2, indices[,1]^2+indices[,2]^2, indices[,1]^2+(indices[,2]+info$colwidth)^2, (indices[,1]-info$rowwidth)^2+(indices[,2]+info$colwidth)^2)
    c.pts <- sqrt(c.pts)

    ## expand out the indices so each one corresponds to one integral,
    ## then stick them back together again
    extract.Limits <- function(vals) as.numeric(levels(as.factor(vals)))
    c.pts <- apply(c.pts,1,extract.Limits)
    lengths <- unlist(lapply(c.pts,length))  ## number of limits for each integral
    indices <- cbind(indices,lengths)

    ## expand the current indices to make numerical integration better
    indicesNEW <- NULL
    ptm <- proc.time()
    for(count.i in 1:length(c.pts))
      {
        c.pts[[count.i]] <- matrix(sort(c(c.pts[[count.i]],c.pts[[count.i]]))[2:(lengths[count.i]*2-1)],lengths[count.i]-1,2,byrow=T)
        ## indicator to say which rows are sufficient for analytic results
        ## 1 means evaluate analytic result
        ## 0 means do not, instead give value of 0
        c.pts[[count.i]] <- cbind(c(1,rep(0,dim(c.pts[[count.i]])[1]-1)),c.pts[[count.i]])
        indicesNEW <- rbind(indicesNEW,cbind(matrix(rep(indices[count.i,],lengths[count.i]-1),lengths[count.i]-1,dim(indices)[2],byrow=T),c.pts[[count.i]]))
      }
    if(cat.level>=1) cat("expanding the limits of integration",(proc.time()-ptm)[3],"seconds\n")

    info <- list(rowwidth=info$rowwidth,colwidth=info$colwidth,rowsep=info$rowsep,colsep=info$colsep,nrows=info$nrows,ncols=info$ncols,lengths=lengths,rowReps=info$rowReps,locations=info$locations,indices=indicesNEW,aniso=aniso,indices.preLimits=info$indices.preLimits)

    if(cat.level) cat("precompute update takes",(proc.time()-preCompTime)[3],"seconds\n")

    info
  }
f.anal.power <- function(coords, h, info)
  {
    ## created Oct 10th
    ## K <- function(d,params) -d^(2*params[1])/params[1]
    ## analytic results for the power model
    ## see write up and notes in binder for more details.

    ## OCT 13th 2003
    ## include an option for the mode of computing the 2F1
    ## values include NI, SUM and NR for
    ## numerical integration, summation and numerical recipes
    
    a <- info$rowwidth
    b <- info$colwidth
    ax <- coords[1]
    bx <- coords[2]
    i <- coords[3]
    j <- coords[4]
    ## 5
    ## 6
    ## 7 - which rows to evaluate - evalFactor
    ## 8 - lengths
    ## 9 - indicator for which rows to evaluate analytic results
    evaluate <- coords[9]
    ## 10 - lower limit for numerical integration
    ## 11 - upper limit for numerical integration

    ret <- 0
    if(evaluate) {
      h <- 2*h
      
      ## most common situation first
      if(i != 1 && j != 1) {
        ## different version, which takes care of possible problems like 0/0, or log(0)
        alphas <- c(ax, ax-a, ax, ax+a)
        betas <- c(bx,bx-b,bx,bx+b)
        asbs <- cbind(as.vector(matrix(alphas,4,4)),as.vector(matrix(1:4,4,4)),as.vector(t(matrix(betas,4,4))),as.vector(t(matrix(1:4,4,4))))
        ret <- 0
        
        ## function is defined below
        ret <- sum(apply(asbs,1,sepRowsepCol.power.res,h=h))
        ret <- ret/(4.*a*a*b*b)
      } else {
        
        if(i!=1 && j==1) {
          ## simplified same row result using mathematica
          ## problems when ax=a with log(ax-a) dealt with below
          alphas <- c(ax, ax-a, ax, ax+a)
          ret <- 0
          
          ## function is defined below
          ret <- sum(apply(cbind(alphas,1:4),1,sameRow.power.res,b,h))
          ret <- ret/(2.*a*a*b*b)
        } else {
          if(i==1 & j!=1) {
            ## SAME COLUMN
            alphas <- c(bx, bx-b, bx, bx+b)
            ret <- 0
            
            ## function is defined below
            ret <- sum(apply(cbind(alphas,1:4),1,sameRow.power.res,b=a,h=h))
            ret <- ret/(2.*a*a*b*b)
          } else {
            if(i==1 && j==1) {
              ## simplified diagaonal result using mathematica
              asqr <- a*a
              bsqr <- b*b
              dsqr <- asqr+bsqr
              ret <- -(2/(h+2)+2/(h+4))*(dsqr)^((h+4)/2)
              ret <- ret + (2/(h+2)-4/(h+3)+2/(h+4))*(a^(h+4)+b^(h+4))
              ret <- ret + 4/3*(dsqr)^(h/2)*(bsqr*bsqr*Hypergeometric2F1(-h/2,1,2.5,bsqr/(dsqr))+asqr*asqr*Hypergeometric2F1(-h/2,1,2.5,asqr/(dsqr)))
              ret <- ret + 4*asqr*bsqr/((h+2)*sqrt(dsqr))*(a^(1+h)*Hypergeometric2F1((3+h)/2,0.5,1.5,bsqr/(dsqr)) + b^(1+h)*Hypergeometric2F1((3+h)/2,0.5,1.5,asqr/(dsqr)))
              ret <- ret /(asqr*bsqr)
            }
          }
        }
      }
      ret <- -ret/(h/2)
    }
    ret
  }

sepRowsepCol.power.res <- function(abVal,h)
  {
    aVal <- abVal[1]
    ii <- abVal[2]
    bVal <- abVal[3]
    jj <- abVal[4]
    asqr <- aVal*aVal
    bsqr <- bVal*bVal
    dsqr <- asqr+bsqr
    sign <- (-1)^(ii+jj+1)
    ab3h <- c(aVal,bVal)^(3+h)
    dsqrh2 <- dsqr^(h/2)
    
    numSteps <- 2-((aVal==0) + (bVal==0))
    ret <- 0
    ret <- ret + sign*(2/(h+2)+2/(h+4))*(dsqr*dsqr*dsqrh2)
    if(is.nan(ret)) ret <- 0
    if(numSteps) {
      ret <- ret - sign*4/3*(dsqrh2)*(asqr*asqr*Hypergeometric2F1(-h/2,1,2.5,asqr/(dsqr))+bsqr*bsqr*Hypergeometric2F1(-h/2,1,2.5,bsqr/(dsqr)))
      ret <- ret - sign*4/((h+2)*sqrt(dsqr))*(ab3h[1]*bsqr*Hypergeometric2F1((3+h)/2,0.5,1.5,bsqr/(dsqr)) + ab3h[2]*asqr*Hypergeometric2F1((3+h)/2,0.5,1.5,asqr/(dsqr)))
    }
    ret
  }

sameRow.power.res <- function(abVal,b,h)
  {
    aVal <- abVal[1]
    ii <- abVal[2]
    asqr <- aVal*aVal
    bsqr <- b*b
    dsqr <- bsqr + asqr
    ret <- 0
    sign <- (-1)^(ii)
    ab1h <- c(aVal,b)^(1+h)
    dsqrh2 <- (dsqr)^(h/2)
    
    ret <- ret - sign*(2/(2+h)+2/(4+h))*(dsqr*dsqr*dsqrh2)
    ret <- ret + sign*(2/(2+h)-4/(3+h)+2/(4+h))*aVal^(4+h)
    
    ret <- ret + sign*4/3*(dsqrh2)*(bsqr*bsqr*Hypergeometric2F1(-h/2,1,2.5,bsqr/(dsqr))+asqr*asqr*Hypergeometric2F1(-h/2,1,2.5,asqr/(dsqr)))
    
    if(aVal) ret <- ret + sign*4*asqr*bsqr/((h+2)*sqrt(dsqr))*(ab1h[1]*Hypergeometric2F1((3+h)/2,0.5,1.5,bsqr/(dsqr)) + ab1h[2]*Hypergeometric2F1((3+h)/2,0.5,1.5,asqr/(dsqr)))            
    ret
  }
