corrections in Turning Bands (still internal) and in anisotropy modelling

DESCRIPTION

 Package: gmGeostats
 Version: 0.11.0-9002
-Date: 2021-12-08
+Date: 2021-12-14
 Title: Geostatistics for Compositional Analysis
 Authors@R: c(person(given = "Raimon", 
         family = "Tolosana-Delgado", 

NEWS.md

 # gmGeostats 0.11.0-9002
 
+* (2021-12-14) bugs in turning bands corrected: getUnitVec was producing a wrong sequence of directions, and bands for exponential and Gaussian variograms did not correct for the difference between parametric and effective range in the right way; 
 * (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()`
  

R/gmAnisotropy.R

@@ -181,7 +181,7 @@ as.AnisotropyRangeMatrix.AnisotropyScaling <- function(x) AnisotropyRangeMatrix(
   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){
+  }else if( (length(ratios)==2) & length(angles)==3){
     anis3D_par2A(ratios=ratios, angles=angles, inv=inv)
   }else{
     stop("anis.GSLIBpar2A: error, ratios and angles length incompatible for both 2D and 3D")

src/gmGeostats.c

@@ -212,7 +212,7 @@ const double *extra /* extra parameter (unused), for consistency with other cova
    * projs points. Return that wave */	
   int i;
   double phase,amp,omega,d1;   
-  omega = norm_rand()  * M_SQRT2 / range; 
+  omega = norm_rand()  * M_SQRT2 / range * M_SQRT_3; /* sqrt(3) needed to produce effective range*/
   phase = unif_rand() * M_2PI;
   amp  = M_SQRT2; /* Lantuejoul (2002), page 191 */
   for(i=0;i<n;i++) {
@@ -283,13 +283,12 @@ void fbandExp(int n, /* number of locations */
     const double *extra /* extra parameter (unused), for consistency with other covariances */
 ){
   int i,j,ns;
-  double d1,x0,x1,sign,effrange;
+  double d1,x0,x1,sign;
   /* Lantuejoul (2002), page 196 */ 
   /* Find the smallest proj; select a random exponential point "x0" left 
    * from it with lambda "2range" away. Domain = (x0, max(projs)) */
   sign = unif_rand()>0.5? 1 : -1; /* start + or - randomly */
-  effrange = range/3.0;  /*  ATTENTION: 3*range, but range is inverted ??? */
-  x0 = projs[0];  /* does min exist?? */
+  x0 = projs[0];  
   x1 = projs[0];
   for(i=1;i<n;i++) {
     if( projs[i]>x1 )
@@ -297,7 +296,7 @@ void fbandExp(int n, /* number of locations */
     else if( projs[i]<x0 )
       x0=projs[i];
   }  
-  x0 -= 2*effrange*exp_rand();
+  x0 -= 2*range*exp_rand();
   
   /* Partition the domain with a Poisson point process of lambda=2range. */
   ns = 0;
@@ -305,7 +304,7 @@ void fbandExp(int n, /* number of locations */
   while( doubleBuf[ns]<x1 ){
     if( ns>=maxIntervals )
       error("fbandExp: too small range; merge with nugget?");
-    doubleBuf[ns+1] = doubleBuf[ns] + 2*effrange*exp_rand();
+    doubleBuf[ns+1] = doubleBuf[ns] + 2*range*exp_rand();
     ns ++;
   }
   /* Assign values*/

vignettes/SimulatingMicrostructures.Rmd

@@ -161,8 +161,8 @@ We produced these matrices in the first part, where the number of phases was def
 This anisotropy range has to be downscaled, therefore a factor is applied prior to the the anisotropy matrices:
 
 ```{r}
-v1_qu = setCgram(type = vg.Gau, sill = s1, anisRanges = 0.005*A1_qu)
-v2_fsp = setCgram(type = vg.Gau, sill = s2, anisRanges = 0.005*A2_fsp)
+v1_qu = setCgram(type = vg.Gau, sill = s1, anisRanges = 0.0005*A1_qu)
+v2_fsp = setCgram(type = vg.Gau, sill = s2, anisRanges = 0.0005*A2_fsp)
 v3_bi = setCgram(type = vg.Exp, sill = s3, anisRanges = 0.005*A3_bi)
 v4_frac = setCgram(type = vg.Exp, sill = s4, anisRanges = 0.05*A4_frac)
 v5_joints = setCgram(type = vg.Exp, sill = s5, anisRanges = 0.05*A5_joints)