subroutine splcl(n,x,y,cox1,coy1,hs1,hs2) C******************************************************************************* C C splcl compute the cubic spline interpolation for closed curve in 2D C given n points (x(i),y(i)) on the curve. We use the length C hs1(1) = 0, hs1(i+1) = hs1(i) + dsqrt( (x(i+1)-x(i))^2 C + (y(i+1)-y(i))^2 ) C as the parameter. The spline is determined by the moments x''(s(i)) C and y''(s(i)). The computed spline in the interval [hs1(i), hs1(i+1)] C for x is: C C x = x(i) + a1 (s-hs1(i)) + a2 (s-hs1(i))^2 + a3 (s-hs1(i))^3 C C where C C a1 = (x(i+1)-x(i))/hs2(i+1) - (2 cox1(i) + cox1(i+1)) hs2(i+1)/6 C a2 = cox1(i)/2 a3 = (cox1(i+1)-cox1(i))/(6 hs2(i+1) ) C C and we have similar formula for y. The computation at some point s C hs1(i) <= s <= hs1(n+1) is carried in the subroutine splval. C This routine calls subroutine splint and uses the cmlib. C C n: On entry, The number of the knots (x_i, y_i). C x: On entry, x(n), The x coordinate of the knots. C y: On entry, y(n), The y coordinate of the knots. C cox1(n+1): Output, The moments of the spline in the x direction as C explained above, C coy1(n+1): Output, The moments of the spline in y direction. C hs1(n+1), hs2(n+1): Output, hs1(i) is the parameter representation of C the knots. hs2(i+1) =hs1(i+1) - hs2(i). C d,c,e,work1, work2, work3: are all n dimensinal array, working space. C C Reference: J. Stoer & R. Bulirsch: Introduction to Numerical Analysis C Springer-Verlag, 1983 C Author: Zhilin Li, December, 1992 C C******************************************************************************* implicit double precision (a-h,o-z) parameter(nn = 20000) dimension x(n+1),y(n+1),cox1(n+1),coy1(n+1),d(nn),c(nn),e(nn) 1 ,work1(nn),work2(nn), work3(nn),hs1(n+1),hs2(n+1) hs1(1) = 0.0d0 do 10 i=1,n-1 hs1(i+1) = hs1(i) + dsqrt( (x(i+1)-x(i))*(x(i+1)-x(i)) + 1 (y(i+1)-y(i))*(y(i+1)-y(i)) ) hs2(i+1) = hs1(i+1) - hs1(i) 10 continue hs1(n+1) = hs1(n) + dsqrt( (x(1)-x(n))*(x(1)-x(n)) + 1 (y(1)-y(n))*(y(1)-y(n)) ) hs2(n+1) = hs1(n+1) - hs1(n) call splint(n,x,hs2,cox1) call splint(n,y,hs2,coy1) return end C************************************************************************** subroutine splint(n,x,hs,cox1) implicit double precision (a-h,o-z) parameter(nn = 20000) dimension x(n+1),cox1(n+1),d(nn),c(nn),e(nn), 1 work1(nn),work2(nn), work3(nn),hs(n+1) C First get the coefficients of the linear system for the moments. do 10 i=2,n-1 d(i) = 2.0d0 temp = 1.0d0/(hs(i) + hs(i+1)) cox1(i) = 6.0d0*temp*( (x(i+1)-x(i))/hs(i+1) - 1 (x(i)-x(i-1))/hs(i) ) c(i) = hs(i) * temp e(i) = hs(i+1)* temp 10 continue d(n) = 2.0d0 temp = 1.0d0/(hs(n) + hs(n+1)) cox1(i) = 6.0d0*temp*( (x(1)-x(n))/hs(n+1) - 1 (x(n)-x(n-1))/hs(n) ) c(n) = hs(n) * temp e(n) = hs(n+1)* temp d(1) = 2.0d0 temp = 1.0d0/(hs(2) + hs(n+1)) cox1(1) = 6.0d0*temp*( (x(2)-x(1))/hs(2) - 1 (x(1)-x(n))/hs(n+1) ) e(1) = hs(2) * temp c(1) = hs(n+1)* temp C Use self-written program to solve the special system. call dgltsl(n,c,d,e,cox1,work1,work2,work3) cox1(n+1) = cox1(1) return end C************************************************************************** subroutine splval(n,s,x,y,cox1,coy1,hs1,hs2,x1,y1,info) C******************************************************************************* C Subroutine splval compute the value of the spline got from the splcl C n, x, y, cox1,coy1,hs1,hs2: On entry produced from splcl. C C s: On entry, The point where the spline is going to be calculated. C x1,y1: Output: The computed coordinates (x1,y1). C info: Output: info = 0, s is in the range of the definition. x1,y1 C are returned. If info = 1, s is out of the range of the C definition. No x1,y1 are returned. C C******************************************************************************* implicit double precision (a-h,o-z) dimension x(n+1),y(n+1),cox1(n+1),coy1(n+1), 1 hs1(n+1),hs2(n+1) do 10 i=1,n if(s .ge. hs1(i) .and. s .le. hs1(i+1)) then xt = x(i+1) a1 = (xt-x(i))/hs2(i+1) - (2.0d0*cox1(i) + 1 cox1(i+1) )*hs2(i+1)/6.0d0 a2 = cox1(i)/2.0d0 a3 = (cox1(i+1) - cox1(i) )/(hs2(i+1)*6.0d0) t = s - hs1(i) x1 = x(i) + t*(a1 + t*(a2 + a3*t)) yt = y(i+1) a1 = (yt-y(i))/hs2(i+1) - (2.0d0*coy1(i) + 1 coy1(i+1) )*hs2(i+1)/6.0d0 a2 = coy1(i)/2.0d0 a3 = (coy1(i+1) - coy1(i) )/(hs2(i+1)*6.0d0) y1 = y(i) + t*(a1 + t*(a2 + a3*t)) info = 0 return endif 10 continue info = 1 return end C************************************************************************** subroutine splval2(s,x1,x2,cox1,cox2,h1,hs2,x3) implicit double precision (a-h,o-z) a1 = (x2-x1)/hs2- (2.0d0*cox1+cox2)*hs2/6.0d0 a2 = cox1/2.0d0 a3 = (cox2 - cox1)/(hs2*6.0d0) t = s - h1 x3 = x1 + t*(a1 + t*(a2 + a3*t)) return end C C************************************************************************** subroutine splval3(s,x1,x2,cox1,cox2,h1,hs2,dxs,dxss) implicit double precision (a-h,o-z) a1 = (x2-x1)/hs2- (2.0d0*cox1+cox2)*hs2/6.0d0 a2 = cox1/2.0d0 a3 = (cox2 - cox1)/(hs2*6.0d0) t = s - h1 dxs = a1 + t*(2.0d0*a2 + 3.0d0*a3*t) dxss = a2 + 6.0d0*t*a3 return end C************************************************************************** subroutine splval4(s,x1,x2,cox1,cox2,h1,hs2,x3,dxs,dxss) implicit double precision (a-h,o-z) a1 = (x2-x1)/hs2- (2.0d0*cox1+cox2)*hs2/6.0d0 a2 = cox1/2.0d0 a3 = (cox2 - cox1)/(hs2*6.0d0) t = s - h1 x3 = x1 + t*(a1 + t*(a2 + a3*t)) dxs = a1 + t*(2.0d0*a2 + 3.0d0*a3*t) dxss = 2.0d0*a2 + 6.0d0*t*a3 return end C-------------------------------------------------------------------------------