From fa191e167322c2c3b0720d07fd7fef3f70c660f0 Mon Sep 17 00:00:00 2001
From: Raimon Tolosana-Delgado <tolosa53@fwg206.ad.fz-rossendorf.de>
Date: Fri, 17 Dec 2021 13:33:55 +0100
Subject: [PATCH] corrections in Turning Bands (still internal) and in
anisotropy modelling
---
DESCRIPTION | 2 +-
NEWS.md | 1 +
R/gmAnisotropy.R | 2 +-
src/gmGeostats.c | 11 +++++------
vignettes/SimulatingMicrostructures.Rmd | 4 ++--
5 files changed, 10 insertions(+), 10 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 71c3655..09f062a 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 0b87a47..4e94965 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 98202c7..3b12f34 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 046f544..7a90c72 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 cadb809..2702fc6 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)
--
GitLab