1 /* surfnm.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__2 = 2;
11
12 /* $Procedure SURFNM ( Surface normal vector on an ellipsoid ) */
surfnm_(doublereal * a,doublereal * b,doublereal * c__,doublereal * point,doublereal * normal)13 /* Subroutine */ int surfnm_(doublereal *a, doublereal *b, doublereal *c__,
14 doublereal *point, doublereal *normal)
15 {
16 /* Initialized data */
17
18 static char mssg[32*7] = "Axis A was nonpositive. " "Axis B was "
19 "nonpositive. " "Axes A and B were nonpositive. " "Axis "
20 "C was nonpositive. " "Axes A and C were nonpositive. "
21 "Axes B and C were nonpositive. " "All three axes were nonposit"
22 "ive.";
23
24 /* System generated locals */
25 address a__1[2];
26 integer i__1, i__2[2];
27 doublereal d__1;
28 char ch__1[35];
29
30 /* Builtin functions */
31 integer s_rnge(char *, integer, char *, integer);
32 /* Subroutine */ int s_cat(char *, char **, integer *, integer *, ftnlen);
33
34 /* Local variables */
35 static doublereal m;
36 extern /* Subroutine */ int chkin_(char *, ftnlen), errch_(char *, char *,
37 ftnlen, ftnlen), errdp_(char *, doublereal *, ftnlen);
38 static doublereal a1, b1, c1;
39 extern /* Subroutine */ int sigerr_(char *, ftnlen), vhatip_(doublereal *)
40 , chkout_(char *, ftnlen), setmsg_(char *, ftnlen);
41 extern logical return_(void);
42 static integer bad;
43
44 /* $ Abstract */
45
46 /* This routine computes the outward-pointing, unit normal vector */
47 /* from a point on the surface of an ellipsoid. */
48
49 /* $ Disclaimer */
50
51 /* THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE */
52 /* CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. */
53 /* GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE */
54 /* ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE */
55 /* PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" */
56 /* TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY */
57 /* WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A */
58 /* PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC */
59 /* SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE */
60 /* SOFTWARE AND RELATED MATERIALS, HOWEVER USED. */
61
62 /* IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA */
63 /* BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT */
64 /* LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, */
65 /* INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, */
66 /* REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE */
67 /* REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. */
68
69 /* RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF */
70 /* THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY */
71 /* CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE */
72 /* ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. */
73
74 /* $ Required_Reading */
75
76 /* None. */
77
78 /* $ Keywords */
79
80 /* ELLIPSOID, GEOMETRY */
81
82 /* $ Declarations */
83 /* $ Brief_I/O */
84
85 /* VARIABLE I/O DESCRIPTION */
86 /* -------- --- -------------------------------------------------- */
87 /* A I Length of the ellipsoid semi-axis along the x-axis. */
88 /* B I Length of the ellipsoid semi-axis along the y-axis. */
89 /* C I Length of the ellipsoid semi-axis along the z-axis. */
90 /* POINT I Body-fixed coordinates of a point on the ellipsoid. */
91 /* NORMAL O Outward pointing unit normal to ellipsoid at POINT. */
92
93 /* $ Detailed_Input */
94
95 /* A This is the length of the semi-axis of the ellipsoid */
96 /* that is parallel to the x-axis of the body-fixed */
97 /* coordinate system. */
98
99 /* B This is the length of the semi-axis of the ellipsoid */
100 /* that is parallel to the y-axis of the body-fixed */
101 /* coordinate system. */
102
103 /* C This is the length of the semi-axis of the ellipsoid */
104 /* that is parallel to the z-axis of the body-fixed */
105 /* coordinate system. */
106
107 /* POINT This is a 3-vector giving the bodyfixed coordinates */
108 /* of a point on the ellipsoid. In bodyfixed coordinates, */
109 /* the semi-axes of the ellipsoid are aligned with the */
110 /* x, y, and z-axes of the coordinate system. */
111
112 /* $ Detailed_Output */
113
114 /* NORMAL A unit vector pointing away from the ellipsoid and */
115 /* normal to the ellipsoid at POINT. */
116
117 /* $ Parameters */
118
119 /* None. */
120
121 /* $ Exceptions */
122
123 /* 1) If any of the axes are non-positive, the error */
124 /* SPICE(BADAXISLENGTH) will be signaled. */
125
126 /* $ Files */
127
128 /* None. */
129
130 /* $ Particulars */
131
132 /* This routine computes the outward pointing unit normal vector to */
133 /* the ellipsoid having semi-axes of length A, B, and C from the */
134 /* point POINT. */
135
136 /* $ Examples */
137
138 /* A typical use of SURFNM would be to find the angle of incidence */
139 /* of the light from the sun at a point on the surface of an */
140 /* ellipsoid. */
141
142 /* Let Q be a 3-vector representing the rectangular body-fixed */
143 /* coordinates of a point on the ellipsoid (we are assuming that */
144 /* the axes of the ellipsoid are aligned with the axes of the */
145 /* body fixed frame.) Let V be the vector from Q to the sun in */
146 /* bodyfixed coordinates. Then the following code fragment could */
147 /* be used to compute angle of incidence of sunlight at Q. */
148
149 /* CALL SURFNM ( A, B, C, Q, NRML ) */
150
151 /* INCIDN = VSEP ( V, NRML ) */
152
153
154 /* $ Restrictions */
155
156 /* It is assumed that the input POINT is indeed on the ellipsoid. */
157 /* No checking for this is done. */
158
159 /* $ Literature_References */
160
161 /* None. */
162
163 /* $ Author_and_Institution */
164
165 /* N.J. Bachman (JPL) */
166 /* B.V. Semenov (JPL) */
167 /* W.L. Taber (JPL) */
168
169 /* $ Version */
170
171 /* - SPICELIB Version 1.3.2, 23-FEB-2016 (NJB) */
172
173 /* Corrected some typos in the header. */
174
175 /* - SPICELIB Version 1.3.1, 18-MAY-2010 (BVS) */
176
177 /* Removed "C$" marker from text in the header. */
178
179 /* - SPICELIB Version 1.3.0, 02-SEP-2005 (NJB) */
180
181 /* Updated to remove non-standard use of duplicate arguments */
182 /* in VHAT call. */
183
184 /* - SPICELIB Version 1.2.0, 07-AUG-1996 (WLT) */
185
186 /* Added a SAVE statement so that the error message will */
187 /* not be lost between separate invocations of the routine. */
188
189 /* - SPICELIB Version 1.1.0, 21-JUL-1995 (WLT) */
190
191 /* A typo in the Examples section was corrected */
192
193 /* - SPICELIB Version 1.0.1, 10-MAR-1992 (WLT) */
194
195 /* Comment section for permuted index source lines was added */
196 /* following the header. */
197
198 /* - SPICELIB Version 1.0.0, 31-JAN-1990 (WLT) */
199
200 /* -& */
201 /* $ Index_Entries */
202
203 /* surface normal vector on an ellipsoid */
204
205 /* -& */
206 /* $ Revisions */
207
208 /* - SPICELIB Version 1.3.0, 02-SEP-2005 (NJB) */
209
210 /* Updated to remove non-standard use of duplicate arguments */
211 /* in VHAT call. */
212
213 /* - Beta Version 2.0.0, 9-JAN-1989 (WLT) */
214
215 /* Error handling added. */
216
217 /* The algorithm was modified from the initial obvious routine */
218 /* to one that is immune to numerical catastrophes (multiplication */
219 /* or division overflows). */
220
221 /* -& */
222
223 /* SPICELIB Functions */
224
225
226 /* Local Variables */
227
228
229 /* Initial values */
230
231
232 /* Standard SPICE error handling. */
233
234 if (return_()) {
235 return 0;
236 } else {
237 chkin_("SURFNM", (ftnlen)6);
238 }
239
240 /* Check the axes to make sure that none of them is less than or */
241 /* equal to zero. If one is, signal an error and return. */
242
243 bad = 0;
244 if (*a <= 0.) {
245 ++bad;
246 }
247 if (*b <= 0.) {
248 bad += 2;
249 }
250 if (*c__ <= 0.) {
251 bad += 4;
252 }
253 if (bad > 0) {
254 /* Writing concatenation */
255 i__2[0] = 32, a__1[0] = mssg + (((i__1 = bad - 1) < 7 && 0 <= i__1 ?
256 i__1 : s_rnge("mssg", i__1, "surfnm_", (ftnlen)251)) << 5);
257 i__2[1] = 3, a__1[1] = " ? ";
258 s_cat(ch__1, a__1, i__2, &c__2, (ftnlen)35);
259 setmsg_(ch__1, (ftnlen)35);
260 errch_(" ? ", "The A,B, and C axes were #, #, and # respectively.", (
261 ftnlen)3, (ftnlen)50);
262 errdp_("#", a, (ftnlen)1);
263 errdp_("#", b, (ftnlen)1);
264 errdp_("#", c__, (ftnlen)1);
265 sigerr_("SPICE(BADAXISLENGTH)", (ftnlen)20);
266 chkout_("SURFNM", (ftnlen)6);
267 return 0;
268 }
269
270 /* Mathematically we want to compute (Px/a**2, Py/b**2, Pz/c**2) */
271 /* and then convert this to a unit vector. However, computationally */
272 /* this can blow up in our faces. But note that only the ratios */
273 /* a/b, b/c and a/c are important in computing the unit normal. */
274 /* We can use the trick below to avoid the unpleasantness of */
275 /* multiplication and division overflows. */
276
277 /* Computing MIN */
278 d__1 = min(*a,*b);
279 m = min(d__1,*c__);
280
281 /* M can be divided by A,B or C without fear of an overflow */
282 /* occurring. */
283
284 a1 = m / *a;
285 b1 = m / *b;
286 c1 = m / *c__;
287
288 /* All of the terms A1,B1,C1 are less than 1. Thus no overflows */
289 /* can occur. */
290
291 normal[0] = point[0] * (a1 * a1);
292 normal[1] = point[1] * (b1 * b1);
293 normal[2] = point[2] * (c1 * c1);
294 vhatip_(normal);
295 chkout_("SURFNM", (ftnlen)6);
296 return 0;
297 } /* surfnm_ */
298
299