1 /* zzspklt0.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 /* Table of constant values */
9
10 static integer c__0 = 0;
11 static integer c__6 = 6;
12 static doublereal c_b19 = -1.;
13
14 /* $Procedure ZZSPKLT0 ( S/P Kernel, light time corrected state ) */
zzspklt0_(integer * targ,doublereal * et,char * ref,char * abcorr,doublereal * stobs,doublereal * starg,doublereal * lt,doublereal * dlt,ftnlen ref_len,ftnlen abcorr_len)15 /* Subroutine */ int zzspklt0_(integer *targ, doublereal *et, char *ref, char
16 *abcorr, doublereal *stobs, doublereal *starg, doublereal *lt,
17 doublereal *dlt, ftnlen ref_len, ftnlen abcorr_len)
18 {
19 /* Initialized data */
20
21 static logical pass1 = TRUE_;
22 static char prvcor[5] = " ";
23
24 /* System generated locals */
25 doublereal d__1, d__2, d__3, d__4;
26
27 /* Builtin functions */
28 integer s_cmp(char *, char *, ftnlen, ftnlen);
29 /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
30
31 /* Local variables */
32 doublereal dist;
33 extern /* Subroutine */ int zzspkgo0_(integer *, doublereal *, char *,
34 integer *, doublereal *, doublereal *, ftnlen);
35 extern doublereal vdot_(doublereal *, doublereal *);
36 static logical xmit;
37 extern /* Subroutine */ int zzvalcor_(char *, logical *, ftnlen);
38 doublereal a, b, c__;
39 integer i__, refid;
40 extern /* Subroutine */ int chkin_(char *, ftnlen);
41 doublereal epoch;
42 extern /* Subroutine */ int errch_(char *, char *, ftnlen, ftnlen);
43 static logical usecn;
44 extern /* Subroutine */ int vlcom_(doublereal *, doublereal *, doublereal
45 *, doublereal *, doublereal *), vsubg_(doublereal *, doublereal *,
46 integer *, doublereal *);
47 doublereal ssblt, lterr;
48 static logical uselt;
49 extern doublereal vnorm_(doublereal *);
50 doublereal prvlt;
51 extern logical failed_(void);
52 extern doublereal clight_(void);
53 logical attblk[15];
54 extern doublereal touchd_(doublereal *);
55 extern /* Subroutine */ int sigerr_(char *, ftnlen), chkout_(char *,
56 ftnlen);
57 integer ltsign;
58 extern /* Subroutine */ int irfnum_(char *, integer *, ftnlen), setmsg_(
59 char *, ftnlen);
60 doublereal ssbtrg[6];
61 integer numitr;
62 extern logical return_(void);
63 logical usestl;
64
65 /* $ Abstract */
66
67 /* Return the state (position and velocity) of a target body */
68 /* relative to an observer, optionally corrected for light time, */
69 /* expressed relative to an inertial reference frame. */
70
71 /* $ Disclaimer */
72
73 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
74 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
75 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
76 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
77 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
78 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
79 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
80 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
81 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
82 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
83
84 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
85 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
86 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
87 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
88 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
89 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
90
91 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
92 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
93 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
94 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
95
96 /* $ Required_Reading */
97
98 /* SPK */
99
100 /* $ Keywords */
101
102 /* EPHEMERIS */
103
104 /* $ Declarations */
105 /* $ Abstract */
106
107 /* Include file zzabcorr.inc */
108
109 /* SPICE private file intended solely for the support of SPICE */
110 /* routines. Users should not include this file directly due */
111 /* to the volatile nature of this file */
112
113 /* The parameters below define the structure of an aberration */
114 /* correction attribute block. */
115
116 /* $ Disclaimer */
117
118 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
119 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
120 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
121 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
122 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
123 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
124 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
125 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
126 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
127 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
128
129 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
130 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
131 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
132 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
133 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
134 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
135
136 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
137 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
138 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
139 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
140
141 /* $ Parameters */
142
143 /* An aberration correction attribute block is an array of logical */
144 /* flags indicating the attributes of the aberration correction */
145 /* specified by an aberration correction string. The attributes */
146 /* are: */
147
148 /* - Is the correction "geometric"? */
149
150 /* - Is light time correction indicated? */
151
152 /* - Is stellar aberration correction indicated? */
153
154 /* - Is the light time correction of the "converged */
155 /* Newtonian" variety? */
156
157 /* - Is the correction for the transmission case? */
158
159 /* - Is the correction relativistic? */
160
161 /* The parameters defining the structure of the block are as */
162 /* follows: */
163
164 /* NABCOR Number of aberration correction choices. */
165
166 /* ABATSZ Number of elements in the aberration correction */
167 /* block. */
168
169 /* GEOIDX Index in block of geometric correction flag. */
170
171 /* LTIDX Index of light time flag. */
172
173 /* STLIDX Index of stellar aberration flag. */
174
175 /* CNVIDX Index of converged Newtonian flag. */
176
177 /* XMTIDX Index of transmission flag. */
178
179 /* RELIDX Index of relativistic flag. */
180
181 /* The following parameter is not required to define the block */
182 /* structure, but it is convenient to include it here: */
183
184 /* CORLEN The maximum string length required by any aberration */
185 /* correction string */
186
187 /* $ Author_and_Institution */
188
189 /* N.J. Bachman (JPL) */
190
191 /* $ Literature_References */
192
193 /* None. */
194
195 /* $ Version */
196
197 /* - SPICELIB Version 1.0.0, 18-DEC-2004 (NJB) */
198
199 /* -& */
200 /* Number of aberration correction choices: */
201
202
203 /* Aberration correction attribute block size */
204 /* (number of aberration correction attributes): */
205
206
207 /* Indices of attributes within an aberration correction */
208 /* attribute block: */
209
210
211 /* Maximum length of an aberration correction string: */
212
213
214 /* End of include file zzabcorr.inc */
215
216 /* $ Brief_I/O */
217
218 /* Variable I/O Description */
219 /* -------- --- -------------------------------------------------- */
220 /* TARG I Target body. */
221 /* ET I Observer epoch. */
222 /* REF I Inertial reference frame of output state. */
223 /* ABCORR I Aberration correction flag. */
224 /* STOBS I State of the observer relative to the SSB. */
225 /* STARG O State of target. */
226 /* LT O One way light time between observer and target. */
227 /* DLT O Derivative of light time with respect to time. */
228
229 /* $ Detailed_Input */
230
231 /* TARG is the NAIF ID code for a target body. The target */
232 /* and observer define a state vector whose position */
233 /* component points from the observer to the target. */
234
235 /* ET is the ephemeris time, expressed as seconds past */
236 /* J2000 TDB, at which the state of the target body */
237 /* relative to the observer is to be computed. ET */
238 /* refers to time at the observer's location. */
239
240 /* REF is the inertial reference frame with respect to which */
241 /* the input state STOBS and the output state STARG are */
242 /* expressed. REF must be recognized by the SPICE */
243 /* Toolkit. The acceptable frames are listed in the */
244 /* Frames Required Reading, as well as in the SPICELIB */
245 /* routine CHGIRF. */
246
247 /* Case and blanks are not significant in the string */
248 /* REF. */
249
250
251 /* ABCORR indicates the aberration corrections to be applied to */
252 /* the state of the target body to account for one-way */
253 /* light time. See the discussion in the Particulars */
254 /* section for recommendations on how to choose */
255 /* aberration corrections. */
256
257 /* If ABCORR includes the stellar aberration correction */
258 /* symbol '+S', this flag is simply ignored. Aside from */
259 /* the possible presence of this symbol, ABCORR may be */
260 /* any of the following: */
261
262 /* 'NONE' Apply no correction. Return the */
263 /* geometric state of the target body */
264 /* relative to the observer. */
265
266 /* The following values of ABCORR apply to the */
267 /* "reception" case in which photons depart from the */
268 /* target's location at the light-time corrected epoch */
269 /* ET-LT and *arrive* at the observer's location at ET: */
270
271 /* 'LT' Correct for one-way light time (also */
272 /* called "planetary aberration") using a */
273 /* Newtonian formulation. This correction */
274 /* yields the state of the target at the */
275 /* moment it emitted photons arriving at */
276 /* the observer at ET. */
277
278 /* The light time correction involves */
279 /* iterative solution of the light time */
280 /* equation (see Particulars for details). */
281 /* The solution invoked by the 'LT' option */
282 /* uses one iteration. */
283
284 /* 'CN' Converged Newtonian light time */
285 /* correction. In solving the light time */
286 /* equation, the 'CN' correction iterates */
287 /* until the solution converges (three */
288 /* iterations on all supported platforms). */
289 /* Whether the 'CN+S' solution is */
290 /* substantially more accurate than the */
291 /* 'LT' solution depends on the geometry */
292 /* of the participating objects and on the */
293 /* accuracy of the input data. In all */
294 /* cases this routine will execute more */
295 /* slowly when a converged solution is */
296 /* computed. See the Particulars section of */
297 /* SPKEZR for a discussion of precision of */
298 /* light time corrections. */
299
300 /* The following values of ABCORR apply to the */
301 /* "transmission" case in which photons *depart* from */
302 /* the observer's location at ET and arrive at the */
303 /* target's location at the light-time corrected epoch */
304 /* ET+LT: */
305
306 /* 'XLT' "Transmission" case: correct for */
307 /* one-way light time using a Newtonian */
308 /* formulation. This correction yields the */
309 /* state of the target at the moment it */
310 /* receives photons emitted from the */
311 /* observer's location at ET. */
312
313 /* 'XCN' "Transmission" case: converged */
314 /* Newtonian light time correction. */
315
316
317 /* Neither special nor general relativistic effects are */
318 /* accounted for in the aberration corrections applied */
319 /* by this routine. */
320
321 /* Case and blanks are not significant in the string */
322 /* ABCORR. */
323
324
325 /* STOBS is the geometric (uncorrected) state of the observer */
326 /* relative to the solar system barycenter at epoch ET. */
327 /* STOBS is a 6-vector: the first three components of */
328 /* STOBS represent a Cartesian position vector; the last */
329 /* three components represent the corresponding velocity */
330 /* vector. STOBS is expressed relative to the inertial */
331 /* reference frame designated by REF. */
332
333 /* Units are always km and km/sec. */
334
335 /* $ Detailed_Output */
336
337 /* STARG is a Cartesian state vector representing the position */
338 /* and velocity of the target body relative to the */
339 /* specified observer. STARG is corrected for the */
340 /* specified aberration, and is expressed with respect */
341 /* to the specified inertial reference frame. The first */
342 /* three components of STARG represent the x-, y- and */
343 /* z-components of the target's position; last three */
344 /* components form the corresponding velocity vector. */
345
346 /* The position component of STARG points from the */
347 /* observer's location at ET to the aberration-corrected */
348 /* location of the target. Note that the sense of the */
349 /* position vector is independent of the direction of */
350 /* radiation travel implied by the aberration */
351 /* correction. */
352
353 /* Units are always km and km/sec. */
354
355 /* LT is the one-way light time between the observer and */
356 /* target in seconds. If the target state is corrected */
357 /* for light time, then LT is the one-way light time */
358 /* between the observer and the light time-corrected */
359 /* target location. */
360
361 /* DLT is the derivative with respect to barycentric */
362 /* dynamical time of the one way light time between */
363 /* target and observer: */
364
365 /* DLT = d(LT)/d(ET) */
366
367 /* DLT can also be described as the rate of change of */
368 /* one way light time. DLT is unitless, since LT and */
369 /* ET both have units of TDB seconds. */
370
371 /* If the observer and target are at the same position, */
372 /* then DLT is set to zero. */
373
374 /* $ Parameters */
375
376 /* None. */
377
378 /* $ Exceptions */
379
380 /* 1) For the convenience of the caller, the input aberration */
381 /* correction flag can call for stellar aberration correction via */
382 /* inclusion of the '+S' suffix. This portion of the aberration */
383 /* correction flag is ignored if present. */
384
385 /* 2) If the value of ABCORR is not recognized, the error */
386 /* is diagnosed by a routine in the call tree of this */
387 /* routine. */
388
389 /* 3) If the reference frame requested is not a recognized */
390 /* inertial reference frame, the error SPICE(BADFRAME) */
391 /* is signaled. */
392
393 /* 4) If the state of the target relative to the solar system */
394 /* barycenter cannot be computed, the error will be diagnosed */
395 /* by routines in the call tree of this routine. */
396
397 /* 5) If the observer and target are at the same position, */
398 /* then DLT is set to zero. This situation could arise, */
399 /* for example, when the observer is Mars and the target */
400 /* is the Mars barycenter. */
401
402 /* 6) If a division by zero error would occur in the computation */
403 /* of DLT, the error SPICE(DIVIDEBYZERO) is signaled. */
404
405 /* $ Files */
406
407 /* This routine computes states using SPK files that have been */
408 /* loaded into the SPICE system, normally via the kernel loading */
409 /* interface routine FURNSH. Application programs typically load */
410 /* kernels once before this routine is called, for example during */
411 /* program initialization; kernels need not be loaded repeatedly. */
412 /* See the routine FURNSH and the SPK and KERNEL Required Reading */
413 /* for further information on loading (and unloading) kernels. */
414
415 /* If any of the ephemeris data used to compute STARG are expressed */
416 /* relative to a non-inertial frame in the SPK files providing those */
417 /* data, additional kernels may be needed to enable the reference */
418 /* frame transformations required to compute the state. Normally */
419 /* these additional kernels are PCK files or frame kernels. Any */
420 /* such kernels must already be loaded at the time this routine is */
421 /* called. */
422
423 /* $ Particulars */
424
425 /* This routine supports higher-level SPK API routines that can */
426 /* perform both light time and stellar aberration corrections. */
427 /* User applications normally will not need to call this routine */
428 /* directly. */
429
430 /* See the header of the routine SPKEZR for a detailed discussion */
431 /* of aberration corrections. */
432
433 /* $ Examples */
434
435 /* The numerical results shown for this example may differ across */
436 /* platforms. The results depend on the SPICE kernels used as */
437 /* input, the compiler and supporting libraries, and the machine */
438 /* specific arithmetic implementation. */
439
440 /* 1) Look up a sequence of states of the Moon as seen from the */
441 /* Earth. Use light time corrections. Compute the first state for */
442 /* the epoch 2000 JAN 1 12:00:00 TDB; compute subsequent states at */
443 /* intervals of 1 hour. For each epoch, display the states, the */
444 /* one way light time between target and observer, and the rate of */
445 /* change of the one way light time. */
446
447 /* Use the following meta-kernel to specify the kernels to */
448 /* load: */
449
450 /* KPL/MK */
451
452 /* File name: spkltc.tm */
453
454 /* This meta-kernel is intended to support operation of SPICE */
455 /* example programs. The kernels shown here should not be */
456 /* assumed to contain adequate or correct versions of data */
457 /* required by SPICE-based user applications. */
458
459 /* In order for an application to use this meta-kernel, the */
460 /* kernels referenced here must be present in the user's */
461 /* current working directory. */
462
463
464 /* \begindata */
465
466 /* KERNELS_TO_LOAD = ( 'de421.bsp', */
467 /* 'pck00010.tpc', */
468 /* 'naif0010.tls' ) */
469
470 /* \begintext */
471
472
473 /* The code example follows: */
474
475 /* PROGRAM EX1 */
476 /* IMPLICIT NONE */
477 /* C */
478 /* C Local constants */
479 /* C */
480 /* C The meta-kernel name shown here refers to a file whose */
481 /* C contents are those shown above. This file and the kernels */
482 /* C it references must exist in your current working directory. */
483 /* C */
484 /* CHARACTER*(*) META */
485 /* PARAMETER ( META = 'spkltc.tm' ) */
486 /* C */
487 /* C Use a time step of 1 hour; look up 5 states. */
488 /* C */
489 /* DOUBLE PRECISION STEP */
490 /* PARAMETER ( STEP = 3600.0D0 ) */
491
492 /* INTEGER MAXITR */
493 /* PARAMETER ( MAXITR = 5 ) */
494 /* C */
495 /* C Local variables */
496 /* C */
497 /* DOUBLE PRECISION DLT */
498 /* DOUBLE PRECISION ET */
499 /* DOUBLE PRECISION ET0 */
500 /* DOUBLE PRECISION LT */
501 /* DOUBLE PRECISION STATE ( 6 ) */
502 /* DOUBLE PRECISION STOBS ( 6 ) */
503 /* INTEGER I */
504
505 /* C */
506 /* C Load the SPK and LSK kernels via the meta-kernel. */
507 /* C */
508 /* CALL FURNSH ( META ) */
509 /* C */
510 /* C Convert the start time to seconds past J2000 TDB. */
511 /* C */
512 /* CALL STR2ET ( '2000 JAN 1 12:00:00 TDB', ET0 ) */
513 /* C */
514 /* C Step through a series of epochs, looking up a */
515 /* C state vector at each one. */
516 /* C */
517 /* DO I = 1, MAXITR */
518
519 /* ET = ET0 + (I-1)*STEP */
520
521 /* C */
522 /* C Look up a state vector at epoch ET using the */
523 /* C following inputs: */
524 /* C */
525 /* C Target: Moon (NAIF ID code 301) */
526 /* C Reference frame: J2000 */
527 /* C Aberration correction: Light time ('LT') */
528 /* C Observer: Earth (NAIF ID code 399) */
529 /* C */
530 /* C Before we can execute this computation, we'll need the */
531 /* C geometric state of the observer relative to the solar */
532 /* C system barycenter at ET, expressed relative to the */
533 /* C J2000 reference frame: */
534 /* C */
535 /* CALL SPKSSB ( 399, ET, 'J2000', STOBS ) */
536 /* C */
537 /* C Now compute the desired state vector: */
538 /* C */
539 /* CALL SPKLTC ( 301, ET, 'J2000', 'LT', */
540 /* . STOBS, STATE, LT, DLT ) */
541
542 /* WRITE (*,*) 'ET = ', ET */
543 /* WRITE (*,*) 'J2000 x-position (km): ', STATE(1) */
544 /* WRITE (*,*) 'J2000 y-position (km): ', STATE(2) */
545 /* WRITE (*,*) 'J2000 z-position (km): ', STATE(3) */
546 /* WRITE (*,*) 'J2000 x-velocity (km/s): ', STATE(4) */
547 /* WRITE (*,*) 'J2000 y-velocity (km/s): ', STATE(5) */
548 /* WRITE (*,*) 'J2000 z-velocity (km/s): ', STATE(6) */
549 /* WRITE (*,*) 'One-way light time (s): ', LT */
550 /* WRITE (*,*) 'Light time rate: ', DLT */
551 /* WRITE (*,*) ' ' */
552
553 /* END DO */
554
555 /* END */
556
557
558 /* On a PC/Linux/gfortran platform, the following output was */
559 /* produced: */
560
561
562 /* ET = 0.0000000000000000 */
563 /* J2000 x-position (km): -291569.26541282982 */
564 /* J2000 y-position (km): -266709.18647825718 */
565 /* J2000 z-position (km): -76099.155118763447 */
566 /* J2000 x-velocity (km/s): 0.64353061322177041 */
567 /* J2000 y-velocity (km/s): -0.66608181700820079 */
568 /* J2000 z-velocity (km/s): -0.30132283179625752 */
569 /* One-way light time (s): 1.3423106103251679 */
570 /* Light time rate: 1.07316908698977495E-007 */
571
572 /* ET = 3600.0000000000000 */
573 /* J2000 x-position (km): -289240.78128184378 */
574 /* J2000 y-position (km): -269096.44087958336 */
575 /* J2000 z-position (km): -77180.899725757539 */
576 /* J2000 x-velocity (km/s): 0.65006211520087476 */
577 /* J2000 y-velocity (km/s): -0.66016273921695667 */
578 /* J2000 z-velocity (km/s): -0.29964267390571342 */
579 /* One-way light time (s): 1.3426939548635302 */
580 /* Light time rate: 1.05652598952224259E-007 */
581
582 /* ET = 7200.0000000000000 */
583 /* J2000 x-position (km): -286888.88736709207 */
584 /* J2000 y-position (km): -271462.30170547962 */
585 /* J2000 z-position (km): -78256.555682137609 */
586 /* J2000 x-velocity (km/s): 0.65653599154284592 */
587 /* J2000 y-velocity (km/s): -0.65419657680401588 */
588 /* J2000 z-velocity (km/s): -0.29794027307420823 */
589 /* One-way light time (s): 1.3430713117337547 */
590 /* Light time rate: 1.03990456898758609E-007 */
591
592 /* ET = 10800.000000000000 */
593 /* J2000 x-position (km): -284513.79173691198 */
594 /* J2000 y-position (km): -273806.60031034052 */
595 /* J2000 z-position (km): -79326.043183274567 */
596 /* J2000 x-velocity (km/s): 0.66295190054599118 */
597 /* J2000 y-velocity (km/s): -0.64818380709706158 */
598 /* J2000 z-velocity (km/s): -0.29621577937090349 */
599 /* One-way light time (s): 1.3434426890693671 */
600 /* Light time rate: 1.02330665243423737E-007 */
601
602 /* ET = 14400.000000000000 */
603 /* J2000 x-position (km): -282115.70368389413 */
604 /* J2000 y-position (km): -276129.16976799071 */
605 /* J2000 z-position (km): -80389.282965712249 */
606 /* J2000 x-velocity (km/s): 0.66930950377548726 */
607 /* J2000 y-velocity (km/s): -0.64212490805688027 */
608 /* J2000 z-velocity (km/s): -0.29446934336246899 */
609 /* One-way light time (s): 1.3438080956559786 */
610 /* Light time rate: 1.00673403630050830E-007 */
611
612
613 /* $ Restrictions */
614
615 /* 1) The routine SPKGEO should be used instead of this routine */
616 /* to compute geometric states. SPKGEO introduces less */
617 /* round-off error when the observer and target have common */
618 /* center that is closer to both objects than is the solar */
619 /* system barycenter. */
620
621 /* 2) The kernel files to be used by SPKLTC must be loaded */
622 /* (normally by the SPICELIB kernel loader FURNSH) before */
623 /* this routine is called. */
624
625 /* 3) Unlike most other SPK state computation routines, this */
626 /* routine requires that the output state be relative to an */
627 /* inertial reference frame. */
628
629 /* $ Literature_References */
630
631 /* None. */
632
633 /* $ Author_and_Institution */
634
635 /* N.J. Bachman (JPL) */
636
637 /* $ Version */
638
639 /* - SPICELIB Version 2.0.0, 04-JUL-2014 (NJB) */
640
641 /* Discussion of light time corrections was updated. Assertions */
642 /* that converged light time corrections are unlikely to be */
643 /* useful were removed. */
644
645 /* Last update was 02-MAY-2012 (NJB) */
646
647 /* Updated to ensure convergence when CN or XCN light time */
648 /* corrections are used. The new algorithm also terminates early */
649 /* (after fewer than three iterations) when convergence is */
650 /* attained. */
651
652 /* Call to ZZPRSCOR was replaced by a call to ZZVALCOR. */
653
654 /* - SPICELIB Version 1.0.0, 11-JAN-2008 (NJB) */
655
656 /* -& */
657 /* $ Index_Entries */
658
659 /* low-level light time correction */
660 /* light-time corrected state from spk file */
661 /* get light-time corrected state */
662
663 /* -& */
664 /* $ Revisions */
665
666 /* None. */
667
668 /* -& */
669
670 /* SPICELIB functions */
671
672
673 /* Local parameters */
674
675
676 /* TOL is the tolerance used for a division-by-zero test */
677 /* performed prior to computation of DLT. */
678
679
680 /* Convergence limit: */
681
682
683 /* Maximum number of light time iterations for any */
684 /* aberration correction: */
685
686
687 /* Local variables */
688
689
690 /* Saved variables */
691
692
693 /* Initial values */
694
695
696 /* Standard SPICE error handling. */
697
698 if (return_()) {
699 return 0;
700 } else {
701 chkin_("ZZSPKLT0", (ftnlen)8);
702 }
703 if (pass1 || s_cmp(abcorr, prvcor, abcorr_len, (ftnlen)5) != 0) {
704
705 /* The aberration correction flag differs from the value it */
706 /* had on the previous call, if any. Analyze the new flag. */
707
708 zzvalcor_(abcorr, attblk, abcorr_len);
709 if (failed_()) {
710 chkout_("ZZSPKLT0", (ftnlen)8);
711 return 0;
712 }
713
714 /* The aberration correction flag is recognized; save it. */
715
716 s_copy(prvcor, abcorr, (ftnlen)5, abcorr_len);
717
718 /* Set logical flags indicating the attributes of the requested */
719 /* correction: */
720
721 /* XMIT is .TRUE. when the correction is for transmitted */
722 /* radiation. */
723
724 /* USELT is .TRUE. when any type of light time correction */
725 /* (normal or converged Newtonian) is specified. */
726
727 /* USECN indicates converged Newtonian light time correction. */
728
729 /* The above definitions are consistent with those used by */
730 /* ZZVALCOR. */
731
732 xmit = attblk[4];
733 uselt = attblk[1];
734 usecn = attblk[3];
735 usestl = attblk[2];
736 pass1 = FALSE_;
737 }
738
739 /* See if the reference frame is a recognized inertial frame. */
740
741 irfnum_(ref, &refid, ref_len);
742 if (refid == 0) {
743 setmsg_("The requested frame '#' is not a recognized inertial frame. "
744 , (ftnlen)60);
745 errch_("#", ref, (ftnlen)1, ref_len);
746 sigerr_("SPICE(BADFRAME)", (ftnlen)15);
747 chkout_("ZZSPKLT0", (ftnlen)8);
748 return 0;
749 }
750
751 /* Find the geometric state of the target body with respect to */
752 /* the solar system barycenter. Subtract the state of the */
753 /* observer to get the relative state. Use this to compute the */
754 /* one-way light time. */
755
756 zzspkgo0_(targ, et, ref, &c__0, ssbtrg, &ssblt, ref_len);
757 if (failed_()) {
758 chkout_("ZZSPKLT0", (ftnlen)8);
759 return 0;
760 }
761 vsubg_(ssbtrg, stobs, &c__6, starg);
762 dist = vnorm_(starg);
763 *lt = dist / clight_();
764 if (*lt == 0.) {
765
766 /* This can happen only if the observer and target are at the */
767 /* same position. We don't consider this an error, but we're not */
768 /* going to compute the light time derivative. */
769
770 *dlt = 0.;
771 chkout_("ZZSPKLT0", (ftnlen)8);
772 return 0;
773 }
774 if (! uselt) {
775
776 /* This is a special case: we're not using light time */
777 /* corrections, so the derivative */
778 /* of light time is just */
779
780 /* (1/c) * d(VNORM(STARG))/dt */
781
782 *dlt = vdot_(starg, &starg[3]) / (dist * clight_());
783
784 /* LT and DLT are both set, so we can return. */
785
786 chkout_("ZZSPKLT0", (ftnlen)8);
787 return 0;
788 }
789
790 /* To correct for light time, find the state of the target body */
791 /* at the current epoch minus the one-way light time. Note that */
792 /* the observer remains where it is. */
793
794 /* Determine the sign of the light time offset. */
795
796 if (xmit) {
797 ltsign = 1;
798 } else {
799 ltsign = -1;
800 }
801
802 /* Let NUMITR be the number of iterations we'll perform to */
803 /* compute the light time. */
804
805 if (usecn) {
806 numitr = 5;
807 } else {
808 numitr = 1;
809 }
810 i__ = 0;
811 lterr = 1.;
812 while(i__ < numitr && lterr > 1e-17) {
813
814 /* LT was set either prior to this loop or */
815 /* during the previous loop iteration. */
816
817 epoch = *et + ltsign * *lt;
818 zzspkgo0_(targ, &epoch, ref, &c__0, ssbtrg, &ssblt, ref_len);
819 if (failed_()) {
820 chkout_("ZZSPKLT0", (ftnlen)8);
821 return 0;
822 }
823 vsubg_(ssbtrg, stobs, &c__6, starg);
824 prvlt = *lt;
825 d__1 = vnorm_(starg) / clight_();
826 *lt = touchd_(&d__1);
827 /* LTERR is the magnitude of the change between the current */
828 /* estimate of light time and the previous estimate, relative to */
829 /* the previous light time corrected epoch. */
830
831 /* Computing MAX */
832 d__3 = 1., d__4 = abs(epoch);
833 d__2 = (d__1 = *lt - prvlt, abs(d__1)) / max(d__3,d__4);
834 lterr = touchd_(&d__2);
835 ++i__;
836 }
837
838 /* At this point, STARG contains the light time corrected */
839 /* state of the target relative to the observer. */
840
841 /* Compute the derivative of light time with respect */
842 /* to time: dLT/dt. Below we derive the formula for */
843 /* this quantity for the reception case. Let */
844
845 /* POBS be the position of the observer relative to the */
846 /* solar system barycenter. */
847
848 /* VOBS be the velocity of the observer relative to the */
849 /* solar system barycenter. */
850
851 /* PTARG be the position of the target relative to the */
852 /* solar system barycenter. */
853
854 /* VTARG be the velocity of the target relative to the */
855 /* solar system barycenter. */
856
857 /* S be the sign of the light time correction. S is */
858 /* negative for the reception case. */
859
860 /* The light-time corrected position of the target relative to */
861 /* the observer at observation time ET, given the one-way */
862 /* light time LT is: */
863
864 /* PTARG(ET+S*LT) - POBS(ET) */
865
866 /* The light-time corrected velocity of the target relative to */
867 /* the observer at observation time ET is */
868
869 /* VTARG(ET+S*LT)*( 1 + S*d(LT)/d(ET) ) - VOBS(ET) */
870
871 /* We need to compute dLT/dt. Below, we use the facts that, */
872 /* for a time-dependent vector X(t), */
873
874 /* ||X|| = <X,X> ** (1/2) */
875
876 /* d(||X||)/dt = (1/2)<X,X>**(-1/2) * 2 * <X,dX/dt> */
877
878 /* = <X,X>**(-1/2) * <X,dX/dt> */
879
880 /* = <X,dX/dt> / ||X|| */
881
882 /* Newtonian light time equation: */
883
884 /* LT = (1/c) * || PTARG(ET+S*LT) - POBS(ET)|| */
885
886 /* Differentiate both sides: */
887
888 /* dLT/dt = (1/c) * ( 1 / || PTARG(ET+S*LT) - POBS(ET) || ) */
889
890 /* * < PTARG(ET+S*LT) - POBS(ET), */
891 /* VTARG(ET+S*LT)*(1+S*d(LT)/d(ET)) - VOBS(ET) > */
892
893
894 /* = (1/c) * ( 1 / || PTARG(ET+S*LT) - POBS(ET) || ) */
895
896 /* * ( < PTARG(ET+S*LT) - POBS(ET), */
897 /* VTARG(ET+S*LT) - VOBS(ET) > */
898
899 /* + < PTARG(ET+S*LT) - POBS(ET), */
900 /* VTARG(ET+S*LT) > * (S*d(LT)/d(ET)) ) */
901
902 /* Let */
903
904 /* A = (1/c) * ( 1 / || PTARG(ET+S*LT) - POBS(ET) || ) */
905
906 /* B = < PTARG(ET+S*LT) - POBS(ET), VTARG(ET+S*LT) - VOBS(ET) > */
907
908 /* C = < PTARG(ET+S*LT) - POBS(ET), VTARG(ET+S*LT) > */
909
910 /* Then */
911
912 /* d(LT)/d(ET) = A * ( B + C * S*d(LT)/d(ET) ) */
913
914 /* which implies */
915
916 /* d(LT)/d(ET) = A*B / ( 1 - S*C*A ) */
917
918
919
920 a = 1. / (clight_() * vnorm_(starg));
921 b = vdot_(starg, &starg[3]);
922 c__ = vdot_(starg, &ssbtrg[3]);
923
924 /* For physically realistic target velocities, S*C*A cannot equal 1. */
925 /* We'll check for this case anyway. */
926
927 if (ltsign * c__ * a > .99999999989999999) {
928 setmsg_("Target range rate magnitude is approximately the speed of l"
929 "ight. The light time derivative cannot be computed.", (ftnlen)
930 110);
931 sigerr_("SPICE(DIVIDEBYZERO)", (ftnlen)19);
932 chkout_("ZZSPKLT0", (ftnlen)8);
933 return 0;
934 }
935
936 /* Compute DLT: the rate of change of light time. */
937
938 *dlt = a * b / (1. - ltsign * c__ * a);
939
940 /* Overwrite the velocity portion of the output state */
941 /* with the light-time corrected velocity. */
942
943 d__1 = ltsign * *dlt + 1.;
944 vlcom_(&d__1, &ssbtrg[3], &c_b19, &stobs[3], &starg[3]);
945 chkout_("ZZSPKLT0", (ftnlen)8);
946 return 0;
947 } /* zzspklt0_ */
948
949