1 /* vsep.f -- translated by f2c (version 19980913).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "f2c.h"
7
8 /* $Procedure VSEP ( Angular separation of vectors, 3 dimensions ) */
vsep_(doublereal * v1,doublereal * v2)9 doublereal vsep_(doublereal *v1, doublereal *v2)
10 {
11 /* System generated locals */
12 doublereal ret_val;
13
14 /* Builtin functions */
15 double asin(doublereal);
16
17 /* Local variables */
18 extern doublereal vdot_(doublereal *, doublereal *);
19 doublereal dmag1, dmag2, vtemp[3];
20 extern /* Subroutine */ int unorm_(doublereal *, doublereal *, doublereal
21 *);
22 extern doublereal vnorm_(doublereal *);
23 doublereal u1[3], u2[3];
24 extern doublereal pi_(void);
25
26 /* $ Abstract */
27
28 /* Find the separation angle in radians between two double */
29 /* precision, 3-dimensional vectors. This angle is defined as zero */
30 /* if either vector is zero. */
31
32 /* $ Disclaimer */
33
34 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
35 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
36 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
37 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
38 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
39 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
40 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
41 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
42 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
43 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
44
45 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
46 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
47 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
48 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
49 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
50 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
51
52 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
53 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
54 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
55 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
56
57 /* $ Required_Reading */
58
59 /* None. */
60
61 /* $ Keywords */
62
63 /* ANGLE, VECTOR */
64
65 /* $ Declarations */
66 /* $ Brief_I/O */
67
68 /* VARIABLE I/O DESCRIPTION */
69 /* -------- --- -------------------------------------------------- */
70 /* V1 I First vector. */
71 /* V2 I Second vector. */
72
73
74 /* $ Detailed_Input */
75
76 /* V1 is an arbitrary double precision, 3-dimensional vector. */
77 /* V2 is also an arbitrary double precision, 3-dimensional */
78 /* vector. V1 or V2 or both may be the zero vector. */
79
80 /* $ Detailed_Output */
81
82 /* VSEP is the angle between V1 and V2 expressed in radians. */
83 /* VSEP is strictly non-negative. If either V1 or V2 is */
84 /* the zero vector, then VSEP is defined to be 0 radians. */
85
86 /* $ Parameters */
87
88 /* None. */
89
90 /* $ Particulars */
91
92 /* In the plane, it is a simple matter to calculate the angle */
93 /* between two vectors once the two vectors have been made to be */
94 /* unit length. Then, since the two vectors form the two equal */
95 /* sides of an isosceles triangle, the length of the third side */
96 /* is given by the expression */
97
98 /* LENGTH = 2.0 * SINE ( VSEP/2.0 ) */
99
100 /* The length is given by the magnitude of the difference of the */
101 /* two unit vectors */
102
103 /* LENGTH = NORM ( U1 - U2 ) */
104
105 /* Once the length is found, the value of VSEP may be calculated */
106 /* by inverting the first expression given above as */
107
108 /* VSEP = 2.0 * ARCSINE ( LENGTH/2.0 ) */
109
110 /* This expression becomes increasingly unstable when VSEP gets */
111 /* larger than PI/2 or 90 degrees. In this situation (which is */
112 /* easily detected by determining the sign of the dot product of */
113 /* V1 and V2) the supplementary angle is calculated first and */
114 /* then VSEP is given by */
115
116 /* VSEP = PI - SUPPLEMENTARY_ANGLE */
117
118 /* $ Examples */
119
120 /* The following table gives sample values for V1, V2 and VSEP */
121 /* implied by the inputs. */
122
123 /* V1 V2 VSEP */
124 /* ---------------------------------------------------------------- */
125 /* (1, 0, 0) (1, 0, 0) 0.0D0 */
126 /* (1, 0, 0) (0, 1, 0) PI/2 (=1.571...) */
127
128
129 /* $ Restrictions */
130
131 /* The user is required to insure that the input vectors will not */
132 /* cause floating point overflow upon calculation of the vector */
133 /* dot product since no error detection or correction code is */
134 /* implemented. In practice, this is not a significant */
135 /* restriction. */
136
137 /* $ Exceptions */
138
139 /* Error free. */
140
141 /* $ Files */
142
143 /* None */
144
145 /* $ Author_and_Institution */
146
147 /* K.R. Gehringer (JPL) */
148 /* W.M. Owen (JPL) */
149 /* W.L. Taber (JPL) */
150
151 /* $ Literature_References */
152
153 /* None */
154
155 /* $ Version */
156
157 /* - SPICELIB Version 1.1.1, 17-APR-2006 (EDW) */
158
159 /* Typo correction to the value of PI/2 in the Examples */
160 /* section, 1.571 instead of 1.71. */
161
162 /* - SPICELIB Version 1.1.0, 29-FEB-1996 (KRG) */
163
164 /* The declaration for the SPICELIB function PI is now */
165 /* preceded by an EXTERNAL statement declaring PI to be an */
166 /* external function. This removes a conflict with any */
167 /* compilers that have a PI intrinsic function. */
168
169 /* - SPICELIB Version 1.0.1, 10-MAR-1992 (WLT) */
170
171 /* Comment section for permuted index source lines was added */
172 /* following the header. */
173
174 /* - SPICELIB Version 1.0.0, 31-JAN-1990 (WMO) */
175
176 /* -& */
177 /* $ Index_Entries */
178
179 /* angular separation of 3-dimensional vectors */
180
181 /* -& */
182 /* $ Revisions */
183
184 /* - SPICELIB Version 1.1.0, 29-FEB-1996 (KRG) */
185
186 /* The declaration for the SPICELIB function PI is now */
187 /* preceded by an EXTERNAL statement declaring PI to be an */
188 /* external function. This removes a conflict with any */
189 /* compilers that have a PI intrinsic function. */
190
191 /* - SPICELIB Version 1.0.1, 10-MAR-1992 (WLT) */
192
193 /* Comment section for permuted index source lines was added */
194 /* following the header. */
195
196 /* - SPICELIB Version 1.0.0, 31-JAN-1990 (WMO) */
197
198 /* - Beta Version 1.0.1, 10-JAN-1989 (WLT) */
199
200 /* Error free specification added. */
201
202 /* -& */
203
204 /* SPICELIB functions */
205
206
207 /* Local Variables */
208
209 /* The following declarations represent, respectively: */
210 /* Magnitudes of V1, V2 */
211 /* Either of the difference vectors: V1-V2 or V1-(-V2) */
212 /* Unit vectors parallel to V1 and V2 */
213
214
215 /* Calculate the magnitudes of V1 and V2; if either is 0, VSEP = 0 */
216
217 unorm_(v1, u1, &dmag1);
218 if (dmag1 == 0.) {
219 ret_val = 0.;
220 return ret_val;
221 }
222 unorm_(v2, u2, &dmag2);
223 if (dmag2 == 0.) {
224 ret_val = 0.;
225 return ret_val;
226 }
227 if (vdot_(u1, u2) > 0.) {
228 vtemp[0] = u1[0] - u2[0];
229 vtemp[1] = u1[1] - u2[1];
230 vtemp[2] = u1[2] - u2[2];
231 ret_val = asin(vnorm_(vtemp) * .5) * 2.;
232 } else if (vdot_(u1, u2) < 0.) {
233 vtemp[0] = u1[0] + u2[0];
234 vtemp[1] = u1[1] + u2[1];
235 vtemp[2] = u1[2] + u2[2];
236 ret_val = pi_() - asin(vnorm_(vtemp) * .5) * 2.;
237 } else {
238 ret_val = pi_() / 2.;
239 }
240 return ret_val;
241 } /* vsep_ */
242
243