From a58ce452fca67a09c33b1153b6fce01c83f68826 Mon Sep 17 00:00:00 2001
From: Raimon Tolosana-Delgado <tolosa53@fwg206.ad.fz-rossendorf.de>
Date: Tue, 14 Dec 2021 11:16:20 +0100
Subject: [PATCH] turning bands bug in getUnitVec corrected
---
R/geostats.R | 6 +++++-
src/gmGeostats.c | 22 +++++++++++++++++++---
2 files changed, 24 insertions(+), 4 deletions(-)
diff --git a/R/geostats.R b/R/geostats.R
index 46bed0e..481d333 100644
--- a/R/geostats.R
+++ b/R/geostats.R
@@ -70,7 +70,11 @@ gsi.calcCgram <- function(X,Y,vgram,ijEqual=FALSE) {
structure(erg$C,dim=c(d*nX,d*nY))
}
-
+gsiGetUnitVec <- function(dimX,ip) {
+ erg = .C("getUnitvecR",dimX=checkInt(dimX,1),ip=checkInt(ip,1),unitvec=numeric(dimX))
+ unitvec = erg$unitvec
+ unitvec
+ }
diff --git a/src/gmGeostats.c b/src/gmGeostats.c
index c2c4cf3..046f544 100644
--- a/src/gmGeostats.c
+++ b/src/gmGeostats.c
@@ -52,7 +52,7 @@ double invBitExp2(int i) {
}
-/* invBitExp2
+/* invBitExp
inverts the sequence of the digits in b-adict representation.
this is used for the generation of almost equally space directions in 3D
@@ -288,7 +288,7 @@ void fbandExp(int n, /* number of locations */
/* 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; /* ATTENTION: 3*range, but range is inverted ??? */
+ effrange = range/3.0; /* ATTENTION: 3*range, but range is inverted ??? */
x0 = projs[0]; /* does min exist?? */
x1 = projs[0];
for(i=1;i<n;i++) {
@@ -334,7 +334,7 @@ void getUnitvec(
if(dimX>3)
error("no expression for unit vectors in dimension larger than 3");
if( dimX==3) {
- d1=invBitExp2(ip)*M_2_PI;
+ d1=invBitExp2(ip)*2*M_PI;
d2=invBitExp(ip,3);
d3 = sqrt(1-d2*d2);
unitvec[2] = d2;
@@ -352,6 +352,20 @@ void getUnitvec(
+static R_NativePrimitiveArgType getUnitvecR_t[] = { /* INTSXP,REALSXP */
+ /* dimX, iP, unitvec */
+ INTSXP,INTSXP,REALSXP
+};
+
+
+void getUnitvecR(int * dimX,
+ int * ip,
+ double *unitvec
+ ) {
+ getUnitvec(*dimX,*ip,unitvec);
+}
+
+
static R_NativePrimitiveArgType CMVTurningBands_t[] = { /* INTSXP,REALSXP */
/* dimX, X, dimZ Z nBands sqrtNug nCgram typeCgr A sqrtSill moreCgr */
INTSXP,REALSXP,INTSXP,REALSXP,INTSXP,REALSXP,INTSXP,INTSXP,REALSXP,REALSXP,REALSXP
@@ -872,9 +886,11 @@ static R_CMethodDef cMethods[] = {
{"CCondSim", (DL_FUNC) & CCondSim, 18, CCondSim_t},
{"anaForwardC", (DL_FUNC) &anaForwardC, 8, anaForwardC_t},
{"anaBackwardC", (DL_FUNC) &anaBackwardC, 8, anaBackwardC_t},
+ {"getUnitvecR", (DL_FUNC) &getUnitvecR, 3, getUnitvecR_t},
{NULL, NULL, 0}
};
+
//{"CMVTurningBands2", (DL_FUNC) &CMVTurningBands2, 11, CMVTurningBands2_t},
--
GitLab