
## PrintLevels
## 0 : no message
## 1 : important error messages
## 2 : warnings
## 3 : minium debugging information
## 5 : extended debugging information

mleRF <-
function(coord,data,model,param,
         lower.kappa=NULL,upper.kappa=NULL,sill=NA,
         ...) {
  mleRF.default(coord=coord,data=data,model=model,param=param,
                lower.kappa=lower.kappa,upper.kappa=upper.kappa,sill=sill,...)
}
  
mleRF.default <-
function(coord,data,model,param,
         lower.kappa=NULL,upper.kappa=NULL,sill=NA,
         use.naturalscaling=TRUE,
         PrintLevel=0, trace.optim=0,
         bins=20,distance.factor=0.5,
         upperbound.scale.factor=10,
         lowerbound.scale.factor=20,
         lowerbound.scale.LS.factor=5,
         upperbound.var.factor=10,
         lowerbound.var.factor=100,
         lowerbound.sill=1E-10,
         scale.max.relative.factor=1000,
         minbounddistance=0.001, minboundreldist=0.02,
         approximate.functioncalls=50,
         pch="*"
) {
  index <- is.na(param)
  if (!any(index)) return(param)
  
  ENVIR <- environment()
  ip<-.C("GetParameterIndices", MEAN= integer(1), VARIANCE= integer(1),
          NUGGET= integer(1), SCALE= integer(1), 
          KAPPA= integer(1), LASTKAPPA= integer(1),integer(1),
          SILL= integer(1), DUP = FALSE)
  MEAN <- ip$MEAN+1
  VARIANCE <- ip$VARIANCE+1
  NUGGET <- ip$NUGGET+1
  SCALE <- ip$SCALE+1
  KAPPA <- ip$KAPPA+1
  LASTKAPPA <- ip$LASTKAPPA+1
  SILL <- ip$SILL+1

  stopifnot(is.character(model),
            all(is.na(param)) || is.numeric(param[!is.na(param)]),
            all(is.finite(as.matrix(coord))),
            all(is.finite(as.matrix(data))))
  covnr <- .C("GetModelNr", as.character(model), nr=integer(1))$nr
  if ((covnr==.C("GetModelNr", as.character("nugget"), nr=integer(1))$nr) &&
      (any(index[-c(MEAN,NUGGET)]) || (index[VARIANCE]!=0.0)))
    stop("if model=nugget effect then the variance must be zero, and only mean and/or nugget can be NA")   
  if (covnr<0) stop("model not ok")
  storage.mode(covnr) <- "integer"
  Xcoord <- coord[,1]
  storage.mode(Xcoord) <-"double"
  Ycoord <- coord[,2]
  storage.mode(Ycoord) <-"double"
  if ( (dim <- ncol(coord))!=2)
    stop("Dimension differs from 2. Sorry, not programmed yet.")
  storage.mode(dim) <- "integer"
  if (length(param)>LASTKAPPA) stop("parameter vector too long!")
  lower <- -Inf * seq(1, 1, l=length(param))
  lower[NUGGET] <- 0 
  if (!is.null(lower.kappa)) {
    lower[KAPPA:(KAPPA+length(lower.kappa)-1)]<-lower.kappa
  }
  upper <- Inf * seq(1,1,l=length(param)) 
  if (!is.null(upper.kappa)) {
    upper[KAPPA:(KAPPA+length(lower.kappa)-1)]<-upper.kappa
  }
  if (length(param)>=KAPPA) {
    l.kappa <- lower[KAPPA:length(param)]
    u.kappa <- upper[KAPPA:length(param)]
    nice.kappa <- (l.kappa+u.kappa)/2
    if (!all(is.finite(nice.kappa)))
      warning("limits for kappa are not all finite")
    nice.kappa[is.na(nice.kappa)] <- 1.0
    nice.kappa[nice.kappa==Inf] <- l.kappa + 1.0
    nice.kappa[nice.kappa==-Inf] <- u.kappa - 1.0
  } else nice.kappa <- NULL
  lc<- nrow(coord)
  storage.mode(lc) <- "integer"
  lpar <- length(param)
  storage.mode(lpar) <- "integer"
  distances <- as.double(dist(coord))
  storage.mode(distances) <- "double"
  mindistances <- min(distances[distances!=0])
  maxdistances <- max(distances)
  vardata <- var(data)
  variables <- c(0,vardata,0,maxdistances,nice.kappa)
  PARAM <- param ## just to make clear in MLEtarget and LStarget what is global
                 ## and not overwrite param

  step<- as.double(distance.factor * maxdistances / (bins-1)) 
  storage.mode(bins) <- "integer"
  if (is.na.mean <- is.na(PARAM[MEAN])) { PARAM[MEAN] <- mean(data) }
  MLEtargetV <- data-PARAM[MEAN] ## needed in the  next line
  ##                                and in MLE if !is.na.mean

  dummy <- .C("binnedvariogram",
              Xcoord,Ycoord,
              as.double(MLEtargetV),  ## change it if CoVariates!!!
              ##                         see all the others PARAM[MEAN], too
              lc,
              step,
              binned.variogram=double(bins),
              binned.n=integer(bins),
              bins, DUP=FALSE)
  binned.variogram <- dummy$binned.variogram
  binned.n <- dummy$binned.n
  dummy <- NULL
  bin.centers <- as.double(c(0,step/2 + (1:(bins-1)) * step))
  max.bin.vario <- max(binned.variogram)

  varnugNA <- FALSE 
  zeronugget <- FALSE
  sillbounded <- !is.na(sill)
  if (sum(index[-c(MEAN,VARIANCE)])==0) {
    if (index[VARIANCE]) PARAM[VARIANCE] <- 1
    cov.matrix <- chol(matrix(.C("CovarianceMatrixNatSc",
                                 distances,
                                 lc,
                                 covnr,
                                 as.double(PARAM),
                                 lpar,
                                 dim,
                                 cov.matrix=double(lc * lc),
                                 as.integer(0),
                                 DUP=FALSE)$cov.matrix     
                              ,ncol=lc))
    stopifnot(all(diag(cov.matrix)>=0))
    cov.matrix<- chol2inv(cov.matrix)
    cc <- crossprod(rep(1,lc), cov.matrix)
    PARAM[MEAN] <- (cc %*% data) / sum(cc)
    if (!index[VARIANCE]) return(PARAM)
    else {
      MLEtargetV <- data-PARAM[MEAN] 
      PARAM[VARIANCE] <-
        sqrt((crossprod(MLEtargetV, cov.matrix) %*% MLEtargetV)/lc)
      return(PARAM)
    }
  }
  if (sillbounded) {
    ## only VARIANCE need to be optimised
    ## NUGGET = SILL - VARIANCE
    if (xor(is.na(PARAM[VARIANCE]),is.na(PARAM[NUGGET])))
      stop("Sill fixed. Then variance and nugget should be given or unknown simultaneously.")
    if (!is.na(PARAM[VARIANCE]) && (PARAM[VARIANCE]+PARAM[NUGGET]!=sill))
      stop("sill !=  variance + nugget") 
    index[NUGGET] <- FALSE;
    PARAM[NUGGET] <- 0; ## just to avoid trouble with NAs later on (*)
    variables[VARIANCE]<-sill/2
    upper[NUGGET] <- upper[VARIANCE]<-sill
    lower[VARIANCE] <- lower[NUGGET] <- 0
  } else {
    if (is.na(PARAM[VARIANCE])) {
      if (is.na(PARAM[NUGGET])) {
         ## both NUGGET and VARIANCE have to be optimised
         ## instead optimise the SILL (which can be done by
         ## a formula) and estimate the part of NUGGET
         ## (which will be stored in NUGGET) !
         ##  consequences:
         ## * estimtation space increased only by 1, not 2
         ## * real nugget: NUGGET * SILL
         ## *  variance  : (1-NUGGET) * SILL
         varnugNA <- TRUE
         lower[NUGGET] <- 0
         upper[NUGGET] <- 1
	 variables[NUGGET] <- 0.5
         index[VARIANCE] <- FALSE
         PARAM[VARIANCE] <- 0 ## does not have any meaning, just to get rid
         ##                       of NA since later (*) checked whether VARIANCE
         ##                        in [0,oo] 
      } else {       
	if (PARAM[NUGGET]==0) {
          ## here SILL==VARIANCE and therefore, variance
          ## can be estimated by a formula without increasing the dimension
          
          ## more than only the variance has to be estimated
          zeronugget <- TRUE
          index[VARIANCE] <- FALSE
          PARAM[VARIANCE] <- 0 ## dito
        } else {
          upper[VARIANCE] <- upperbound.var.factor * max(binned.variogram)
          lower[VARIANCE] <- (vardata-PARAM[NUGGET])/lowerbound.var.factor
          if (lower[VARIANCE]<lowerbound.sill) {
            if (PrintLevel>0)
              cat("low.var=",lower[VARIANCE]," low.sill",lowerbound.sill,"\n")
            warning("param[NUGGET] might not be given correctly.")
            lower[VARIANCE] <- lowerbound.sill
          }
        }
      }
    }
    else {
      if (PARAM[VARIANCE]==0.0) {
        
      }
      if (is.na(PARAM[NUGGET])) { ## and !is.na(param[VARIANCE])
        lower[NUGGET]<-(var(data)-PARAM[VARIANCE])/lowerbound.var.factor
        if (lower[NUGGET]<lowerbound.sill) {
          if (PrintLevel>0)
            cat("\nlow.nug=",lower[NUGGET]," low.sill",lowerbound.sill,"\n")
          warning("param[VARIANCE] might not be given correctly!")
          lower[NUGGET]<-lowerbound.sill
        }
        upper[NUGGET] <- upperbound.var.factor * max.bin.vario
      }
    }
  }  

  if (index[SCALE]) {
    lower[SCALE]<- mindistances / lowerbound.scale.LS.factor
    upper[SCALE]<- maxdistances * upperbound.scale.factor
    if (use.naturalscaling) {
      ##  11: exact or numeric
      ##  13: MLE or numeric
      scalingmethod <- as.integer(11)  ## check with RFsimu.h, NATSCALEMLE!
      ##          or   3 ?
      if (.C("GetNaturalScaling",covnr,as.double(variables),
             as.integer(length(variables)),
             scalingmethod,double(1),error=integer(1),DUP=FALSE)$error)
        {
          scalingmethod <- as.integer(0)
          if (PrintLevel>0) cat("No natural scaling.")
        }
    } else {
      scalingmethod <- as.integer(0)## check with RFsimu.h, NATSCALEMLE!
    }
  } else {
    scalingmethod <- as.integer(0)## check with RFsimu.h, NATSCALEMLE!
  }

  stopifnot(length(index)==length(variables))
  if (any((PARAM[!index]>upper[!index]) | (PARAM[!index]<lower[!index]))){
    ## (*) : this part makes trouble if index has been reset (TRUE->FALSE)
    ##       and PARAM is not changed accordingly
    stop("fixed parameters out of range")
  }  
  index[MEAN] <- FALSE ## directly estimated
  variab <- variables[index]
  
  LB  <- lower[index]
  UB  <- upper[index]
  if (PrintLevel>2) {print("LB,UB");print(cbind(LB,UB))}

  weights <- cbind(sqrt(bins:1 * binned.n), bins:1, sqrt(bins:1), sqrt(binned.n))

  show <- function(text,col,param) {
    if (PrintLevel>3) {
      cat(text,format(param,dig=3),"\n")
      if (PrintLevel>4) {
        plot(bin.centers,binned.variogram)
        model.values <- .C("CovarianceNatSc",bin.centers,bins,
                          covnr,as.double(param),lpar,
                          dim,model.values=double(bins),
                          scalingmethod,DUP=FALSE)$model.values
        lines(bin.centers,param[VARIANCE]-model.values,col=col)
      }
    }
  }
  
 
  CoVariate <- matrix(1,ncol=1,nrow=lc)
  ## note "global" variables MLEMIN, MLEPARAM
  MLEtarget<-function(variab) {
    if (PrintLevel>2) print(variab,dig=10)
    if (any((variab<LB) | (variab>UB))) {
      if (PrintLevel>1) cat("WARNING! forbidden variab values !\n")  
      penalty <- variab 
      variab<-pmax(LB,pmin(UB,variab)) 
      penalty <- sum(variab-penalty)^2
      res <- MLEtarget(variab)
      if (res<=0) return(penalty + res * (1-penalty))
      if (PrintLevel>3) {
        cat("penalty",format(c(variab,penalty + res * (1+penalty))),"\n")
      }
      return(penalty + res * (1+penalty))
    }
    param <- PARAM
    param[index]<-variab  
    if (sillbounded) param[NUGGET]<-sill - param[VARIANCE]
    if (varnugNA) {
      param[VARIANCE] <- 1.0 - param[NUGGET]
    }
    if (zeronugget) param[VARIANCE] <- 1.0
    
    options(show.error.messages = FALSE)
    cov.matrix <- try(chol(matrix(.C("CovarianceMatrixNatSc",
                                     distances,
                                     lc,
                                     covnr,
                                     as.double(param),
                                     lpar,
                                     dim,
                                     cov.matrix=double(lc * lc),
                                     scalingmethod,
                                     DUP=FALSE)$cov.matrix     
                                 ,ncol=lc)))
    options(show.error.messages = TRUE)
    if (is.numeric(cov.matrix)) {
      if (any(diag(cov.matrix)<0)) {stop("chol det<0!")}
      logdet <- 2 * sum(log(diag(cov.matrix))) 
      cov.matrix<- chol2inv(cov.matrix)
    } else {
      return(1E300)
    }
    if (!is.finite(logdet)) logdet<-1E300 ## the best ?! 
    
    if (is.na.mean)  {
      dummy <- crossprod(CoVariate,cov.matrix)
      m <- solve(dummy %*% CoVariate, dummy %*% data)
      MLEtargetV  <- data - CoVariate %*% m
    }
    quadratic  <- crossprod(MLEtargetV, cov.matrix) %*% MLEtargetV
    
    if (varnugNA || zeronugget) {
      res <- logdet + lc * log(quadratic) 
    } else {
      res <- logdet + quadratic
    }   
    if (res<MLEMIN) {
      if (is.na.mean) param[MEAN] <- m
      if (sillbounded) param[NUGGET]<-sill-param[VARIANCE]
      if (varnugNA) {
        sill <- quadratic / lc
        param[VARIANCE] <- sill * (1.0 - param[NUGGET])
        param[NUGGET] <- sill * param[NUGGET]
      }
      if (zeronugget) {
        param[VARIANCE] <- quadratic / lc
      }
      assign("MLEMIN",res,envir=ENVIR)
      assign("MLEPARAM",param,envir=ENVIR)
    }
    if (PrintLevel>2) cat("result MLE",res,"\n")
    return(res)
  }
    
  ## Note
  ## creating "global" variables LSMIN, LSPARAM
  ## using WEIGHTS, BINNEDSQUARE
  LStarget <- function(variab) {
    if (PrintLevel>2) print(variab,dig=20)
    if (any((variab<LB) | (variab>UB))) {
      if (PrintLevel>1) cat("WARNING! forbidden variab values !\n")
      penalty <- variab 
      variab<-pmax(LB,pmin(UB,variab)) 
      penalty <- sum(variab-penalty)^2
      res <- LStarget(variab)
      if (res<=0) return(penalty + res * (1-penalty))
      return(penalty + res * (1+penalty))
    }
    param <- PARAM
    param[index]<-variab
    if (sillbounded) param[NUGGET]<-sill-param[VARIANCE]
    if (varnugNA) param[VARIANCE] <- 1.0 - param[NUGGET]
    if (zeronugget) param[VARIANCE] <- 1.0
    
    ## if CoVariates
    ##binned.variogram <- .C("binnedvariogram",Xcoord,Ycoord,
    ##     as.double(data-f(coord,parameter)),   ###
    ##     lc,step,binned.variogram=double(bins),integer(bins),bins,
    ##     DUP=FALSE)$binned.variogram
    model.values <-
      .C("VariogramNatSc", bin.centers, bins, covnr, as.double(param), lpar,
         dim,model.values=double(bins), scalingmethod, DUP=FALSE)$model.values
    if (varnugNA || zeronugget) {
      bgw <- sum(binned.variogram * model.values * WEIGHTS)
      g2w <- sum(model.values^2 * WEIGHTS)
      ergb <- BINNEDSQUARE - bgw^2/g2w
    } else {
      ergb<-sum((binned.variogram - model.values)^2 * WEIGHTS) 
    }
    if (ergb<LSMIN) {
      if (sillbounded) param[NUGGET]<-sill-param[VARIANCE]
      if (varnugNA) {
        sill <- bgw/g2w
        param[VARIANCE] <- sill * (1.0 - param[NUGGET])
        param[NUGGET] <- sill * param[NUGGET]
      }
      if (zeronugget) {
        param[VARIANCE] <- bgw/g2w
      }
      assign("LSMIN",ergb,envir=ENVIR)
      assign("LSPARAM",param,envir=ENVIR)
      show("LSX ",1,LSPARAM)
    }
    return(ergb)
  }

  ## Note
  ## creating "global" variable WEIGHTS, LSMIN, BINNEDSQUARE, LSpMIN,LSpMAX
  ## using LSPARAM, LSMIN
  LSoptimize <- function(index,variab) {
    minimum <- Inf
    LSpmin <- rep(Inf,length(variab))
    LSpmax <- rep(-Inf,length(variab))
    for (W in 1:ncol(weights)) {
      assign("WEIGHTS",weights[,W],envir=ENVIR)
      if (varnugNA || zeronugget) {
        assign("BINNEDSQUARE",sum(binned.variogram^2*WEIGHTS),envir=ENVIR)
      }
      assign("LSMIN", Inf,envir=ENVIR)
      
      options(show.error.messages = FALSE) ##
      neuvariab <- try(optim(variab, LStarget, method ="L-BFGS-B", lower = LB,
                            upper = UB, control=list(trace=trace.optim))$par)
      ## side effect: minimum so far is in LSMIN and LSPARAM
      ## even if the algorithm finally fails
      options(show.error.messages = TRUE) ##
      
      if (LSMIN>1E299) next ### ==1E300 or Inf,,see LStarget,
      ##                         and assign(LSMIN) above

      neuvariab <- LSPARAM[index]
      LSpmin <- pmin(LSpmin,neuvariab)
      LSpmax <- pmax(LSpmax,neuvariab)

      mlet <- MLEtarget(neuvariab) ## necessary since LStarget is optimised!

      if (mlet<minimum) {
        minimum.param <- LSPARAM
        minimum.variab <- LSPARAM[index]
        minimum <- mlet
        minimum.index <- W## only for debugging reasons.
        ##                   Which LS weights have been the best?
      }
    }    
    if (!is.finite(minimum)) {stop("Minimum not found")}
    if (PrintLevel>3)
      cat("LS. . . ",format(minimum.param,dig=3),"(",format(minimum.index),")")
    return(param=minimum.param,variab=minimum.variab,minimum=minimum,
           LSpmin=LSpmin,LSpmax=LSpmax)
  }

  ## find a good initial value for MLE using  weighted least squares
  ## and binned variogram
  ##
  ## background: if the number of observations (and the observation
  ## field) tends to infinity then any least square algorithm should
  ## yield the same result as MLE
  ## so the hope is that for a finite number of points the least squares
  ## find an acceptable initial values  
  assign("MLEMIN",Inf,envir=ENVIR) ## necessary here, as LSoptim calls MLEtarget!
  cat(pch)
  LSopt <- LSoptimize(index,variab)
  PARAM <- LSopt$param    ## could use also LSPARAM directly...
  variab <- LSopt$variab
  LSpmin <- LSopt$LSpmin
  LSpmax <- LSopt$LSpmax
   
  ## if rather close to nugget and nugget==0 (fixed), do exception handling.
  ## This part is active only if
  ## scale.max.relative.factor < lowerbound.scale.LS.factor
  ## By default it is not, i.e., the following if-condition
  ## will "always" be FALSE.
  if ((PARAM[SCALE] < mindistances/scale.max.relative.factor) &&
      (!index[NUGGET]))
    {
      if (PrintLevel>1) cat("Looks like if a nugget effect is included...\n")
      if (PARAM[NUGGET]!=0) stop("Chosen model does not seem to be appropriate!")
      param[NUGGET] <- NA;
      return(mleRF(coord=coord, data=data, model=model, param=param,
                   lower.kappa=lower.kappa, upper.kappa=upperkappa, sill=sill,
                   use.naturalscaling=use.naturalscaling,
                   PrintLevel=PrintLevel, trace.optim=trace.optim,
                   bins=bin, distance.factor=distance.factor,
                   upperbound.scale.factor=upperbound.scale.factor,
                   lowerbound.scale.factor=lowerbound.scale.factor,
                   lowerbound.scale.LS.factor=lowerbound.scale.LS.factor,
                   upperbound.var.factor=upperbound.var.factor,
                   lowerbound.var.factor=lowerbound.var.factor,
                   lowerbound.sill=lowerbound.sill,
                   scale.max.relative.factor=scale.max.relative.factor,
                   minbounddistance=minbounddistance,
                   minboundreldist=minboundreldist,
                   approximate.functioncalls=approximate.functioncalls,
                   pch=pch
                   )
             )
    }

  show("Final",3,PARAM)
  
  ## now MLE
  assign("MLEMIN",LSopt$minimum,envir=ENVIR)
  assign("MLEPARAM",PARAM,envir=ENVIR)

  ## lowerbound.scale.LS.factor <  lowerbound.scale.factor, usually
  ## LS optimisation should not run to a boundary (what often happens
  ## for the scale) since a boundary value is usually a bad initial
  ## value for MLE (heuristic statement). Therefore a small
  ## lowerbound.scale.LS.factor is used for LS optimisation.
  ## For MLE estimation we should include the true value of the scale;
  ## so the bounds must be larger. Here lower[SCALE] is corrected
  ## to be suitable for MLE estimation
  lower[SCALE] <-
    lower[SCALE] * lowerbound.scale.LS.factor / lowerbound.scale.factor
  LB  <- lower[index]

  ##save(file="debug.save",MLEMIN,MLEPARAM,variab,LB,UB,trace.optim,covnr,index,
  ##     PrintLevel,param,sillbounded,varnugNA,zeronugget,MEAN,NUGGET,VARIANCE,
  ##     is.na.mean,SCALE,lc,distances,lpar,dim,scalingmethod,CoVariate,data)

  distances <- as.double(distances)
  options(show.error.messages = FALSE) ##

  cat(pch)
  variab <- try(optim(variab, MLEtarget, method="L-BFGS-B", lower = LB,
                      upper = UB, control=list(trace=trace.optim))$par)
  
  options(show.error.messages = TRUE) ##
  
  if (!is.numeric(variab) && (PrintLevel>0)) cat("MLEtarget I failed.\n")

  PARAM <- MLEPARAM
  variab <- PARAM[index]
 
  mindistance <- pmax(minbounddistance, minboundreldist * abs(variab))
  onborderline <- any(abs(variab-LB) <
                      pmax(mindistance,              ## absolute difference
                           minboundreldist * abs(LB) ## relative difference
                           )) ||
                  any(abs(variab-UB) <
                      pmax(mindistance, minboundreldist * abs(UB))) 
  ## if the MLE result is close to the border, it usually means that
  ## the algorithm has failed, especially because of a bad starting
  ## value (least squares do not always give a good starting point, helas)
  ## so the brutal method:
  ## calculate the MLE values on a grid and start the optimization with
  ## the best grid point. Again, there is the believe that the
  ## least square give at least a hint what a good grid is
  if ( onborderline ) {
    if (PrintLevel>3) cat("MLE a ", format(c(MLEPARAM,MLEMIN),dig=4),"\n")
    gridlength <- round(approximate.functioncalls^(1/sum(index)))
    ## grid is given by the extremes of the LS results
    ## so, therefore we should examine above at least 4 different sets
    ## of weights 
    step <- (LSpmax-LSpmin) / (gridlength-2) # grid starts bit outside

    LSpmax <- pmin(LSpmax + step/2,UB)       # the extremes of LS
    LSpmin <- pmax(LSpmin - step/2,LB)
    step <- (LSpmax-LSpmin) / (gridlength-1)
    i <- 1
    zk <-  paste("LSpmin[",i,"] + step[",i,"] * (0:",gridlength-1,")")
    if (length(step)>1)
      for (i in 2:length(step))
        zk <- paste(zk,",LSpmin[",i,"] + step[",i,"] * (0:",gridlength-1,")")
    zk <- paste("expand.grid(",zk,")")
    startingvalues <- eval(parse(text=zk))
    MLEMIN.old <- MLEMIN
    MLEPARAM.old <- MLEPARAM
    MLEMIN <- Inf
    cat(pch)
    apply(startingvalues,1,MLEtarget) ## side effect: Minimum is in
    ##                                   MLEMIN !
    PARAM <- MLEPARAM
    variab <- PARAM[index]
    options(show.error.messages = FALSE) ##

    cat(pch)
    variab <- try(optim(variab, MLEtarget,
                        method ="L-BFGS-B",
                        lower = LB, upper = UB,
                        control=list(trace=trace.optim))$par)
    options(show.error.messages = TRUE) ##
    if (!is.numeric(variab) && (PrintLevel>0)) cat("MLEtarget II failed.\n")
    ## do not check anymore whether there had been convergence or not.
    ## just take the best of the two strategies (initial value given by
    ## LS, initial value given by a grid), and be happy.
    if (MLEMIN<MLEMIN.old)  PARAM <- MLEPARAM
    else PARAM <- MLEPARAM.old
    variab <- PARAM[index]        
  }
  show("MLE",col=2,MLEPARAM)

  ## if the covariance functions use natural scaling, just
  ## correct the final output by GNS$natscale
  ## (GNS$natscale==1 if no rescaling was used)
  GNS <- .C("GetNaturalScaling",
            covnr,
            as.double(PARAM),
            as.integer(length(PARAM)),
            scalingmethod,
            natscale=double(1),
            error=integer(1),DUP=FALSE)
  if (GNS$error)
    stop(paste("Error",error,"occured whilst rescaling"))
  PARAM[SCALE] <- PARAM[SCALE] * GNS$natscale
  if (PrintLevel>2)
    cat("MLE(end)",format(c(PARAM,,MLEMIN),dig=4),"\n")
  if (pch!="") cat("\n")
  return(PARAM)
}


ShowModels <- function(covx=ifelse(is.null(empirical),4,max(empirical$c)),
                       x=NULL, y=NULL,
                       fixed.rs=FALSE,
                       method=NULL,
                       empirical=NULL,
                       model=NULL,
                       param=NULL,
                       all.param=NULL,
                       PracticalRange=FALSE,   
                       legend = TRUE,
                       register=0,
                       ...){

  stopifnot(all(covx>=0))
  if (length(covx)>1) {
    covx <- covx[covx>0]
    covx <- c(0, sort(c(max(covx)/10000, covx)))
  } else {
    covx <- c(0,seq(covx/10000, covx, l=100))
  }
  
  scr <- split.screen(figs=rbind(
                        c(0.01,0.60,0.01,0.49),
                        c(0.62,0.99,0.01,0.49),
                        c(0.01,0.60,0.51,0.99),
                        c(0.62,0.99,0.51,0.99)))
  if (isnullX <- is.null(x)) {
    if (!is.null(y)) stop("x is null, but not y")
  }

  if (!is.null(y)) {
    ll <-  list(...)
    opts <- lapply(names(ll), pmatch,c("col","zlim"),no=0)
    if (length(opts)==0) opts <- 0 else opts <- as.integer(opts)

    if (sum(index <- (opts==1))==1) col <- ll[index][[1]]
    else {
      if (exists("image.default")) ## R-1.3.0 onwards
        text <- paste("col <-",
                        paste(as.character(as.list(args(image.default))$col),
                              collapse="("),")")
     else  ## R-1.2.3
       text <- paste("col <-",
                    paste(as.character(as.list(args(image))$col),
                          collapse="("),")")         
      eval(parse(text=text))
  }
    if (!(missing.zlim <- (sum(index <- (opts==2))!=1))) zlim <- ll[index][[1]]
    cn <- length(col)
    if (cn==1) stop ("more than one colour must be given!")
    lwd <- 1
    y.i <-  0.03
    if (big.legend <- (cn>50)) {
      filler <- vector("character",length=(cn-3)/2)
      if (cn>80) y.i <-  2.4 / cn
    } else {
      filler <- vector("character",length=cn-2)
      if (cn==2) {
        lwd <- 5
        y.i <- 1
      }
      else if (cn<30) {
        y.i <- 0.9 / (cn-2)
        lwd <- 1 + 10/(cn-2)
      }
    }
    lu.x <- min(x)
    lu.y <- max(y)
  }
  runif(1)
  rs <- .Random.seed
  screen(scr[4])
  maxrow<- 14
  namen <- GetModelNames()
  n <- length(namen)
  par(mar=c(0,0,2,0))
  ncol <- 1 + as.integer( (n-1) / maxrow)
  top <- min(maxrow+1,n,na.rm=TRUE)
  plot(-1,-1,xlim=c(0,ncol),ylim=c(1,top+1),axes=FALSE)
  maintitle <- "Models"
  title(main=maintitle)
  for (i in 1:n) {
    text(as.integer((i-1) / maxrow), top - (i-1) %% maxrow,
         labels=paste(namen[i]),adj=c(0,0))
  }
  p <- .C("GetParameterIndices",integer(1),integer(1),integer(1),
          scale=integer(1),
          kappa=integer(1),lastkappa=integer(1),integer(1),integer(1))
  p <- lapply(p,function(k) k+1)
  nkappa <- p$lastkappa - p$kappa + 1
  
  globals <- rbind(c("Variogram","CovarianceFct"),
                   c("practical range","math. def")) ##"" to note numericals
  shortglobals <- rbind(c("Vario","CovFct"),
                   c("pract","math")) ##"" to note numericals
  globmin <- c(1,1)
  globmax <- c(2,2)
  globminsteps <- c(1,1)
  globmaxsteps <- c(1,1)
  if (is.null(empirical)) glob <- c(2,2-PracticalRange)
  else {
    glob <- c(1,2-PracticalRange)
    screen(scr[3])
    par(mar=c(2,2,0,0))
    ylim <- range(empirical$e, na.rm=TRUE)
    plot(empirical$c, empirical$e, pch="*",xlim=range(covx, na.rm=TRUE),
         ylim=ylim) 
  } 
  nglob <- length(glob)
  mins <- c(-Inf,0,0,0.001,rep(-Inf,nkappa))
  minsteps <- c(0.0001,0.0001,0.0001,0.0001,rep(0.0001,nkappa))
  maxsteps <- c(1,1,1,1,rep(1,nkappa))
  maxs <- c(Inf,Inf,Inf,Inf,rep(Inf,nkappa))
  pnames <- c("mean","variance","nugget","scale","a","b","c","d","e")
  
  currentparam <- list()
  if (unknown.par <- is.null(all.param)) all.param <- c(0,1,0,1)
  else stopifnot(length(all.param)==4)

  
  for (i in 1:n) 
    currentparam[[i]] <-  c(all.param,rep(1,.C("GetNrParameters",
                                             as.integer(i-1),
                                             k=integer(1))$k))
  i <- .C("GetModelNr","cone",nr=integer(1))$nr + 1
  currentparam[[i]] <- c(all.param,rep(0.5,.C("GetNrParameters",
                                          as.integer(i-1),
                                          k=integer(1))$k))
  currentparam[[.C("GetModelNr","cauchy",nr=integer(1))$nr+1]] <-
    c(all.param,2)
  currentparam[[.C("GetModelNr","cauchytbm",nr=integer(1))$nr+1]] <-
    c(all.param,1,1,2)
  currentparam[[.C("GetModelNr","power",nr=integer(1))$nr+1]] <-
    c(all.param,2)
  
  if (unknown.par) { all.param[p$scale] <- 0.3 }                     
  currentparam[[.C("GetModelNr","bessel",nr=integer(1))$nr+1]] <-
    c(all.param,1)
  currentparam[[.C("GetModelNr","wave",nr=integer(1))$nr+1]] <- c(all.param)
  currentparam[[.C("GetModelNr","gengneiting",nr=integer(1))$nr+1]] <-
    c(all.param,1,3)

  special <- rep(0,n)
  paramAdiscrete1 <- 1
  special[.C("GetModelNr","gengneiting",nr=integer(1))$nr+1] <- paramAdiscrete1
  
  covlist <- c("bessel","cauchy","cauchytbm","circular",
               "cone","cubic","exponential","gaussian",
               "gencauchy","gengneiting","gneiting","gneitingdiff",
               "holeeffect",
               "hyperbolic","nugget","pentamodel","power",
               "qexponential","spherical","stable","wave",
               "whittlematern")
  exprlist <- c(expression(2^a *Gamma(a+1)*x^{-a}* J[a](x)),
                expression((1+x^2)^{-a}),
                expression((1+(1-b/c)*x^a)*(1+x^a)^{-b/a-1}),
                expression((1-2/pi*(x *sqrt(1-x^2)+asin(x))) * (x<1)),
                
                expression("C(x) not available"),
                expression((1-7*x^2 + 8.75 *x^3 - 3.5* x^4+0.75 *x^6)*(x<1)),
                expression(e^{-x}),
                expression(exp(-x^2)),
                
                expression((1+x^a)^{-b/a}),
                expression("formula see help page for CovarianceFct"),
                expression((1 + 8 *s *x + 25* s^2* x^2 + 32*
                    s^3 *x^3)*(1-s* x)^8 * (x<1)),
                expression((1 + 8 *x/b + 25 *(x/b)^2 + 32* (x/b)^3)*(1-x/b)^8
                    * 2^{1-a} *Gamma(a)^{-1}* x^a* K[a](x)*(x<1) ),

                expression(e^{-a* x}* cos(x)),
                
                expression(const[a][b][c] * (c^2 +x^2)^{b/2} *
                    K[b](sqrt(a(c^2 + x^2)))),
                expression((x==0)),
                expression(( 1 - 22*x^2/3  +33 *x^4  - 77*x^5/2  + 33* x^7/2
                    - 11* x^9/2 + 5 *x^11 /6)*(x<1)),
                expression((1-x)^a  * (x<1)),
                
                expression((2*e^{-x}-a*e^{-2*x})/(2-a)),
                expression((1 - 1.5* x + 0.5* x^3)*(x<1)),
                expression(exp(-x^a)),
                expression(sin(x)/x ),
                    
                expression(2^{1-a}* Gamma(a)^{-1}* x^a * K[a](x))
                )
  expr <- NULL
  for (i in 1:n) {
    j <- pmatch(namen[i],covlist)
    if (!is.na(j)) expr[i] <- exprlist[j] else expr[j] <- "C(x) unknown"
  }
  DeleteRegister(register)
  
  oldPracticalRange <-  RFparameters()$PracticalRange
  RFparameters(PracticalRange=2-glob[2])

  first <- TRUE
  repeat {
    #screen(scr[4]);title(main=maintitle,col="red")
    screen(scr[4],new=FALSE)

    if (first && (!is.null(model))) {
      first <- FALSE
      covnr <- as.integer(.C("GetModelNr", as.character(model),
                             nr = integer(1))$nr) + 1
     
      if ((covnr<=0) || (covnr>n)) {
        cat("warning: unknown covariance model")
        next;
      } else {
        if (!is.null(param)) currentparam[[covnr]] <- param
      }
    } else {      
      if (length(loc <-locator(1))==0) break;
      loc <- lapply(loc,floor)
      covnr <- (1+top-loc$y) + maxrow * loc$x
      if ((covnr<=0) || (covnr>n) || (loc$y<0) || (loc$y>top))  next;   
    }
    cov <- namen[covnr];

    cp <- length(currentparam[[covnr]])
    currentp <- cp  + nglob
    ##screen(scr[4]);title(main=maintitle,col="black")
                   
    textcol2 <- "red"
    locy <- nglob
    repeat {
      ##erase.screen(scr[2])
      screen(scr[2])
      par(mar=c(0,0,2,0))
      plot(-1,-1,xlim=c(0,1),ylim=c(0,currentp+1.5),axes=FALSE)
      if (textcol2=="black")  title(main=cov,col=textcol2)
      else title(paste("*",main=cov,"*"),col=textcol2)
      lines(c(0.5,0.5),c(-0.1,currentp+0.5),col="lightgrey")
     
      if (locy<nglob) {
        locy <- locy+1
        glob[locy] <- min(globmax[locy],max(globmin[locy],
              glob[locy] + 2 * (loc$x-0.5) * globmaxsteps[locy] 
                         + sign(loc$x-0.5) * globminsteps[locy]))
        if (locy==2) {
          DeleteRegister(register)
          RFparameters(PracticalRange=2-glob[locy])
        }
      } else if (locy>nglob) {
        locy <- locy - nglob
        
        if (spec <- (locy>=p$kappa) && (special[covnr]!=0)) {
          spec <- FALSE
          switch(special[covnr],
                 { ##paramAdiscrete1
                   if (locy==p$kappa) {
                     spec <- TRUE
                     currentparam[[covnr]][locy] <-
                       min(maxs[locy], max(mins[locy],
                                           currentparam[[covnr]][locy] +
                                           sign(loc$x-0.5)))
                   }
                 }
                 )
        }
        if (!spec)
          currentparam[[covnr]][locy] <- min(maxs[locy],max(mins[locy],
               currentparam[[covnr]][locy] + 2 * (loc$x-0.5) * maxsteps[locy]
                                     + sign(loc$x-0.5) * minsteps[locy]))
      }
      text(0,currentp+1,lab=expr[covnr],adj=c(0,0))
      positions <- ( 1:(currentp+1) )[-c(nglob+1)] - 1
      labels <- rep("-----",length(positions))
      labels[(1:nglob)[globals[,2]!=""]] <- shortglobals[globals[,2]!="",1]
      text(0,positions,labels=labels,adj=c(0,0),col="blue")
      labels <- rep("+++++",length(positions))
      labels[(1:nglob)[globals[,2]!=""]] <- shortglobals[globals[,2]!="",2]
      text(1,positions,labels=labels,adj=c(1,0),col="red")
      txt <- globals[cbind(c(1,2),glob)]
      txt[globals[,2]==""] <- globals[globals[,2]=="",1]
      text(0.5, (1:nglob)-1, labels=txt, adj=c(0.5,0))
      text(0.5, nglob+(1:cp),
           labels=paste(pnames[1:cp]," (",
             format(currentparam[[covnr]],d=3),")",sep=""),
           adj=c(0.5,0))

      if (textcol2=="black") break;

      ##erase.screen(scr[3])
      screen(scr[3])
      par(mar=c(2,2,0,0))

      eval(parse(text=paste("covvalue <- ",globals[1,glob[1]],
                   "(covx,model=cov,param=currentparam[[covnr]],dim=2-is.null(y))")))
     if (is.na(covvalue[1])) {
        plot(0,0,col=0,axes=FALSE)
        text(0,0,label="plot not available",adj=c(0.5,0.5))
      }
      else {
        if (!is.null(empirical) && (glob[1]==1))
          ylim <- range(covvalue,empirical$e, na.rm=TRUE)
        else ylim <- range(covvalue, na.rm=TRUE)
        plot(covx[1], covvalue[1], type="p", xlim=range(covx, na.rm=TRUE),
             ylim=ylim)
        lines(covx[-1],covvalue[-1])
        if (!is.null(empirical) && (glob[1]==1)) {
          points(empirical$c, empirical$e, pch="*")
        }
      }
      
      if (!isnullX) {
        screen(scr[1],FALSE)
        if (!is.null(y)) lines(range(x),range(y),col="grey")
        if (fixed.rs) set.seed(rs)
        z <- GaussRF(x,y,model=cov,param=currentparam[[covnr]],grid=TRUE,
                     register=register,method=method)

        ##erase.screen(scr[1])
        screen(scr[1])
        par(mar=c(2,2,0,0))
        if (is.null(z)) {
          plot(0,0,col=0,axes=FALSE)
          text(0,0,label="image not available",adj=c(0.5,0.5))
        } else {          
          if (is.null(y)) { plot(x,z,...) }
          else {
            image(x,y,z,...)
            if (legend) {
              if (missing.zlim) zlim <- range(z)
              if (big.legend) ml <- c(format(mean(zlim),d=2),filler)
              else ml <- NULL
              legend(lu.x, lu.y, y.i=y.i, x.i=0.1, 
                     legend=c(format(zlim[2],d=2),filler,
                       ml, format(zlim[1],d=2)),
                     lty=1, lwd=lwd, col=col[length(col):1])
            }
          }
        } # is.null(z)
      }

      repeat {
        screen(scr[2],FALSE)
        if (length(loc <-locator(1))==0) {
          locy <- nglob; textcol2 <- "black"; break; }
        locy <- floor(loc$y)
        if (locy>currentp+1)  {locy <- nglob; textcol2 <- "black"; break;}
        if ((locy>=0) && (locy<=currentp)&& (locy!=nglob) && (loc$x>=0)) break;
      }
    }
  }
 
  RFparameters(PracticalRange=oldPracticalRange)
  close.screen(scr)
  if (exists("covnr"))
    return(list(model=cov,param=currentparam[[covnr]],
                PracticalRange=as.logical(2-glob[2])))
  else return(NULL)
}
  
"InitMaxStableRF" <-
function(x, y = NULL, z = NULL, grid, model, param, maxstable,
         method = NULL, register = 0, gridtriple = FALSE) 
{
  MaxStableList <- c("extremalGauss","BooleanFunction")
  MaxStableNr <- pmatch(maxstable,MaxStableList)
  if (is.na(MaxStableNr)) stop(paste("Unknown max-stable random field",
                                     "\nPossible values for `maxstable': \"",
                                     paste(MaxStableList,collapse="\", \""),
                                     "\"",sep=""))
  if (MaxStableNr==2) {
    if (is.null(method)) method <- "max.MPP"
    else{
      if (!is.character(method)) stop("Method must be a string.")
      if (.C("GetMethodNr",method, nr = integer(1))$nr !=
          .C("GetMethodNr","max.MPP", nr = integer(1))$nr) { 
        warning("Method does not match max-stable random field definition. Set to `max.MPP'.")  
        method <- "max.MPP"
      }
    }
  }
  return (InitSimulateRF(x=x, y=y, z=z, grid=grid, model=model, param=param,
                         method=method, register=register, gridtriple=gridtriple,
                         distribution="MaxStable")
          )
}

"MaxStableRF" <-
function (x, y = NULL, z = NULL, grid, model, param,  maxstable,
          method = NULL, n = 1, register = 0, gridtriple = FALSE) 
{
    if (InitSimulateRF(x=x, y=y, z=z, grid=grid, model=model, param=param,
                       method=method, register=register, gridtriple=gridtriple,
                       distribution="MaxStable") <= 0) {
        return(DoSimulateRF(n=n, reg=register))
    }
    else {
        return(NULL)
    }
}

## to.do replace pure nugget effect model as error.model by arbitrarz model

Kriging <- function(krige.method, x, y=NULL, z=NULL, grid, gridtriple=FALSE,
                    model, param, given, data, pch=".") {
  krige.methodlist <- c("S","O")
  if (is.na(krige.method.nr <- pmatch(krige.method,krige.methodlist)))
    stop(paste("kriging method not identifiable from the list.",
               "Possible values are",paste(krige.methodlist,collapse=",")))
  ip<-.C("GetParameterIndices", MEAN= integer(1), VARIANCE= integer(1),
          NUGGET= integer(1), SCALE= integer(1), 
          KAPPA= integer(1), LASTKAPPA= integer(1),integer(1),
          SILL= integer(1), DUP = FALSE)
  ip <- lapply(ip,function(i){i+1})

  x <- CheckXYZ(x, y, z, grid, gridtriple)
  x <- cbind(x$x, x$y, x$z)
  y <- z <- NULL
  ncol.data <- ncol(data)
  is.matrix.data <- is.matrix(data) && (ncol.data>1)
  given <- as.matrix(given)
  nd <- nrow(given)
  covmatrix <- diag(nd) * (CovarianceFct(0,model,param) / 2)
  covmatrix[lower.tri(covmatrix)] <- CovarianceFct(dist(given),model,param)
  covmatrix <- covmatrix + t(covmatrix)

  tgiven <-  t(as.matrix(given))
  given <- NULL

  if (ncol(x)!=nrow(tgiven))
    stop("dimension of the kriged points and the given points do not match")
  if (grid) {
    dimension <- 1 + floor((x[2,]-x[1,])/x[3,])
    l <- length(dimension)
    if (l==1) x <- matrix(seq(x[1],x[2],x[3]),ncol=1)
    else {
     text <- paste("x <- expand.grid(",
                      paste("seq(x[1,",1:l,
                            "],x[2,",1:l,
                            "],x[3,",1:l,"])",collapse=","),
                      ")")
     eval(parse(text=text))
   }
  } 

  nd.step <-  ceiling(nrow(x) / 79)
  env <- environment()
  assign("enu",0,envir=env)
  
  covnr <- as.integer(.C("GetModelNr", as.character(model), nr = integer(1))$nr)
  pp <- param
  if (RFparameters()$PracticalRange) {
    pp[ip$SCALE] <- pp[ip$SCALE] / GetPracticalRange(model,pp[-1:(-ip$KAPPA+1)])
  }
  
  storage.mode(pp) <- "double"
  nn <- as.integer(ncol(tgiven))
  np <- as.integer(length(pp))

  switch(krige.method.nr,
         { ## simple kriging
           data <- as.matrix(solve(covmatrix,data-param[ip$MEAN]))
           res <- apply(x, 1,
                        function(z){
                          if ((enu %% nd.step) == 0) cat(pch);
                          assign("enu",enu+1,envir=env);
                          ## to.do replace the following by unchecked !
                          ## CovFct has been check already !
                          .C("UncheckedCovFct",
                             as.double(.C("onePdist",as.double(tgiven),
                                          as.integer(nrow(tgiven)),
                                          as.integer(ncol(tgiven)),
                                          as.double(z),
                                          dx=double(ncol(tgiven)),DUP=FALSE)$dx),
                             nn,covnr,pp,np,cov=double(nn),DUP=FALSE)$cov %*%
                               data}
                        ) + pp[ip$MEAN]
         }, {
           ## ordinary kriging
           covmatrix <- rbind(cbind(covmatrix,1), c(rep(1,nd),0)) 
           data <- solve(covmatrix,rbind(as.matrix(data),0))
           res <-
             apply(x, 1,
                   function(z){
                     if ((enu %% nd.step) == 0) cat(pch);
                     assign("enu",enu+1,envir=env);
                     c(.C("UncheckedCovFct",
                          as.double(.C("onePdist",as.double(tgiven),
                                       as.integer(nrow(tgiven)),
                                       as.integer(ncol(tgiven)),
                                       as.double(z),
                                       dx=double(ncol(tgiven)),DP=FALSE)$dx),
                          nn,covnr,pp,np,cov=double(nn),DUP=FALSE)$cov,
                       1) %*% data 
                   })
         }
         ) # switch

  if (pch!="") cat("\n")
  x <- data <- NULL
  if (is.matrix.data) {
    if (grid) return(array(t(res),dim=c(dimension,ncol.data)))
    else return(t(res))
  } else {
    if (grid) return(array(res,dim=dimension))
    else return(res)
  }
}



CondSimu <- function(krige.method, x, y=NULL, z=NULL, grid, model, param, 
                     method=NULL, n=1, register=0, gridtriple=FALSE,
                     err.model=NULL, err.param=NULL, err.method=NULL,
                     err.register=1,
                     given, data, tol=1E-5, pch=".") {
  if (is.character(method) && (is.na(pmatch(method, c("S","O")))))
    stop("Sorry. The parameters of the function `CondSimu' as been redefined. Use `krige.method' instead of `method'. See help(CondSimu) for details")

  ip<-.C("GetParameterIndices", MEAN=integer(1), VARIANCE=integer(1),
          NUGGET=integer(1), SCALE=integer(1), 
          KAPPA=integer(1), LASTKAPPA=integer(1),integer(1),
          SILL=integer(1), DUP=FALSE)
  ip <- lapply(ip,function(i){i+1})
  
  if (is.null(err.model)) {
    krige.model <- model
    krige.param <- param
  } else {
    if (.C("GetModelNr", as.character(err.model), nr=integer(1))$nr !=
        .C("GetModelNr", as.character("nugget"), nr=integer(1))$nr)
      stop("currently only pure nugget effect allowed as error.model")
    krige.model <- model  ## this has to be replaced by an additive model ->to.do
    krige.param <- param  ## dito
    krige.param[ip$NUGGET] <- krige.param[ip$NUGGET] +
      err.param[ip$VARIANCE] + err.param[ip$NUGGET]
  }
  
  x <- CheckXYZ(x, y, z, grid, gridtriple)
  y <- NULL   
  total <- x$total
  x <-  cbind(x$x, x$y, x$z)
  
  simu.grid <- grid
  given <- as.matrix(given)

  if (grid) {
    ind <- 1 + (t(given) - x[1,]) / x[3,]
    index <-  round(ind)
    dimension<- 1 + floor((x[2,]-x[1,])/x[3,])
    outside.grid <- apply((abs(ind-index)>tol) | (index<1) |
                          (index>dimension),2,any)   
    if (any(outside.grid))
      {
      ## at least some data points are not on the grid
      ## simulate as if there is not grid
      simu.grid <- FALSE
      l <- length(dimension)
      if (l==1) xx <- matrix(seq(x[1],x[2],x[3]),nrow=1)
      else  eval(parse(text=paste("xx <-  t(as.matrix(expand.grid(",
                         paste("seq(x[1,",1:l,
                               "],x[2,",1:l,
                               "],x[3,",1:l,"])",collapse=","),
                         ")))")))
      eval(parse(text=paste("ll <- c(",
                   paste("length(seq(x[1,",1:l,"],x[2,",1:l,"],x[3,",1:l,"]))",
                         collapse=","),
                   ")")))
     # print(ll)
     # print(outside.grid)
     # print(index[,!outside.grid,drop=FALSE])
      new.index <- rep(0,ncol(index))
      if (!all(outside.grid)) {
        new.index[!outside.grid] <- 1 +
          apply((index[,!outside.grid,drop=FALSE]-1) *
                cumprod(c(1,ll[-length(ll)])), 2, sum)
      }
                                        # print(new.index)
      index <- new.index
      new.index <- NULL
    } else {
      ## data points are all lying on the grid
      z <- GaussRF(x, grid=TRUE, model=model, param=param,
                   method=method, n=n, register=register,
                   gridtriple=TRUE)
      index <- t(index)
    }
  } else {
    xx <- t(as.matrix(x)) ## not a grid
    ## the next step can be pretty time consuming!!!
    ## to.do: programme it in C
    one2ncol.xx <- 1:ncol(xx)
    index <- apply(as.matrix(given),1,function(z){
      i <- one2ncol.xx[apply(abs(xx - z),2,sum) < tol]
      if (length(i)==0) return(0)
      if (length(i)==1) return(i)
      return(NA)
    })
  }

  if (!simu.grid) {
    tol <- tol * nrow(xx)
    if (any(is.na(index)))
      stop("identification of the given data points is not unique -- `tol' too large?")
    if (any(notfound <- (index==0))) {
      index[notfound] <- (ncol(xx)+1):(ncol(xx)+sum(notfound))
    }

    xx <- rbind(t(xx), given[notfound,,drop=FALSE])
    z <- GaussRF(xx, grid=FALSE, model=model, param=param,
                 method=method, n=n, register=register)

    xx <- NULL
  }
  if (is.null(z)) stop("random fields simulation failed")

  if (n==1) {
    ## z values at the `given' locations
    zgiven <- z[index]
    z <- z[1:total]
  } else {
    ## this is a bit more complicated since index is either a vector or
    ## a matrix of dimension dim(z)-1
    zgiven <- apply(z,length(dim(z)),function(m) m[index])
    z <- as.vector(apply(z,length(dim(z)),function(m) m[1:total]))
  }

  if (!is.null(err.model)) {
     error <- GaussRF(given, grid=FALSE, model=err.model, param=err.param,
                     method=err.method, n=n, register=err.register)  
     if (is.null(error)) stop("error field simulation failed")
     zgiven <- zgiven + as.vector(error)
     error <- NULL
   }
 
  param[1] <- 0
  z + Kriging(krige.method=krige.method,
              x=x, grid=grid,
              model=krige.model,param=krige.param,
              given=given,data=as.vector(data)-zgiven,
              gridtriple=TRUE, pch=pch)
}  
  
"CovarianceFct" <-
function (x, model, param, dim = 1) 
{
  if (!is.finite(param[1])) {
        param[1] <- 0
    }
  stopifnot(length(x)!=0,
            all(is.finite(x)),
            is.character(model),
            all(is.finite(param)),
            length(dim) == 1,
            is.finite(dim)
            )
  
  covnr <- as.integer(.C("GetModelNr", as.character(model), nr = integer(1))$nr)
    if (covnr < 0) {
        .C("PrintModelList")
        stop("given model cannot be (uniquely) identified from the above list")
    }
    kappas <- as.integer(.C("GetNrParameters", covnr, k = integer(1),DUP=FALSE)$k)
    if (length(param) < 4 + kappas) 
        stop("not enough parameters for the covariance function")
    storage.mode(x) <- "double"
    storage.mode(param) <- "double"
    storage.mode(dim) <- "integer"
    return(.C("Covariance", x, as.integer(length(x)), covnr, param, 
              as.integer(length(param)), dim, res=double(length(x)), DUP = FALSE)$res
           )
}
"DeleteAllRegisters" <-
function () 
{
    .C("DeleteAllKeys")
    return(NULL)
}
"DeleteRegister" <-
function (nr = 0) 
{
   stopifnot( length(nr) == 1, is.finite(nr) ) 
    .C("DeleteKey", as.integer(nr))
    return(NULL)
}
"DoSimulateRF" <-
function (n = 1, register = 0) 
{
  stopifnot(length(n) == 1,
            is.finite(n),
            length(register) == 1,
            is.finite(register)
            )
    DoNameList <- c("DoSimulateRF","DoSimulateRF","DoMaxStableRF")    
    MAXDIM <- 3
    dimensions <- .C("GetKeyInfo", as.integer(register), 
        total = integer(1), len = integer(MAXDIM), dim = integer(1), 
        grid = integer(1), distr = integer(1), DUP=FALSE)
    param <- RFparameters()
    if (dimensions$total <= 0) {
        stop(paste("register", register, "does not look initialized"))
    }
    if (n == 1) {
      x <- .C(DoNameList[1+dimensions$distr], as.integer(register),
                 result=double(dimensions$total), error=integer(1),DUP=FALSE)
      if (x$error) 
        return(NULL)
      else if (dimensions$grid) 
        return(array(x$result, dim = dimensions$len[1:dimensions$dim])) 
      else return(x$result)
    }
    else {
        if (!param$Storing) 
            stop("Use option `RFparameters(Storing=TRUE)' if n>1")

        if (dimensions$grid) {
            result <- array(NA, dim = c(dimensions$len[1:dimensions$dim], 
                n))
            s <- paste(rep(",",dimensions$dim),collapse="")
          }
        else {
          result <- matrix(nrow = dimensions$total, ncol = n)
          s <- ","
        }
        s <- paste("result[", s, "i] <- dummy$res")
        pch <- param$pch
        for (i in 1:n) {
           cat(pch)
           dummy <- .C(DoNameList[1+dimensions$distr], as.integer(register),
                       res=double(dimensions$total),
                       error=integer(1), DUP = FALSE, NAOK = TRUE)
            if (dummy$error) 
                break
            eval(parse(text = s))
        }
        if (pch!="") cat("\n")
        if (dummy$error) 
            return(NULL)
        else return(result)
    }
}
"EmpiricalVariogram" <-
function (x, y = NULL, z = NULL, data, grid, bin, gridtriple = FALSE) 
{
    stopifnot(all(is.finite(data)),
            length(bin)>=2,
            all(is.finite(bin)),
              bin[1] <= 0
            )
  
    new <- CheckXYZ(x, y, z, grid, gridtriple)
    
    repet <- as.integer(length(data)/new$total)
    if (length(data) != new$total * repet) 
        stop("number of data does not match coordinates")
    
    centers <- pmax(0,(bin[-1] + bin[-length(bin)])/2)
    emp.vario <- double(length(bin) - 1)
    .C("empiricalvariogram",
       as.double(new$x), as.double(new$y), as.double(new$z),
       as.integer(new$dim), as.integer(new$l), 
       as.double(data), as.integer(repet), as.integer(grid), 
       as.double(bin), as.integer(length(bin) - 1), as.integer(0), 
       emp.vario, DUP = FALSE)
    emp.vario[is.na(emp.vario) & (centers==0)] <- 0
    return(list(centers=centers, emp.vario=emp.vario))
}
".First.lib" <- function (lib, pkg) 
{
    library.dynam("RandomFields", pkg, lib)
    library(mva)
}
"InitGaussRF" <-
  function(x, y = NULL, z = NULL, grid, model, param, method = NULL, 
    register = 0, gridtriple = FALSE) 
{
  return(InitSimulateRF(x=x, y=y, z=z, grid=grid, model=model,
                        param=param, method=method, register=register,
                        gridtriple=gridtriple, distribution="Gauss"))
}

CheckXYZ <- function(x, y, z, grid, gridtriple){
  stopifnot(length(x) != 0,
            is.logical(grid),
            is.logical(gridtriple)
            )
  if (is.data.frame(x)) {
    if (ncol(x)==1) x <- as.vector(x) else x <- as.matrix(x)
  }
  stopifnot(all(is.finite(x))) 
  if (is.matrix(x)) {
        if (!is.null(y) || !is.null(z)) {
            stop("If x is a matrix, then y and z may not be given")
        }
        dim <- ncol(x)
        if (dim > 3) {
            dim <- 3
            warning("only the first three columns are considered as coordinates")
        }
        if (dim > 1) {
            y <- x[, 2]
            if (dim > 2) 
                z <- x[, 3]
        }
        x <- x[, 1]
        l <- length(x)
    }
    else {
        l <- length(x)
        storage.mode(x) <- "double"
        if (is.null(y)) {
            dim <- 1
            if (!is.null(z)) {
                stop("y is not given, but z")
            }
        }
        else {
            if ((l != length(y)) && (!grid || gridtriple)) {
                stop("x and y differ in length")
            }
            stopifnot(all(is.finite(y)))
            storage.mode(y) <- "double"
            if (is.null(z)) {
                dim <- 2
            }
            else {
                dim <- 3
                if (l != length(z) && (!grid || gridtriple)) {
                  stop("x and z differ in length")
                }
                stopifnot(all(is.finite(z))) 
                storage.mode(z) <- "double"
            }
        }
    }
    if (l == 1) 
        stop("Use grid=FALSE if only a single point is simulated")
    else if (grid) {
        if (gridtriple) {
            if (l != 3) {
                stop("In case of simulating a grid with option gridtriple, exactly 3 numbers are needed for each direction")
            }
         lx <- length(seq(x[1],x[2],x[3]))
         ly <- lz <- 1
         x[2] <- x[1] + (lx - 0.999) * x[3]
         if (dim > 1) {
             ly <- length(seq(y[1],y[2],y[3]))
             y[2] <- y[1] + (ly - 0.999) * y[3]
             if (dim > 2) {
               lz <- length(seq(z[1],z[2],z[3]))
               z[2] <- z[1] + (lz - 0.999) * z[3]
              }
           }
         total <- lx * ly * lz
       if (total==0) stop("incorrect grid specification (one or no points)")
        
        }
        else {
            eqdist <- function(x) {
                step <- diff(x)                
                if (sum(abs(step - step[1]))/step[1] > 1e-10) {
                  stop("Grid must have equal distances in each direction.")
                }
                return(c(x[1], x[length(x)]+0.001*step[1], step[1]))
            }
            total <- length(x)
            x <- eqdist(x)            
            if (dim > 1) {
                total <- total * length(y)
                y <- eqdist(y)
                if (dim > 2) {
                  total <- total * length(z)
                  z <- eqdist(z)
                }
            }
            l <- 3
        }
    }
    else {
        total <- length(x)
    }
  return(list(x=x, y=y, z=z, total=total, l=l, dim=dim))
}

"InitSimulateRF" <-
function (x, y = NULL, z = NULL, grid, model, param, method = NULL, 
          register = 0, gridtriple = FALSE,distribution=NA)  
{
   distributionList <- c("Gauss","Poisson","MaxStable")
   InitNameList <- c("InitSimulateRF","InitSimulateRF","InitMaxStableRF")
   if (is.na(distribution)) {
     stop("This function is an internal function.\nUse `GaussRF', `InitGaussRF', `MaxStableRF', etc., instead.\n")    
   }
   distrNr <- pmatch(distribution,distributionList)   
   if (is.na(distrNr)) stop("Unknown distribution -- do not use this internal function directly")
   else {
     if ((distrNr==2) && !exists("PoissonRF")) stop("Sorry. Not programmed yet.")
     InitName <- InitNameList[distrNr]
     distrNr <- distrNr - 1;
   }
   stopifnot(is.character(model)) 
    covnr <- as.integer(.C("GetModelNr", as.character(model), 
        nr = integer(1))$nr)
    if (covnr < 0) {
        .C("PrintModelList")
        stop("given model cannot be (uniquely) identified from the above list")
    }
    if (is.null(method)) {
        method <- -1
    }
    else {
        if (!is.character(method)) {
            stop("method must be NULL or a string")
        }
        method <- as.integer(.C("GetMethodNr", method, nr = integer(1))$nr)
        if (method < 0) {
          .C("PrintMethods")
          stop("given method cannot be (uniquely) identified from the above list")
        }
    }
    kappas <- as.integer(.C("GetNrParameters", covnr, k = integer(1))$k)
    if (length(param) < 4 + kappas) 
        stop("not enough parameters for the covariance function")
    storage.mode(param) <- "double"

    new <- CheckXYZ(x, y, z, grid, gridtriple)
   
    return(
           max(0,.C(InitName, as.double(new$x), as.double(new$y),
                    as.double(new$z), as.integer(new$dim), as.integer(new$l), 
              as.integer(grid), covnr, param, as.integer(length(param)), 
              as.integer(method), as.integer(distrNr),
                    as.integer(register), error=integer(1),
              DUP = FALSE)$error)
    )
}
"PrintModelList" <-
function () 
{
    .C("PrintModelList")
    return(NULL)
}
"PrintMethodList" <-
function () 
{
    .C("PrintMethods")
    return(NULL)
}
"RFparameters" <-function (...) {
  RFparameters.default(...)
}
"RFparameters.default" <-
function (Storing = storing, PrintLevel = printlevel,
          PracticalRange = practicalrange, 
    CE.force = ce.force, CE.mmin = ce.mmin, CE.tolRe = ce.tolRe, 
    CE.tolIm = ce.tolIm, CE.trials = ce.trials,
    direct.checkprecision = directcheckprecision, 
    direct.maxvariables = directmaxvariables,
    direct.method = directmethod, 
    direct.requiredprecision = directrequiredprecision,
    spectral.lines = spectrallines, 
    spectral.grid = spectralgrid, TBMCE.force = tbmceforce, TBMCE.mmin = tbmcemmin, 
    TBMCE.tolRe = tbmcetolre, TBMCE.tolIm = tbmcetolim, TBMCE.trials = tbmcetrials, 
    TBM2.lines = tbm2lines, TBM2.linesimufactor = tbm2linesimufactor, 
    TBM2.linesimustep = tbm2linesimustep, TBM3D2.lines = tbm3D2lines, 
    TBM3D2.linesimufactor = tbm3D2linesimufactor, TBM3D2.linesimustep = tbm3D2linesimustep, 
    TBM3D3.lines = tbm3D3lines, TBM3D3.linesimufactor = tbm3D3linesimufactor, 
    TBM3D3.linesimustep = tbm3D3linesimustep,
    MPP.approxzero=mppapproxzero, add.MPP.realisations=addmpprealisations,
    MPP.radius=mppradius,
    maxstable.maxGauss=maxstablemaxGauss,
    pch = pchx
 ) {
    x <- .C("SetParam", as.integer(1), storing = integer(1), 
            printlevel = integer(1), practicalrange = integer(1),
            pch="  ")
    storing <- x$storing
    printlevel <- x$printlevel
    practicalrange <- as.logical(x$practicalrange)
    if (!is.finite(Storing + PrintLevel + PracticalRange)) 
        stop("some parameters are not finite")
    pchx <- x$pch
    if (!is.character(pch)) stop("pch is not a character")
   
    ## do not allow integer values for users!
    PracticalRange <- as.logical(PracticalRange)
    .C("SetParam", as.integer(0), as.integer(Storing), as.integer(PrintLevel), 
        as.integer(PracticalRange), pch)
    
    x <- .C("SetParamCircEmbed", as.integer(1), force = integer(1), 
        tolRe = double(1), tolIm = double(1), trials = integer(1), 
        mmin = integer(1))
    ce.force <- x$force
    ce.tolRe <- x$tolRe
    ce.tolIm <- x$tolIm
    ce.trials <- x$trials
    ce.mmin <- x$mmin
    if (!is.finite(CE.force + CE.mmin + CE.tolRe + CE.tolIm + 
        CE.trials)) 
        stop("some parameters are not finite")
    .C("SetParamCircEmbed", as.integer(0), as.integer(CE.force), 
        as.double(CE.tolRe), as.double(CE.tolIm), as.integer(CE.trials), 
        as.integer(CE.mmin))
    x <- .C("SetParamTBMCE", as.integer(1), tbmceforce = integer(1), 
        tbmcetolre = double(1), tbmcetolim = double(1), tbmcetrials = integer(1), 
        tbmcemmin = integer(1))
    tbmceforce <- x$tbmceforce
    tbmcetolre <- x$tbmcetolre
    tbmcetolim <- x$tbmcetolim
    tbmcetrials <- x$tbmcetrials
    tbmcemmin <- x$tbmcemmin
    if (!is.finite(TBMCE.force + TBMCE.mmin + TBMCE.tolRe + TBMCE.tolIm + 
        TBMCE.trials)) 
        stop("some parameters are not finite")
    .C("SetParamTBMCE", as.integer(0), as.integer(TBMCE.force), 
        as.double(TBMCE.tolRe), as.double(TBMCE.tolIm), as.integer(TBMCE.trials), 
        as.integer(TBMCE.mmin))
    x <- .C("SetParamTBM2", as.integer(1), tbm2lines = integer(1), 
        tbm2linesimufactor = double(1), tbm2linesimustep = double(1))
    tbm2lines <- x$tbm2lines
    tbm2linesimufactor <- x$tbm2linesimufactor
    tbm2linesimustep <- x$tbm2linesimustep
    if (!is.finite(TBM2.lines + TBM2.linesimufactor + TBM2.linesimustep)) 
        stop("some parameters are not finite")
    .C("SetParamTBM2", as.integer(0), as.integer(TBM2.lines), 
        as.double(TBM2.linesimufactor), as.double(TBM2.linesimustep))
    x <- .C("SetParamTBM3D2", as.integer(1), tbm3D2lines = integer(1), 
        tbm3D2linesimufactor = double(1), tbm3D2linesimustep = double(1))
    tbm3D2lines <- x$tbm3D2lines
    tbm3D2linesimufactor <- x$tbm3D2linesimufactor
    tbm3D2linesimustep <- x$tbm3D2linesimustep
    if (!is.finite(TBM3D2.lines + TBM3D2.linesimufactor + TBM3D2.linesimustep)) 
        stop("some parameters are not finite")
    .C("SetParamTBM3D2", as.integer(0), as.integer(TBM3D2.lines), 
        as.double(TBM3D2.linesimufactor), as.double(TBM3D2.linesimustep))
    x <- .C("SetParamTBM3D3", as.integer(1), tbm3D3lines = integer(1), 
        tbm3D3linesimufactor = double(1), tbm3D3linesimustep = double(1))
    tbm3D3lines <- x$tbm3D3lines
    tbm3D3linesimufactor <- x$tbm3D3linesimufactor
    tbm3D3linesimustep <- x$tbm3D3linesimustep
    if (!is.finite(TBM3D3.lines + TBM3D3.linesimufactor + TBM3D3.linesimustep)) 
        stop("some parameters are not finite")
    .C("SetParamTBM3D3", as.integer(0), as.integer(TBM3D3.lines), 
        as.double(TBM3D3.linesimufactor), as.double(TBM3D3.linesimustep))
    x <- .C("SetParamSpectral", as.integer(1), spectrallines = integer(1), 
        spectralgrid = integer(1))
    spectrallines <- x$spectrallines
    spectralgrid <- x$spectralgrid
    if (!is.finite(spectral.lines + spectral.grid)) 
        stop("some parameters are not finite")
    .C("SetParamSpectral", as.integer(0), as.integer(spectral.lines), 
        as.integer(spectral.grid))
    x <- .C("SetParamDirectGauss", as.integer(1), directmethod = integer(1), 
        directcheckprecision = integer(1), directrequiredprecision = double(1), 
        directmaxvariables = integer(1))
    directmethod <- x$directmethod
    directcheckprecision <- x$directcheckprecision
    directrequiredprecision <- x$directrequiredprecision
    directmaxvariables <- x$directmaxvariables
    if (!is.finite(direct.checkprecision + direct.maxvariables + 
        direct.method + direct.requiredprecision)) 
        stop("some parameters are not finite")
    .C("SetParamDirectGauss", as.integer(0), as.integer(direct.method), 
        as.integer(direct.checkprecision), as.double(direct.requiredprecision), 
        as.integer(direct.maxvariables))

    x <- .C("SetMPP", as.integer(1),
            mppapproxzero=double(1), addmpprealisations=double(1),
            mppradius=double(1))
    mppapproxzero <- x$mppapproxzero;
    addmpprealisations <- x$addmpprealisations;
    mppradius <- x$mppradius;
    if (!is.finite(MPP.approxzero+add.MPP.realisations+MPP.radius)) 
        stop("some parameters are not finite")
    .C("SetMPP", as.integer(0),
       as.double(MPP.approxzero),as.double(add.MPP.realisations),
       as.double(MPP.radius))
  
    
    x <- .C("SetExtremes", as.integer(1),
            maxstablemaxGauss=double(1))
    maxstablemaxGauss <- x$maxstablemaxGauss;
    if (!is.finite(maxstable.maxGauss)) 
        stop("some parameters are not finite")
    .C("SetExtremes", as.integer(0),
       as.double(maxstable.maxGauss))
     
    if (length(as.list(match.call()))>1) return(NULL)
    else return(list(Storing = as.logical(Storing), PrintLevel = PrintLevel, 
        PracticalRange = as.logical(PracticalRange),
        CE.force = as.logical(CE.force), CE.mmin = CE.mmin,
        CE.tolRe = CE.tolRe, CE.tolIm = CE.tolIm, CE.trials = CE.trials, 
        direct.checkprecision = as.logical(direct.checkprecision), 
        direct.maxvariables = direct.maxvariables,
        direct.method = direct.method,
        direct.requiredprecision = direct.requiredprecision, 
         spectral.lines = spectral.lines, 
        spectral.grid = as.logical(spectral.grid),
        TBMCE.force = as.logical(TBMCE.force), TBMCE.mmin = TBMCE.mmin, 
        TBMCE.tolRe = TBMCE.tolRe, TBMCE.tolIm = TBMCE.tolIm, 
        TBMCE.trials = TBMCE.trials, 
        TBM2.lines = TBM2.lines, TBM2.linesimufactor = TBM2.linesimufactor, 
        TBM2.linesimustep = TBM2.linesimustep, TBM3D2.lines = TBM3D2.lines, 
        TBM3D2.linesimufactor = TBM3D2.linesimufactor, TBM3D2.linesimustep = TBM3D2.linesimustep, 
        TBM3D3.lines = TBM3D3.lines, TBM3D3.linesimufactor = TBM3D3.linesimufactor, 
        TBM3D3.linesimustep = TBM3D3.linesimustep,
        MPP.approxzero=MPP.approxzero, add.MPP.realisations=add.MPP.realisations,
        MPP.radius=MPP.radius,
        maxstable.maxGauss=maxstable.maxGauss,
        pch=pch
        ))
}
"GaussRF" <-
function (x, y = NULL, z = NULL, grid, model, param, method = NULL, 
    n = 1, register = 0, gridtriple = FALSE) 
{
    if (InitSimulateRF(x=x, y=y, z=z, grid=grid, model=model,
                       param=param, method=method, register=register,
                       gridtriple=gridtriple, distribution="Gauss") <= 0) {
        return(DoSimulateRF(n=n, reg=register))
    }
    else {
        return(NULL)
    }
}
"SimulateRF" <-
function (x, y = NULL, z = NULL, grid, model, param, method = NULL, 
    n = 1, register = 0, gridtriple = FALSE, distribution=NA) 
{
    if (InitSimulateRF(x=x, y=y, z=z, grid=grid, model=model,
                       param=param, method=method, 
                       reg=register, gridtriple=gridtriple,
                       distr=distribution) <= 0) {
        return(DoSimulateRF(n=n, reg=register))
    }
    else {
        return(NULL)
    }
}
"Variogram" <-
function (x, model, param, dim = 1) 
{
  if (!is.finite(param[1])) param[1] <- 0
  stopifnot(length(x) != 0,
            all(is.finite(x)),
            is.character(model),
            all(is.finite(param)),
            length(dim) == 1,
            is.finite(dim)) 
    covnr <- as.integer(.C("GetModelNr", as.character(model), nr = integer(1))$nr)
    if (covnr < 0) {
        .C("PrintModelList")
        stop("given model cannot be (uniquely) identified from the above list")
    }
    kappas <- as.integer(.C("GetNrParameters", covnr, k = integer(1))$k)
    if (length(param) < 4 + kappas) 
        stop("not enough parameters for the variogram model")
    storage.mode(x) <- "double"
    storage.mode(param) <- "double"
    storage.mode(dim) <- "integer"
    return(
           .C("Variogram", x, as.integer(length(x)), covnr, param, as.integer(length(param)), 
              dim, res=double(length(x)), DUP = FALSE)$res
           )
}

GetModelNames <- function() {
  p <- .C("GetrfParameters",covmaxchar=integer(1),methodmaxchar=integer(1),
          covnr=integer(1),methodnr=integer(1))
  l <- character(p$covnr)
  for (i in 1:p$covnr) {
    l[i] <- .C("GetModelName",as.integer(i-1),
               n=paste(rep(" ",p$covmaxchar),collapse=""))$n
  }
  return (l)
}


GetMethodNames <- function() {
  p <- .C("GetrfParameters",covmaxchar=integer(1),methodmaxchar=integer(1),
          covnr=integer(1),methodnr=integer(1))
  l <- character(p$methodnr)
  for (i in 1:p$methodnr) {
    l[i] <- .C("GetMethodName",as.integer(i-1),
               n=paste(rep(" ",p$methodmaxchar),collapse=""))$n
  }
  return(l)
}


GetPracticalRange <- function(model,kappas=NULL) {
  covnr <- as.integer(.C("GetModelNr", as.character(model),
                         nr = integer(1))$nr)
  if (covnr < 0) {
    .C("PrintModelList")
    stop("given model cannot be (uniquely) identified from the above list")
  }
  if (length(kappas)!=.C("GetNrParameters", covnr, k = integer(1),DUP=FALSE)$k)
    stop("incorrect number of parameters!")
  nat.scl <- double(1)
  error <- integer(1)
  .C("GetNaturalScaling",
     as.integer(covnr),
     as.double(c(0,1,0,1,kappas)),         ## not stable w.r.t. to changings !!
     as.integer(4+ length(kappas)),        ## dito
     as.integer(11),
     nat.scl,
     error,DUP=FALSE)
  if (error) stop("natural scaling could not be obtained")
  return(1.0 / nat.scl)
}


