link graphics
procedure main()
    WOpen("size=800,200",
	  "label=spline demo; left to set next point, right to run")
    points := []
    Fg("red")
    while e := Event() do {
	if e === &lpress then {
	    put(points, Point(&x,&y))
	    FillCircle(&x,&y,3)
	}
	else if e === &rrelease then {
	    break
	}
    }
    Fg("black")
    WAttrib("label=are you ready?")
    delay(1000)
    WAttrib("label=are you ready? ...3")
    delay(1000)
    WAttrib("label=are you ready? ...2")
    delay(1000)
    WAttrib("label=are you ready? ...1")
    delay(1000)
    gencurve(points)
    Event()
end

record Point(x,y)
procedure gencurve(p)

#   float  ax, ay, bx, b_y, stepsize, stepsize2, stepsize3;
#   float  x, dx, d2x, d3x, y, dy, d2y, d3y;

   if (p[1].x == p[-1].x) & (p[1].y == p[-1].y) then {
      closed := 1
      push(p, copy(p[-2]))
      put(p, copy(p[2]))
      }
   else {
      push(p,copy(p[1]))
      put(p,copy(p[-1]))
   }

   every i := 4 to *p do {
      #
      # build the coefficients ax, ay, bx and b_y, using:
      #                             _              _   _    _
      #   i                 i    1 | -1   3  -3   1 | | Pi-3 |
      #  Q (t) = T * M   * G   = - |  2  -5   4  -1 | | Pi-2 |
      #               CR    Bs   2 | -1   0   1   0 | | Pi-1 |
      #                            |_ 0   2   0   0_| |_Pi  _|
      #

      ax := p[i].x - 3 * p[i-1].x + 3 * p[i-2].x - p[i-3].x
      ay := p[i].y - 3 * p[i-1].y + 3 * p[i-2].y - p[i-3].y
      bx := 2 * p[i-3].x - 5 * p[i-2].x + 4 * p[i-1].x - p[i].x
      b_y := 2 * p[i-3].y - 5 * p[i-2].y + 4 * p[i-1].y - p[i].y
            #
      # calculate the forward differences for the function using
      # intervals of size 0.1
      #
      steps := max(abs(p[i-1].x - p[i-2].x), abs(p[i-1].y - p[i-2].y)) + 10

      thepoints := [ ]

      stepsize := 1.0 / steps
      stepsize2 := stepsize * stepsize
      stepsize3 := stepsize * stepsize2

      x := p[i-2].x
      y := p[i-2].y
      put(thepoints, x, y)
      dx := (stepsize3*0.5)*ax + (stepsize2*0.5)*bx +
           (stepsize*0.5)*(p[i-1].x-p[i-3].x)
      dy := (stepsize3*0.5)*ay + (stepsize2*0.5)*b_y +
           (stepsize*0.5)*(p[i-1].y-p[i-3].y)
      d2x := (stepsize3*3) * ax + stepsize2 * bx
      d2y := (stepsize3*3) * ay + stepsize2 * b_y
      d3x := (stepsize3*3) * ax
      d3y := (stepsize3*3) * ay

      # calculate the points for drawing (this piece of) the curve

      every 1 to steps do {
         x   +:= dx
         y   +:= dy
         dx  +:= d2x
         dy  +:= d2y
         d2x +:= d3x
         d2y +:= d3y
         put(thepoints, x, y)
	 DrawPoint(x,y)
	 WAttrib("label=x="||integer(x)||", y="||integer(y)||
		 ", dx="||left(dx,5)||", dy="||left(dy,5)||
		", d2x="||left(d2x,5)||", d2y="||left(d2y,5))
	 delay(500)
         }
      DrawLine ! thepoints
      }
end