SUBROUTINE AKM_CSPL (n,x,f,c) * DESCRIPTION * Subroutine akm_cspl is based on the method of interpolation described * by Hiroshi Akima (1970 Journal of the Association for Computing Machinery, * Vol. 17, pp. 580-602) and also by Carl de Boor (1978 Applied Mathematical * Sciences, Vol. 27, "A Practical Guide to Splines"). * This routine computes the Akima spline interpolant with the "not-a-knot * end conditions which are useful when there is no information about the * derivatives at the endpoints. In this case P(1)=P(2) and P(n-2)=P(n-1); * in other words, x(2) and x(n-1) are not knots even though they define * the endpoints of a polynomial piece --- the spline must pass through * these two points. Thus, it is a requirement that f''' be continuous * across x(2) and x(n-1). * DIMENSIONS * The internal work arrays dd, w, and s are restricted to 2000 elements * by the parameter nmx; consequently the data arrays x and f are also * limited to 2000 elements. The spline coefficients are contained in * the array c, a 4x2000 2 dimensional array. * CALL PARAMETERS * n - number of data pairs * x(n) - array, independent variable * f(n) - array, dependent variable * c(4,n) - array, spline coefficients implicit double precision (a-h, o-z) parameter (nmx = 2000) parameter (zero = 0.0d0, two = 2.0d0, three = 3.0d0) dimension x(n),f(n),c(4,n) dimension dd(nmx),w(nmx),s(nmx) nm1 = n - 1 nm2 = n - 2 nm3 = n - 3 * Create two additional points outside the region of the spline by fitting * a quadratic to the three adjacent data points. This is a special feature * of the Akima spline. call akm_ends (x(1),x(2),x(3),f(1),f(2),f(3),f0,fm1) call akm_ends (x(n),x(nm1),x(nm2),f(n),f(nm1),f(nm2),fnp1,fnp2) * Compute the divided differences ddm1 = (f0 - fm1)/(x(2) - x(1)) dd0 = (f(1) - f0)/(x(3) - x(2)) do i = 1, nm1 ip1 = i + 1 dd(i) = (f(ip1) - f(i))/(x(ip1) - x(i)) end do dd(n) = (f(n) - fnp1)/(x(nm2) - x(nm1)) ddnp1 = (fnp1 - fnp2)/(x(nm1) - x(n)) * Compute the Akima weights w0 = abs (dd0 - ddm1) w(1) = abs (dd(1) - dd0) do i = 2, nm1 im1 = i - 1 w(i) = abs(dd(i) - dd(im1)) end do w(n) = abs (dd(n) - dd(nm1)) wnp1 = abs (ddnp1 - dd(n)) * Compute Akima slopes at interior knots if (w(2).eq.zero.and.w0.eq.zero) then s(1) = 5.0d-1*(dd(1) + dd0) else s(1) = (w(2)*dd0 + w0*dd(1))/(w0 + w(2)) end if do i = 2, nm1 im1 = i - 1 ip1 = i + 1 if (w(ip1).eq.zero.and.w(im1).eq.zero) then s(i) = 5.0d-1*(dd(i) + dd(im1)) else s(i) = (w(ip1)*dd(im1) + w(im1)*dd(i))/(w(im1) + w(ip1)) end if end do if (wnp1.eq.zero.and.w(nm1).eq.zero) then s(n) = 5.0d-1*(dd(n) + dd(nm1)) else s(n) = (wnp1*dd(nm1) + w(nm1)*dd(n))/(w(nm1) + wnp1) end if * Spline coefficients for all the polynomial pieces do i = 1, nm1 ip1 = i + 1 dx = x(ip1) - x(i) c(1,i) = f(i) c(2,i) = s(i) c(3,i) = (three*dd(i) - two*s(i) - s(ip1))/dx c(4,i) = (s(ip1) + s(i) - two*dd(i))/dx/dx end do return end SUBROUTINE AKM_ENDS(x1,x2,x3,f1,f2,f3,f00,f01) * DESCRIPTION * Subroutine akm_ends is based on the method of interpolation described * by Hiroshi Akima (1970 Journal of the Association for Computing Machinery, * Vol. 17, pp. 580-602). * This routine is required by the routine AKM_CSPL. It sets up the two * additional external points required by the Akima method of constructing * a spline. * CALL PARAMETERS * (x1,f1), (x2,f2) and (x3,f3) are the known data pairs at the ends of the * region of interpolation. f00 and f01 are the computed values of the * interpolating polynomial external to the region of interpolation. implicit double precision (a-h,o-z) g0 = f1 df21 = f2 - f1 df31 = f3 - f1 dx21 = x2 - x1 dx31 = x3 - x1 dx32 = x3 - x2 den = dx21*dx32*dx31 g1 = (df21*dx31*dx31 - df31*dx21*dx21)/den g2 = (df31*dx21 - df21*dx31)/den f00 = g0 - g1*dx32 + g2*dx32*dx32 f01 = g0 - g1*dx31 + g2*dx31*dx31 return end SUBROUTINE AKM_EVAL (n,x,c,xp,fp,dfp) * DESCRIPTION * This routine is used to evaluate the spline defined by the * routine AKM_CSPL and its derivative at an interpolating point. * CALL PARAMETERS * n - number of data pairs * x(n) - array, independent variable * c(4,n) - array, spline coefficients * xp - interpolating point * fp - value of the spline at xp * dfp - first derivative of the spline at xp implicit double precision (a-h, o-z) dimension x(n), c(4,n) logical more * Statement function defines cubic spline interpolation csval(dx,a,b,e,d) = a +dx*(b + dx*(e + dx*d)) * Statement function defines cubic spline derivative csder(dx,b,e,d) = b + dx*(2.0d0*e + dx*3.0d0*d) * Check direction of spline if (x(n).gt.x(1)) then * Check to see that x(1) < xp < x(n) if (xp.lt.x(1)) then write(*,'('' ***Test point is less than x(1)***'')') stop else if (xp.gt.x(n)) then write(*,'('' ***Test point is greater than x(n)***'')') stop end if * Find the polynomial piece containing the test point i = 2 more = .true. do while (more) if (xp.lt.x(i).or.i.eq.n) then more = .false. nl = i - 1 else i = i + 1 end if end do else * The spline is backwards * Check to see that x(n) < xp < x(1) if (xp.gt.x(1)) then write(*,'('' ***Test point is greater than x(1)***'')') stop else if (xp.lt.x(n)) then write(*,'('' ***Test point is less than x(n)***'')') stop end if * Find the polynomial piece containing the test point i = 2 more = .true. do while (more) if (xp.gt.x(i).or.i.eq.n) then more = .false. nl = i - 1 else i = i + 1 end if end do end if * Evaluate spline at the test point xp dx = xp - x(nl) fp = csval(dx,c(1,nl),c(2,nl),c(3,nl),c(4,nl)) * Evaluate the derivative at the test point dfp = csder(dx,c(2,nl),c(3,nl),c(4,nl)) return end SUBROUTINE AKM_INT (a,b,n,x,c,sum) * DESCRIPTION * This routine is used to integrate the spline defined by the * routine AKM_CSPL. * CALL PARAMETERS * a - lower bound * b - upper bound * n - number of data pairs * x(n) - array, independent variable * c(4,n) - array, spline coefficients * sum - value of the integral over the interval [a,b] implicit double precision (a-h, o-z) dimension x(n),c(4,n) * Statement function defined by cubic spline integration csint(dx,a,b,g,d) = abs(dx*(a + dx*(b/2.0d0 + dx*(g/3.0d0 1 + dx*d/4.0d0)))) * Check direction of integration if (x(1).lt.x(n)) then * Forward integration * Check limits of integration if (a.lt.x(1)) then write(*,'('' ***Lower bound is less than x(1)***'')') stop else if (b.gt.x(n)) then write(*,'('' ***Upper bound is greater than x(n)***'')') stop end if * Find the polynomial piece containing the lower bound i = 2 do while (a.gt.x(i)) i = i + 1 end do nl = i - 1 * Find the polynomial piece containing the upper bound i = n - 1 do while (b.lt.x(i)) i = i - 1 end do nu = i else * Backward integration * Check limits of integration if (a.gt.x(1)) then write(*,'('' ***Upper bound is greater than x(1)***'')') stop else if (b.lt.x(n)) then write(*,'('' ***Lower bound is less than x(n)***'')') stop end if * Find the polynomial piece containing the upper bound i = 2 do while (a.lt.x(i)) i = i + 1 end do nl = i - 1 * Find the polynomial piece containing the lower bound i = n - 1 do while (b.gt.x(i)) i = i - 1 end do nu = i end if * Initialize the sum sum = 0.0d0 * Subtract the portion of the 1st polynomial piece outside the lower bound dx = a - x(nl) if (dx.ne.0.0d0) then sum = sum - csint(dx,c(1,nl),c(2,nl),c(3,nl),c(4,nl)) end if * Integrate from x(nl) to x(nu) if (nu.gt.nl) then nt = nu - 1 do i = nl, nt dx = x(i+1) - x(i) sum = sum + csint(dx,c(1,i),c(2,i),c(3,i),c(4,i)) end do end if * Add the portion of the last polynomial piece within the upper bound dx = b - x(nu) if (dx.ne.0.0d0) then sum = sum + csint(dx,c(1,nu),c(2,nu),c(3,nu),c(4,nu)) end if return end SUBROUTINE AKM_BISECTION(n,c,x,x1,x2,froot,xroot) * DESCRIPTION * This routine is used to invert the spline defined by the * routine AKM_CSPL in order to find a root by the bisection method. * The routine AKM_EVAL is called within this routine. * CALL PARAMETERS * n - number of data pairs * c(4,n) - array, spline coefficients * x(n) - array, independent variable * x1 - estimated lower limit of region to be searched * x2 - estimated upper limit of region to be searched * froot - value of the spline at the location of the root * xroot - calculated value of the root implicit double precision (a-h, o-z) dimension x(*), c(4,n) logical more * This routine must be preceded by a call to akm_cspl(n,x,f,c) to set up * the spline which is used to find the root. xa = x1 xb = x2 call akm_eval(n,x,c,xa,fa,dfdx) fa = froot - fa call akm_eval(n,x,c,xb,fb,dfdx) fb = froot - fb kl = 0 do while (fa*fb.ge.0.0d0.and.kl.lt.20) kl = kl + 1 if (abs(fa).lt.abs(fb)) then xa = xa + 1.6d0*(xa - xb) xa = max(xa,x(1)) call akm_eval(n,x,c,xa,fa,dfdx) fa = froot - fa else xb = xb + 1.6d0*(xb - xa) xb = min(xb,x(n)) call akm_eval(n,x,c,xb,fb,dfdx) fb = froot - fb end if write(*,'(i4,1pd14.7,0pf7.3,1pd14.7,0pf7.3)') kl,xa,fa,xb,fb end do if (kl.eq.20) then write(*,'('' *** akm_bisection failed to find the region'')') write(*,'('' containing the root after 20 iterations.'')') stop end if diff = abs(xb - xa) more = .true. kl = 0 do while (more) kl = kl + 1 xab = 5.0d-1*(xa + xb) call akm_eval(n,x,c,xab,fb,dfdx) fb = froot - fb if (fa*fb.le.0.0d0) then xb = xab else xa = xab fa = fb end if diff = abs(xb - xa)/abs(5.0d-1*(xa + xb)) more = (diff.gt.1.0d-8.and.kl.lt.40) end do if (kl.eq.40) then write(*,'('' ***akm_bisection root did not converge within'')') write(*,'('' one in 1.0d-8 after 40 iterations.'')') write(*,'('' Converged within '',1pd12.5)') diff end if xroot = 5.0d-1*(xa + xb) return end