diff --git a/DESCRIPTION b/DESCRIPTION index 71c365593e0eba075c53944dd4075dc27997e1c4..09f062af514461cdf5d6cd0a6607481f82bba540 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ 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", diff --git a/NEWS.md b/NEWS.md index 0b87a47263b4a8f7efbdec1cec8368627902fafc..4e94965c3a0ecfed07b82d28af8eed057db82949 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,6 @@ # 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()` diff --git a/R/gmAnisotropy.R b/R/gmAnisotropy.R index 98202c706416ffeb65a87be612a9382423e5a09d..3b12f34e5d414b72a6ca5f809e101abe114b1195 100644 --- a/R/gmAnisotropy.R +++ b/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") diff --git a/src/gmGeostats.c b/src/gmGeostats.c index 046f54492c6606bda985b2fcc50f90efe4d616bc..7a90c723f8f1cc082400dd8980d181aec8873c0d 100644 --- a/src/gmGeostats.c +++ b/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*/ diff --git a/vignettes/SimulatingMicrostructures.Rmd b/vignettes/SimulatingMicrostructures.Rmd index cadb809954f2b9b0467f4f946506791e88a29285..2702fc6e4ed6d9b601dca33e2da6667030369070 100644 --- a/vignettes/SimulatingMicrostructures.Rmd +++ b/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)