diff --git a/Mesh/Interpolation.cpp b/Mesh/Interpolation.cpp
index 30583f26f34e49c1afa59ed95b89f25d0c4dd905..ec1b13c04dbfb76df39b21bc5f70e99bb486efcc 100644
--- a/Mesh/Interpolation.cpp
+++ b/Mesh/Interpolation.cpp
@@ -1,4 +1,4 @@
-// $Id: Interpolation.cpp,v 1.11 2001-08-12 20:45:02 geuzaine Exp $
+// $Id: Interpolation.cpp,v 1.12 2001-08-17 07:41:58 geuzaine Exp $
 
 #include "Gmsh.h"
 #include "Numeric.h"
@@ -50,7 +50,7 @@ Vertex InterpolateCurve (Curve * Curve, double u, int derivee){
     List_Read (Curve->Control_Points, i, &v[1]);
     List_Read (Curve->Control_Points, i + 1, &v[2]);
     
-    V.lc = t * v[2]->lc + (1. - t) * v[1]->lc;
+    V.lc = t * v[2]->lc + (1. - t) * v[1]->lc; // ?????
     V.Pos.X = v[1]->Pos.X + t * (v[2]->Pos.X - v[1]->Pos.X);
     V.Pos.Y = v[1]->Pos.Y + t * (v[2]->Pos.Y - v[1]->Pos.Y);
     V.Pos.Z = v[1]->Pos.Z + t * (v[2]->Pos.Z - v[1]->Pos.Z);
@@ -61,7 +61,8 @@ Vertex InterpolateCurve (Curve * Curve, double u, int derivee){
     V.Pos.X = evaluate_scalarfunction ("t", u, Curve->functu);
     V.Pos.Y = evaluate_scalarfunction ("t", u, Curve->functv);
     V.Pos.Z = evaluate_scalarfunction ("t", u, Curve->functw);
-    V.lc = (u * Curve->beg->lc + (1. - u) * Curve->end->lc);
+
+    V.lc = (u * Curve->end->lc + (1. - u) * Curve->beg->lc);// ?????
     V.w = (u * Curve->beg->w + (1. - u) * Curve->end->w);
     return V;
     
@@ -90,22 +91,17 @@ Vertex InterpolateCurve (Curve * Curve, double u, int derivee){
     V.Pos.Y += Curve->Circle.v[2]->Pos.Y;
     V.Pos.Z += Curve->Circle.v[2]->Pos.Z;
     V.w = (u * Curve->beg->w + (1. - u) * Curve->end->w);
-
-    // ?????
-    V.lc = (u * Curve->end->lc + (1. - u) * Curve->beg->lc);
+    V.lc = (u * Curve->end->lc + (1. - u) * Curve->beg->lc);// ?????
     return V;
     
   case MSH_SEGM_BSPLN:
   case MSH_SEGM_BEZIER:
-    V.lc = (u * Curve->beg->lc + (1. - u) * Curve->end->lc);
     return InterpolateUBS (Curve, u, derivee);
-
+    
   case MSH_SEGM_NURBS:
-    V.lc = (u * Curve->beg->lc + (1. - u) * Curve->end->lc);
     return InterpolateNurbs (Curve, u, derivee);
 
   case MSH_SEGM_SPLN:
-    V.lc = (u * Curve->beg->lc + (1. - u) * Curve->end->lc);
     N = List_Nbr (Curve->Control_Points);
 
     /* 
@@ -138,7 +134,7 @@ Vertex InterpolateCurve (Curve * Curve, double u, int derivee){
     List_Read (Curve->Control_Points, i, &v[1]);
     List_Read (Curve->Control_Points, i + 1, &v[2]);
     
-    V.lc = t * v[1]->lc + (1. - t) * v[2]->lc;
+    V.lc = t * v[2]->lc + (1. - t) * v[1]->lc; //?????
     
     if (!i){
       v[0] = &temp1;
diff --git a/Mesh/Nurbs.cpp b/Mesh/Nurbs.cpp
index 196a1d08bf92f28e9822a2e3fec1f9d44b32ca9c..ced535206d890b3d9a42773abb3a0c180935e86a 100644
--- a/Mesh/Nurbs.cpp
+++ b/Mesh/Nurbs.cpp
@@ -1,4 +1,4 @@
-// $Id: Nurbs.cpp,v 1.5 2001-01-08 08:05:46 geuzaine Exp $
+// $Id: Nurbs.cpp,v 1.6 2001-08-17 07:41:58 geuzaine Exp $
 
 #include "Gmsh.h"
 #include "Mesh.h"
@@ -10,7 +10,8 @@ Vertex InterpolateCubicSpline (Vertex * v[4], double t, double mat[4][4],
   double vec[4], T[4];
 
   V.Pos.X = V.Pos.Y = V.Pos.Z = 0.0;
-  V.lc = t * v[1]->lc + (1. - t) * v[2]->lc;
+  //V.lc = t * v[1]->lc + (1. - t) * v[2]->lc;
+  V.lc = (1-t) * v[1]->lc + t * v[2]->lc; // ???????
 
   if (derivee){
     T[3] = 0.;
diff --git a/demos/splines.geo b/demos/splines.geo
new file mode 100644
index 0000000000000000000000000000000000000000..db9b753474a4359cba5838fc7c81d730ae2dfdaa
--- /dev/null
+++ b/demos/splines.geo
@@ -0,0 +1,49 @@
+
+lc = 1e-1 ;
+xx = 0;
+
+// Lines
+p = newp;
+Point(p) = {xx,  0,  0, 0.3*lc} ;
+Point(p+1) = {xx+.5, 0,  0, lc} ;
+Point(p+2) = {xx+.5, 1, 0, 0.1*lc} ;
+Point(p+3) = {xx, 1, 0, lc} ;
+l = newreg;
+Line(l) = {p+3,p+2};
+Line(l+1) = {p+2,p+1};
+Line(l+2) = {p+1,p};
+Line(l+3) = {p,p+3};
+s = newreg;
+Line Loop(s) = {-l,-(l+1),-(l+2),-(l+3)};
+Plane Surface(s+1) = {s};
+
+xx += 1;
+
+// B-Splines
+p = newp;
+Point(p) = {xx,  0,  0, 0.3*lc} ;
+Point(p+1) = {xx+.5, 0,  0, lc} ;
+Point(p+2) = {xx+.5, 1, 0, 0.1*lc} ;
+Point(p+3) = {xx, 1, 0, lc} ;
+l = newreg;
+BSpline(l) = {p+3,p+3,p+3,p+2,p+1,p+1,p,p,p};
+BSpline(l+1) = {p,p,p,p+3,p+3,p+3};
+s = newreg;
+Line Loop(s) = {-l,-(l+1)};
+Plane Surface(s+1) = s;
+
+xx += 1;
+
+// Splines (CatmullRom)
+p = newp;
+Point(p) = {xx,  0,  0, 0.3*lc} ;
+Point(p+1) = {xx+.5, 0,  0, lc} ;
+Point(p+2) = {xx+.5, 1, 0, 0.1*lc} ;
+Point(p+3) = {xx, 1, 0, lc} ;
+l = newreg;
+Spline(l) = {p+3,p+2,p+1,p};
+Spline(l+1) = {p,p+3};
+s = newreg;
+Line Loop(s) = {-l,-(l+1)};
+Plane Surface(s+1) = s;
+