1      real*8 function dblint(tx,nx,ty,ny,c,kx,ky,xb,xe,yb,ye,wrk)
2c  function dblint calculates the double integral
3c         / xe  / ye
4c        |     |      s(x,y) dx dy
5c    xb /  yb /
6c  with s(x,y) a bivariate spline of degrees kx and ky, given in the
7c  b-spline representation.
8c
9c  calling sequence:
10c     aint = dblint(tx,nx,ty,ny,c,kx,ky,xb,xe,yb,ye,wrk)
11c
12c  input parameters:
13c   tx    : real array, length nx, which contains the position of the
14c           knots in the x-direction.
15c   nx    : integer, giving the total number of knots in the x-direction
16c   ty    : real array, length ny, which contains the position of the
17c           knots in the y-direction.
18c   ny    : integer, giving the total number of knots in the y-direction
19c   c     : real array, length (nx-kx-1)*(ny-ky-1), which contains the
20c           b-spline coefficients.
21c   kx,ky : integer values, giving the degrees of the spline.
22c   xb,xe : real values, containing the boundaries of the integration
23c   yb,ye   domain. s(x,y) is considered to be identically zero out-
24c           side the rectangle (tx(kx+1),tx(nx-kx))*(ty(ky+1),ty(ny-ky))
25c
26c  output parameters:
27c   aint  : real , containing the double integral of s(x,y).
28c   wrk   : real array of dimension at least (nx+ny-kx-ky-2).
29c           used as working space.
30c           on exit, wrk(i) will contain the integral
31c                / xe
32c               | ni,kx+1(x) dx , i=1,2,...,nx-kx-1
33c           xb /
34c           with ni,kx+1(x) the normalized b-spline defined on
35c           the knots tx(i),...,tx(i+kx+1)
36c           wrk(j+nx-kx-1) will contain the integral
37c                / ye
38c               | nj,ky+1(y) dy , j=1,2,...,ny-ky-1
39c           yb /
40c           with nj,ky+1(y) the normalized b-spline defined on
41c           the knots ty(j),...,ty(j+ky+1)
42c
43c  other subroutines required: fpintb
44c
45c  references :
46c    gaffney p.w. : the calculation of indefinite integrals of b-splines
47c                   j. inst. maths applics 17 (1976) 37-41.
48c    dierckx p. : curve and surface fitting with splines, monographs on
49c                 numerical analysis, oxford university press, 1993.
50c
51c  author :
52c    p.dierckx
53c    dept. computer science, k.u.leuven
54c    celestijnenlaan 200a, b-3001 heverlee, belgium.
55c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
56c
57c  latest update : march 1989
58c
59c  ..scalar arguments..
60      integer nx,ny,kx,ky
61      real*8 xb,xe,yb,ye
62c  ..array arguments..
63      real*8 tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),wrk(nx+ny-kx-ky-2)
64c  ..local scalars..
65      integer i,j,l,m,nkx1,nky1
66      real*8 res
67c  ..
68      nkx1 = nx-kx-1
69      nky1 = ny-ky-1
70c  we calculate the integrals of the normalized b-splines ni,kx+1(x)
71      call fpintb(tx,nx,wrk,nkx1,xb,xe)
72c  we calculate the integrals of the normalized b-splines nj,ky+1(y)
73      call fpintb(ty,ny,wrk(nkx1+1),nky1,yb,ye)
74c  calculate the integral of s(x,y)
75      dblint = 0.
76      do 200 i=1,nkx1
77        res = wrk(i)
78        if(res.eq.0.) go to 200
79        m = (i-1)*nky1
80        l = nkx1
81        do 100 j=1,nky1
82          m = m+1
83          l = l+1
84          dblint = dblint+res*wrk(l)*c(m)
85 100    continue
86 200  continue
87      return
88      end
89