diff --git a/pub/CHANGELOG b/pub/CHANGELOG
index 1b3b1e0a6ca1f7c9f534aa95122baa2478f17dc7..5adb31512d59f3c7904cbaf967be7581c6227a7f 100644
--- a/pub/CHANGELOG
+++ b/pub/CHANGELOG
@@ -1,10 +1,12 @@
-Release
+Release 2.3.0d of 1jul15:
 
 - New functionality:
   - 'mpe-' to remove error bars
   - on upstart, test whether Gnuplot supports X11
 - Bugfixes:
   - 'msj' can detect incongruent error bar status
+- Improved behavior:
+  - modified parametrization of ornuhl
 
 Release 2.3.0c of 28may15:
 
diff --git a/pub/lib/fbase.cpp b/pub/lib/fbase.cpp
index 7d826ee8ed6621d969b01cf19e6b53f297428e09..10fba78f011deacffcec8bf0aabc2c984397f9d2 100644
--- a/pub/lib/fbase.cpp
+++ b/pub/lib/fbase.cpp
@@ -353,15 +353,17 @@ double func_rotdiff( double w, double tau, double qb )
 }
 
 double func_ornuhl( double _w, double _a, double _b )
-// Ornstein-Uhlenbeck spectrum (cf kww branch parfour1)
+// Ornstein-Uhlenbeck spectrum FT[exp(-|b|t-a(1-exp(-t))),
+//    computed for w->0 by expanding the outer exponential in e^(a e^(-x/w))
+//    which turns out to converge rapidly for all w
 {
     const double ornuhl_delta=2.2e-16, ornuhl_eps=5.5e-20;
     long double w = _w;
     long double a = _a;
-    long double b = _b;
+    long double b = fabs(_b);
     if ( a==0 )
-        return 1/(1+SQR(w));
-    double alog = log(fabs(a));
+        return b/(SQR(b)+SQR(w));
+    long double a_log = logl( a );
     long double S=0;      // summed series
     long double h;        // exponent
     long double f;        // Cauchy-Lorentz function
@@ -369,15 +371,18 @@ double func_ornuhl( double _w, double _a, double _b )
     long double u_next=0; // - next value [initialized to avoid warning]
 
     // sum the expansion
-    for ( int i=0; i<1000; ++i ) {
+    for ( int k=0; k<1000; ++k ) {
         u = u_next;
-        h = i*alog - lgammal((long double)(i+1));
+        h = k*a_log - lgammal((long double)(k+1));
         if ( h>DBL_MAX_EXP )
             return -3; // exponent overflow
-        f = (1+i*b)/(SQR(1+i*b)+SQR(w));
+        if ( h<DBL_MIN_EXP )
+            return -4; // exponent underflow
+        f = (k+b)/(SQR(k+b)+SQR(w));
         u_next = expl( h ) * f;
         S += u;
-        //// ornuhl_num_of_terms = i;
+        if( k==0 )
+            continue;
         // termination criteria
         if ( ornuhl_eps*S+u_next <= ornuhl_delta*S )
             return exp(-a)*S; // reached required precision
@@ -632,7 +637,7 @@ void fbase_initialize()
     G->register_fct_d_ddd( "zorn2", func_zorn_gauss );
     G->register_fct_meta ( "rrdm", 3, "(w*t0,EA_mean/T,EA_stdv/EA_mean: rotational rate distribution model" );
     G->register_fct_d_ddd( "rrdm", func_rrdm );
-    G->register_fct_meta ( "ornuhl", 3, "Ornstein-Uhlenbeck spectrum FT[exp(-t-a(1-exp(-bt)))]" );
+    G->register_fct_meta ( "ornuhl", 3, "(w,a,b) Ornstein-Uhlenbeck spectrum FT[exp(-|b|t-a(1-exp(-t)))]" );
     G->register_fct_d_ddd( "ornuhl", func_ornuhl );
     G->register_fct_meta ( "rotdiff", 3, "(w,tau,qb: rotational diffusion spectrum)" );
     G->register_fct_d_ddd( "rotdiff", func_rotdiff );