1 /* dvhat.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 DVHAT ( Derivative and unit vector "V-hat" of a state) */
dvhat_(doublereal * s1,doublereal * sout)9 /* Subroutine */ int dvhat_(doublereal *s1, doublereal *sout)
10 {
11     /* System generated locals */
12     doublereal d__1;
13 
14     /* Local variables */
15     extern /* Subroutine */ int vequ_(doublereal *, doublereal *), vperp_(
16 	    doublereal *, doublereal *, doublereal *), unorm_(doublereal *,
17 	    doublereal *, doublereal *);
18     doublereal length;
19     extern /* Subroutine */ int vsclip_(doublereal *, doublereal *);
20 
21 /* $ Abstract */
22 
23 /*     Find the unit vector corresponding to a state vector and the */
24 /*     derivative of the unit vector. */
25 
26 /* $ Disclaimer */
27 
28 /*     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
29 /*     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
30 /*     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
31 /*     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
32 /*     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
33 /*     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
34 /*     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
35 /*     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
36 /*     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
37 /*     SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
38 
39 /*     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
40 /*     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
41 /*     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
42 /*     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
43 /*     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
44 /*     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
45 
46 /*     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
47 /*     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
48 /*     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
49 /*     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
50 
51 /* $ Required_Reading */
52 
53 /*     None. */
54 
55 /* $ Keywords */
56 
57 /*     VECTOR */
58 /*     DERIVATIVE */
59 /*     MATH */
60 
61 /* $ Declarations */
62 /* $ Brief_I/O */
63 
64 /*     VARIABLE  I/O  DESCRIPTION */
65 /*     --------  ---  -------------------------------------------------- */
66 /*     S1        I     State to be normalized. */
67 /*     SOUT      O     Unit vector S1 / |S1|, and its time derivative. */
68 
69 /* $ Detailed_Input */
70 
71 /*     S1       This is any double precision state. If the position */
72 /*              component of the state is the zero vector, this routine */
73 /*              will detect it and will not attempt to divide by zero. */
74 
75 /* $ Detailed_Output */
76 
77 /*     SOUT     SOUT is a state containing the unit vector pointing in */
78 /*              the direction of position component of S1 and the */
79 /*              derivative of the unit vector with respect to time. */
80 
81 /*              SOUT may overwrite S1. */
82 
83 /* $ Parameters */
84 
85 /*     None. */
86 
87 /* $ Exceptions */
88 
89 /*     Error free. */
90 
91 /*     1) If S1 represents the zero vector, then the position */
92 /*        component of SOUT will also be the zero vector.  The */
93 /*        velocity component will be the velocity component */
94 /*        of S1. */
95 
96 /* $ Files */
97 
98 /*     None. */
99 
100 /* $ Particulars */
101 
102 /*     Let S1 be a state vector with position and velocity components P */
103 /*     and V respectively.  From these components one can compute the */
104 /*     unit vector parallel to P, call it U and the derivative of U */
105 /*     with respect to time, DU.  This pair (U,DU) is the state returned */
106 /*     by this routine in SOUT. */
107 
108 /* $ Examples */
109 
110 /*   Any numerical results shown for this example may differ between */
111 /*   platforms as the results depend on the SPICE kernels used as input */
112 /*   and the machine specific arithmetic implementation. */
113 
114 /*   Suppose that STATE gives the apparent state of a body with */
115 /*   respect to an observer.  This routine can be used to compute the */
116 /*   instantaneous angular rate of the object across the sky as seen */
117 /*   from the observers vantage. */
118 
119 /*           PROGRAM DVHAT_T */
120 /*           IMPLICIT NONE */
121 
122 /*           DOUBLE PRECISION      ET */
123 /*           DOUBLE PRECISION      LT */
124 /*           DOUBLE PRECISION      OMEGA */
125 /*           DOUBLE PRECISION      STATE  (6) */
126 /*           DOUBLE PRECISION      USTATE (6) */
127 
128 /*           DOUBLE PRECISION      VNORM */
129 
130 /*           CHARACTER*(32)        EPOCH */
131 /*           CHARACTER*(32)        TARGET */
132 /*           CHARACTER*(32)        FRAME */
133 /*           CHARACTER*(32)        ABCORR */
134 /*           CHARACTER*(32)        OBSRVR */
135 
136 /*     C */
137 /*     C     Load SPK, PCK, and LSK kernels, use a meta kernel for */
138 /*     C     convenience. */
139 /*     C */
140 /*           CALL FURNSH ( 'standard.tm' ) */
141 
142 /*     C */
143 /*     C     Define an arbitrary epoch, convert the epoch to ephemeris */
144 /*     C     time. */
145 /*     C */
146 /*           EPOCH = 'Jan 1 2009' */
147 /*           CALL STR2ET ( EPOCH, ET ) */
148 
149 /*     C */
150 /*     C     Calculate the state of the moon with respect to the */
151 /*     C     earth-moon barycenter in J2000, corrected for light time */
152 /*     C     and stellar aberration at ET. */
153 /*     C */
154 /*           TARGET   = 'MOON' */
155 /*           FRAME    = 'J2000' */
156 /*           ABCORR   = 'LT+S' */
157 /*           OBSRVR   = 'EARTH BARYCENTER' */
158 
159 /*           CALL SPKEZR ( TARGET, ET, FRAME, ABCORR, OBSRVR, STATE, LT ) */
160 
161 /*     C */
162 /*     C     Calculate the unit vector of STATE and the derivative of the */
163 /*     C     unit vector. */
164 /*     C */
165 /*           CALL DVHAT ( STATE, USTATE ) */
166 
167 /*     C */
168 /*     C     Calculate the instantaneous angular velocity from the */
169 /*     C     magnitude of the derivative of the unit vector. */
170 /*     C */
171 /*     C          v = r x omega */
172 /*     C */
173 /*     C          ||omega|| = ||v||  for  r . v = 0 */
174 /*     C                      ----- */
175 /*     C                      ||r|| */
176 /*     C */
177 /*     C          ||omega|| = ||v||  for  ||r|| = 1 */
178 /*     C */
179 /*           OMEGA = VNORM( USTATE(4) ) */
180 
181 /*           WRITE(*,*) 'Instantaneous angular velocity, rad/sec', OMEGA */
182 
183 /*           END */
184 
185 /*   The program outputs: */
186 
187 /*       Instantaneous angular velocity, rad/sec  2.48106658E-06 */
188 
189 /* $ Restrictions */
190 
191 /*     None. */
192 
193 /* $ Literature_References */
194 
195 /*     None. */
196 
197 /* $ Author_and_Institution */
198 
199 /*     N.J. Bachman    (JPL) */
200 /*     W.L. Taber      (JPL) */
201 
202 /* $ Version */
203 
204 /* -    SPICELIB Version 1.1.1, 06-MAY-2010 (EDW) */
205 
206 /*        Expanded the code example into a complete program. */
207 
208 /*        Reordered header sections to proper NAIF convention. */
209 /*        Removed Revision section, it listed a duplication of a */
210 /*        Version section entry. */
211 
212 /* -    SPICELIB Version 1.1.0, 02-SEP-2005 (NJB) */
213 
214 /*        Updated to remove non-standard use of duplicate arguments */
215 /*        in VPERP and VSCL calls. */
216 
217 /* -    SPICELIB Version 1.0.0, 15-JUN-1995 (WLT) */
218 
219 /* -& */
220 /* $ Index_Entries */
221 
222 /*     State of a unit vector parallel to a state vector */
223 
224 /* -& */
225 
226 /*     Get the position portion of the output state and the length of */
227 /*     the input position. */
228 
229     unorm_(s1, sout, &length);
230     if (length == 0.) {
231 
232 /*        If the length of the input position is zero, just copy */
233 /*        the input velocity to the output velocity. */
234 
235 	vequ_(&s1[3], &sout[3]);
236     } else {
237 
238 /*        Otherwise the derivative of the unit vector is just the */
239 /*        component of the input velocity perpendicular to the input */
240 /*        position, scaled by the reciprocal of the length of the */
241 /*        input position. */
242 
243 	vperp_(&s1[3], sout, &sout[3]);
244 	d__1 = 1. / length;
245 	vsclip_(&d__1, &sout[3]);
246     }
247     return 0;
248 } /* dvhat_ */
249 
250