Geometry/shape_point_dist_2d.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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
subroutine shape_point_dist_2d ( pc, p1, side_num, p, dist )

!*****************************************************************************80
!
!! SHAPE_POINT_DIST_2D: distance ( regular shape, point ) in 2D.
!
!  Discussion:
!
!    The "regular shape" is assumed to be an equilateral and equiangular
!    polygon, such as the standard square, pentagon, hexagon, and so on.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    04 January 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, double precision PC(2), the center of the shape.
!
!    Input, double precision P1(2), the first vertex of the shape.
!
!    Input, integer SIDE_NUM, the number of sides.
!
!    Input, double precision P(2), the point to be checked.
!
!    Output, double precision DIST, the distance from the point to the shape.
!
  implicit none

  integer, parameter :: dim_num = 2

  double precision angle
  double precision angle_deg_2d
  double precision angle2
  double precision degrees_to_radians
  double precision dist
  double precision p(dim_num)
  double precision p1(dim_num)
  double precision pa(dim_num)
  double precision pb(dim_num)
  double precision pc(dim_num)
  double precision, parameter :: r8_pi = 3.141592653589793D+00
  double precision radius
  double precision sector_angle
  integer sector_index
  integer side_num
!
!  Determine the angle subtended by a single side.
!
  sector_angle = 360.0D+00 / dble(side_num)
!
!  How long is the half-diagonal?
!
  radius = sqrt ( sum ( ( p1(1:dim_num) - pc(1:dim_num) )**2 ) )
!
!  If the radius is zero, then the shape is a point and the computation is easy.
!
  if ( radius == 0.0D+00 ) then
    dist = sqrt ( sum ( ( p(1:dim_num) - pc(1:dim_num) )**2 ) )
    return
  end if
!
!  If the test point is at the pc, then the computation is easy.
!  The angle subtended by any side is ( 2 * PI / SIDE_NUM ) and the
!  nearest distance is the midpoint of any such side.
!
  if ( all ( p(1:dim_num) == pc(1:dim_num) ) ) then
    dist = radius * cos ( r8_pi / dble(side_num) )
    return
  end if
!
!  Determine the angle between the ray to the first corner,
!  and the ray to the test point.
!
  angle = angle_deg_2d ( p1(1:2), pc(1:2), p(1:2) )
!
!  Determine the sector of the point.
!
  sector_index = int ( angle / sector_angle ) + 1
!
!  Generate the two corner points that terminate the SECTOR-th side.
!
  angle2 = dble(sector_index - 1) * sector_angle
  angle2 = degrees_to_radians ( angle2 )

  call vector_rotate_base_2d ( p1, pc, angle2, pa )

  angle2 = dble(sector_index) * sector_angle
  angle2 = degrees_to_radians ( angle2 )

  call vector_rotate_base_2d ( p1, pc, angle2, pb )
!
!  Determine the distance from the test point to the line segment that
!  is the SECTOR-th side.
!
  call segment_point_dist_2d ( pa, pb, p, dist )

  return
end