Geometry/ball_unit_sample_3d.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
subroutine ball_unit_sample_3d ( seed, p )

!*****************************************************************************80
!
!! BALL_UNIT_SAMPLE_3D picks a random point in the unit ball in 3D.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    26 August 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer SEED, a seed for the random
!    number generator.
!
!    Output, double precision P(3), the sample point.
!
  implicit none

  integer, parameter :: dim_num = 3

  double precision p(dim_num)
  double precision phi
  double precision r
  double precision r8_acos
  double precision, parameter :: r8_pi = 3.141592653589793D+00
  integer seed
  double precision theta
  double precision u(dim_num)
  double precision vdot

  call r8vec_uniform_01 ( dim_num, seed, u )
!
!  Pick a uniformly random VDOT, which must be between -1 and 1.
!  This represents the dot product of the random vector with the Z unit vector.
!
!  Note: this works because the surface area of the sphere between
!  Z and Z + dZ is independent of Z.  So choosing Z uniformly chooses
!  a patch of area uniformly.
!
  vdot = 2.0D+00 * u(1) - 1.0D+00

  phi = r8_acos ( vdot )
!
!  Pick a uniformly random rotation between 0 and 2 Pi around the
!  axis of the Z vector.
!
  theta = 2.0D+00 * r8_pi * u(2)
!
!  Pick a random radius R.
!
  r = u(3) ** ( 1.0D+00 / 3.0D+00 )

  p(1) = r * cos ( theta ) * sin ( phi )
  p(2) = r * sin ( theta ) * sin ( phi )
  p(3) = r * cos ( phi )

  return
end