Repatch (from linux) of OSX IRAF

1 files changed, 138 insertions, 0 deletions

diff --git a/math/deboor/bvalue.f b/math/deboor/bvalue.f
new file mode 100644
index 00000000..6b038939
--- /dev/null
+++ b/math/deboor/bvalue.f
@@ -0,0 +1,138 @@
+      real function bvalue ( t, bcoef, n, k, x, jderiv )
+c  from  * a practical guide to splines *  by c. de boor
+calls  interv
+c
+calculates value at  x	of  jderiv-th derivative of spline from b-repr.
+c  the spline is taken to be continuous from the right.
+c
+c******  i n p u t ******
+c  t, bcoef, n, k......forms the b-representation of the spline  f  to
+c	 be evaluated. specifically,
+c  t.....knot sequence, of length  n+k, assumed nondecreasing.
+c  bcoef.....b-coefficient sequence, of length	n .
+c  n.....length of  bcoef  and dimension of spline(k,t),
+c	 a s s u m e d	positive .
+c  k.....order of the spline .
+c
+c  w a r n i n g . . .	 the restriction  k .le. kmax (=20)  is imposed
+c	 arbitrarily by the dimension statement for  aj, dl, dr  below,
+c	 but is  n o w h e r e	c h e c k e d  for.
+c
+c  x.....the point at which to evaluate .
+c  jderiv.....integer giving the order of the derivative to be evaluated
+c	 a s s u m e d	to be zero or positive.
+c
+c******  o u t p u t  ******
+c  bvalue.....the value of the (jderiv)-th derivative of  f  at  x .
+c
+c******  m e t h o d  ******
+c     the nontrivial knot interval  (t(i),t(i+1))  containing  x  is lo-
+c  cated with the aid of  interv . the	k  b-coeffs of	f  relevant for
+c  this interval are then obtained from  bcoef (or taken to be zero if
+c  not explicitly available) and are then differenced  jderiv  times to
+c  obtain the b-coeffs of  (d**jderiv)f  relevant for that interval.
+c  precisely, with  j = jderiv, we have from x.(12) of the text that
+c
+c     (d**j)f  =  sum ( bcoef(.,j)*b(.,k-j,t) )
+c
+c  where
+c		    / bcoef(.), 		    ,  j .eq. 0
+c		    /
+c    bcoef(.,j)  =  / bcoef(.,j-1) - bcoef(.-1,j-1)
+c		    / ----------------------------- ,  j .gt. 0
+c		    /	 (t(.+k-j) - t(.))/(k-j)
+c
+c     then, we use repeatedly the fact that
+c
+c    sum ( a(.)*b(.,m,t)(x) )  =  sum ( a(.,x)*b(.,m-1,t)(x) )
+c  with
+c		  (x - t(.))*a(.) + (t(.+m-1) - x)*a(.-1)
+c    a(.,x)  =	  ---------------------------------------
+c		  (x - t(.))	  + (t(.+m-1) - x)
+c
+c  to write  (d**j)f(x)  eventually as a linear combination of b-splines
+c  of order  1 , and the coefficient for  b(i,1,t)(x)  must then be the
+c  desired number  (d**j)f(x). (see x.(17)-(19) of text).
+c
+c     parameter kmax = 20
+      integer jderiv,k,n,   i,ilo,imk,j,jc,jcmin,jcmax,jj,kmj,km1,mflag
+     *			   ,nmi,jdrvp1
+      real bcoef(n),t(1),x,   aj(20),dl(20),dr(20),fkmj
+c     real bcoef(n),t(1),x,   aj(kmax),dl(kmax),dr(kmax),fkmj
+c     dimension t(n+k)
+current fortran standard makes it impossible to specify the length of  t
+c  precisely without the introduction of otherwise superfluous addition-
+c  al arguments.
+      bvalue = 0.
+      if (jderiv .ge. k)		go to 99
+c
+c  *** find  i	 s.t.	1 .le. i .lt. n+k   and   t(i) .lt. t(i+1)   and
+c      t(i) .le. x .lt. t(i+1) . if no such i can be found,  x	lies
+c      outside the support of  the spline  f  and  bvalue = 0.
+c      (the asymmetry in this choice of  i  makes  f  rightcontinuous)
+      call interv ( t, n+k, x, i, mflag )
+      if (mflag .ne. 0) 		go to 99
+c  *** if k = 1 (and jderiv = 0), bvalue = bcoef(i).
+      km1 = k - 1
+      if (km1 .gt. 0)			go to 1
+      bvalue = bcoef(i)
+					go to 99
+c
+c  *** store the k b-spline coefficients relevant for the knot interval
+c     (t(i),t(i+1)) in aj(1),...,aj(k) and compute dl(j) = x - t(i+1-j),
+c     dr(j) = t(i+j) - x, j=1,...,k-1 . set any of the aj not obtainable
+c     from input to zero. set any t.s not obtainable equal to t(1) or
+c     to t(n+k) appropriately.
+    1 jcmin = 1
+      imk = i - k
+      if (imk .ge. 0)			go to 8
+      jcmin = 1 - imk
+      do 5 j=1,i
+    5	 dl(j) = x - t(i+1-j)
+      do 6 j=i,km1
+	 aj(k-j) = 0.
+    6	 dl(j) = dl(i)
+					go to 10
+    8 do 9 j=1,km1
+    9	 dl(j) = x - t(i+1-j)
+c
+   10 jcmax = k
+      nmi = n - i
+      if (nmi .ge. 0)			go to 18
+      jcmax = k + nmi
+      do 15 j=1,jcmax
+   15	 dr(j) = t(i+j) - x
+      do 16 j=jcmax,km1
+	 aj(j+1) = 0.
+   16	 dr(j) = dr(jcmax)
+					go to 20
+   18 do 19 j=1,km1
+   19	 dr(j) = t(i+j) - x
+c
+   20 do 21 jc=jcmin,jcmax
+   21	 aj(jc) = bcoef(imk + jc)
+c
+c		*** difference the coefficients  jderiv  times.
+      if (jderiv .eq. 0)		go to 30
+      do 23 j=1,jderiv
+	 kmj = k-j
+	 fkmj = float(kmj)
+	 ilo = kmj
+	 do 23 jj=1,kmj
+	    aj(jj) = ((aj(jj+1) - aj(jj))/(dl(ilo) + dr(jj)))*fkmj
+   23	    ilo = ilo - 1
+c
+c  *** compute value at  x  in (t(i),t(i+1)) of jderiv-th derivative,
+c     given its relevant b-spline coeffs in aj(1),...,aj(k-jderiv).
+   30 if (jderiv .eq. km1)		go to 39
+      jdrvp1 = jderiv + 1
+      do 33 j=jdrvp1,km1
+	 kmj = k-j
+	 ilo = kmj
+	 do 33 jj=1,kmj
+	    aj(jj) = (aj(jj+1)*dl(ilo) + aj(jj)*dr(jj))/(dl(ilo)+dr(jj))
+   33	    ilo = ilo - 1
+   39 bvalue = aj(1)
+c
+   99					return
+      end