~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/interpolate/fitpack/insert.f

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
ImportĀ upstreamĀ versionĀ 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
      subroutine insert(iopt,t,n,c,k,x,tt,nn,cc,nest,ier)
 
2
c  subroutine insert inserts a new knot x into a spline function s(x)
 
3
c  of degree k and calculates the b-spline representation of s(x) with
 
4
c  respect to the new set of knots. in addition, if iopt.ne.0, s(x)
 
5
c  will be considered as a periodic spline with period per=t(n-k)-t(k+1)
 
6
c  satisfying the boundary constraints
 
7
c       t(i+n-2*k-1) = t(i)+per  ,i=1,2,...,2*k+1
 
8
c       c(i+n-2*k-1) = c(i)      ,i=1,2,...,k
 
9
c  in that case, the knots and b-spline coefficients returned will also
 
10
c  satisfy these boundary constraints, i.e.
 
11
c       tt(i+nn-2*k-1) = tt(i)+per  ,i=1,2,...,2*k+1
 
12
c       cc(i+nn-2*k-1) = cc(i)      ,i=1,2,...,k
 
13
c
 
14
c  calling sequence:
 
15
c     call insert(iopt,t,n,c,k,x,tt,nn,cc,nest,ier)
 
16
c
 
17
c  input parameters:
 
18
c    iopt : integer flag, specifying whether (iopt.ne.0) or not (iopt=0)
 
19
c           the given spline must be considered as being periodic.
 
20
c    t    : array,length nest, which contains the position of the knots.
 
21
c    n    : integer, giving the total number of knots of s(x).
 
22
c    c    : array,length nest, which contains the b-spline coefficients.
 
23
c    k    : integer, giving the degree of s(x).
 
24
c    x    : real, which gives the location of the knot to be inserted.
 
25
c    nest : integer specifying the dimension of the arrays t,c,tt and cc
 
26
c           nest > n.
 
27
c
 
28
c  output parameters:
 
29
c    tt   : array,length nest, which contains the position of the knots
 
30
c           after insertion.
 
31
c    nn   : integer, giving the total number of knots after insertion
 
32
c    cc   : array,length nest, which contains the b-spline coefficients
 
33
c           of s(x) with respect to the new set of knots.
 
34
c    ier  : error flag
 
35
c      ier = 0 : normal return
 
36
c      ier =10 : invalid input data (see restrictions)
 
37
c
 
38
c  restrictions:
 
39
c    nest > n
 
40
c    t(k+1) <= x <= t(n-k)
 
41
c    in case of a periodic spline (iopt.ne.0) there must be
 
42
c       either at least k interior knots t(j) satisfying t(k+1)<t(j)<=x
 
43
c       or at least k interior knots t(j) satisfying x<=t(j)<t(n-k)
 
44
c
 
45
c  other subroutines required: fpinst.
 
46
c
 
47
c  further comments:
 
48
c   subroutine insert may be called as follows
 
49
c        call insert(iopt,t,n,c,k,x,t,n,c,nest,ier)
 
50
c   in which case the new representation will simply replace the old one
 
51
c
 
52
c  references :
 
53
c    boehm w : inserting new knots into b-spline curves. computer aided
 
54
c              design 12 (1980) 199-201.
 
55
c   dierckx p. : curve and surface fitting with splines, monographs on
 
56
c                numerical analysis, oxford university press, 1993.
 
57
c
 
58
c  author :
 
59
c    p.dierckx
 
60
c    dept. computer science, k.u.leuven
 
61
c    celestijnenlaan 200a, b-3001 heverlee, belgium.
 
62
c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
 
63
c
 
64
c  latest update : march 1987
 
65
c
 
66
c  ..scalar arguments..
 
67
      integer iopt,n,k,nn,nest,ier
 
68
      real*8 x
 
69
c  ..array arguments..
 
70
      real*8 t(nest),c(nest),tt(nest),cc(nest)
 
71
c  ..local scalars..
 
72
      integer kk,k1,l,nk,nk1
 
73
c  ..
 
74
c  before starting computations a data check is made. if the input data
 
75
c  are invalid control is immediately repassed to the calling program.
 
76
      ier = 10
 
77
      if(nest.le.n) go to 40
 
78
      k1 = k+1
 
79
      nk = n-k
 
80
      if(x.lt.t(k1) .or. x.gt.t(nk)) go to 40
 
81
c  search for knot interval t(l) <= x < t(l+1).
 
82
      nk1 = nk-1
 
83
      l = k1
 
84
  10  if(x.lt.t(l+1) .or. l.eq.nk1) go to 20
 
85
      l = l+1
 
86
      go to 10
 
87
  20  if(t(l).ge.t(l+1)) go to 40
 
88
      if(iopt.eq.0) go to 30
 
89
      kk = 2*k
 
90
      if(l.le.kk .and. l.ge.(n-kk)) go to 40
 
91
  30  ier = 0
 
92
c  insert the new knot.
 
93
      call fpinst(iopt,t,n,c,k,x,l,tt,nn,cc,nest)
 
94
  40  return
 
95
      end