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