diff -Naur gnuplot-3.7/fnproto.h gnuplot-3.7_p/fnproto.h
--- gnuplot-3.7/fnproto.h	Wed Apr 15 21:22:24 1998
+++ gnuplot-3.7_p/fnproto.h	Wed Oct 27 13:52:28 1999
@@ -129,6 +129,10 @@
 void f_inverse_normal __PROTO((void));
 void f_inverse_erf __PROTO((void));
 
+/* new function in "specfun.c" 
+ * Lamberts W-function */
+void f_omega __PROTO((void));
+
 /* prototypes from file "datafile.c" */
 
 void f_dollars __PROTO((union argument *x));
diff -Naur gnuplot-3.7/plot.c gnuplot-3.7_p/plot.c
--- gnuplot-3.7/plot.c	Wed Dec  9 16:24:14 1998
+++ gnuplot-3.7_p/plot.c	Wed Oct 27 13:52:21 1999
@@ -208,6 +208,8 @@
     {"tm_wday", (FUNC_PTR) f_tmwday},	/* for timeseries */
     {"tm_yday", (FUNC_PTR) f_tmyday},	/* for timeseries */
 
+    {"omega", (FUNC_PTR) f_omega},      /* self added omega */
+
     {NULL, NULL}
 };
 
diff -Naur gnuplot-3.7/specfun.c gnuplot-3.7_p/specfun.c
--- gnuplot-3.7/specfun.c	Thu Dec 10 19:30:35 1998
+++ gnuplot-3.7_p/specfun.c	Wed Oct 27 13:52:17 1999
@@ -113,6 +113,8 @@
 static double ranf __PROTO((double init));
 static double inverse_normal_func __PROTO((double p));
 static double inverse_error_func __PROTO((double p));
+/* added by defining omega */
+static double omega __PROTO((double x));
 
 #ifndef GAMMA
 /* Provide GAMMA function for those who do not already have one */
@@ -867,4 +869,69 @@
     }
     y = sigma * f;
     return (y);
+}
+
+/* Implementation of Lamberts W-function which is defined as
+ * w(x)*e^(w(x))=x
+ * Implementation by Gunter Kuhnle, gk@uni-leipzig.de 
+ * Algorith originally developed by
+ * KEITH BRIGGS, DEPARTMENT OF PLANT SCIENCES,
+ * e-mail:kmb28@cam.ac.uk
+ * http://epidem13.plantsci.cam.ac.uk/~kbriggs/W-ology.html */
+
+
+
+static double omega(x)
+double x;
+{
+    double p,e,t,w,eps;
+    int i;
+
+    
+	
+    eps = MACHEPS;
+
+    if (x>=-exp(-1))
+    {
+        if (x==0) 
+	{
+	    return 0.0;
+	} else
+	{
+	    if (x < 1)
+	    {
+	       	p=sqrt(2.0*(exp(1.0)*x+1.0)); 
+	      	w=-1.0+p-p*p/3.0+11.0/72.0*p*p*p;
+	    } else
+	    {
+	        w = log(x);
+	    }
+
+	    if (x > 3)
+	    {
+	 	w = w-log(w);
+	    }
+
+	    for (i=0;i<20;i++)
+	    {
+		e=exp(w);
+		t=w*e-x;
+		t=t/(e*(w+1.0)-0.5*(w+2.0)*t/(w+1.0));
+		w=w-t;
+		if (fabs(t)<eps*(1.0+fabs(w))) return w;
+	    } 
+	    return -1;
+       } 
+       return -1;
+    }
+}
+
+void f_omega()
+{
+    struct value a;
+    double x;
+
+    x=real(pop(&a));
+
+    push(Gcomplex(&a,omega(x),0.0));
 }