Geometry/radec_distance_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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
subroutine radec_distance_3d ( ra1, dec1, ra2, dec2, theta )

!*****************************************************************************80
!
!! RADEC_DISTANCE_3D - angular distance, astronomical units, sphere in 3D.
!
!  Discussion:
!
!    Right ascension is measured in hours, between 0 and 24, and
!    essentially measures longitude.
!
!    Declination measures the angle from the equator towards the north pole,
!    and ranges from -90 (South Pole) to 90 (North Pole).
!
!    On the unit sphere, the angular separation between two points is
!    equal to their geodesic or great circle distance.  On any other
!    sphere, multiply the angular separation by the radius of the
!    sphere to get the geodesic or great circle distance.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    02 December 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, double precision RA1, DEC1, RA2, DEC2, the right ascension and
!    declination of the two points.
!
!    Output, double precision THETA, the angular separation between the points,
!    in radians.
!
  implicit none

  integer, parameter :: dim_num = 3

  double precision cos_theta
  double precision dec1
  double precision dec2
  double precision degrees_to_radians
  double precision norm_v1
  double precision norm_v2
  double precision phi1
  double precision phi2
  double precision r8_acos
  double precision ra1
  double precision ra2
  double precision theta
  double precision theta1
  double precision theta2
  double precision v1(dim_num)
  double precision v2(dim_num)

  theta1 = degrees_to_radians ( 15.0D+00 * ra1 )
  phi1 = degrees_to_radians ( dec1 )

  v1(1:dim_num) = (/ cos ( theta1 ) * cos ( phi1 ), &
                     sin ( theta1 ) * cos ( phi1 ), &
                                      sin ( phi1 ) /)

  norm_v1 = sqrt ( sum ( v1(1:dim_num)**2 ) )

  theta2 = degrees_to_radians ( 15.0D+00 * ra2 )
  phi2 = degrees_to_radians ( dec2 )

  v2(1:dim_num) = (/ cos ( theta2 ) * cos ( phi2 ), &
                     sin ( theta2 ) * cos ( phi2 ), &
                                      sin ( phi2 ) /)

  norm_v2 = sqrt ( sum ( v2(1:dim_num)**2 ) )

  cos_theta = dot_product ( v1(1:dim_num), v2(1:dim_num) ) &
    / ( norm_v1 * norm_v2 )

  theta = r8_acos ( cos_theta )

  return
end