cr_spline_interpol.F File Reference

#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	cr_spline_interpol (pts, knots, t, c, c_d, c_dd)

Function/Subroutine Documentation

cr_spline_interpol()

subroutine cr_spline_interpol	(	dimension(0:3,3), intent(in)	pts,
		dimension(0:3), intent(in)	knots,
		intent(in)	t,
		dimension(3), intent(inout)	c,
		dimension(3), intent(inout)	c_d,
		dimension(3), intent(inout)	c_dd )

Definition at line 29 of file cr_spline_interpol.F.

C-----------------------------------------------
        !FOR A GIVEN POSITION T in [0,1] which is position on [P1,P2], GET VALUE ON SPLINE P0-P1-P2-P3
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   D e s c r i p t i o n
C-----------------------------------------------
        ! INPUT  - PTS(0:3,:) are 4 control points
        ! INPUT  - T is in [0,1]  : T=0 => VAL = P(2,:)
        ! OUTPUT - C is interpolated value.
        ! OUTPUT - C_D is 1st derivative
        ! OUTPUT - C_DD is 2nd derivatice
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
        my_real, INTENT(IN)    :: pts(0:3,3),knots(0:3),t
        my_real, INTENT(INOUT) :: c(3),c_d(3),c_dd(3)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
        my_real :: a1(3),a2(3),a3(3)
        my_real :: b1(3),b2(3)
        my_real :: a1p(3),a2p(3),a3p(3), b1p(3),b2p(3),b1pp(3),b2pp(3)
        my_real :: tt  !in [T1,T2]=[KNOT(1),KNOT(2)]
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
        tt=knots(1)+t*(knots(2)-knots(1))
        a1(1:3)   = ((knots(1)-tt)*pts(0,1:3)+(tt-knots(0))*pts(1,1:3)) /(knots(1)-knots(0))
        a2(1:3)   = ((knots(2)-tt)*pts(1,1:3)+(tt-knots(1))*pts(2,1:3)) /(knots(2)-knots(1))
        a3(1:3)   = ((knots(3)-tt)*pts(2,1:3)+(tt-knots(2))*pts(3,1:3)) /(knots(3)-knots(2))
        b1(1:3)   = ((knots(2)-tt)*a1(1:3) + (tt-knots(0))*a2(1:3)) / (knots(2)-knots(0))
        b2(1:3)   = ((knots(3)-tt)*a2(1:3) + (tt-knots(1))*a3(1:3)) / (knots(3)-knots(1))
        c(1:3)    = ((knots(2)-tt)*b1(1:3) + (tt-knots(1))*b2(1:3)) / (knots(2)-knots(1))
        a1p(1:3)  = (pts(1,1:3)-pts(0,1:3))/(knots(1)-knots(0))
        a2p(1:3)  = (pts(2,1:3)-pts(1,1:3))/(knots(2)-knots(1))
        a3p(1:3)  = (pts(3,1:3)-pts(2,1:3))/(knots(3)-knots(2))
        b1p(1:3)  =  (a2(1:3)-a1(1:3))/(knots(2)-knots(0))
     .            + (knots(2)-tt)/(knots(2)-knots(0))*a1p(1:3)
     .            + (tt-knots(0))/(knots(2)-knots(0))*a2p(1:3)
        b2p(1:3)  =  (a3(1:3)-a2(1:3))/(knots(3)-knots(1))
     .            + (knots(3)-tt)/(knots(3)-knots(1))*a2p(1:3)
     .            + (tt-knots(1))/(knots(3)-knots(1))*a3p(1:3)
        c_d(1:3)  =  (b2(1:3)-b1(1:3))/(knots(2)-knots(1))
     .            + (knots(2)-tt)/(knots(2)-knots(1))*b1p(1:3)
     .            + (tt-knots(1))/(knots(2)-knots(1))*b2p(1:3)
        b1pp(1:3) =  (2.*a2p(1:3)-2.*a1p(1:3))/(knots(2)-knots(0))
        b2pp(1:3) =  (2.*a3p(1:3)-2.*a2p(1:3))/(knots(3)-knots(1))
        c_dd(1:3) = (2.*b2p(1:3)-2.*b1p(1:3)+(knots(2)-tt)*b1pp(1:3)+(tt-knots(1))*b2pp(1:3))/(knots(2)-knots(1))
        RETURN