test_splines.F File Reference

#include "implicit_f.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	test_spline (npt0, nsub, xf, yf, xx, yy)

Function/Subroutine Documentation

test_spline()

subroutine test_spline	(	integer, intent(in)	npt0,
		integer, intent(in)	nsub,
		intent(in)	xf,
		intent(in)	yf,
		intent(out)	xx,
		intent(out)	yy )

Definition at line 28 of file test_splines.F.

C-----------------------------------------------
C   D e s c r i p t i o n
C-----------------------------------------------
       !compute SPLINE length with third order method (SIMPSON)
       !  INPUT  - LOCAL_PT : 4 control points P0,P1,P2,P3
       !  INPUT  - ALPHA : CCR SPLINE PARAMETER
       !  INPUT  - T : position [0,1] on SPLINE [T1,T2].  1.0 means full spline length
       !  OUTPUT - LEN : length of the curve parametrised with t in [0,T]  T<=1.0
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER  ,INTENT(IN)    :: NPT0   ! number of points of input function
      INTEGER  ,INTENT(IN)    :: NSUB   ! 
      my_real  ,DIMENSION(NPT0)   ,INTENT(IN)  :: xf,yf   ! initial curve coordinates
      my_real  ,DIMENSION((NPT0-1)*NSUB+1) ,INTENT(OUT) :: xx,yy   ! curve coordinates build with splines
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I,J,K,NPTS,NSEG
      my_real :: tt,dx,dy,nx,ny
      my_real ,DIMENSION(:,:)  ,ALLOCATABLE :: spline_knots
      my_real ,DIMENSION(:)    ,ALLOCATABLE :: ctrl_ptx,ctrl_pty
      my_real ,DIMENSION(4)  :: ptx,pty,knots
      my_real ,PARAMETER  :: alpha = 0.5
C-----------------------------------------------
C   S o u r c e   L i n e s
c=======================================================================
      npts = npt0 + 2
      nseg = npt0 - 1
 
      ALLOCATE (ctrl_ptx(npts))
      ALLOCATE (ctrl_pty(npts))
      ALLOCATE (spline_knots(nseg,4))
 
c     calculate spline control points
      ctrl_ptx(2:npts-1) = xf(1:npt0)
      ctrl_pty(2:npts-1) = yf(1:npt0)
      ! Add start point - minimum of bending energy            
      ctrl_ptx(1) = (half*ctrl_ptx(2) - four*ctrl_ptx(3) + ctrl_ptx(4)) * half                  
      ctrl_pty(1) = (half*ctrl_pty(2) - four*ctrl_pty(3) + ctrl_pty(4)) * half                  
      ! Add end point - minimum of bending energy               
      ctrl_ptx(npts) = (ctrl_ptx(npts-3) - four*ctrl_ptx(npts-2) + five*ctrl_ptx(npts-1)) * half                  
      ctrl_pty(npts) = (ctrl_pty(npts-3) - four*ctrl_pty(npts-2) + five*ctrl_pty(npts-1)) * half                  
c
      k = 0
      DO i = 1,nseg                         
        ptx(1) = ctrl_ptx(i)                  
        ptx(2) = ctrl_ptx(i+1)                  
        ptx(3) = ctrl_ptx(i+2)                  
        ptx(4) = ctrl_ptx(i+3)
        pty(1) = ctrl_pty(i)                  
        pty(2) = ctrl_pty(i+1)                  
        pty(3) = ctrl_pty(i+2)                  
        pty(4) = ctrl_pty(i+3)
c
        knots(1) = zero
        dx = ptx(2) - ptx(1)
        dy = pty(2) - pty(1)
        knots(2) = spline_knots(i,1) + exp(alpha*log(sqrt(dx**2 + dy**2)))
        dx = ptx(3) - ptx(2)
        dy = pty(3) - pty(2)
        knots(3) = knots(2) + exp(alpha*log(sqrt(dx**2 + dy**2)))
        dx = ptx(4) - ptx(3)
        dy = pty(4) - pty(3)
        knots(4) = knots(3) + exp(alpha*log(sqrt(dx**2 + dy**2)))
        spline_knots(i,1:4) = knots(1:4) 
c
        DO j = 1,nsub                                      
          tt = (j-one) / nsub
          k = k + 1
          CALL spline_interpol_2d(ptx, pty ,knots, tt,  nx  ,ny  )
          xx(k) = nx
          yy(k) = ny
        ENDDO                                            
        
      ENDDO 
      ! last point                                           
      tt = one
      k = k + 1
      xx(k) = xf(npt0)
      yy(k) = yf(npt0)
 
c      CALL SPLINE_INTERPOL_2D(PTX, PTY ,KNOTS, TT,  NX  ,NY  )
c      XX(K) = NX
c      YY(K) = NY
c
      DEALLOCATE (spline_knots)
      DEALLOCATE (ctrl_pty)
      DEALLOCATE (ctrl_ptx)
c-----------
      RETURN