1 /*
2  *  - - - - - - - - -
3  *   g a l _ p v 2 s
4  *  - - - - - - - - -
5  *
6  *  This routine is part of the General Astrodynamics Library
7  *
8  *  Description:
9  *
10  *  Convert position/velocity from Cartesian to spherical coordinates.
11  *
12  *  This routine is an independent translation of a FORTRAN routine
13  *  that is part of IAU's SOFA software collection.
14  *
15  *  Status:
16  *
17  *     vector/matrix support routine.
18  *
19  *  Given:
20  *
21  *     pv               d[2][3]     pv-vector
22  *
23  *  Returned:
24  *
25  *     *theta              d        longitude angle (radians)
26  *     *phi                d        latitude angle (radians)
27  *     *r                  d        radial distance
28  *     *td                 d        rate of change of THETA
29  *     *pd                 d        rate of change of PHI
30  *     *rd                 d        rate of change of R
31  *
32  *  Notes:
33  *
34  *  1) If the position part of pv is null, theta, phi, td and pd
35  *     are indeterminate.  This is handled by extrapolating the
36  *     position through unit time by using the velocity part of
37  *     pv.  This moves the origin without changing the direction
38  *     of the velocity component.  If the position and velocity
39  *     components of pv are both null, zeroes are returned for all
40  *     six results.
41  *
42  *  2) If the position is a pole, theta, td and pd are indeterminate.
43  *     In such cases zeroes are returned for theta, td and pd.
44  *
45  *  This revision:
46  *
47  *     2001 January 4 ( c version 2008 February 4 )
48  *
49  *
50  *  Copyright (C) 2008 Paul C. L. Willmott. See notes at end.
51  *
52  *-----------------------------------------------------------------------
53  */
54 
55 #include <math.h>
56 #include "gal_pv2s.h"
57 
58 void
gal_pv2s(double pv[2][3],double * theta,double * phi,double * r,double * td,double * pd,double * rd)59 gal_pv2s
60  (
61     double pv[2][3],
62     double *theta,
63     double *phi,
64     double *r,
65     double *td,
66     double *pd,
67     double *rd
68  )
69 {
70 
71     double x, y, z, xd, yd, zd, rxy2, rxy, r2, rtrue, rw, xyp ;
72 
73 /*
74  * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75  */
76 
77 /*
78  * Components of position/velocity vector.
79  */
80 
81     x  = pv[0][0] ;
82     y  = pv[0][1] ;
83     z  = pv[0][2] ;
84     xd = pv[1][0] ;
85     yd = pv[1][1] ;
86     zd = pv[1][2] ;
87 
88 /*
89  * Component of R in XY plane squared.
90  */
91 
92     rxy2 = x * x + y * y ;
93 
94 /*
95  * Modulus squared.
96  */
97 
98     r2 = rxy2 + z * z ;
99 
100 /*
101  * Modulus.
102  */
103 
104     rtrue = sqrt ( r2 ) ;
105 
106 /*
107  * If null vector, move the origin along the direction of movement.
108  */
109 
110     rw = rtrue ;
111     if ( rtrue == 0.0 ) {
112         x = xd ;
113         y = yd ;
114         z = zd ;
115         rxy2 = x * x + y * y ;
116         r2 = rxy2 + z * z ;
117         rw = sqrt ( r2 ) ;
118     }
119 
120 /*
121  * Position and velocity in spherical coordinates.
122  */
123 
124     rxy = sqrt ( rxy2 ) ;
125     xyp = x * xd + y * yd ;
126     if ( rxy2 != 0.0 ) {
127         *theta = atan2 ( y, x ) ;
128         *phi = atan2 ( z, rxy ) ;
129         *td = ( x * yd - y * xd ) / rxy2 ;
130         *pd = ( zd * rxy2 - z * xyp ) / ( r2 * rxy ) ;
131     }
132     else {
133         *theta = 0.0 ;
134         if ( z != 0.0 ) {
135             *phi = atan2 ( z, rxy ) ;
136         }
137         else {
138             *phi = 0.0 ;
139         }
140         *td = 0.0 ;
141         *pd = 0.0 ;
142     }
143     *r = rtrue ;
144     if ( rw != 0.0 ) {
145         *rd = ( xyp + z * zd ) / rw ;
146     }
147     else {
148         *rd = 0.0 ;
149     }
150 
151 /*
152  * Finished.
153  */
154 
155 }
156 
157 /*
158  *  gal - General Astrodynamics Library
159  *  Copyright (C) 2008 Paul C. L. Willmott
160  *
161  *  This program is free software; you can redistribute it and/or modify
162  *  it under the terms of the GNU General Public License as published by
163  *  the Free Software Foundation; either version 2 of the License, or
164  *  (at your option) any later version.
165  *
166  *  This program is distributed in the hope that it will be useful,
167  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
168  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
169  *  GNU General Public License for more details.
170  *
171  *  You should have received a copy of the GNU General Public License along
172  *  with this program; if not, write to the Free Software Foundation, Inc.,
173  *  51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
174  *
175  *  Contact:
176  *
177  *  Paul Willmott
178  *  vp9mu@amsat.org
179  */
180