diff --git a/DESCRIPTION b/DESCRIPTION
index 0fe1f73599dccb001f322855749bb51eb9735491..71c365593e0eba075c53944dd4075dc27997e1c4 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: gmGeostats
-Version: 0.11.0-9001
-Date: 2021-11-05
+Version: 0.11.0-9002
+Date: 2021-12-08
 Title: Geostatistics for Compositional Analysis
 Authors@R: c(person(given = "Raimon", 
         family = "Tolosana-Delgado", 
@@ -60,7 +60,6 @@ VignetteBuilder: knitr
 LazyData: true
 Collate: 
     'Anamorphosis.R'
-    '__preliminary_3D_Rosie.R'
     'preparations.R'
     'gstatCompatibility.R'
     'gmAnisotropy.R'
diff --git a/NAMESPACE b/NAMESPACE
index b17e337f0b52dc811374b4ff9566c10d147ed65d..091250b1c8f4cbffceafa762627eb5cd6cf44505 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -24,8 +24,6 @@ S3method(UWEDGE,default)
 S3method(UWEDGE,rcomp)
 S3method(accuracy,DataFrameStack)
 S3method(accuracy,data.frame)
-S3method(as.AnisotropyScaling,AnisotropyScaling)
-S3method(as.AnisotropyScaling,numeric)
 S3method(as.CompLinModCoReg,CompLinModCoReg)
 S3method(as.CompLinModCoReg,LMCAnisCompo)
 S3method(as.DataFrameStack,array)
@@ -136,6 +134,8 @@ S3method(xvErrorMeasures,DataFrameStack)
 S3method(xvErrorMeasures,data.frame)
 S3method(xvErrorMeasures,default)
 export("stackDim<-")
+export(AnisotropyRangeMatrix)
+export(AnisotropyScaling)
 export(CholeskyDecomposition)
 export(DSpars)
 export(DataFrameStack)
@@ -153,7 +153,10 @@ export(accuracy)
 export(ana)
 export(anaBackward)
 export(anaForward)
-export(anis2D.par2A)
+export(anis2D_par2A)
+export(anis3D_par2A)
+export(anis_GSLIBpar2A)
+export(as.AnisotropyRangeMatrix)
 export(as.AnisotropyScaling)
 export(as.CompLinModCoReg)
 export(as.DataFrameStack)
diff --git a/NEWS.md b/NEWS.md
index eae19e3d28e65f74b283f966717cc7ab588b5e22..0b87a47263b4a8f7efbdec1cec8368627902fafc 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,8 @@
+# gmGeostats 0.11.0-9002
+
+* (2021-12-08) anisotropy objects complete, restructured and everywhere correctly used: if `class(A)=="AnisotropyScaling"` and `class(M)=="AnisotropyRangeMatrix"` (the two possible classes), then `A %*% t(A) == solve(M)`; with them, `u = sqrt(h %*% solve(M, h))` can be used in a scalar variogram (i.e. `M` is a matrix of ranges and orientations), and `v = t(A) %*% h` an isotropic variogram (i.e. `u=sqrt(sum(v^2))`,  and `A` is a scaling matrix that makes the space isotropic)
+* (2021-12-08) bugs corrected in `gsi.EVario3D()`
+ 
 # gmGeostats 0.11.0-9001
 
 * (2021-11-05) bugs corrected in internal 2D empirical variogram function `gsi.EVario2D()` . First version of 3D internal variogram function `gsi.EVario3D()` available. Usage of these functions is only for specialists foreseen. In the near future a user-friendlier wrapper will be provided.
diff --git a/R/compositionsCompatibility.R b/R/compositionsCompatibility.R
index 86e9b81a734321c2d40460d848878a02a6b7477b..b8fd5314b7ffb4d2f9191732d952bf2aef56b606 100644
--- a/R/compositionsCompatibility.R
+++ b/R/compositionsCompatibility.R
@@ -553,6 +553,8 @@ fit_lmc.logratioVariogramAnisotropy <- function(v, g, model, ...){
 
 #################################################
 ## anisotropic variogram models and accesory functions 
+# h : matrix of nx3 (or nx2) lag vectors
+# A : isotropizing matrix
 anish2Dist <- function (h, A = NULL){
   if (is.data.frame(h)) 
     h <- as.matrix(h)
@@ -606,7 +608,8 @@ anvgram.nugget = function (h, nugget = 0, sill = 0, range = 1, A = NULL, ...){
 #'
 #' @param Z compositional data set, used to derive the compositional dimension and colnames
 #' @param models string (or vector of strings) specifying which reference model(s) to use 
-#' @param azimuths typically a vector providing, for each model, the direction of maximal continuity
+#' @param azimuths typically a vector providing, for each model, the direction 
+#' of maximal continuity (measured from North clockwise)
 #' @param ranges typically a `G`-column matrix providing the minimal and maximal ranges, with one row per model
 #' (with `G` specified below)
 #' @param sillarray array of sills for each model. It can be null (to be estimated in the future). If specified, 
@@ -634,6 +637,12 @@ LMCAnisCompo = function(Z, models=c("nugget", "sph", "sph"),
   DD = D*(D+1)/2
   dd = d*D/2
   K = length(models)
+  # check dimension of azimuths
+  if(length(dim(azimuths))==0){
+    dim(azimuths) = c(length(azimuths),1)
+  }else if( !(ncol(azimuths) %in% c(1,3))){
+    stop("LMCAnisCompo: `azimuths` should be a (column)vector in 2D, or a 3-column matrix in 3D")
+  }
   
   # consider the sill matrices and try to interpret its shapes
   if(is.null(sillarray)){
@@ -676,7 +685,7 @@ LMCAnisCompo = function(Z, models=c("nugget", "sph", "sph"),
     sill = darstellung(sillarray, i)
     range = ranges[i,2]
     # if(G==2)
-    A = anis2D.par2A(ratio=ranges[i,2]/ranges[i,1], angle=azimuths[i])
+    A = anis_GSLIBpar2A(ratios=ranges[i,-1]/ranges[i,1], angles=azimuths[i,] )
     # elseif(G==3) and so on...
     rs = list(model=model, range=range, A=A, sill=sill)
     class(rs) = "variostructure"
@@ -825,7 +834,7 @@ as.CompLinModCoReg.LMCAnisCompo = function(v, ...){
 #     vvgg = rep(par[1], min(length(hh), nrow(hh)))
 #     for(istruc in 1:length(variomodelpars$models)){
 #       fn = paste("anvgram", variomodelpars$models[istruc], sep=".")
-#       A = anis2D.par2A(par[istruc*4+2], par[istruc*4+3])
+#       A = anis2D_par2A(par[istruc*4+2], par[istruc*4+3])
 #       g = eval(call(fn, h=hh, range=par[istruc*4], sill=par[istruc*4+1], A=A))
 #       # g = eval(call(fn, h=hh, range=par[istruc*4], sill=par[istruc*4+1]))
 #       vvgg = vvgg + g
@@ -1289,7 +1298,7 @@ as.gmCgram.LMCAnisCompo = function(m, V=attr(m,"contrasts"),
   }
   # put anisotropy L matrix
   inonugget = which(m["model",]!="nugget")
-  Ms = sapply(inonugget, function(i) gsi.powM(m["A",i][[1]], -2))
+  Ms = sapply(inonugget, function(i) as.AnisotropyRangeMatrix.AnisotropyScaling(m["A",i][[1]]) )
   dim(Ms) = c(dim(m["A",1][[1]]),length(inonugget))
   output$M = aperm(Ms, c(3,1,2))
   # put sills
diff --git a/R/gmAnisotropy.R b/R/gmAnisotropy.R
index 83eef9b96e2e4ec293b2e7e044206e392ebc7731..3e2ad18cde33c6be7cf371f82928866931562231 100644
--- a/R/gmAnisotropy.R
+++ b/R/gmAnisotropy.R
@@ -1,20 +1,6 @@
+
+
 #### Anisotropy ------------------
-#' Check for any anisotropy class
-#' 
-#' Check that an object contains a valid specification of anisotropy, in any form 
-#'
-#' @param x object to check
-#'
-#' @return a logical, TRUE if the object is an anisotropy specification; FALSE otherwise 
-#' @export
-#'
-#' @examples
-#' a =  anis2D.par2A(0.5, 30)
-#' a
-#' is.anisotropySpecification(a)
-is.anisotropySpecification = function(x){
-  "Anisotropy" %in% class(x)  |  inherits(x, "Anisotropy")
-} 
 
 
 #' Convert to anisotropy scaling matrix
@@ -24,7 +10,7 @@ is.anisotropySpecification = function(x){
 #' @param x an object convertible to an anisotropy scaling matrix; see details
 #'
 #' @return A matrix \eqn{A} such that for any lag vector \eqn{h}, the variogram model turns 
-#'  isotropic in terms of \eqn{u=A\cdot h}. 
+#'  isotropic in terms of \eqn{u'=h'\cdot A}. 
 #'
 #' @details Method `as.AnisotropyScaling.numeric()` expects a vector of two numbers in 2D, 
 #' or a vector of 5 numbers in 3D. These are in 2D, the azimuth of maximum continuity (in 
@@ -32,66 +18,228 @@ is.anisotropySpecification = function(x){
 #' these are: 1,2) the azimuth and the dip of the direction of maximal continuity; 3) the 
 #' angle of rotation around the axis of the first direction; 4,5) the anisotropy ratios of 
 #' the ranges of the second/first and third/first directions of maximal continuity. All angles
-#' are given in degrees, all ratios must be smaller or equal to 1. 
+#' are given in degrees, all ratios must be smaller or equal to 1. This follows gstat (and hence
+#' GSlib) conventions; see gstat::vgm() for details.
 #'    
 #' @export
 #' @examples
-#' as.AnisotropyScaling(c(30,0.5))
-as.AnisotropyScaling <- function(x) UseMethod("as.AnisotropyScaling", x)
-
+#' ( l = anis_GSLIBpar2A(angles=30, ratios=0.5)) 
+#' ( ll = unclass(l) )
+#' AnisotropyScaling(l)
+AnisotropyScaling = function(x){
+  class(x) = "AnisotropyScaling"
+  return(x)
+}
 
-#' @describeIn as.AnisotropyScaling Convert to anisotropy scaling matrix
-#' @method as.AnisotropyScaling AnisotropyScaling
+#' @describeIn AnisotropyScaling coerce function to anisotropy scaling matrix
 #' @export
-as.AnisotropyScaling.AnisotropyScaling = function(x) x
+as.AnisotropyScaling <- function(x){ UseMethod("as.AnisotropyScaling", x) }
+
+# @describeIn AnisotropyScaling identity
+# @method as.AnisotropyScaling AnisotropyScaling
+# @export
+as.AnisotropyScaling.AnisotropyScaling = function(x){
+  x
+} 
+
+
+# @describeIn AnisotropyScaling  from a vector of numbers
+# @method as.AnisotropyScaling numeric
+# @export
+as.AnisotropyScaling.double <-as.AnisotropyScaling.numeric <- function(x){
+  if(length(dim(x))==0){
+    if(length(x)==2) return(anis2D_par2A(ratio=x[2], angle=x[1], inv=TRUE))
+    if(length(x)==5) return(anis3D_par2A(ratios=x[4:5], angles=x[1:3], inv=TRUE))
+    stop("as.AnisotropyScaling.numeric: only works for length=2 1.angle+1.ratio, or lenth=5 3.angles+2.ratios")
+  }else{
+    warning("as.AnisotropyScaling.numeric: expected vector, but matrix given; attemting an interpretation as AnisotropyScaling or as AnisotropyRangeMatrix")
+    if( all(abs(x-t(x))<1e-12) ){
+      return(as.AnisotropyScaling.AnisotropyRangeMatrix(x))
+    }else{
+      return(AnisotropyScaling(x))
+    }
+  }
+}
+
+# @describeIn AnisotropyScaling  from an AnisotropicRangeMatrix
+# @method as.AnisotropyScaling AnisotropyRangeMatrix
+# @export
+as.AnisotropyScaling.AnisotropyRangeMatrix = function(x){
+  AnisotropyScaling(solve(chol(x)))
+}
+
 
 
-#' @describeIn as.AnisotropyScaling  Convert to anisotropy scaling matrix
-#' @method as.AnisotropyScaling numeric
+#' Force a matrix to be anisotropy range matrix,
+#' 
+#' Force a matrix M to be considered an anisotropy range matrix, i.e 
+#' with ranges and orientations, 
+#' such that \eqn{u = sqrt(h' * M^{-1} * h)} allows to use an isotropic
+#' variogram.
+#'
+#' @param x  matrix simmetric positive definite (i.e. M above)
+#' @param checkValidity boolean, should validity be checked?
+#'
+#' @return the same matrix with a class attribute
 #' @export
-as.AnisotropyScaling.numeric = function(x){
-  if(length(x)==2) return(anis2D.par2A(ratio=x[2], angle=x[1]))
-  if(length(x)==5){
-    if(sum( (x[4:5]-1)^2)<1e-12){
-      return(diag(3))
+#'
+AnisotropyRangeMatrix = function(x, checkValidity=TRUE){
+  if(checkValidity){
+    if(length(dim(x))==3){
+      odim = dim(x)
+      x = t(sapply(1:(odim[1]), function(i){ AnisotropyRangeMatrix(x[i,,])} ))
+      dim(x) = odim
+    }else{
+      if(nrow(x)!=ncol(x)) stop("AnisotropyRangeMatrix: square matrix needed")
+      if(any(abs(x-t(x))>1e-12 )) stop("AnisotropyRangeMatrix: symmetric matrix needed")
+      if(any(eigen(x, only.values = TRUE)$values<(-1e-12))) stop("AnisotropyRangeMatrix: positive definite matrix needed")
     }
-    stop("as.AnisotropyScaling: 3D from arbitrary 5-vector values not yet implemented") # return(anis3D.par2A(ratios=x[4:5], angles=x[1:3]))
   }
-  stop("as.AnisotropyScaling.numeric: only works for length=2 1.angle+1.ratio, or lenth=5 3.angles+2.ratios")
+  class(x) = "AnisotropyRangeMatrix"
+  return(x)
 }
 
+#' @describeIn AnisotropyRangeMatrix  Convert to anisotropy range matrix
+#' @export
+as.AnisotropyRangeMatrix <- function(x){ UseMethod("as.AnisotropyRangeMatrix", x) }
+
+# @describeIn AnisotropyRangeMatrix  Convert to anisotropy range matrix
+# @method as.AnisotropyRangeMatrix default
+# @export
+as.AnisotropyRangeMatrix.default <- function(x) as.AnisotropyRangeMatrix(as.AnisotropyScaling(x))
+
+
+# @describeIn AnisotropyRangeMatrix  Convert to anisotropy range matrix
+# @method as.AnisotropyRangeMatrix AnisotropyRangeMatrix
+# @export
+as.AnisotropyRangeMatrix.AnisotropyRangeMatrix <- function(x) x
+
+# @describeIn AnisotropyRangeMatrix  Convert to anisotropy range matrix
+# @method as.AnisotropyRangeMatrix AnisotropyScaling
+# @export
+as.AnisotropyRangeMatrix.AnisotropyScaling <- function(x) AnisotropyRangeMatrix(solve(tcrossprod(x)), checkValidity=FALSE)
 
-#' Produce anisotropy matrix from angle and anisotropy ratios
+
+
+
+
+
+#' Produce anisotropy scaling matrix from angle and anisotropy ratios
 #' 
-#' Produce anisotropy matrix (in Cholesky decomposition) from angle and anisotropy ratios
+#' Produce anisotropy matrix (as the transposed of the Cholesky 
+#' decomposition) from angle and anisotropy ratios
 #' 
-#' @param ratio an anisotropy ratio (min/max range)
-#' @param angle direction of maximum range, i.e. largest spatial continuity, measured 
-#' counterclockwise from East
+#' @param ratios vector of two values between 0 and 1 giving the anisotropy ratios of 
+#' medium/largest smallest/largest ranges
+#' @param angles, as defined in gstat::vgm (and indeed GSLIB)
+
+#' @param inv boolean or integer, see `return` for details
 #' 
-#' @return a 3x3 matrix of anisotropy, in Cholesky form
+#' @return a 3x3 matrix of anisotropy. 
+#' 
+#' If `inv=TRUE` (or 1) the output is a matrix `A` such that `norm(h %*% A)` 
+#' allows to use isotropic variograms, being `h = c(hx, hy, hz)` the lag vectprs. 
+#' 
+#' If `inv=FALSE` (or 0) the output is a matrix `A` such that `norm(h %*% solve(A))` 
+#' allows to use isotropic variograms.  
+#' 
+#' Other values are meaningless.
 #'
 #' @export
-anis2D.par2A = function(ratio, angle){
-  ## Cholesky decomp of the rotation matrix from angle of maximum
+anis_GSLIBpar2A = function(ratios, angles, inv=FALSE){
+  if( (length(ratios)==1) & length(angles)==1){
+    anis2D_par2A(ratio=ratios, angle=angles, inv=inv)
+  }else if( (length(ratios)==1) & length(angles)==1){
+    anis3D_par2A(ratios=ratios, angles=angles, inv=inv)
+  }else{
+    stop("anis.GSLIBpar2A: error, ratios and angles length incompatible for both 2D and 3D")
+  }
+}
+
+
+
+
+
+#' @describeIn anis_GSLIBpar2A 2D case
+#' @param ratio an anisotropy ratio (min/max range)
+#' @param angle direction of maximum range, i.e. largest spatial continuity, measured 
+#' clockwise from North
+#' @export
+#' @examples
+#' ## ratio=0.5, azimuth 30?? (i.e. direction 60??) 
+#' A = anis2D_par2A(1, 30)
+#' A
+#' AAt = A %*% t(A)
+#'  #  project the bisector 1:1 (i.e. 45??)
+#' (k = c(1,1,0) %*% A)
+#' atan2(k[2], k[1]) * 180/pi  # should be 15
+#' sqrt(sum(k^2))
+#' sqrt( c(1,1,0) %*% AAt %*% c(1,1,0) )
+#' A = anis2D_par2A(0.5, 60)
+#' rd = 60 * pi/180
+#' A
+#' A %*% t(A)
+#' c(cos(rd), sin(rd),0) %*% A #  should be 1
+#' c(-sin(rd), cos(rd),0) %*% A #  should be +/- sqrt(2)
+anis2D_par2A = function(ratio, angle, inv=FALSE){
+  ## transposed Cholesky decomp of the rotation matrix from angle of maximum
   #     spatial continuity + ratio min/max ranges
   a = angle * pi/180
-  R = matrix(c(cos(a), sin(a), -sin(a), cos(a)), ncol=2)
-  S = diag(c(1, sqrt(ratio)))
-  A = t(R) %*% S
+  
+  # matrix of coordinates of the new basis (in columns)=
+  #    =matrix of transformation from new to old coordinate system
+  R = matrix(c(sin(a), cos(a), cos(a), -sin(a)), ncol=2)
+  # matrix of transformation from old to new coordinate system
+  # Rinv = t(R)
+  # scaling:
+  S = diag(c(1, ratio^((1-2*inv)/2) ))
+  # rotation to new coordinates and scaling:
+  A = R %*% S # t(Rinv) %*% S # because t(Rinv) = R in orthogonal matrices
   A = rbind(cbind(A,0),c(0,0,1))
   class(A) = c("AnisotropyScaling","Anisotropy")
-    return(A)
-}  ### NOT SURE OF IT
-#  anis2D.par2A(1,0)
-#  anis2D.par2A(0.5,0)
-#  A = anis2D.par2A(0.5, 30)
-#    h = rbind(c(0,0), c(2,1), c(1,-2))
-#      anish2Dist(h, A = A)
-#aa = seq(from=0, to=2*pi, length.out=360)
-#A = anis2D.par2A(ratio=1/2, angle=45)
-#x = cbind(sin(aa), cos(aa)) %*% A
-#plot(x, type="l", asp=1)
-
+  return(A)
+}
 
 
+#' @describeIn anis_GSLIBpar2A 3D case
+#' @export
+#' @examples 
+#' c60 = cos(60*pi/180)
+#' s60 = sin(60*pi/180)
+#' c30 = cos(30*pi/180)
+#' s30 = sin(30*pi/180)
+#' #  in the new coordinates, 60cwN is (0,1,0)
+#' R60p = anis3D_par2A(ratios=c(1,1), angles=c(60,0,0))
+#' c(s60, c60, 0) %*% R60p 
+#' R6030 = anis3D_par2A(ratios=c(1,1), angles=c(60,30,0))
+#' # the original X axis is positive on newX and newY, but negative on newZ
+#' c(1,0,0) %*% R6030 
+#' # rotate first direction 60 degrees azimuth, then dip 30degrees upwards
+#' c( c60*c30, -s60*c30, s30) %*% R6030 
+#' (Ranis = anis3D_par2A(ratios=c(0.5,0.25), angles=c(60,30,0)) )
+anis3D_par2A = function(ratios, angles, inv=FALSE){
+  # A should be a matrix such that for h=[hx hy, hz] lag vectpr
+  #   hn = norm( h %*% A[1:ncol(h),]) allows to use the 
+  #   isotropic correlogram
+  ## transposed Cholesky decomp of the rotation matrix from angle of maximum
+  #     spatial continuity + ratio min/max ranges
+  a = angles * pi/180
+  # matrix of coordinates of the new basis (in columns)=
+  #    =matrix of transformation from new to old coordinate system
+  r = function(a, i=3, cw=TRUE){
+    if(i==2) cw=!cw
+    m = matrix( c(cos(a), (-1)^cw*sin(a), 0, -(-1)^cw*sin(a), cos(a), 0, 0,0,1), ncol=3)
+    ii = 1:3
+    jj = order(c(ii[-i],ii[i]))
+    return(m[jj,jj])
+  }
+  R = r(a[1],i=3,cw=T) %*% r(a[2],i=2,cw=T) %*% r(a[3],i=1,cw=F)
+  # matrix of transformation from old to new coordinate system
+  # Rinv = t(R)
+  # scaling:
+  S = diag(c(1, ratios^((1-2*inv)/2) ))
+  # rotation to new coordinates and scaling:
+  A = R %*% S # t(Rinv) %*% S # because t(Rinv) = R in orthogonal matrices
+  class(A) = c("AnisotropyScaling","Anisotropy")
+  return(A)
+}
diff --git a/R/gstatCompatibility.R b/R/gstatCompatibility.R
index 8f141edee7f75ddf23824464a1195cbc429626a9..725e19a15046e466dffb398938c43b04b2dba2d5 100644
--- a/R/gstatCompatibility.R
+++ b/R/gstatCompatibility.R
@@ -573,7 +573,6 @@ as.variogramModel.LMCAnisCompo <- function(m, V=NULL, prefix=NULL, ensurePSD=TRU
   eqtabmodels = factor(c(nugget="Nug", exp="Exp", sph="Sph", gau="Gau"), levels=levels(vgm()[,1]))
   models = eqtabmodels[sapply(1:ncol(m), function(j) m[,j]$model)]
   # express all sill matrices in the desired logratio 
-  # recode A in azimuth, range and range ratio
   for(j in 1:ncol(m)){
     aux = -0.5 * t(V) %*% m[,j]$sill  %*% V
     colnames(aux) <- rownames(aux) <- noms
@@ -586,6 +585,8 @@ as.variogramModel.LMCAnisCompo <- function(m, V=NULL, prefix=NULL, ensurePSD=TRU
     }
     m[,j]$sill = aux
   }
+  # recode A in azimuth, range and range ratio
+  # TODO: correct and extend to 3D (check consistency with `anis2D.par2A()` and `anish2Dist()`)
   anis = matrix(0, nrow=ncol(m), ncol=3)
   colnames(anis) = c("range", "ang1", "anis1")
   for(j in 1:ncol(m)){
@@ -691,7 +692,9 @@ as.LMCAnisCompo.variogramModelList <-
       model = eqtabmodels[ m[[1]]$model[i] ]
       sill = clrvar2variation(t(W) %*% darstellung(m, i)  %*% W)
       range = m[[1]]$range[i] * m[[1]]$anis1[i]
-      A = anis2D.par2A(ratio=m[[1]]$anis1[i], angle=m[[1]]$ang1[i])
+      # A = anis2D_par2A(ratio=m[[1]]$anis1[i], angle=m[[1]]$ang1[i], inv=FALSE)
+      A = anis_GSLIBpar2A(ratios=c(m[[1]]$anis1[i],m[[1]]$anis2[i]), 
+                          angles=c(m[[1]]$ang1[i],m[[1]]$ang2[i],m[[1]]$ang3[i]) )
       rs = list(model=model, range=range, A=A, sill=sill)
       class(rs) = "variostructure"
       return(rs)
diff --git a/R/variograms.R b/R/variograms.R
index 256bf6ff54f875deaeb7eb8107c54ac09c5e75dd..01502ace7235d15f473bb3dad2037ca77ec88155 100644
--- a/R/variograms.R
+++ b/R/variograms.R
@@ -50,7 +50,7 @@ setCgram = function(type, nugget=sill*0, sill, anisRanges, extraPar=0){
                 data=extraPar,
                 nugget=nugget,
                 sill=sill, #(nstru, nvar, nvar)
-                M=anisRanges     # these are "classical" ranges (i.e. distances)
+                M=AnisotropyRangeMatrix(anisRanges)     # these are "classical" ranges (i.e. distances)
   )
   class(vgout) = "gmCgram"
   return(vgout)
@@ -453,6 +453,7 @@ plot.gmCgram = function(x, xlim.up=NULL, xlim.lo=NULL, vdir.up= NULL, vdir.lo= N
 
 
 
+
 #' Check for anisotropy of a theoretical variogram
 #'
 #' Checks for anisotropy of a theoretical variogram or covariance function model
@@ -473,10 +474,10 @@ is.isotropic.default = function(x, tol=1e-10, ...) NA
 #' @method is.isotropic gmCgram
 #' @export
 is.isotropic.gmCgram = function(x, tol=1e-10, ...){
-   all(apply(x$M, 1, function(y){
-     ev = eigen(y, only.values=TRUE)[[1]]
-     all(abs(ev-ev[1])<tol)
-   }))
+  all(apply(x$M, 1, function(y){
+    ev = eigen(y, only.values=TRUE)[[1]]
+    all(abs(ev-ev[1])<tol)
+  }))
 }
 
 #' @method is.isotropic variogramModel
@@ -503,6 +504,25 @@ is.isotropic.LMCAnisCompo = function(x, tol=1e-10, ...){
 
 
 
+#' Check for any anisotropy class
+#' 
+#' Check that an object contains a valid specification of anisotropy, in any form 
+#'
+#' @param x object to check
+#'
+#' @return a logical, TRUE if the object is an anisotropy specification; FALSE otherwise 
+#' @export
+#'
+#' @examples
+#' a =  anis2D_par2A(0.5, 30)
+#' a
+#' is.anisotropySpecification(a)
+is.anisotropySpecification = function(x){
+  "Anisotropy" %in% class(x)  |  inherits(x, "Anisotropy")
+} 
+
+
+
 
 
 
@@ -662,7 +682,7 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)),
   azs = azimuths * pi/180
   res = sapply(1:Na, function(i){
     tk_a = (azs[i,1]<=XXaz) & (azs[i,2]>=XXaz)
-    zz = ZZ[tk_a,,]
+    zz = ZZ[tk_a,]
     xxabs = XXabs[tk_a]
     xxaz = XXaz[tk_a]
     tk_h = outer(xxabs, lags[,1],">=") & outer(xxabs, lags[,2],"<=")
@@ -674,7 +694,7 @@ gsi.EVario2D = function(X,Z,Ff=rep(1, nrow(X)),
     for(j in 1:Nh){
       if(!is.na(n[j,i]) && n[j,i]>minpairs){
         if(cov){ # covariance function
-          vg[j,,,i] = gmApply((ZZ[tk_a,,][tk_h[,j],,]), c(2,3),"sum")/(n[j,i])
+          vg[j,,,i] = gmApply((ZZ[tk_a,][tk_h[,j],,]), c(2,3),"sum")/(n[j,i])
         }else{   # semi-variogram
           aux = rowSums(
             gmApply(X=ZZ[tk_a,][tk_h[,j],], 1, function(x)outer(x,x,"*"))
@@ -834,11 +854,11 @@ gsi.EVario3D = function(X,Z,Ff=rep(1, nrow(X)),
         vg[j,,,i]=NA
       }
     }
-    return(list(gamma=vg[,,,i], lags=gmGeostats:::gsi.lagClass(lags), npairs =n[,i]))
+    return(list(gamma=vg[,,,i], lags=gsi.lagClass(lags), npairs =n[,i]))
   })
 
   # output
-  attr(res, "directions") = gmGeostats:::gsi.directorVector(dirvecs[,1:3])
+  attr(res, "directions") = gsi.directorVector(dirvecs[,1:3])
   # attr(res, "lags") = gsi.lagClass(lags)
   attr(res, "type") = ifelse(cov, "covariance","semivariogram")
   class(res) = "gmEVario"
diff --git a/vignettes/SimulatingMicrostructures.Rmd b/vignettes/SimulatingMicrostructures.Rmd
index 508917c13d18fc5726a162aa80ab6e6731342efa..a25d0dd9af691b727aa2da7f9e29628fd71ac82a 100644
--- a/vignettes/SimulatingMicrostructures.Rmd
+++ b/vignettes/SimulatingMicrostructures.Rmd
@@ -83,7 +83,7 @@ s1
 
 ### Characteristics of phases
 
-The anisotropy of the elongation for each phase can be obtained by the function `anis2D.par2A`.
+The anisotropy of the elongation for each phase can be obtained by the function `anis2D_._par2A`.
 The function parameter `ratio` provides the ratio of shortest to longest range, and should therefore be smaller than 1.
 E.g. `ratio = 1` is a sphere (isotrop), while `ratio = 0.5` would describe a mineral which grains are typicall double as long as thick.
 The function parameter `angle` gives the orientation of this anisotropy, measured in degrees counterclockwise from the east (E).
@@ -91,15 +91,15 @@ E.g. `angle = 45` with `ratio = 0.5` would result into NE-SW elongated crystals.
 
 For Quartz, let's assume we have mostly roundish/quadratic grains and only a minor preferred orientation of the crystals, so that they are slightly oriented along NE-SW:
 ```{r p-qu}
-p1_qu = anis2D.par2A(ratio = .9, angle = 45)
+p1_qu = anis2D_par2A(ratio = .9, angle = 45)
 ```
 
 The Feldspars have a more pronounced elongated shape, also with preferred orientation along the same direction as the Quartz grains:
 ```{r p-fsp}
-p2_fsp = anis2D.par2A(ratio = .6, angle = 45)
+p2_fsp = anis2D_par2A(ratio = .6, angle = 45)
 ```
 
-Function `anis2D.par2A` returns the Cholesky decomposition of the anisotropy matrix. So, if the anistropy of the mineral phases had been generated by this function, then we need to calculate the cross product of the matrices to obtain the anistropy matrices `A`:
+Function `anis2D_par2A` returns the Cholesky decomposition of the anisotropy matrix. So, if the anistropy of the mineral phases had been generated by this function, then we need to calculate the cross product of the matrices to obtain the anistropy matrices `A`:
 ```{r}
 A1_qu = tcrossprod(p1_qu)   
 A2_fsp = tcrossprod(p2_fsp)
@@ -287,7 +287,7 @@ For example, in the simulation above, the geologist would expect, that the fract
 
 Changing the orientation of the fractures:
 ```{r}
-p4_frac = anis2D.par2A(ratio = .01, angle = 95)
+p4_frac = anis2D_par2A(ratio = .01, angle = 95)
 A4_frac = tcrossprod(p4_frac)
 A4_frac
 ```