1 /* illumf.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__100 = 100;
11 static integer c__3 = 3;
12 static doublereal c_b56 = 1.;
13 
14 /* $Procedure ILLUMF ( Illumination angles, general source, return flags ) */
illumf_(char * method,char * target,char * ilusrc,doublereal * et,char * fixref,char * abcorr,char * obsrvr,doublereal * spoint,doublereal * trgepc,doublereal * srfvec,doublereal * phase,doublereal * incdnc,doublereal * emissn,logical * visibl,logical * lit,ftnlen method_len,ftnlen target_len,ftnlen ilusrc_len,ftnlen fixref_len,ftnlen abcorr_len,ftnlen obsrvr_len)15 /* Subroutine */ int illumf_(char *method, char *target, char *ilusrc,
16 	doublereal *et, char *fixref, char *abcorr, char *obsrvr, doublereal *
17 	spoint, doublereal *trgepc, doublereal *srfvec, doublereal *phase,
18 	doublereal *incdnc, doublereal *emissn, logical *visibl, logical *lit,
19 	 ftnlen method_len, ftnlen target_len, ftnlen ilusrc_len, ftnlen
20 	fixref_len, ftnlen abcorr_len, ftnlen obsrvr_len)
21 {
22     /* Initialized data */
23 
24     static logical first = TRUE_;
25     static logical pri = FALSE_;
26     static integer nsurf = 0;
27     static char prvcor[5] = "     ";
28     static char prvmth[500] = "                                             "
29 	    "                                                                "
30 	    "                                                                "
31 	    "                                                                "
32 	    "                                                                "
33 	    "                                                                "
34 	    "                                                                "
35 	    "                                                                "
36 	    "       ";
37     static integer shape = 0;
38 
39     /* Builtin functions */
40     integer s_cmp(char *, char *, ftnlen, ftnlen);
41     /* Subroutine */ int s_copy(char *, char *, ftnlen, ftnlen);
42 
43     /* Local variables */
44     extern /* Subroutine */ int zzbods2c_(integer *, char *, integer *,
45 	    logical *, char *, integer *, logical *, ftnlen, ftnlen);
46     extern doublereal vsep_(doublereal *, doublereal *);
47     extern /* Subroutine */ int vequ_(doublereal *, doublereal *);
48     integer type__;
49     static logical xmit;
50     extern /* Subroutine */ int zzmaxrad_(doublereal *), zznamfrm_(integer *,
51 	    char *, integer *, char *, integer *, ftnlen, ftnlen), zzvalcor_(
52 	    char *, logical *, ftnlen), zzsbfnrm_(integer *, integer *,
53 	    integer *, doublereal *, integer *, doublereal *, doublereal *),
54 	    zzsudski_(integer *, integer *, integer *, integer *), zzctruin_(
55 	    integer *), zzprsmet_(integer *, char *, integer *, char *, char *
56 	    , logical *, integer *, integer *, char *, char *, ftnlen, ftnlen,
57 	     ftnlen, ftnlen, ftnlen);
58     integer n;
59     extern /* Subroutine */ int zzsrftrk_(integer *, logical *), zzraysfx_(
60 	    doublereal *, doublereal *, doublereal *, doublereal *, logical *)
61 	    ;
62     doublereal s, scale, radii[3];
63     extern /* Subroutine */ int chkin_(char *, ftnlen), errch_(char *, char *,
64 	     ftnlen, ftnlen);
65     logical found;
66     extern /* Subroutine */ int vlcom_(doublereal *, doublereal *, doublereal
67 	    *, doublereal *, doublereal *);
68     extern logical eqstr_(char *, char *, ftnlen, ftnlen);
69     static logical uselt, svfnd1, svfnd2;
70     static integer svctr1[2], svctr2[2];
71     extern logical failed_(void);
72     static integer svctr3[2], svctr4[2];
73     doublereal lt;
74     integer obscde;
75     extern doublereal halfpi_(void);
76     extern /* Subroutine */ int bodvcd_(integer *, char *, integer *, integer
77 	    *, doublereal *, ftnlen);
78     integer fixfid;
79     static integer trgcde;
80     doublereal maxrad;
81     extern logical return_(void);
82     char pntdef[20], shpstr[9], subtyp[20], trmstr[20];
83     doublereal normal[3], obspos[3], opstat[6], tistat[6], vertex[3];
84     static integer center, srflst[100];
85     integer svnsrf, typeid;
86     logical attblk[15], surfup;
87     static char svtarg[36];
88     logical fnd;
89     static integer svtcde;
90     static char svobsr[36];
91     static integer svobsc;
92     static char svfref[32];
93     static integer svrefc;
94     extern /* Subroutine */ int chkout_(char *, ftnlen), setmsg_(char *,
95 	    ftnlen), sigerr_(char *, ftnlen), frinfo_(integer *, integer *,
96 	    integer *, integer *, logical *), errint_(char *, integer *,
97 	    ftnlen), spkcpt_(doublereal *, char *, char *, doublereal *, char
98 	    *, char *, char *, char *, doublereal *, doublereal *, ftnlen,
99 	    ftnlen, ftnlen, ftnlen, ftnlen, ftnlen), spkcpo_(char *,
100 	    doublereal *, char *, char *, char *, doublereal *, char *, char *
101 	    , doublereal *, doublereal *, ftnlen, ftnlen, ftnlen, ftnlen,
102 	    ftnlen, ftnlen), vminus_(doublereal *, doublereal *), surfnm_(
103 	    doublereal *, doublereal *, doublereal *, doublereal *,
104 	    doublereal *), vhatip_(doublereal *);
105     doublereal lti, xpt[3];
106 
107 /* $ Abstract */
108 
109 /*     Compute the illumination angles---phase, incidence, and */
110 /*     emission---at a specified point on a target body. Return logical */
111 /*     flags indicating whether the surface point is visible from */
112 /*     the observer's position and whether the surface point is */
113 /*     illuminated. */
114 
115 /*     The target body's surface is represented using topographic data */
116 /*     provided by DSK files, or by a reference ellipsoid. */
117 
118 /*     The illumination source is a specified ephemeris object. */
119 
120 /* $ Disclaimer */
121 
122 /*     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
123 /*     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
124 /*     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
125 /*     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
126 /*     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
127 /*     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
128 /*     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
129 /*     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
130 /*     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
131 /*     SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
132 
133 /*     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
134 /*     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
135 /*     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
136 /*     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
137 /*     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
138 /*     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
139 
140 /*     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
141 /*     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
142 /*     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
143 /*     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
144 
145 /* $ Required_Reading */
146 
147 /*     DSK */
148 /*     FRAMES */
149 /*     NAIF_IDS */
150 /*     PCK */
151 /*     SPK */
152 /*     TIME */
153 
154 /* $ Keywords */
155 
156 /*     ANGLES */
157 /*     GEOMETRY */
158 /*     ILLUMINATION */
159 
160 /* $ Declarations */
161 
162 /*     File: dsk.inc */
163 
164 
165 /*     Version 1.0.0 05-FEB-2016 (NJB) */
166 
167 /*     Maximum size of surface ID list. */
168 
169 
170 /*     End of include file dsk.inc */
171 
172 /* $ Abstract */
173 
174 /*     This file contains public, global parameter declarations */
175 /*     for the SPICELIB Geometry Finder (GF) subsystem. */
176 
177 /* $ Disclaimer */
178 
179 /*     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
180 /*     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
181 /*     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
182 /*     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
183 /*     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
184 /*     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
185 /*     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
186 /*     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
187 /*     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
188 /*     SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
189 
190 /*     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
191 /*     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
192 /*     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
193 /*     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
194 /*     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
195 /*     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
196 
197 /*     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
198 /*     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
199 /*     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
200 /*     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
201 
202 /* $ Required_Reading */
203 
204 /*     GF */
205 
206 /* $ Keywords */
207 
208 /*     GEOMETRY */
209 /*     ROOT */
210 
211 /* $ Restrictions */
212 
213 /*     None. */
214 
215 /* $ Author_and_Institution */
216 
217 /*     N.J. Bachman      (JPL) */
218 /*     L.E. Elson        (JPL) */
219 /*     E.D. Wright       (JPL) */
220 
221 /* $ Literature_References */
222 
223 /*     None. */
224 
225 /* $ Version */
226 
227 /* -    SPICELIB Version 2.0.0  29-NOV-2016 (NJB) */
228 
229 /*        Upgraded to support surfaces represented by DSKs. */
230 
231 /*        Bug fix: removed declaration of NVRMAX parameter. */
232 
233 /* -    SPICELIB Version 1.3.0, 01-OCT-2011 (NJB) */
234 
235 /*       Added NWILUM parameter. */
236 
237 /* -    SPICELIB Version 1.2.0, 14-SEP-2010 (EDW) */
238 
239 /*       Added NWPA parameter. */
240 
241 /* -    SPICELIB Version 1.1.0, 08-SEP-2009 (EDW) */
242 
243 /*       Added NWRR parameter. */
244 /*       Added NWUDS parameter. */
245 
246 /* -    SPICELIB Version 1.0.0, 21-FEB-2009 (NJB) (LSE) (EDW) */
247 
248 /* -& */
249 
250 /*     Root finding parameters: */
251 
252 /*     CNVTOL is the default convergence tolerance used by the */
253 /*     high-level GF search API routines. This tolerance is */
254 /*     used to terminate searches for binary state transitions: */
255 /*     when the time at which a transition occurs is bracketed */
256 /*     by two times that differ by no more than CNVTOL, the */
257 /*     transition time is considered to have been found. */
258 
259 /*     Units are TDB seconds. */
260 
261 
262 /*     NWMAX is the maximum number of windows allowed for user-defined */
263 /*     workspace array. */
264 
265 /*        DOUBLE PRECISION      WORK   ( LBCELL : MW, NWMAX ) */
266 
267 /*     Currently no more than twelve windows are required; the three */
268 /*     extra windows are spares. */
269 
270 /*     Callers of GFEVNT can include this file and use the parameter */
271 /*     NWMAX to declare the second dimension of the workspace array */
272 /*     if necessary. */
273 
274 
275 /*     Callers of GFIDST should declare their workspace window */
276 /*     count using NWDIST. */
277 
278 
279 /*     Callers of GFSEP should declare their workspace window */
280 /*     count using NWSEP. */
281 
282 
283 /*     Callers of GFRR should declare their workspace window */
284 /*     count using NWRR. */
285 
286 
287 /*     Callers of GFUDS should declare their workspace window */
288 /*     count using NWUDS. */
289 
290 
291 /*     Callers of GFPA should declare their workspace window */
292 /*     count using NWPA. */
293 
294 
295 /*     Callers of GFILUM should declare their workspace window */
296 /*     count using NWILUM. */
297 
298 
299 /*     ADDWIN is a parameter used to expand each interval of the search */
300 /*     (confinement) window by a small amount at both ends in order to */
301 /*     accommodate searches using equality constraints. The loaded */
302 /*     kernel files must accommodate these expanded time intervals. */
303 
304 
305 /*     FRMNLN is a string length for frame names. */
306 
307 
308 /*     FOVTLN -- maximum length for FOV string. */
309 
310 
311 /*     Specify the character strings that are allowed in the */
312 /*     specification of field of view shapes. */
313 
314 
315 /*     Character strings that are allowed in the */
316 /*     specification of occultation types: */
317 
318 
319 /*     Occultation target shape specifications: */
320 
321 
322 /*     Specify the number of supported occultation types and occultation */
323 /*     type string length: */
324 
325 
326 /*     Instrument field-of-view (FOV) parameters */
327 
328 /*     Maximum number of FOV boundary vectors: */
329 
330 
331 /*     FOV shape parameters: */
332 
333 /*        circle */
334 /*        ellipse */
335 /*        polygon */
336 /*        rectangle */
337 
338 
339 /*     End of file gf.inc. */
340 
341 /* $ Abstract */
342 
343 /*     Include file zzabcorr.inc */
344 
345 /*     SPICE private file intended solely for the support of SPICE */
346 /*     routines.  Users should not include this file directly due */
347 /*     to the volatile nature of this file */
348 
349 /*     The parameters below define the structure of an aberration */
350 /*     correction attribute block. */
351 
352 /* $ Disclaimer */
353 
354 /*     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
355 /*     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
356 /*     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
357 /*     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
358 /*     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
359 /*     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
360 /*     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
361 /*     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
362 /*     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
363 /*     SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
364 
365 /*     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
366 /*     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
367 /*     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
368 /*     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
369 /*     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
370 /*     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
371 
372 /*     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
373 /*     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
374 /*     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
375 /*     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
376 
377 /* $ Parameters */
378 
379 /*     An aberration correction attribute block is an array of logical */
380 /*     flags indicating the attributes of the aberration correction */
381 /*     specified by an aberration correction string.  The attributes */
382 /*     are: */
383 
384 /*        - Is the correction "geometric"? */
385 
386 /*        - Is light time correction indicated? */
387 
388 /*        - Is stellar aberration correction indicated? */
389 
390 /*        - Is the light time correction of the "converged */
391 /*          Newtonian" variety? */
392 
393 /*        - Is the correction for the transmission case? */
394 
395 /*        - Is the correction relativistic? */
396 
397 /*    The parameters defining the structure of the block are as */
398 /*    follows: */
399 
400 /*       NABCOR    Number of aberration correction choices. */
401 
402 /*       ABATSZ    Number of elements in the aberration correction */
403 /*                 block. */
404 
405 /*       GEOIDX    Index in block of geometric correction flag. */
406 
407 /*       LTIDX     Index of light time flag. */
408 
409 /*       STLIDX    Index of stellar aberration flag. */
410 
411 /*       CNVIDX    Index of converged Newtonian flag. */
412 
413 /*       XMTIDX    Index of transmission flag. */
414 
415 /*       RELIDX    Index of relativistic flag. */
416 
417 /*    The following parameter is not required to define the block */
418 /*    structure, but it is convenient to include it here: */
419 
420 /*       CORLEN    The maximum string length required by any aberration */
421 /*                 correction string */
422 
423 /* $ Author_and_Institution */
424 
425 /*     N.J. Bachman    (JPL) */
426 
427 /* $ Literature_References */
428 
429 /*     None. */
430 
431 /* $ Version */
432 
433 /* -    SPICELIB Version 1.0.0, 18-DEC-2004 (NJB) */
434 
435 /* -& */
436 /*     Number of aberration correction choices: */
437 
438 
439 /*     Aberration correction attribute block size */
440 /*     (number of aberration correction attributes): */
441 
442 
443 /*     Indices of attributes within an aberration correction */
444 /*     attribute block: */
445 
446 
447 /*     Maximum length of an aberration correction string: */
448 
449 
450 /*     End of include file zzabcorr.inc */
451 
452 /* $ Abstract */
453 
454 /*     This include file defines the dimension of the counter */
455 /*     array used by various SPICE subsystems to uniquely identify */
456 /*     changes in their states. */
457 
458 /* $ Disclaimer */
459 
460 /*     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
461 /*     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
462 /*     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
463 /*     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
464 /*     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
465 /*     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
466 /*     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
467 /*     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
468 /*     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
469 /*     SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
470 
471 /*     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
472 /*     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
473 /*     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
474 /*     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
475 /*     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
476 /*     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
477 
478 /*     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
479 /*     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
480 /*     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
481 /*     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
482 
483 /* $ Parameters */
484 
485 /*     CTRSIZ      is the dimension of the counter array used by */
486 /*                 various SPICE subsystems to uniquely identify */
487 /*                 changes in their states. */
488 
489 /* $ Author_and_Institution */
490 
491 /*     B.V. Semenov    (JPL) */
492 
493 /* $ Literature_References */
494 
495 /*     None. */
496 
497 /* $ Version */
498 
499 /* -    SPICELIB Version 1.0.0, 29-JUL-2013 (BVS) */
500 
501 /* -& */
502 
503 /*     End of include file. */
504 
505 
506 /*     File: zzdsk.inc */
507 
508 
509 /*     Version 4.0.0 13-NOV-2015 (NJB) */
510 
511 /*        Changed parameter LBTLEN to CVTLEN. */
512 /*        Added parameter LMBCRV. */
513 
514 /*     Version 3.0.0 05-NOV-2015 (NJB) */
515 
516 /*        Added parameters */
517 
518 /*           CTRCOR */
519 /*           ELLCOR */
520 /*           GUIDED */
521 /*           LBTLEN */
522 /*           PNMBRL */
523 /*           TANGNT */
524 /*           TMTLEN */
525 /*           UMBRAL */
526 
527 /*     Version 2.0.0 04-MAR-2015 (NJB) */
528 
529 /*        Removed declaration of parameter SHPLEN. */
530 /*        This name is already in use in the include */
531 /*        file gf.inc. */
532 
533 /*     Version 1.0.0 26-JAN-2015 (NJB) */
534 
535 
536 /*     Parameters supporting METHOD string parsing: */
537 
538 
539 /*     Local method length. */
540 
541 
542 /*     Length of sub-point type string. */
543 
544 
545 /*     Length of curve type string. */
546 
547 
548 /*     Limb type parameter codes. */
549 
550 
551 /*     Length of terminator type string. */
552 
553 
554 /*     Terminator type and limb parameter codes. */
555 
556 
557 /*     Length of aberration correction locus string. */
558 
559 
560 /*     Aberration correction locus codes. */
561 
562 
563 /*     End of include file zzdsk.inc */
564 
565 /* $ Brief_I/O */
566 
567 /*     Variable  I/O  Description */
568 /*     --------  ---  -------------------------------------------------- */
569 /*     METHOD     I   Computation method. */
570 /*     TARGET     I   Name of target body. */
571 /*     ILUSRC     I   Name of illumination source. */
572 /*     ET         I   Epoch in ephemeris seconds past J2000 TDB. */
573 /*     FIXREF     I   Body-fixed, body-centered target body frame. */
574 /*     ABCORR     I   Desired aberration correction. */
575 /*     OBSRVR     I   Name of observing body. */
576 /*     SPOINT     I   Body-fixed coordinates of a target surface point. */
577 /*     TRGEPC     O   Target surface point epoch. */
578 /*     SRFVEC     O   Vector from observer to target surface point. */
579 /*     PHASE      O   Phase angle at the surface point. */
580 /*     INCDNC     O   Source incidence angle at the surface point. */
581 /*     EMISSN     O   Emission angle at the surface point. */
582 /*     VISIBL     O   Visibility flag (.TRUE. = visible). */
583 /*     LIT        O   Illumination flag (.TRUE. = illuminated). */
584 
585 /* $ Detailed_Input */
586 
587 /*      METHOD     is a short string providing parameters defining */
588 /*                 the computation method to be used. In the syntax */
589 /*                 descriptions below, items delimited by brackets */
590 /*                 are optional. */
591 
592 /*                 METHOD may be assigned the following values: */
593 
594 /*                    'ELLIPSOID' */
595 
596 /*                       The illumination angle computation uses a */
597 /*                       triaxial ellipsoid to model the surface of the */
598 /*                       target body. The ellipsoid's radii must be */
599 /*                       available in the kernel pool. */
600 
601 
602 /*                    'DSK/UNPRIORITIZED[/SURFACES = <surface list>]' */
603 
604 /*                       The illumination angle computation uses */
605 /*                       topographic data to model the surface of the */
606 /*                       target body. These data must be provided by */
607 /*                       loaded DSK files. */
608 
609 /*                       The surface list specification is optional. The */
610 /*                       syntax of the list is */
611 
612 /*                          <surface 1> [, <surface 2>...] */
613 
614 /*                       If present, it indicates that data only for the */
615 /*                       listed surfaces are to be used; however, data */
616 /*                       need not be available for all surfaces in the */
617 /*                       list. If absent, loaded DSK data for any surface */
618 /*                       associated with the target body are used. */
619 
620 /*                       The surface list may contain surface names or */
621 /*                       surface ID codes. Names containing blanks must */
622 /*                       be delimited by double quotes, for example */
623 
624 /*                          SURFACES = "Mars MEGDR 128 PIXEL/DEG" */
625 
626 /*                       If multiple surfaces are specified, their names */
627 /*                       or IDs must be separated by commas. */
628 
629 /*                       See the Particulars section below for details */
630 /*                       concerning use of DSK data. */
631 
632 
633 /*                 Neither case nor white space are significant in */
634 /*                 METHOD, except within double-quoted strings. For */
635 /*                 example, the string ' eLLipsoid ' is valid. */
636 
637 /*                 Within double-quoted strings, blank characters are */
638 /*                 significant, but multiple consecutive blanks are */
639 /*                 considered equivalent to a single blank. Case is */
640 /*                 not significant. So */
641 
642 /*                    "Mars MEGDR 128 PIXEL/DEG" */
643 
644 /*                 is equivalent to */
645 
646 /*                    " mars megdr  128  pixel/deg " */
647 
648 /*                 but not to */
649 
650 /*                    "MARS MEGDR128PIXEL/DEG" */
651 
652 
653 /*     TARGET      is the name of the target body. TARGET is */
654 /*                 case-insensitive, and leading and trailing blanks in */
655 /*                 TARGET are not significant. Optionally, you may */
656 /*                 supply a string containing the integer ID code for */
657 /*                 the object. For example both 'MOON' and '301' are */
658 /*                 legitimate strings that indicate the Moon is the */
659 /*                 target body. */
660 
661 
662 /*     ILUSRC      is the name of the illumination source. This source */
663 /*                 may be any ephemeris object. Case, blanks, and */
664 /*                 numeric values are treated in the same way as for the */
665 /*                 input TARGET. */
666 
667 
668 /*     ET          is the epoch, expressed as seconds past J2000 TDB, */
669 /*                 for which the apparent illumination angles at the */
670 /*                 specified surface point on the target body, as seen */
671 /*                 from the observing body, are to be computed. */
672 
673 
674 /*     FIXREF      is the name of a body-fixed reference frame centered */
675 /*                 on the target body. FIXREF may be any such frame */
676 /*                 supported by the SPICE system, including built-in */
677 /*                 frames (documented in the Frames Required Reading) */
678 /*                 and frames defined by a loaded frame kernel (FK). The */
679 /*                 string FIXREF is case-insensitive, and leading and */
680 /*                 trailing blanks in FIXREF are not significant. */
681 
682 /*                 The input surface point SPOINT and the output vector */
683 /*                 SRFVEC are expressed relative to this reference */
684 /*                 frame. */
685 
686 
687 /*     ABCORR      is the aberration correction to be used in computing */
688 /*                 the position and orientation of the target body and */
689 /*                 the location of the illumination source. */
690 
691 /*                 For remote sensing applications, where the apparent */
692 /*                 illumination angles seen by the observer are desired, */
693 /*                 normally either of the corrections */
694 
695 /*                    'LT+S' */
696 /*                    'CN+S' */
697 
698 /*                 should be used. These and the other supported options */
699 /*                 are described below. ABCORR may be any of the */
700 /*                 following: */
701 
702 /*                    'NONE'     No aberration correction. */
703 
704 /*                 Let LT represent the one-way light time between the */
705 /*                 observer and SPOINT (note: NOT between the observer */
706 /*                 and the target body's center). The following values */
707 /*                 of ABCORR apply to the "reception" case in which */
708 /*                 photons depart from SPOINT at the light-time */
709 /*                 corrected epoch ET-LT and *arrive* at the observer's */
710 /*                 location at ET: */
711 
712 /*                    'LT'       Correct both the position of SPOINT as */
713 /*                               seen by the observer, and the position */
714 /*                               of the illumination source as seen by */
715 /*                               the target, for light time. Correct the */
716 /*                               orientation of the target for light */
717 /*                               time. */
718 
719 /*                    'LT+S'     Correct both the position of SPOINT as */
720 /*                               seen by the observer, and the position */
721 /*                               of the illumination source as seen by */
722 /*                               the target, for light time and stellar */
723 /*                               aberration. Correct the orientation of */
724 /*                               the target for light time. */
725 
726 /*                    'CN'       Converged Newtonian light time */
727 /*                               correction. In solving the light time */
728 /*                               equations for target and the */
729 /*                               illumination source, the 'CN' */
730 /*                               correction iterates until the solution */
731 /*                               converges. */
732 
733 /*                    'CN+S'     Converged Newtonian light time and */
734 /*                               stellar aberration corrections. This */
735 /*                               option produces a solution that is at */
736 /*                               least as accurate at that obtainable */
737 /*                               with the 'LT+S' option. Whether the */
738 /*                               'CN+S' solution is substantially more */
739 /*                               accurate depends on the geometry of the */
740 /*                               participating objects and on the */
741 /*                               accuracy of the input data. In all */
742 /*                               cases this routine will execute more */
743 /*                               slowly when a converged solution is */
744 /*                               computed. */
745 
746 
747 /*                 The following values of ABCORR apply to the */
748 /*                 "transmission" case in which photons *arrive* at */
749 /*                 SPOINT at the light-time corrected epoch ET+LT and */
750 /*                 *depart* from the observer's location at ET: */
751 
752 /*                    'XLT'      "Transmission" case: correct for */
753 /*                               one-way light time using a Newtonian */
754 /*                               formulation. This correction yields the */
755 /*                               illumination angles at the moment that */
756 /*                               SPOINT receives photons emitted from the */
757 /*                               observer's location at ET. */
758 
759 /*                               The light time correction uses an */
760 /*                               iterative solution of the light time */
761 /*                               equation. The solution invoked by the */
762 /*                               'XLT' option uses one iteration. */
763 
764 /*                               Both the target position as seen by the */
765 /*                               observer, and rotation of the target */
766 /*                               body, are corrected for light time. */
767 
768 /*                    'XLT+S'    "Transmission" case: correct for */
769 /*                               one-way light time and stellar */
770 /*                               aberration using a Newtonian */
771 /*                               formulation  This option modifies the */
772 /*                               angles obtained with the 'XLT' option */
773 /*                               to account for the observer's and */
774 /*                               target's velocities relative to the */
775 /*                               solar system barycenter (the latter */
776 /*                               velocity is used in computing the */
777 /*                               direction to the apparent illumination */
778 /*                               source). */
779 
780 /*                    'XCN'      Converged Newtonian light time */
781 /*                               correction. This is the same as XLT */
782 /*                               correction but with further iterations */
783 /*                               to a converged Newtonian light time */
784 /*                               solution. */
785 
786 /*                    'XCN+S'    "Transmission" case: converged */
787 /*                               Newtonian light time and stellar */
788 /*                               aberration corrections. This option */
789 /*                               produces a solution that is at least as */
790 /*                               accurate at that obtainable with the */
791 /*                               'XLT+S' option. Whether the 'XCN+S' */
792 /*                               solution is substantially more accurate */
793 /*                               depends on the geometry of the */
794 /*                               participating objects and on the */
795 /*                               accuracy of the input data. In all */
796 /*                               cases this routine will execute more */
797 /*                               slowly when a converged solution is */
798 /*                               computed. */
799 
800 /*                 Neither case nor white space are significant in */
801 /*                 ABCORR. For example, the string */
802 
803 /*                   'Lt + s' */
804 
805 /*                 is valid. */
806 
807 
808 /*     OBSRVR      is the name of the observing body. The observing body */
809 /*                 is an ephemeris object: it typically is a spacecraft */
810 /*                 or a surface point on an extended body such as a */
811 /*                 planet. OBSRVR is case-insensitive, and leading and */
812 /*                 trailing blanks in OBSRVR are not significant. */
813 /*                 Optionally, you may supply a string containing the */
814 /*                 integer ID code for the object. For example both */
815 /*                 'MOON' and '301' are legitimate strings that indicate */
816 /*                 the Moon is the observer. */
817 
818 /*                 OBSRVR may be not be identical to TARGET. */
819 
820 
821 /*     SPOINT      is a surface point on the target body, expressed in */
822 /*                 Cartesian coordinates, relative to the body-fixed */
823 /*                 target frame designated by FIXREF. */
824 
825 /*                 SPOINT need not be visible from the observer's */
826 /*                 location at the epoch ET. */
827 
828 /*                 The components of SPOINT have units of km. */
829 
830 
831 /* $ Detailed_Output */
832 
833 
834 /*     TRGEPC      is the "target surface point epoch." TRGEPC is */
835 /*                 defined as follows: letting LT be the one-way light */
836 /*                 time between the observer and the input surface point */
837 /*                 SPOINT, TRGEPC is either the epoch ET-LT or ET */
838 /*                 depending on whether the requested aberration */
839 /*                 correction is, respectively, for received radiation */
840 /*                 or omitted. LT is computed using the method indicated */
841 /*                 by ABCORR. */
842 
843 /*                 TRGEPC is expressed as seconds past J2000 TDB. */
844 
845 
846 /*     SRFVEC      is the vector from the observer's position at ET to */
847 /*                 the aberration-corrected (or optionally, geometric) */
848 /*                 position of SPOINT, where the aberration corrections */
849 /*                 are specified by ABCORR. SRFVEC is expressed in the */
850 /*                 target body-fixed reference frame designated by */
851 /*                 FIXREF, evaluated at TRGEPC. */
852 
853 /*                 The components of SRFVEC are given in units of km. */
854 
855 /*                 One can use the SPICELIB function VNORM to obtain the */
856 /*                 distance between the observer and SPOINT: */
857 
858 /*                    DIST = VNORM ( SRFVEC ) */
859 
860 /*                 The observer's position OBSPOS, relative to the */
861 /*                 target body's center, where the center's position is */
862 /*                 corrected for aberration effects as indicated by */
863 /*                 ABCORR, can be computed via the call: */
864 
865 /*                    CALL VSUB ( SPOINT, SRFVEC, OBSPOS ) */
866 
867 /*                 To transform the vector SRFVEC to a time-dependent */
868 /*                 reference frame REF at ET, a sequence of calls is */
869 /*                 required. For example, let XFORM be 3x3 matrix */
870 /*                 describing the transformation between the target */
871 /*                 body-fixed frame at TRGEPC to the time-dependent */
872 /*                 frame REF at ET. Then SRFVEC can be transformed to */
873 /*                 the result REFVEC as follows: */
874 
875 /*                     CALL PXFRM2 ( FIXREF, REF,    TRGEPC, ET, XFORM ) */
876 /*                     CALL MXV    ( XFORM,  SRFVEC, REFVEC ) */
877 
878 
879 
880 /*     The following outputs depend on the existence of a well-defined */
881 /*     outward normal vector to the surface at SPOINT. See restriction 1. */
882 
883 
884 /*     PHASE       is the phase angle at SPOINT, as seen from OBSRVR at */
885 /*                 time ET. This is the angle between the negative of */
886 /*                 the vector SRFVEC and the SPOINT-source vector at */
887 /*                 TRGEPC. Units are radians. The range of PHASE is */
888 /*                 [0, pi]. */
889 
890 
891 /*     INCDNC      is the illumination source incidence angle at SPOINT, */
892 /*                 as seen from OBSRVR at time ET. This is the angle */
893 /*                 between the surface normal vector at SPOINT and the */
894 /*                 SPOINT-source vector at TRGEPC. Units are radians. */
895 /*                 The range of INCDNC is [0, pi]. */
896 
897 
898 /*     EMISSN      is the emission angle at SPOINT, as seen from OBSRVR */
899 /*                 at time ET. This is the angle between the surface */
900 /*                 normal vector at SPOINT and the negative of the */
901 /*                 vector SRFVEC. Units are radians. The range of EMISSN */
902 /*                 is [0, pi]. */
903 
904 
905 /*     VISIBL      is a logical flag indicating whether the surface */
906 /*                 point is visible to the observer. VISIBL takes into */
907 /*                 account whether the target surface occults SPOINT, */
908 /*                 regardless of the emission angle at SPOINT. VISIBL is */
909 /*                 returned with the value .TRUE. if SPOINT is visible; */
910 /*                 otherwise it is .FALSE. */
911 
912 
913 /*     LIT         is a logical flag indicating whether the surface */
914 /*                 point is illuminated; the point is considered to be */
915 /*                 illuminated if the vector from the point to the */
916 /*                 center of the illumination source doesn't intersect */
917 /*                 the target surface. LIT takes into account whether */
918 /*                 the target surface casts a shadow on SPOINT, */
919 /*                 regardless of the incidence angle at SPOINT. LIT is */
920 /*                 returned with the value .TRUE. if SPOINT is */
921 /*                 illuminated; otherwise it is .FALSE. */
922 
923 /* $ Parameters */
924 
925 /*     None. */
926 
927 /* $ Exceptions */
928 
929 
930 /*     1)  If the specified aberration correction is relativistic or */
931 /*         calls for stellar aberration but not light time correction, */
932 /*         the error SPICE(NOTSUPPORTED) is signaled. If the specified */
933 /*         aberration correction is any other unrecognized value, the */
934 /*         error will be diagnosed and signaled by a routine in the call */
935 /*         tree of this routine. */
936 
937 /*     2)  If any of the target, observer, or illumination source */
938 /*         input strings cannot be converted to an integer ID code, the */
939 /*         error SPICE(IDCODENOTFOUND) is signaled. */
940 
941 /*     3)  If OBSRVR and TARGET map to the same NAIF integer ID code, */
942 /*         the error SPICE(BODIESNOTDISTINCT) is signaled. */
943 
944 /*     4)  If the input target body-fixed frame FIXREF is not */
945 /*         recognized, the error SPICE(NOFRAME) is signaled. A frame */
946 /*         name may fail to be recognized because a required frame */
947 /*         specification kernel has not been loaded; another cause is a */
948 /*         misspelling of the frame name. */
949 
950 /*     5)  If the input frame FIXREF is not centered at the target body, */
951 /*         the error SPICE(INVALIDFRAME) is signaled. */
952 
953 /*     6)  If the input argument METHOD cannot be parsed, the error */
954 /*         is signaled either by this routine or by a routine in the */
955 /*         call tree of this routine. */
956 
957 /*     7)  If insufficient ephemeris data have been loaded prior to */
958 /*         calling ILLUMF, the error will be diagnosed and signaled by a */
959 /*         routine in the call tree of this routine. Note that when */
960 /*         light time correction is used, sufficient ephemeris data must */
961 /*         be available to propagate the states of observer, target, and */
962 /*         the illumination source to the solar system barycenter. */
963 
964 /*     8)  If the computation method specifies an ellipsoidal target */
965 /*         shape and triaxial radii of the target body have not been */
966 /*         loaded into the kernel pool prior to calling ILLUMF, the */
967 /*         error will be diagnosed and signaled by a routine in the call */
968 /*         tree of this routine. */
969 
970 /*     9)  The target must be an extended body: if any of the radii of */
971 /*         the target body are non-positive, the error will be */
972 /*         diagnosed and signaled by routines in the call tree of this */
973 /*         routine. */
974 
975 /*     10) If PCK data specifying the target body-fixed frame */
976 /*         orientation have not been loaded prior to calling ILLUMF, */
977 /*         the error will be diagnosed and signaled by a routine in the */
978 /*         call tree of this routine. */
979 
980 /*     11) If METHOD specifies that the target surface is represented by */
981 /*         DSK data, and no DSK files are loaded for the specified */
982 /*         target, the error is signaled by a routine in the call tree */
983 /*         of this routine. */
984 
985 /*     12) If METHOD specifies that the target surface is represented */
986 /*         by DSK data, and data representing the portion of the surface */
987 /*         on which SPOINT is located are not available, an error will */
988 /*         be signaled by a routine in the call tree of this routine. */
989 
990 /*     13) If METHOD specifies that the target surface is represented */
991 /*         by DSK data, SPOINT must lie on the target surface, not above */
992 /*         or below it. A small tolerance is used to allow for round-off */
993 /*         error in the calculation determining whether SPOINT is on the */
994 /*         surface. If, in the DSK case, SPOINT is too far from the */
995 /*         surface, an error will be signaled by a routine in the call */
996 /*         tree of this routine. */
997 
998 /*         If the surface is represented by a triaxial ellipsoid, SPOINT */
999 /*         is not required to be close to the ellipsoid; however, the */
1000 /*         results computed by this routine will be unreliable if SPOINT */
1001 /*         is too far from the ellipsoid. */
1002 
1003 /* $ Files */
1004 
1005 /*     Appropriate kernels must be loaded by the calling program before */
1006 /*     this routine is called. */
1007 
1008 /*     The following data are required: */
1009 
1010 /*        - SPK data: ephemeris data for target, observer, and the */
1011 /*          illumination source must be loaded. If aberration */
1012 /*          corrections are used, the states of target, observer, and */
1013 /*          the illumination source relative to the solar system */
1014 /*          barycenter must be calculable from the available ephemeris */
1015 /*          data. Typically ephemeris data are made available by loading */
1016 /*          one or more SPK files via FURNSH. */
1017 
1018 /*        - PCK data: rotation data for the target body must be */
1019 /*          loaded. These may be provided in a text or binary PCK file. */
1020 
1021 /*        - Shape data for the target body: */
1022 
1023 /*            PCK data: */
1024 
1025 /*               If the target body shape is modeled as an ellipsoid, */
1026 /*               triaxial radii for the target body must be loaded into */
1027 /*               the kernel pool. Typically this is done by loading a */
1028 /*               text PCK file via FURNSH. */
1029 
1030 /*               Triaxial radii are also needed if the target shape is */
1031 /*               modeled by DSK data, but the DSK NADIR method is */
1032 /*               selected. */
1033 
1034 /*            DSK data: */
1035 
1036 /*               If the target shape is modeled by DSK data, DSK files */
1037 /*               containing topographic data for the target body must be */
1038 /*               loaded. If a surface list is specified, data for at */
1039 /*               least one of the listed surfaces must be loaded. */
1040 
1041 /*     The following data may be required: */
1042 
1043 /*        - Frame data: if a frame definition is required to convert the */
1044 /*          observer and target states to the body-fixed frame of the */
1045 /*          target, that definition must be available in the kernel */
1046 /*          pool. Typically the definition is supplied by loading a */
1047 /*          frame kernel via FURNSH. */
1048 
1049 /*        - Surface name-ID associations: if surface names are specified */
1050 /*          in METHOD, the association of these names with their */
1051 /*          corresponding surface ID codes must be established by */
1052 /*          assignments of the kernel variables */
1053 
1054 /*             NAIF_SURFACE_NAME */
1055 /*             NAIF_SURFACE_CODE */
1056 /*             NAIF_SURFACE_BODY */
1057 
1058 /*          Normally these associations are made by loading a text */
1059 /*          kernel containing the necessary assignments. An example */
1060 /*          of such an assignment is */
1061 
1062 /*             NAIF_SURFACE_NAME += 'Mars MEGDR 128 PIXEL/DEG' */
1063 /*             NAIF_SURFACE_CODE += 1 */
1064 /*             NAIF_SURFACE_BODY += 499 */
1065 
1066 /*     In all cases, kernel data are normally loaded once per program */
1067 /*     run, NOT every time this routine is called. */
1068 
1069 
1070 /* $ Particulars */
1071 
1072 /*     SPICELIB contains four routines that compute illumination angles: */
1073 
1074 /*        ILLUMF (this routine) */
1075 
1076 /*        ILLUMG (same as this routine, except that */
1077 /*                output flags are not returned.) */
1078 
1079 /*        ILUMIN (same as ILLUMG, except that the sun is fixed */
1080 /*                as the illumination source.) */
1081 
1082 /*        ILLUM  (deprecated) */
1083 
1084 /*     This routine is the most capable of the set. */
1085 
1086 
1087 /*     Illumination angles */
1088 /*     =================== */
1089 
1090 /*     The term "illumination angles" refers to the following set of */
1091 /*     angles: */
1092 
1093 
1094 /*        phase angle              Angle between the vectors from the */
1095 /*                                 surface point to the observer and */
1096 /*                                 from the surface point to the */
1097 /*                                 illumination source. */
1098 
1099 /*        incidence angle          Angle between the surface normal at */
1100 /*                                 the specified surface point and the */
1101 /*                                 vector from the surface point to the */
1102 /*                                 illumination source. */
1103 
1104 /*        emission angle           Angle between the surface normal at */
1105 /*                                 the specified surface point and the */
1106 /*                                 vector from the surface point to the */
1107 /*                                 observer. */
1108 
1109 /*     The diagram below illustrates the geometric relationships */
1110 /*     defining these angles. The labels for the incidence, emission, */
1111 /*     and phase angles are "inc.", "e.", and "phase". */
1112 
1113 
1114 /*                                                      * */
1115 /*                                              illumination source */
1116 
1117 /*                    surface normal vector */
1118 /*                              ._                 _. */
1119 /*                              |\                 /|  illumination */
1120 /*                                \    phase      /    source vector */
1121 /*                                 \   .    .    / */
1122 /*                                 .            . */
1123 /*                                   \   ___   / */
1124 /*                              .     \/     \/ */
1125 /*                                    _\ inc./ */
1126 /*                             .    /   \   / */
1127 /*                             .   |  e. \ / */
1128 /*         *             <--------------- *  surface point on */
1129 /*      viewing            vector            target body */
1130 /*      location           to viewing */
1131 /*      (observer)         location */
1132 
1133 
1134 /*     Note that if the target-observer vector, the target normal vector */
1135 /*     at the surface point, and the target-illumination source vector */
1136 /*     are coplanar, then phase is the sum of the incidence and emission */
1137 /*     angles. This rarely occurs; usually */
1138 
1139 /*        phase angle  <  incidence angle + emission angle */
1140 
1141 /*     All of the above angles can be computed using light time */
1142 /*     corrections, light time and stellar aberration corrections, or no */
1143 /*     aberration corrections. In order to describe apparent geometry as */
1144 /*     observed by a remote sensing instrument, both light time and */
1145 /*     stellar aberration corrections should be used. */
1146 
1147 /*     The way aberration corrections are applied by this routine */
1148 /*     is described below. */
1149 
1150 /*        Light time corrections */
1151 /*        ====================== */
1152 
1153 /*           Observer-target surface point vector */
1154 /*           ------------------------------------ */
1155 
1156 /*           Let ET be the epoch at which an observation or remote */
1157 /*           sensing measurement is made, and let ET - LT ("LT" stands */
1158 /*           for "light time") be the epoch at which the photons */
1159 /*           received at ET were emitted from the surface point SPOINT. */
1160 /*           Note that the light time between the surface point and */
1161 /*           observer will generally differ from the light time between */
1162 /*           the target body's center and the observer. */
1163 
1164 
1165 /*           Target body's orientation */
1166 /*           ------------------------- */
1167 
1168 /*           Using the definitions of ET and LT above, the target body's */
1169 /*           orientation at ET - LT is used. The surface normal is */
1170 /*           dependent on the target body's orientation, so the body's */
1171 /*           orientation model must be evaluated for the correct epoch. */
1172 
1173 
1174 /*           Target body -- illumination source vector */
1175 /*           ----------------------------------------- */
1176 
1177 /*           The surface features on the target body near SPOINT will */
1178 /*           appear in a measurement made at ET as they were at ET-LT. */
1179 /*           In particular, lighting on the target body is dependent on */
1180 /*           the apparent location of the illumination source as seen */
1181 /*           from the target body at ET-LT. So, a second light time */
1182 /*           correction is used to compute the position of the */
1183 /*           illumination source relative to the surface point. */
1184 
1185 
1186 /*        Stellar aberration corrections */
1187 /*        ============================== */
1188 
1189 /*        Stellar aberration corrections are applied only if */
1190 /*        light time corrections are applied as well. */
1191 
1192 /*           Observer-target surface point body vector */
1193 /*           ----------------------------------------- */
1194 
1195 /*           When stellar aberration correction is performed, the */
1196 /*           direction vector SRFVEC is adjusted so as to point to the */
1197 /*           apparent position of SPOINT: considering SPOINT to be an */
1198 /*           ephemeris object, SRFVEC points from the observer's */
1199 /*           position at ET to the light time and stellar aberration */
1200 /*           corrected position of SPOINT. */
1201 
1202 /*           Target body-illumination source vector */
1203 /*           -------------------------------------- */
1204 
1205 /*           The target body-illumination source vector is the apparent */
1206 /*           position of the illumination source, corrected for light */
1207 /*           time and stellar aberration, as seen from the target body */
1208 /*           at time ET-LT. */
1209 
1210 
1211 
1212 /*     Using DSK data */
1213 /*     ============== */
1214 
1215 /*        DSK loading and unloading */
1216 /*        ------------------------- */
1217 
1218 /*        DSK files providing data used by this routine are loaded by */
1219 /*        calling FURNSH and can be unloaded by calling UNLOAD or */
1220 /*        KCLEAR. See the documentation of FURNSH for limits on numbers */
1221 /*        of loaded DSK files. */
1222 
1223 /*        For run-time efficiency, it's desirable to avoid frequent */
1224 /*        loading and unloading of DSK files. When there is a reason to */
1225 /*        use multiple versions of data for a given target body---for */
1226 /*        example, if topographic data at varying resolutions are to be */
1227 /*        used---the surface list can be used to select DSK data to be */
1228 /*        used for a given computation. It is not necessary to unload */
1229 /*        the data that are not to be used. This recommendation presumes */
1230 /*        that DSKs containing different versions of surface data for a */
1231 /*        given body have different surface ID codes. */
1232 
1233 
1234 /*        DSK data priority */
1235 /*        ----------------- */
1236 
1237 /*        A DSK coverage overlap occurs when two segments in loaded DSK */
1238 /*        files cover part or all of the same domain---for example, a */
1239 /*        given longitude-latitude rectangle---and when the time */
1240 /*        intervals of the segments overlap as well. */
1241 
1242 /*        When DSK data selection is prioritized, in case of a coverage */
1243 /*        overlap, if the two competing segments are in different DSK */
1244 /*        files, the segment in the DSK file loaded last takes */
1245 /*        precedence. If the two segments are in the same file, the */
1246 /*        segment located closer to the end of the file takes */
1247 /*        precedence. */
1248 
1249 /*        When DSK data selection is unprioritized, data from competing */
1250 /*        segments are combined. For example, if two competing segments */
1251 /*        both represent a surface as sets of triangular plates, the */
1252 /*        union of those sets of plates is considered to represent the */
1253 /*        surface. */
1254 
1255 /*        Currently only unprioritized data selection is supported. */
1256 /*        Because prioritized data selection may be the default behavior */
1257 /*        in a later version of the routine, the UNPRIORITIZED keyword is */
1258 /*        required in the METHOD argument. */
1259 
1260 
1261 /*        Syntax of the METHOD input argument */
1262 /*        ----------------------------------- */
1263 
1264 /*        The keywords and surface list in the METHOD argument */
1265 /*        are called "clauses." The clauses may appear in any */
1266 /*        order, for example */
1267 
1268 /*           DSK/<surface list>/UNPRIORITIZED */
1269 /*           DSK/UNPRIORITIZED/<surface list> */
1270 /*           UNPRIORITIZED/<surface list>/DSK */
1271 
1272 /*        The simplest form of the METHOD argument specifying use of */
1273 /*        DSK data is one that lacks a surface list, for example: */
1274 
1275 /*           'DSK/UNPRIORITIZED' */
1276 
1277 /*        For applications in which all loaded DSK data for the target */
1278 /*        body are for a single surface, and there are no competing */
1279 /*        segments, the above string suffices. This is expected to be */
1280 /*        the usual case. */
1281 
1282 /*        When, for the specified target body, there are loaded DSK */
1283 /*        files providing data for multiple surfaces for that body, the */
1284 /*        surfaces to be used by this routine for a given call must be */
1285 /*        specified in a surface list, unless data from all of the */
1286 /*        surfaces are to be used together. */
1287 
1288 /*        The surface list consists of the string */
1289 
1290 /*           SURFACES = */
1291 
1292 /*        followed by a comma-separated list of one or more surface */
1293 /*        identifiers. The identifiers may be names or integer codes in */
1294 /*        string format. For example, suppose we have the surface */
1295 /*        names and corresponding ID codes shown below: */
1296 
1297 /*           Surface Name                              ID code */
1298 /*           ------------                              ------- */
1299 /*           'Mars MEGDR 128 PIXEL/DEG'                1 */
1300 /*           'Mars MEGDR 64 PIXEL/DEG'                 2 */
1301 /*           'Mars_MRO_HIRISE'                         3 */
1302 
1303 /*        If data for all of the above surfaces are loaded, then */
1304 /*        data for surface 1 can be specified by either */
1305 
1306 /*           'SURFACES = 1' */
1307 
1308 /*        or */
1309 
1310 /*           'SURFACES = "Mars MEGDR 128 PIXEL/DEG"' */
1311 
1312 /*        Double quotes are used to delimit the surface name because */
1313 /*        it contains blank characters. */
1314 
1315 /*        To use data for surfaces 2 and 3 together, any */
1316 /*        of the following surface lists could be used: */
1317 
1318 /*           'SURFACES = 2, 3' */
1319 
1320 /*           'SURFACES = "Mars MEGDR  64 PIXEL/DEG", 3' */
1321 
1322 /*           'SURFACES = 2, Mars_MRO_HIRISE' */
1323 
1324 /*           'SURFACES = "Mars MEGDR 64 PIXEL/DEG", Mars_MRO_HIRISE' */
1325 
1326 /*        An example of a METHOD argument that could be constructed */
1327 /*        using one of the surface lists above is */
1328 
1329 /*        'DSK/UNPRIORITIZED/SURFACES = "Mars MEGDR 64 PIXEL/DEG", 3' */
1330 
1331 
1332 /*        Aberration corrections */
1333 /*        ---------------------- */
1334 
1335 /*        For irregularly shaped target bodies, the distance between the */
1336 /*        observer and the nearest surface intercept need not be a */
1337 /*        continuous function of time; hence the one-way light time */
1338 /*        between the intercept and the observer may be discontinuous as */
1339 /*        well. In such cases, the computed light time, which is found */
1340 /*        using an iterative algorithm, may converge slowly or not at */
1341 /*        all. In all cases, the light time computation will terminate, */
1342 /*        but the result may be less accurate than expected. */
1343 
1344 /* $ Examples */
1345 
1346 /*     The numerical results shown for this example may differ across */
1347 /*     platforms. The results depend on the SPICE kernels used as */
1348 /*     input, the compiler and supporting libraries, and the machine */
1349 /*     specific arithmetic implementation. */
1350 
1351 /*     1) Find the phase, solar incidence, and emission angles at the */
1352 /*        sub-solar and sub-spacecraft points on Mars as seen from the */
1353 /*        Mars Global Surveyor spacecraft at a specified UTC time. Use */
1354 /*        light time and stellar aberration corrections. */
1355 
1356 /*        Use the meta-kernel shown below to load the required SPICE */
1357 /*        kernels. */
1358 
1359 
1360 /*           KPL/MK */
1361 
1362 /*           File: illumf_ex1.tm */
1363 
1364 /*           This meta-kernel is intended to support operation of SPICE */
1365 /*           example programs. The kernels shown here should not be */
1366 /*           assumed to contain adequate or correct versions of data */
1367 /*           required by SPICE-based user applications. */
1368 
1369 /*           In order for an application to use this meta-kernel, the */
1370 /*           kernels referenced here must be present in the user's */
1371 /*           current working directory. */
1372 
1373 /*           The names and contents of the kernels referenced */
1374 /*           by this meta-kernel are as follows: */
1375 
1376 /*              File name                        Contents */
1377 /*              ---------                        -------- */
1378 /*              de430.bsp                        Planetary ephemeris */
1379 /*              mar097.bsp                       Mars satellite ephemeris */
1380 /*              pck00010.tpc                     Planet orientation and */
1381 /*                                               radii */
1382 /*              naif0011.tls                     Leapseconds */
1383 /*              mgs_ext12_ipng_mgs95j.bsp        MGS ephemeris */
1384 /*              megr90n000cb_plate.bds           Plate model based on */
1385 /*                                               MEGDR DEM, resolution */
1386 /*                                               4 pixels/degree. */
1387 
1388 /*           \begindata */
1389 
1390 /*              KERNELS_TO_LOAD = ( 'de430.bsp', */
1391 /*                                  'mar097.bsp', */
1392 /*                                  'pck00010.tpc', */
1393 /*                                  'naif0011.tls', */
1394 /*                                  'mgs_ext12_ipng_mgs95j.bsp', */
1395 /*                                  'megr90n000cb_plate.bds'      ) */
1396 /*           \begintext */
1397 
1398 
1399 
1400 /*        Example code begins here. */
1401 
1402 
1403 /*           PROGRAM EX1 */
1404 /*           IMPLICIT NONE */
1405 /*     C */
1406 /*     C     SPICELIB functions */
1407 /*     C */
1408 /*           DOUBLE PRECISION      DPR */
1409 /*     C */
1410 /*     C     Local parameters */
1411 /*     C */
1412 /*           CHARACTER*(*)         F1 */
1413 /*           PARAMETER           ( F1     = '(A,F15.9)' ) */
1414 
1415 /*           CHARACTER*(*)         F2 */
1416 /*           PARAMETER           ( F2     = '(A)' ) */
1417 
1418 /*           CHARACTER*(*)         F3 */
1419 /*           PARAMETER           ( F3     = '(A,2(2X,L))' ) */
1420 
1421 /*           CHARACTER*(*)         META */
1422 /*           PARAMETER           ( META   = 'illumf_ex1.tm' ) */
1423 
1424 /*           INTEGER               NAMLEN */
1425 /*           PARAMETER           ( NAMLEN = 32 ) */
1426 
1427 /*           INTEGER               TIMLEN */
1428 /*           PARAMETER           ( TIMLEN = 25 ) */
1429 
1430 /*           INTEGER               CORLEN */
1431 /*           PARAMETER           ( CORLEN = 5 ) */
1432 
1433 /*           INTEGER               MTHLEN */
1434 /*           PARAMETER           ( MTHLEN = 50 ) */
1435 
1436 /*           INTEGER               NMETH */
1437 /*           PARAMETER           ( NMETH  = 2 ) */
1438 /*     C */
1439 /*     C     Local variables */
1440 /*     C */
1441 /*           CHARACTER*(CORLEN)    ABCORR */
1442 /*           CHARACTER*(NAMLEN)    FIXREF */
1443 /*           CHARACTER*(MTHLEN)    ILUMTH ( NMETH ) */
1444 /*           CHARACTER*(NAMLEN)    OBSRVR */
1445 /*           CHARACTER*(MTHLEN)    SUBMTH ( NMETH ) */
1446 /*           CHARACTER*(NAMLEN)    TARGET */
1447 /*           CHARACTER*(TIMLEN)    UTC */
1448 
1449 /*           DOUBLE PRECISION      ET */
1450 /*           DOUBLE PRECISION      SRFVEC ( 3 ) */
1451 /*           DOUBLE PRECISION      SSCEMI */
1452 /*           DOUBLE PRECISION      SSCPHS */
1453 /*           DOUBLE PRECISION      SSCPT  ( 3 ) */
1454 /*           DOUBLE PRECISION      SSCSOL */
1455 /*           DOUBLE PRECISION      SSLEMI */
1456 /*           DOUBLE PRECISION      SSLPHS */
1457 /*           DOUBLE PRECISION      SSLSOL */
1458 /*           DOUBLE PRECISION      SSOLPT ( 3 ) */
1459 /*           DOUBLE PRECISION      TRGEPC */
1460 
1461 /*           INTEGER               I */
1462 
1463 /*           LOGICAL               SSCLIT */
1464 /*           LOGICAL               SSCVIS */
1465 /*           LOGICAL               SSLLIT */
1466 /*           LOGICAL               SSLVIS */
1467 
1468 /*     C */
1469 /*     C     Initial values */
1470 /*     C */
1471 /*           DATA                  ILUMTH / 'Ellipsoid', */
1472 /*          .                               'DSK/Unprioritized' / */
1473 
1474 /*           DATA                  SUBMTH / 'Near Point/Ellipsoid', */
1475 /*          .                            'DSK/Nadir/Unprioritized' / */
1476 
1477 /*     C */
1478 /*     C     Load kernel files. */
1479 /*     C */
1480 /*           CALL FURNSH ( META ) */
1481 /*     C */
1482 /*     C     Convert the UTC request time string to seconds past */
1483 /*     C     J2000 TDB. */
1484 /*     C */
1485 /*           UTC = '2003 OCT 13 06:00:00 UTC' */
1486 
1487 /*           CALL UTC2ET ( UTC, ET ) */
1488 
1489 /*           WRITE (*,F2) ' ' */
1490 /*           WRITE (*,F2) 'UTC epoch is '//UTC */
1491 /*     C */
1492 /*     C     Assign observer and target names. The acronym MGS */
1493 /*     C     indicates Mars Global Surveyor. See NAIF_IDS for a */
1494 /*     C     list of names recognized by SPICE. Also set the */
1495 /*     C     aberration correction flag. */
1496 /*     C */
1497 /*           TARGET = 'Mars' */
1498 /*           OBSRVR = 'MGS' */
1499 /*           FIXREF = 'IAU_MARS' */
1500 /*           ABCORR = 'CN+S' */
1501 
1502 /*           DO I = 1, NMETH */
1503 /*     C */
1504 /*     C        Find the sub-solar point on Mars as */
1505 /*     C        seen from the MGS spacecraft at ET. Use the */
1506 /*     C        "near point" style of sub-point definition */
1507 /*     C        when the shape model is an ellipsoid, and use */
1508 /*     C        the "nadir" style when the shape model is */
1509 /*     C        provided by DSK data. This makes it easy to */
1510 /*     C        verify the solar incidence angle when */
1511 /*     C        the target is modeled as an  ellipsoid. */
1512 /*     C */
1513 /*              CALL SUBSLR ( SUBMTH(I),  TARGET,  ET, */
1514 /*          .                 FIXREF,     ABCORR,  OBSRVR, */
1515 /*          .                 SSOLPT,     TRGEPC,  SRFVEC  ) */
1516 /*     C */
1517 /*     C        Now find the sub-spacecraft point. */
1518 /*     C */
1519 /*              CALL SUBPNT ( SUBMTH(I),  TARGET,  ET, */
1520 /*          .                 FIXREF,     ABCORR,  OBSRVR, */
1521 /*          .                 SSCPT,      TRGEPC,  SRFVEC ) */
1522 /*     C */
1523 /*     C        Find the phase, solar incidence, and emission */
1524 /*     C        angles at the sub-solar point on Mars as */
1525 /*     C        seen from MGS at time ET. */
1526 /*     C */
1527 /*              CALL ILLUMF ( ILUMTH(I), TARGET,  'SUN', */
1528 /*          .                 ET,        FIXREF,  ABCORR, */
1529 /*          .                 OBSRVR,    SSOLPT,  TRGEPC, */
1530 /*          .                 SRFVEC,    SSLPHS,  SSLSOL, */
1531 /*          .                 SSLEMI,    SSLVIS,  SSLLIT ) */
1532 /*     C */
1533 /*     C        Do the same for the sub-spacecraft point. */
1534 /*     C */
1535 /*              CALL ILLUMF ( ILUMTH(I), TARGET,  'SUN', */
1536 /*          .                 ET,        FIXREF,  ABCORR, */
1537 /*          .                 OBSRVR,    SSCPT,   TRGEPC, */
1538 /*          .                 SRFVEC,    SSCPHS,  SSCSOL, */
1539 /*          .                 SSCEMI,    SSCVIS,  SSCLIT  ) */
1540 /*     C */
1541 /*     C        Convert the angles to degrees and write them out. */
1542 /*     C */
1543 /*              SSLPHS = DPR() * SSLPHS */
1544 /*              SSLSOL = DPR() * SSLSOL */
1545 /*              SSLEMI = DPR() * SSLEMI */
1546 
1547 /*              SSCPHS = DPR() * SSCPHS */
1548 /*              SSCSOL = DPR() * SSCSOL */
1549 /*              SSCEMI = DPR() * SSCEMI */
1550 
1551 /*              WRITE (*,F2) ' ' */
1552 /*              WRITE (*,F2) '   ILLUMF method: '//ILUMTH(I) */
1553 /*              WRITE (*,F2) '   SUBPNT method: '//SUBMTH(I) */
1554 /*              WRITE (*,F2) '   SUBSLR method: '//SUBMTH(I) */
1555 /*              WRITE (*,F2) ' ' */
1556 /*              WRITE (*,F2) '      Illumination angles at the ' */
1557 /*          .   //           'sub-solar point:' */
1558 /*              WRITE (*,F2) ' ' */
1559 
1560 /*              WRITE (*,F1) '      Phase angle           (deg.): ', */
1561 /*          .                SSLPHS */
1562 /*              WRITE (*,F1) '      Solar incidence angle (deg.): ', */
1563 /*          .                SSLSOL */
1564 /*              WRITE (*,F1) '      Emission angle        (deg.): ', */
1565 /*          .                SSLEMI */
1566 /*              WRITE (*,F3) '      Visible, Lit flags:           ', */
1567 /*          .                SSLVIS, SSLLIT */
1568 /*              WRITE (*,F2) ' ' */
1569 
1570 /*              IF ( I .EQ. 1 ) THEN */
1571 /*                 WRITE (*,F2) '        The solar incidence angle ' */
1572 /*          .      //           'should be 0.' */
1573 /*                 WRITE (*,F2) '        The emission and phase ' */
1574 /*          .      //           'angles should be equal.' */
1575 /*                 WRITE (*,F2) ' ' */
1576 /*              END IF */
1577 
1578 
1579 /*              WRITE (*,F2) '      Illumination angles at the ' */
1580 /*          .   //          'sub-s/c point:' */
1581 /*              WRITE (*,F2) ' ' */
1582 /*              WRITE (*,F1) '      Phase angle           (deg.): ', */
1583 /*          .               SSCPHS */
1584 /*              WRITE (*,F1) '      Solar incidence angle (deg.): ', */
1585 /*          .               SSCSOL */
1586 /*              WRITE (*,F1) '      Emission angle        (deg.): ', */
1587 /*          .               SSCEMI */
1588 /*              WRITE (*,F3) '      Visible, Lit flags:           ', */
1589 /*          .                SSCVIS, SSCLIT */
1590 /*              WRITE (*,F2) ' ' */
1591 
1592 /*              IF ( I .EQ. 1 ) THEN */
1593 /*                 WRITE (*,F2) '        The emission angle ' */
1594 /*          .      //           'should be 0.' */
1595 /*                 WRITE (*,F2) '        The solar incidence ' */
1596 /*          .      //           'and phase angles should be equal.' */
1597 /*              END IF */
1598 
1599 /*           END DO */
1600 
1601 /*           END */
1602 
1603 
1604 /*     When this program was executed on a PC/Linux/gfortran 64-bit */
1605 /*     platform, the output was: */
1606 
1607 
1608 /*        UTC epoch is 2003 OCT 13 06:00:00 UTC */
1609 
1610 /*           ILLUMF method: Ellipsoid */
1611 /*           SUBPNT method: Near Point/Ellipsoid */
1612 /*           SUBSLR method: Near Point/Ellipsoid */
1613 
1614 /*              Illumination angles at the sub-solar point: */
1615 
1616 /*              Phase angle           (deg.):   138.370270685 */
1617 /*              Solar incidence angle (deg.):     0.000000000 */
1618 /*              Emission angle        (deg.):   138.370270685 */
1619 /*              Visible, Lit flags:             F  T */
1620 
1621 /*                The solar incidence angle should be 0. */
1622 /*                The emission and phase angles should be equal. */
1623 
1624 /*              Illumination angles at the sub-s/c point: */
1625 
1626 /*              Phase angle           (deg.):   101.439331040 */
1627 /*              Solar incidence angle (deg.):   101.439331041 */
1628 /*              Emission angle        (deg.):     0.000000002 */
1629 /*              Visible, Lit flags:             T  F */
1630 
1631 /*                The emission angle should be 0. */
1632 /*                The solar incidence and phase angles should be equal. */
1633 
1634 /*           ILLUMF method: DSK/Unprioritized */
1635 /*           SUBPNT method: DSK/Nadir/Unprioritized */
1636 /*           SUBSLR method: DSK/Nadir/Unprioritized */
1637 
1638 /*              Illumination angles at the sub-solar point: */
1639 
1640 /*              Phase angle           (deg.):   138.387071677 */
1641 /*              Solar incidence angle (deg.):     0.967122745 */
1642 /*              Emission angle        (deg.):   137.621480599 */
1643 /*              Visible, Lit flags:             F  T */
1644 
1645 /*              Illumination angles at the sub-s/c point: */
1646 
1647 /*              Phase angle           (deg.):   101.439331359 */
1648 /*              Solar incidence angle (deg.):   101.555993667 */
1649 /*              Emission angle        (deg.):     0.117861156 */
1650 /*              Visible, Lit flags:             T  F */
1651 
1652 
1653 /* $ Restrictions */
1654 
1655 /*     1) Results from this routine are not meaningful if the input point */
1656 /*        lies on a ridge or vertex of a surface represented by DSK data, */
1657 /*        or if for any other reason the direction of the outward normal */
1658 /*        vector at the point is undefined. */
1659 
1660 /*     2) The illumination state indicated by the output argument `lit' */
1661 /*        is computed treating the illumination source as a single */
1662 /*        point. Surface points that are illuminated by part of the */
1663 /*        source are classified as "lit" or not depending on whether the */
1664 /*        center of the source is visible from those points. */
1665 
1666 
1667 /* $ Literature_References */
1668 
1669 /*     None. */
1670 
1671 /* $ Author_and_Institution */
1672 
1673 /*     N.J. Bachman   (JPL) */
1674 /*     B.V. Semenov   (JPL) */
1675 
1676 /* $ Version */
1677 
1678 /* -    SPICELIB Version 2.0.0, 04-APR-2017 (NJB) (BVS) */
1679 
1680 /*       07-APR-2016 (NJB) (BVS) */
1681 
1682 /*        Now uses surface mapping tracking capability. */
1683 
1684 /*        Updated surface ID codes in header comments. */
1685 
1686 /*       30-MAR-2015 (NJB) */
1687 
1688 /*        Now uses illumination angles to determine whether */
1689 /*        self-intersection tests are necessary, for the DSK */
1690 /*        case. Now imports SHPLEN parameter from gf.inc. */
1691 
1692 /*        Original version 09-FEB-2015 (NJB) (BVS) */
1693 
1694 /* -& */
1695 /* $ Index_Entries */
1696 
1697 /*     illumination angles general source with flags */
1698 /*     lighting angles general source with flags */
1699 /*     phase angle general source with flags */
1700 /*     incidence angle general source with flags */
1701 /*     emission angle general source with flags */
1702 
1703 /* -& */
1704 /* $ Revisions */
1705 
1706 /*     None. */
1707 
1708 /* -& */
1709 
1710 /*     SPICELIB functions */
1711 
1712 
1713 /*     Local parameters */
1714 
1715 
1716 /*     Saved body name length. */
1717 
1718 
1719 /*     Saved frame name length. */
1720 
1721 
1722 /*     Local variables */
1723 
1724 
1725 /*     Saved name/ID item declarations. */
1726 
1727 
1728 /*     Saved frame name/ID item declarations. */
1729 
1730 
1731 /*     Saved surface name/ID item declarations. */
1732 
1733 
1734 /*     Saved variables */
1735 
1736 
1737 /*     Saved name/ID items. */
1738 
1739 
1740 /*     Saved frame name/ID items. */
1741 
1742 
1743 /*     Saved surface name/ID items. */
1744 
1745 
1746 /*     Note: XMIT need not be saved, since it's used only */
1747 /*     for error checking when an aberration correction flag */
1748 /*     is parsed. */
1749 
1750 
1751 /*     Initial values */
1752 
1753 
1754 /*     Standard SPICE error handling. */
1755 
1756     if (return_()) {
1757 	return 0;
1758     }
1759     chkin_("ILLUMF", (ftnlen)6);
1760 
1761 /*     Counter initialization is done separately. */
1762 
1763     if (first) {
1764 
1765 /*        Initialize counters. */
1766 
1767 	zzctruin_(svctr1);
1768 	zzctruin_(svctr2);
1769 	zzctruin_(svctr3);
1770     }
1771 
1772 /*     If necessary, parse the aberration correction flag. */
1773 
1774     if (first || s_cmp(abcorr, prvcor, abcorr_len, (ftnlen)5) != 0) {
1775 
1776 /*        Make sure the results of this block won't be reused */
1777 /*        if we bail out due to an error. */
1778 
1779 	s_copy(prvcor, " ", (ftnlen)5, (ftnlen)1);
1780 
1781 /*        The aberration correction flag differs from the value it */
1782 /*        had on the previous call, if any. Analyze the new flag. */
1783 
1784 	zzvalcor_(abcorr, attblk, abcorr_len);
1785 	if (failed_()) {
1786 	    chkout_("ILLUMF", (ftnlen)6);
1787 	    return 0;
1788 	}
1789 /*        Set logical flags indicating the attributes of the requested */
1790 /*        correction: */
1791 
1792 /*           XMIT is .TRUE. when the correction is for transmitted */
1793 /*           radiation. */
1794 
1795 /*           USELT is .TRUE. when any type of light time correction */
1796 /*           (normal or converged Newtonian) is specified. */
1797 
1798 /*           USESTL indicates stellar aberration corrections. */
1799 
1800 
1801 /*        The above definitions are consistent with those used by */
1802 /*        ZZVALCOR. */
1803 
1804 	xmit = attblk[4];
1805 	uselt = attblk[1];
1806 
1807 /*        The aberration correction flag is recognized; save it. */
1808 
1809 	s_copy(prvcor, abcorr, (ftnlen)5, abcorr_len);
1810 
1811 /*        We do NOT set FIRST to .FALSE. here, since we're not */
1812 /*        yet done with it. */
1813 
1814     }
1815 
1816 /*     Get the target ID code here, since it will be needed */
1817 /*     for the initialization calls below. */
1818 
1819     zzbods2c_(svctr1, svtarg, &svtcde, &svfnd1, target, &trgcde, &fnd, (
1820 	    ftnlen)36, target_len);
1821     if (! fnd) {
1822 	setmsg_("The target, '#', is not a recognized name for an ephemeris "
1823 		"object. The cause of this problem may be that you need an up"
1824 		"dated version of the SPICE Toolkit, or that you failed to lo"
1825 		"ad a kernel containing a name-ID mapping for this body.", (
1826 		ftnlen)234);
1827 	errch_("#", target, (ftnlen)1, target_len);
1828 	sigerr_("SPICE(IDCODENOTFOUND)", (ftnlen)21);
1829 	chkout_("ILLUMF", (ftnlen)6);
1830 	return 0;
1831     }
1832 
1833 /*     Check whether the surface name/ID mapping has been updated. */
1834 
1835     zzsrftrk_(svctr4, &surfup);
1836 
1837 /*     If necessary, parse the method specification. PRVMTH */
1838 /*     and the derived flags NEAR and ELIPSD start out with */
1839 /*     valid values. PRVMTH records the last valid value of */
1840 /*     METHOD; ELIPSD is the corresponding shape flag. */
1841 
1842     if (first || surfup || s_cmp(method, prvmth, method_len, (ftnlen)500) !=
1843 	    0) {
1844 
1845 /*        Set the previous method string to an invalid value, so it */
1846 /*        cannot match any future, valid input. This will force this */
1847 /*        routine to parse the input method on the next call if any */
1848 /*        failure occurs in this branch. Once success is assured, we can */
1849 /*        record the current method in the previous method string. */
1850 
1851 	s_copy(prvmth, " ", (ftnlen)500, (ftnlen)1);
1852 
1853 /*        Parse the method string. If the string is valid, the */
1854 /*        outputs SHAPE and SUBTYP will always be be set. However, */
1855 /*        SUBTYP is not used in this routine. */
1856 
1857 /*        For DSK shapes, the surface list array and count will be set */
1858 /*        if the method string contains a surface list. */
1859 
1860 	zzprsmet_(&trgcde, method, &c__100, shpstr, subtyp, &pri, &nsurf,
1861 		srflst, pntdef, trmstr, method_len, (ftnlen)9, (ftnlen)20, (
1862 		ftnlen)20, (ftnlen)20);
1863 	if (failed_()) {
1864 	    chkout_("ILLUMF", (ftnlen)6);
1865 	    return 0;
1866 	}
1867 	if (eqstr_(shpstr, "ELLIPSOID", (ftnlen)9, (ftnlen)9)) {
1868 	    shape = 1;
1869 	} else if (eqstr_(shpstr, "DSK", (ftnlen)9, (ftnlen)3)) {
1870 	    shape = 2;
1871 	} else {
1872 
1873 /*           This is a backstop check. */
1874 
1875 	    setmsg_("Returned shape value from method string was <#>.", (
1876 		    ftnlen)48);
1877 	    errch_("#", shpstr, (ftnlen)1, (ftnlen)9);
1878 	    sigerr_("SPICE(BUG)", (ftnlen)10);
1879 	    chkout_("ILLUMF", (ftnlen)6);
1880 	    return 0;
1881 	}
1882 
1883 /*        There should be no subtype specification in the method */
1884 /*        string. */
1885 
1886 	if (s_cmp(subtyp, " ", (ftnlen)20, (ftnlen)1) != 0) {
1887 	    setmsg_("Spurious sub-observer point type <#> was present in the"
1888 		    " method string #. The sub-observer type is valid in the "
1889 		    "method strings for SUBPNT and SUBSLR, but is not applica"
1890 		    "ble for ILLUMF.", (ftnlen)182);
1891 	    errch_("#", subtyp, (ftnlen)1, (ftnlen)20);
1892 	    errch_("#", method, (ftnlen)1, method_len);
1893 	    sigerr_("SPICE(INVALIDMETHOD)", (ftnlen)20);
1894 	    chkout_("ILLUMF", (ftnlen)6);
1895 	    return 0;
1896 	}
1897 	s_copy(prvmth, method, (ftnlen)500, method_len);
1898     }
1899 
1900 /*     We're done with all tasks that must be executed on the first */
1901 /*     pass. */
1902 
1903     first = FALSE_;
1904 
1905 /*     Obtain integer codes for the observer and */
1906 /*     illumination source. */
1907 
1908     zzbods2c_(svctr2, svobsr, &svobsc, &svfnd2, obsrvr, &obscde, &fnd, (
1909 	    ftnlen)36, obsrvr_len);
1910     if (! fnd) {
1911 	setmsg_("The observer, '#', is not a recognized name for an ephemeri"
1912 		"s object. The cause of this problem may be that you need an "
1913 		"updated version of the SPICE Toolkit, or that you failed to "
1914 		"load a kernel containing a name-ID mapping for this body.", (
1915 		ftnlen)236);
1916 	errch_("#", obsrvr, (ftnlen)1, obsrvr_len);
1917 	sigerr_("SPICE(IDCODENOTFOUND)", (ftnlen)21);
1918 	chkout_("ILLUMF", (ftnlen)6);
1919 	return 0;
1920     }
1921 
1922 /*     Check the observer and target body codes. If they are equal, */
1923 /*     signal an error. */
1924 
1925     if (obscde == trgcde) {
1926 	setmsg_("In computing illumination angles, the observing body and ta"
1927 		"rget body are the same. Both are #.", (ftnlen)94);
1928 	errch_("#", obsrvr, (ftnlen)1, obsrvr_len);
1929 	sigerr_("SPICE(BODIESNOTDISTINCT)", (ftnlen)24);
1930 	chkout_("ILLUMF", (ftnlen)6);
1931 	return 0;
1932     }
1933 
1934 /*     Determine the attributes of the frame designated by FIXREF. */
1935 
1936     zznamfrm_(svctr3, svfref, &svrefc, fixref, &fixfid, (ftnlen)32,
1937 	    fixref_len);
1938     frinfo_(&fixfid, &center, &type__, &typeid, &fnd);
1939     if (failed_()) {
1940 	chkout_("ILLUMF", (ftnlen)6);
1941 	return 0;
1942     }
1943     if (! fnd) {
1944 	setmsg_("Reference frame # is not recognized by the SPICE frame subs"
1945 		"ystem. Possibly a required frame definition kernel has not b"
1946 		"een loaded.", (ftnlen)130);
1947 	errch_("#", fixref, (ftnlen)1, fixref_len);
1948 	sigerr_("SPICE(NOFRAME)", (ftnlen)14);
1949 	chkout_("ILLUMF", (ftnlen)6);
1950 	return 0;
1951     }
1952 
1953 /*     Make sure that FIXREF is centered at the target body's center. */
1954 
1955     if (center != trgcde) {
1956 	setmsg_("Reference frame # is not centered at the target body #. The"
1957 		" ID code of the frame center is #.", (ftnlen)93);
1958 	errch_("#", fixref, (ftnlen)1, fixref_len);
1959 	errch_("#", target, (ftnlen)1, target_len);
1960 	errint_("#", &center, (ftnlen)1);
1961 	sigerr_("SPICE(INVALIDFRAME)", (ftnlen)19);
1962 	chkout_("ILLUMF", (ftnlen)6);
1963 	return 0;
1964     }
1965 
1966 /*     Get the sign S prefixing LT in the expression for TRGEPC. */
1967 /*     When light time correction is not used, setting S = 0 */
1968 /*     allows us to seamlessly set TRGEPC equal to ET. */
1969 
1970     if (uselt) {
1971 	if (xmit) {
1972 	    s = 1.;
1973 	} else {
1974 	    s = -1.;
1975 	}
1976     } else {
1977 	s = 0.;
1978     }
1979 
1980 /*     Look up the state of the surface point relative to the observer. */
1981 /*     The body-fixed frame of the surface point is to be evaluated */
1982 /*     at the epoch of the surface point, not at the epoch of the */
1983 /*     center of the frame; we indicate this by setting the input */
1984 /*     argument REFLOC to 'TARGET'. */
1985 
1986     spkcpt_(spoint, target, fixref, et, fixref, "TARGET", abcorr, obsrvr,
1987 	    opstat, &lt, target_len, fixref_len, fixref_len, (ftnlen)6,
1988 	    abcorr_len, obsrvr_len);
1989     if (failed_()) {
1990 	chkout_("ILLUMF", (ftnlen)6);
1991 	return 0;
1992     }
1993 
1994 /*     TRGEPC is the epoch associated with the surface point. Below, */
1995 /*     since S is set to 0.D0 if no aberration corrections are used, we */
1996 /*     require no logical branch. */
1997 
1998     *trgepc = *et + s * lt;
1999 
2000 /*     Now find the state of the illumination source as seen by */
2001 /*     the surface point at TRGEPC. We want to evaluate the orientation */
2002 /*     of the body-fixed frame of the surface point at the epoch */
2003 /*     associated with the surface point, not at the epoch associated */
2004 /*     with the frame's center;  we indicate this by setting the input */
2005 /*     argument REFLOC to 'OBSERVER'. */
2006 
2007     spkcpo_(ilusrc, trgepc, fixref, "OBSERVER", abcorr, spoint, target,
2008 	    fixref, tistat, &lti, ilusrc_len, fixref_len, (ftnlen)8,
2009 	    abcorr_len, target_len, fixref_len);
2010     if (failed_()) {
2011 	chkout_("ILLUMF", (ftnlen)6);
2012 	return 0;
2013     }
2014 
2015 /*     We'll need the negative of the observer-surface point position */
2016 /*     for the following angle computation. Set the output SRFVEC while */
2017 /*     we're at it. */
2018 
2019     vequ_(opstat, srfvec);
2020     vminus_(srfvec, obspos);
2021 
2022 /*     Find the surface normal at SPOINT. This computation depends */
2023 /*     on target body shape model. */
2024 
2025     if (shape == 1) {
2026 
2027 /*        We'll need the radii of the target body. */
2028 
2029 	bodvcd_(&trgcde, "RADII", &c__3, &n, radii, (ftnlen)5);
2030 	surfnm_(radii, &radii[1], &radii[2], spoint, normal);
2031 	if (failed_()) {
2032 	    chkout_("ILLUMF", (ftnlen)6);
2033 	    return 0;
2034 	}
2035     } else if (shape == 2) {
2036 
2037 /*        Initialize the normal vector algorithm to use a DSK model */
2038 /*        for the surface of the target body. This initialization is */
2039 /*        required to enable later use of ZZRAYSFX and ZZMAXRAD. */
2040 
2041 	svnsrf = 0;
2042 	zzsudski_(&trgcde, &svnsrf, srflst, &fixfid);
2043 
2044 /*        Compute the outward unit normal at SPOINT on the surface */
2045 /*        defined by the designated DSK data. */
2046 
2047 	zzsbfnrm_(&trgcde, &svnsrf, srflst, trgepc, &fixfid, spoint, normal);
2048 	if (failed_()) {
2049 	    chkout_("ILLUMF", (ftnlen)6);
2050 	    return 0;
2051 	}
2052 	vhatip_(normal);
2053     } else {
2054 
2055 /*        We've already checked the computation method input argument, */
2056 /*        so we don't expect to arrive here. This code is present for */
2057 /*        safety. */
2058 
2059 	setmsg_("The computation method # was not recognized. ", (ftnlen)45);
2060 	errch_("#", method, (ftnlen)1, method_len);
2061 	sigerr_("SPICE(INVALIDMETHOD)", (ftnlen)20);
2062 	chkout_("ILLUMF", (ftnlen)6);
2063 	return 0;
2064     }
2065 
2066 /*     Check for errors before calling math routines. */
2067 
2068 
2069 /*     Find the illumination angles. VSEP will give us angular */
2070 /*     separation in radians. */
2071 
2072     *phase = vsep_(obspos, tistat);
2073     *incdnc = vsep_(normal, tistat);
2074     *emissn = vsep_(normal, obspos);
2075 
2076 /*     Set the visibility and illumination flags. */
2077 
2078 
2079 /*     Set default values of VISIBL and LIT. */
2080 
2081     *visibl = *emissn <= halfpi_();
2082     *lit = *incdnc <= halfpi_();
2083     if (shape == 2) {
2084 
2085 /*        There is a possibility that the surface-observer vector */
2086 /*        or the surface-illumination source vector intersects the */
2087 /*        surface. This is possible only if these vectors have */
2088 /*        positive elevation. */
2089 	if (*lit || *visibl) {
2090 
2091 /*           We need to check for self-occultation of at least one of */
2092 /*           the lines-of-sight from surface to observer or illumination */
2093 /*           source. */
2094 
2095 /*           We'll produce a point slightly above the surface point, */
2096 /*           from which the visibility of the observer and light source */
2097 /*           can be determined. We want to avoid detection of spurious */
2098 /*           intersections near the input surface point. */
2099 
2100 /*           Obtain an upper bound on the target body radius. */
2101 
2102 	    zzmaxrad_(&maxrad);
2103 	    if (failed_()) {
2104 		chkout_("ILLUMF", (ftnlen)6);
2105 		return 0;
2106 	    }
2107 
2108 /*           Compute an offset that can be used to produce a point */
2109 /*           slightly above the surface point, and compute the "raised" */
2110 /*           point. NORMAL is a unit vector here. */
2111 
2112 	    scale = maxrad * 1e-10;
2113 	    vlcom_(&c_b56, spoint, &scale, normal, vertex);
2114 	    if (*visibl) {
2115 
2116 /*              Find the surface intercept, if any, of a ray emanating */
2117 /*              from VERTEX and pointing toward the observer. */
2118 
2119 		zzraysfx_(vertex, obspos, trgepc, xpt, &found);
2120 		*visibl = ! found;
2121 	    }
2122 	    if (*lit) {
2123 
2124 /*              Find the surface intercept, if any, of a ray emanating */
2125 /*              from VERTEX and pointing toward the illumination source. */
2126 
2127 		zzraysfx_(vertex, tistat, trgepc, xpt, &found);
2128 		*lit = ! found;
2129 	    }
2130 	}
2131     }
2132 
2133 /*     TRGEPC and SRVFEC were already set before the illumination */
2134 /*     angle computation. */
2135 
2136     chkout_("ILLUMF", (ftnlen)6);
2137     return 0;
2138 } /* illumf_ */
2139 
2140