diff --git a/infer.cc b/infer.cc
index cff380ce262b74619db70e0f053a935be2140f16..bf4ad2410873e7bbfcf0b27c44c9df48c09c8c5f 100644
--- a/infer.cc
+++ b/infer.cc
@@ -28,23 +28,28 @@ Infer::Infer(float udLh, float ut, unsigned int uV) {
   t = ut;
   V = uV;
   if (V<=0.0) V=1.0;
+  // dimensionality
+  D = 3.0;
+  // if (zsize <= 1)  D = 2.0;  // to be done by calling program
 
+  // NB: the (sqr(t) -1) is previous D=3 version (from where??)
   if (fabs(t)<13.0) {
-    Em_ = V * pow(2*M_PI,-2) * dLh * (sqr(t) - 1) * exp(-sqr(t)/2.0);
+    Em_ = V * pow(2*M_PI,-(D+1)/2) * dLh * pow((sqr(t) - 1), (D-1)/2) *
+      exp(-sqr(t)/2.0); 
   } else {
     Em_ = 0.0;  // underflowed exp()
   }
 
   if (fabs(t)<8.0) {
-    B_ = pow((gamma(2.5)*Em_)/(V*(0.5 + 0.5*erf(-t/sqrt(2)))),(2.0/3.0));
+    B_ = pow((gamma(1.0+D/2.0)*Em_)/(V*(0.5 + 0.5*erf(-t/sqrt(2)))),(2.0/D));
   } else {
-    // the large t approximation
-    float a1 = V * dLh / (4 * M_PI * M_PI);
-    float a3 = pow((gamma(2.5)  / V ),(2.0/3.0));
+    // the large t approximation  (see appendix below)
+    float a1 = V * dLh * pow(2*M_PI,-(D+1)/2);
+    float a3 = pow((gamma(1+D/2.0)  / V ),(2.0/D));
     float tsq = t*t;
-    float Em_q = a1 * (2.0 * M_PI) * (tsq - 1.0)*t 
-                     / ( 1.0 - 1.0/tsq + 3.0/(tsq*tsq)) ;
-    B_ = a3 * pow(Em_q,(2.0/3.0));
+    float c = pow(2*M_PI,-1.0/2.0) * t / ( 1.0 - 1.0/tsq + 3.0/(tsq*tsq)) ;
+    float Em_q = a1 * pow((tsq - 1.0),(D-1)/2) * c;
+    B_ = a3 * pow(Em_q,(2.0/D));
   }
   
 
@@ -54,10 +59,10 @@ Infer::Infer(float udLh, float ut, unsigned int uV) {
   
 float Infer::operator() (unsigned int k) {
   // ideally returns the following:
-  //    return 1 - exp(-Em_ * exp(-B_ * pow( k , 2.0 / 3.0)));
+  //    return 1 - exp(-Em_ * exp(-B_ * pow( k , 2.0 / D)));
   // but in practice must be careful about ranges
   // Assumes that exp(+/-87) => 1e+/38 is OK for floats
-  float arg1 = -B_ * pow(k , 2.0 / 3.0);
+  float arg1 = -B_ * pow(k , 2.0 / D);
   if (fabs(arg1)<87.0) {
     float exp1 = exp(arg1);
     float arg2 = -Em_ * exp1;
@@ -66,16 +71,49 @@ float Infer::operator() (unsigned int k) {
       return 1.0 - exp2;
     }
   }
-  // now must invoke the approximation instead
-    float a1 = V * dLh / (4 * M_PI * M_PI);
-    float a3 = pow((gamma(2.5)  / V ),(2.0/3.0));
+  // now must invoke the approximation instead  (see appendix below)
+    float a1 = V * dLh * pow(2*M_PI,-(D+1)/2);
+    float a3 = pow((gamma(1.0+D/2.0)  / V ),(2.0/D));
     float tsq = t*t;
-    float c = 1.0 / ( 1.0 - 1.0/tsq + 3.0/(tsq*tsq));
-    float arg = -tsq*(0.5 + a3*pow(a1,2.0/3.0)*pow(2.0*M_PI,1.0/3.0)*c*pow(k,2.0/3.0));
+    float c =  pow(2*M_PI,-1.0/2.0) * t / ( 1.0 - 1.0/tsq + 3.0/(tsq*tsq));
+    float Em_q = a1 * pow((tsq - 1.0),(D-1)/2) * c;
+    float beta = a3 * pow(Em_q , 2.0/D);
+    float arg = -tsq*0.5 - beta*pow(k,2.0/D));
     float exp1 = 0.0;
     if (arg>-87.0) {
       exp1 = exp(arg);
     }
-    float p = a1*(tsq-1.0)*exp1;
+    float p = a1*pow((tsq-1.0),(D-1)/2)*exp1;
     return p;
 }
+
+
+
+// MATHEMATICAL APPENDIX
+
+/*
+
+The formulas that need to be calculated are:
+(1) E_m = V * dLh * (2*pi)^(-(D+1)/2) * (t^2 -1)^((D-1)/2) * exp(-t^2 /2)
+(2) Beta = (Gamma(D/2+1)/V * E_m / Phi(-t) )^(2/D)
+(3) p = 1 - exp( - E_m * exp(-Beta*k^(2/D)))
+
+where Phi(-t) = Gaussian cumulant = (1/2 + 1/2*erf(-t/sqrt(2)))
+
+These are approximated by:
+
+(2a) Beta = (Gamma(D/2+1)/V)^(2/D) * (Em1)^(2/D) * Ct^(2/D)
+where Em1 = V * dLh * (2*pi)^(-(D+1)/2) * (t^2 -1)^((D-1)/2)
+      Ct = (2*pi)^(-1/2) * t / ( 1.0 - 1.0/t^2 + 3.0/t^4 )
+      which approximates ( exp(-t^2 /2) / Phi(-t) )^(2/D)
+      using 1/2 - 1/2*erf(t/sqrt(2)) = (2*pi)^(1/2) * exp(-t^2 /2) * ...
+                                        (1-1/t^2+3/t^4) / t
+					(this is derived in TR00MJ1)
+
+and
+
+(3a) p = Em1 * exp( -t^2 /2 - beta * k^(2/D))
+using the approximation (2a) for beta
+
+ */
+