Geometry/polyloop_points_nd.f90 at master · lguez/Geometry · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
subroutine polyloop_points_nd ( dim_num, nk, pk, nt, pt )

!*****************************************************************************80
!
!! POLYLOOP_POINTS_ND computes equally spaced points on a polyloop in ND.
!
!  Discussion:
!
!    A polyloop of order NK is the geometric structure consisting of
!    the NK line segments that lie between successive elements of a list
!    of NK points, including a segment from the last point to the first.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    03 March 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer NK, the number of points defining the polyloop.
!
!    Input, double precision PK(DIM_NUM,NK), the points defining the polyloop.
!
!    Input, integer NT, the number of points to be sampled.
!
!    Input, double precision PT(DIM_NUM,NT), equally spaced points
!    on the polyloop.
!
  implicit none

  integer dim_num
  integer nk
  integer nt

  integer it
  integer i4_wrap
  integer j
  integer jp1
  double precision pk(dim_num,nk)
  double precision pt(dim_num,nt)
  double precision sk(nk+1)
  double precision st

  call polyloop_arclength_nd ( dim_num, nk, pk, sk )

  j = 1

  do it = 1,  nt

    st = ( dble(nt - it) * 0.0D+00 + &
           dble(it - 1) * sk(nk+1) ) &
         / dble(nt      - 1)

    do

      if ( sk(j) <= st .and. st <= sk(j+1) ) then
        exit
      end if

      if ( nk <= j ) then
        exit
      end if

      j = j + 1

    end do

    jp1 = i4_wrap ( j + 1, 1, nk )

    pt(1:dim_num,it) = ( ( sk(j+1) - st         ) * pk(1:dim_num,j) &
                       + (           st - sk(j) ) * pk(1:dim_num,jp1) ) &
                       / ( sk(j+1)      - sk(j) )

  end do

  return
end