turning bands bug in getUnitVec corrected

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
+    }
 
 
 

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},