From 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 Mon Sep 17 00:00:00 2001
From: Joe Hunkeler <jhunkeler@gmail.com>
Date: Tue, 11 Aug 2015 16:51:37 -0400
Subject: Repatch (from linux) of OSX IRAF

---
 math/deboor/progs/prog1.f | 45 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 45 insertions(+)
 create mode 100644 math/deboor/progs/prog1.f

(limited to 'math/deboor/progs/prog1.f')

diff --git a/math/deboor/progs/prog1.f b/math/deboor/progs/prog1.f
new file mode 100644
index 00000000..2a0b869c
--- /dev/null
+++ b/math/deboor/progs/prog1.f
@@ -0,0 +1,45 @@
+chapter ii. runge example
+c  from  * a practical guide to splines *  by c. de boor
+      integer i,istep,j,k,n,nmk,nm1
+      real aloger,algerp,d(20),decay,dx,errmax,g,h,pnatx,step,tau(20),x
+      data step, istep /20., 20/
+      g(x) = 1./(1.+(5.*x)**2)
+      print 600
+  600 format(28h  n   max.error   decay exp.//)
+      decay = 0.
+      do 40 n=2,20,2
+c	 choose interpolation points  tau(1), ..., tau(n) , equally
+c	 spaced in (-1,1), and set  d(i) = g(tau(i)), i=1,...,n.
+	 nm1 = n-1
+	 h = 2./float(nm1)
+	 do 10 i=1,n
+	    tau(i) = float(i-1)*h - 1.
+   10	    d(i) = g(tau(i))
+c	 calculate the divided differences for the newton form.
+c
+	 do 20 k=1,nm1
+	    nmk = n-k
+	    do 20 i=1,nmk
+   20	       d(i) = (d(i+1)-d(i))/(tau(i+k)-tau(i))
+c
+c	 estimate max.interpolation error on (-1,1).
+	 errmax = 0.
+	 do 30 i=2,n
+	    dx = (tau(i)-tau(i-1))/step
+	    do 30 j=1,istep
+	       x = tau(i-1) + float(j)*dx
+c	       evaluate interp.pol. by nested multiplication
+c
+	       pnatx = d(1)
+	       do 29 k=2,n
+   29		  pnatx = d(k) + (x-tau(k))*pnatx
+c
+   30	       errmax = amax1(errmax,abs(g(x)-pnatx))
+	 aloger = alog(errmax)
+	 if (n .gt. 2)	decay =
+     *	     (aloger - algerp)/alog(float(n)/float(n-2))
+	 algerp = aloger
+   40	 print 640,n,errmax,decay
+  640 format(i3,e12.4,f11.2)
+					stop
+      end
-- 
cgit