1 /* $Id: ssrfpack.c,v 1.2 2017-10-06 06:41:30 gilles-duvert Exp $
2 * ssrfpack.c: Translated via f2c then massaged so that f2c include and lib
3 * are not required to compile and link the sph supplement.
4 */
5 #ifdef __APPLE__
6 #if MAC_OS_X_VERSION_MIN_REQUIRED >= 1090
7 #define sincos(x, s, c) __sincos(x, s, c)
8 #define sincosf(x, s, c) __sincosf(x, s, c)
9 #else
10 #define sincos(x,s,c) (*s = sin(x), *c = cos(x))
11 #endif
12 #endif
13
14 #ifdef __FreeBSD__
15 # define sincos(x,s,c) (*s = sin(x), *c = cos(x))
16 #endif
17
18 #if defined(__FreeBSD__) || defined(__DragonFly__)
19 #define sincos(x,s,c) *s = sin(x); *c = cos(x)
20 #endif
21
22 #ifdef __cplusplus
23 extern "C" {
24 #endif
25 //#include "f2c.h"
26 typedef DDouble doublereal;
27 typedef DLong integer;
28
29 #define dbg_verbose 0
30
31 namespace ssrfpack {
32
33 /* Table of constant values */
34
35 static doublereal c_b24 = 1.;
36
d_sign(doublereal * a,doublereal * b)37 doublereal d_sign (doublereal *a, doublereal *b)
38 {
39 doublereal x;
40 x = (*a >= 0 ? *a : - *a);
41 return (*b >= 0 ? x : -x);
42 }
43
44 /* Table of constant values */
45
46 static doublereal c_b23 = 1.;
47
aplyr_(doublereal * x,doublereal * y,doublereal * z__,doublereal * cx,doublereal * sx,doublereal * cy,doublereal * sy,doublereal * xp,doublereal * yp,doublereal * zp)48 /* Subroutine */ int aplyr_(doublereal *x, doublereal *y, doublereal *z__,
49 doublereal *cx, doublereal *sx, doublereal *cy, doublereal *sy,
50 doublereal *xp, doublereal *yp, doublereal *zp)
51 {
52 /* Builtin functions */
53 double sqrt(doublereal);
54
55 /* Local variables */
56 static doublereal t;
57
58
59 /* *********************************************************** */
60
61 /* From SSRFPACK */
62 /* Robert J. Renka */
63 /* Dept. of Computer Science */
64 /* Univ. of North Texas */
65 /* renka@cs.unt.edu */
66 /* 05/09/92 */
67
68 /* This subroutine applies the rotation R defined by Sub- */
69 /* routine CONSTR to the unit vector (X Y Z)**T, i,e. (X,Y,Z) */
70 /* is rotated to (XP,YP,ZP). If (XP,YP,ZP) lies in the */
71 /* southern hemisphere (ZP < 0), (XP,YP) are set to the */
72 /* coordinates of the nearest point of the equator, ZP re- */
73 /* maining unchanged. */
74
75 /* On input: */
76
77 /* X,Y,Z = Coordinates of a point on the unit sphere. */
78
79 /* CX,SX,CY,SY = Elements of the rotation defined by */
80 /* Subroutine CONSTR. */
81
82 /* Input parameters are not altered except as noted below. */
83
84 /* On output: */
85
86 /* XP,YP,ZP = Coordinates of the rotated point on the */
87 /* sphere unless ZP < 0, in which case */
88 /* (XP,YP,0) is the closest point of the */
89 /* equator to the rotated point. Storage */
90 /* for XP, YP, and ZP may coincide with */
91 /* storage for X, Y, and Z, respectively, */
92 /* if the latter need not be saved. */
93
94 /* Modules required by APLYR: None */
95
96 /* Intrinsic function called by APLYR: SQRT */
97
98 /* *********************************************************** */
99
100
101 /* Local parameter: */
102
103 /* T = Temporary variable */
104
105 t = *sx * *y + *cx * *z__;
106 *yp = *cx * *y - *sx * *z__;
107 *zp = *sy * *x + *cy * t;
108 *xp = *cy * *x - *sy * t;
109 if (*zp >= 0.) {
110 return 0;
111 }
112
113 /* Move (XP,YP,ZP) to the equator. */
114
115 t = sqrt(*xp * *xp + *yp * *yp);
116 if (t == 0.) {
117 goto L1;
118 }
119 *xp /= t;
120 *yp /= t;
121 return 0;
122
123 /* Move the south pole to an arbitrary point of the equator. */
124
125 L1:
126 *xp = 1.;
127 *yp = 0.;
128 return 0;
129 } /* aplyr_ */
130
aplyrt_(doublereal * g1p,doublereal * g2p,doublereal * cx,doublereal * sx,doublereal * cy,doublereal * sy,doublereal * g)131 /* Subroutine */ int aplyrt_(doublereal *g1p, doublereal *g2p, doublereal *cx,
132 doublereal *sx, doublereal *cy, doublereal *sy, doublereal *g)
133 {
134 static doublereal t;
135
136
137 /* *********************************************************** */
138
139 /* From SSRFPACK */
140 /* Robert J. Renka */
141 /* Dept. of Computer Science */
142 /* Univ. of North Texas */
143 /* renka@cs.unt.edu */
144 /* 05/09/92 */
145
146 /* This subroutine applies the inverse (transpose) of the */
147 /* rotation defined by Subroutine CONSTR to the vector */
148 /* (G1P G2P 0)**T, i.e., the gradient (G1P,G2P,0) in the rot- */
149 /* ated coordinate system is mapped to (G1,G2,G3) in the */
150 /* original coordinate system. */
151
152 /* On input: */
153
154 /* G1P,G2P = X and Y components, respectively, of the */
155 /* gradient in the rotated coordinate system. */
156
157 /* CX,SX,CY,SY = Elements of the rotation R constructed */
158 /* by Subroutine CONSTR. */
159
160 /* Input parameters are not altered by this routine. */
161
162 /* On output: */
163
164 /* G = X, Y, and Z components (in that order) of the */
165 /* inverse rotation applied to (G1P,G2P,0) -- */
166 /* gradient in the original coordinate system. */
167
168 /* Modules required by APLYRT: None */
169
170 /* *********************************************************** */
171
172
173 /* Local parameters: */
174
175 /* T = Temporary variable */
176
177 /* Parameter adjustments */
178 --g;
179
180 /* Function Body */
181 t = *sy * *g1p;
182 g[1] = *cy * *g1p;
183 g[2] = *cx * *g2p - *sx * t;
184 g[3] = -(*sx) * *g2p - *cx * t;
185 return 0;
186 } /* aplyrt_ */
187
arcint_(doublereal * p,doublereal * p1,doublereal * p2,doublereal * f1,doublereal * f2,doublereal * g1,doublereal * g2,doublereal * sigma,doublereal * f,doublereal * g,doublereal * gn)188 /* Subroutine */ int arcint_(doublereal *p, doublereal *p1, doublereal *p2,
189 doublereal *f1, doublereal *f2, doublereal *g1, doublereal *g2,
190 doublereal *sigma, doublereal *f, doublereal *g, doublereal *gn)
191 {
192
193 /* Builtin functions */
194 double sqrt(doublereal), exp(doublereal);
195
196 /* Local variables */
197 static doublereal a, e;
198 static integer i__;
199 static doublereal s, b1, b2, d1, d2, e1, e2, al, cm, gt, sm, tm, un[3],
200 ts, cm2, sb1, sb2, sm2, tm1, tm2, tp1, tp2, cmm, sig, ems, tau1,
201 tau2, sinh__, sinh2, dummy, unorm;
202 extern doublereal arclen_(doublereal *, doublereal *);
203 extern /* Subroutine */ int snhcsh_(doublereal *, doublereal *,
204 doublereal *, doublereal *);
205
206
207 /* *********************************************************** */
208
209 /* From SSRFPACK */
210 /* Robert J. Renka */
211 /* Dept. of Computer Science */
212 /* Univ. of North Texas */
213 /* renka@cs.unt.edu */
214 /* 11/21/96 */
215
216 /* Given 3 points P, P1, and P2 lying on a common geodesic */
217 /* of the unit sphere with P between P1 and P2, along with */
218 /* data values and gradients at P1 and P2, this subroutine */
219 /* computes an interpolated value F and a gradient vector G */
220 /* AT P. F and the tangential component of G are taken to be */
221 /* the value and derivative (with respect to arc-length) of */
222 /* a Hermite interpolatory tension spline defined by the end- */
223 /* point values and tangential gradient components. The nor- */
224 /* mal component of G is obtained by linear interpolation of */
225 /* the normal components of the gradients at P1 and P2. */
226
227 /* On input: */
228
229 /* P = Cartesian coordinates of a point lying on the */
230 /* arc defined by P1 and P2. P(1)**2 + P(2)**2 + */
231 /* P(3)**2 = 1. */
232
233 /* P1,P2 = Coordinates of distinct points on the unit */
234 /* sphere defining an arc with length less than */
235 /* 180 degrees. */
236
237 /* F1,F2 = Data values associated with P1 and P2, */
238 /* respectively. */
239
240 /* G1,G2 = Gradient vectors associated with P1 and P2. */
241 /* G1 and G2 are orthogonal to P1 and P2, */
242 /* respectively. */
243
244 /* SIGMA = Tension factor associated with P1-P2. */
245
246 /* The above parameters are not altered by this routine. */
247
248 /* G = Array of length 3. */
249
250 /* On output: */
251
252 /* F = Interpolated value at P. */
253
254 /* G = Interpolated gradient at P. */
255
256 /* GN = Normal component of G with the direction */
257 /* P1 X P2 taken to be positive. The extrapola- */
258 /* tion procedure requires this component. */
259
260 /* For each vector V, V(1), V(2), and V(3) contain X, Y, */
261 /* and Z components, respectively. */
262
263 /* SSRFPACK modules required by ARCINT: ARCLEN, SNHCSH */
264
265 /* Intrinsic functions called by ARCINT: ABS, EXP, SQRT */
266
267 /* *********************************************************** */
268
269 /* Parameter adjustments */
270 --g;
271 --g2;
272 --g1;
273 --p2;
274 --p1;
275 --p;
276
277 /* Function Body */
278
279 /* Local parameters: */
280
281 /* A = Angle in radians (arc-length) between P1 and */
282 /* P2 */
283 /* AL = Arc-length between P1 and P */
284 /* B1,B2 = Local coordinates of P with respect to P1-P2 */
285 /* CM,CMM = Coshm(SIG) and Coshmm(SIG) -- refer to SNHCSH */
286 /* CM2 = Coshm(SB2) */
287 /* DUMMY = Dummy parameter for SNHCSH */
288 /* D1,D2 = Scaled second differences */
289 /* E = CM**2 - SM*Sinh = SIG*SM - 2*CMM (scaled by */
290 /* 2*EMS if SIG > .5) */
291 /* EMS = Exp(-SIG) */
292 /* E1,E2 = Exp(-SB1), Exp(-SB2) */
293 /* GT = Tangential component of G -- component in the */
294 /* direction UN X P */
295 /* I = DO-loop index */
296 /* LUN = Logical unit for error messages */
297 /* S = Slope: (F2-F1)/A */
298 /* SB1,SB2 = SIG*B1, SIG*B2 */
299 /* SIG = Abs(SIGMA) */
300 /* SINH = Sinh(SIGMA) */
301 /* SINH2 = Sinh(SB2) */
302 /* SM,SM2 = Sinhm(SIG), Sinhm(SB2) */
303 /* TAU1,TAU2 = Tangential derivatives (components of G1,G2) */
304 /* at P1 and P2 */
305 /* TM = 1-EMS */
306 /* TM1,TM2 = 1-E1, 1-E2 */
307 /* TP1,TP2 = 1+E1, 1+E2 */
308 /* TS = TM**2 */
309 /* UN = Unit normal to the plane of P, P1, and P2 */
310 /* UNORM = Euclidean norm of P1 X P2 -- used to normalize */
311 /* UN */
312
313
314 /* Compute unit normal UN. */
315
316 un[0] = p1[2] * p2[3] - p1[3] * p2[2];
317 un[1] = p1[3] * p2[1] - p1[1] * p2[3];
318 un[2] = p1[1] * p2[2] - p1[2] * p2[1];
319 unorm = sqrt(un[0] * un[0] + un[1] * un[1] + un[2] * un[2]);
320 if (unorm == 0.) {
321 goto L2;
322 }
323
324 /* Normalize UN. */
325
326 for (i__ = 1; i__ <= 3; ++i__) {
327 un[i__ - 1] /= unorm;
328 /* L1: */
329 }
330
331 /* Compute tangential derivatives at the endpoints: */
332 /* TAU1 = (G1,UN X P1) = (G1,P2)/UNORM and */
333 /* TAU2 = (G2,UN X P2) = -(G2,P1)/UNORM. */
334
335 tau1 = (g1[1] * p2[1] + g1[2] * p2[2] + g1[3] * p2[3]) / unorm;
336 tau2 = -(g2[1] * p1[1] + g2[2] * p1[2] + g2[3] * p1[3]) / unorm;
337
338 /* Compute arc-lengths A, AL. */
339
340 a = arclen_(&p1[1], &p2[1]);
341 if (a == 0.) {
342 goto L2;
343 }
344 al = arclen_(&p1[1], &p[1]);
345
346 /* Compute local coordinates, slope, and second differences. */
347
348 b2 = al / a;
349 b1 = 1. - b2;
350 s = (*f2 - *f1) / a;
351 d1 = s - tau1;
352 d2 = tau2 - s;
353
354 /* Test the range of SIGMA. */
355
356 sig = fabs(*sigma);
357 if (sig < 1e-9) {
358
359 /* Hermite cubic interpolation. */
360
361 *f = *f1 + al * (tau1 + b2 * (d1 + b1 * (d1 - d2)));
362 gt = tau1 + b2 * (d1 + d2 + b1 * 3. * (d1 - d2));
363 } else if (sig <= .5) {
364
365 /* 0 < SIG .LE. .5. Use approximations designed to avoid */
366 /* cancellation error in the hyperbolic functions. */
367
368 sb2 = sig * b2;
369 snhcsh_(&sig, &sm, &cm, &cmm);
370 snhcsh_(&sb2, &sm2, &cm2, &dummy);
371 sinh__ = sm + sig;
372 sinh2 = sm2 + sb2;
373 e = sig * sm - cmm - cmm;
374 *f = *f1 + al * tau1 + a * ((cm * sm2 - sm * cm2) * (d1 + d2) + sig *
375 (cm * cm2 - sinh__ * sm2) * d1) / (sig * e);
376 gt = tau1 + ((cm * cm2 - sm * sinh2) * (d1 + d2) + sig * (cm * sinh2
377 - sinh__ * cm2) * d1) / e;
378 } else {
379
380 /* SIG > .5. Use negative exponentials in order to avoid */
381 /* overflow. Note that EMS = EXP(-SIG). */
382
383 sb1 = sig * b1;
384 sb2 = sig - sb1;
385 e1 = exp(-sb1);
386 e2 = exp(-sb2);
387 ems = e1 * e2;
388 tm = 1. - ems;
389 ts = tm * tm;
390 tm1 = 1. - e1;
391 tm2 = 1. - e2;
392 e = tm * (sig * (ems + 1.) - tm - tm);
393 *f = *f1 + al * s + a * (tm * tm1 * tm2 * (d1 + d2) + sig * ((e2 *
394 tm1 * tm1 - b1 * ts) * d1 + (e1 * tm2 * tm2 - b2 * ts) * d2))
395 / (sig * e);
396 tp1 = e1 + 1.;
397 tp2 = e2 + 1.;
398 gt = s + (tm1 * (tm * tp2 - sig * e2 * tp1) * d1 - tm2 * (tm * tp1 -
399 sig * e1 * tp2) * d2) / e;
400 }
401
402 /* Compute GN. */
403
404 *gn = b1 * (un[0] * g1[1] + un[1] * g1[2] + un[2] * g1[3]) + b2 * (un[0] *
405 g2[1] + un[1] * g2[2] + un[2] * g2[3]);
406
407 /* Compute G = GT*(UN X P) + GN*UN. */
408
409 g[1] = gt * (un[1] * p[3] - un[2] * p[2]) + *gn * un[0];
410 g[2] = gt * (un[2] * p[1] - un[0] * p[3]) + *gn * un[1];
411 g[3] = gt * (un[0] * p[2] - un[1] * p[1]) + *gn * un[2];
412 return 0;
413
414 /* P1 X P2 = 0. Print an error message and terminate */
415 /* processing. */
416
417 /* 2 WRITE (LUN,100) (P1(I),I=1,3), (P2(I),I=1,3) */
418 /* 100 FORMAT ('1','ERROR IN ARCINT -- P1 = ',2(F9.6,', '), */
419 /* . F9.6/1X,19X,'P2 = ',2(F9.6,', '),F9.6) */
420 if (dbg_verbose) fprintf (stderr, "ERROR IN ARCINT -- P1 = %9.6f %9.6f %9.6f P2 = %9.6f %9.6f %9.6f\n", p1[1], p1[2], p1[3], p2[1], p2[2], p2[3]);
421 L2:
422 i__ = 1;
423 return 0;
424 } /* arcint_ */
425
arclen_(doublereal * p,doublereal * q)426 doublereal arclen_(doublereal *p, doublereal *q)
427 {
428 /* System generated locals */
429 doublereal ret_val, d__1;
430
431 /* Builtin functions */
432 double atan(doublereal), sqrt(doublereal);
433
434 /* Local variables */
435 static doublereal d__;
436 static integer i__;
437
438
439 /* *********************************************************** */
440
441 /* From SSRFPACK */
442 /* Robert J. Renka */
443 /* Dept. of Computer Science */
444 /* Univ. of North Texas */
445 /* renka@cs.unt.edu */
446 /* 05/09/92 */
447
448 /* This function computes the arc-length (angle in radians) */
449 /* between a pair of points on the unit sphere. */
450
451 /* On input: */
452
453 /* P,Q = Arrays of length 3 containing the X, Y, and Z */
454 /* coordinates (in that order) of points on the */
455 /* unit sphere. */
456
457 /* Input parameters are not altered by this function. */
458
459 /* On output: */
460
461 /* ARCLEN = Angle in radians between the unit vectors */
462 /* P and Q. 0 .LE. ARCLEN .LE. PI. */
463
464 /* Modules required by ARCLEN: None */
465
466 /* Intrinsic functions called by ARCLEN: ATAN, SQRT */
467
468 /* *********************************************************** */
469
470
471 /* Local parameters: */
472
473 /* D = Euclidean norm squared of P+Q */
474 /* I = DO-loop index */
475
476 /* Parameter adjustments */
477 --q;
478 --p;
479
480 /* Function Body */
481 d__ = 0.;
482 for (i__ = 1; i__ <= 3; ++i__) {
483 /* Computing 2nd power */
484 d__1 = p[i__] + q[i__];
485 d__ += d__1 * d__1;
486 /* L1: */
487 }
488 if (d__ == 0.) {
489
490 /* P and Q are separated by 180 degrees. */
491
492 ret_val = atan(1.) * 4.;
493 } else if (d__ >= 4.) {
494
495 /* P and Q coincide. */
496
497 ret_val = 0.;
498 } else {
499 ret_val = atan(sqrt((4. - d__) / d__)) * 2.;
500 }
501 return ret_val;
502 } /* arclen_ */
503
constr_(doublereal * xk,doublereal * yk,doublereal * zk,doublereal * cx,doublereal * sx,doublereal * cy,doublereal * sy)504 /* Subroutine */ int constr_(doublereal *xk, doublereal *yk, doublereal *zk,
505 doublereal *cx, doublereal *sx, doublereal *cy, doublereal *sy)
506 {
507 /* Builtin functions */
508 double sqrt(doublereal);
509
510
511 /* *********************************************************** */
512
513 /* From SSRFPACK */
514 /* Robert J. Renka */
515 /* Dept. of Computer Science */
516 /* Univ. of North Texas */
517 /* renka@cs.unt.edu */
518 /* 05/09/92 */
519
520 /* This subroutine constructs the elements of a 3 by 3 */
521 /* orthogonal matrix R which rotates a point (XK,YK,ZK) on */
522 /* the unit sphere to the north pole, i.e., */
523
524 /* (XK) (CY 0 -SY) (1 0 0) (XK) (0) */
525 /* R * (YK) = ( 0 1 0) * (0 CX -SX) * (YK) = (0) */
526 /* (ZK) (SY 0 CY) (0 SX CX) (ZK) (1) */
527
528 /* On input: */
529
530 /* XK,YK,ZK = Components of a unit vector to be */
531 /* rotated to (0,0,1). */
532
533 /* Input parameters are not altered by this routine. */
534
535 /* On output: */
536
537 /* CX,SX,CY,SY = Elements of R: CX,SX define a rota- */
538 /* tion about the X-axis and CY,SY define */
539 /* a rotation about the Y-axis. */
540
541 /* Modules required by CONSTR: None */
542
543 /* Intrinsic function called by CONSTR: SQRT */
544
545 /* *********************************************************** */
546
547 *cy = sqrt(*yk * *yk + *zk * *zk);
548 *sy = *xk;
549 if (*cy != 0.) {
550 *cx = *zk / *cy;
551 *sx = *yk / *cy;
552 } else {
553
554 /* (XK,YK,ZK) lies on the X-axis. */
555
556 *cx = 1.;
557 *sx = 0.;
558 }
559 return 0;
560 } /* constr_ */
561
fval_(doublereal * b1,doublereal * b2,doublereal * b3,doublereal * v1,doublereal * v2,doublereal * v3,doublereal * f1,doublereal * f2,doublereal * f3,doublereal * g1,doublereal * g2,doublereal * g3,doublereal * sig1,doublereal * sig2,doublereal * sig3)562 doublereal fval_(doublereal *b1, doublereal *b2, doublereal *b3, doublereal *
563 v1, doublereal *v2, doublereal *v3, doublereal *f1, doublereal *f2,
564 doublereal *f3, doublereal *g1, doublereal *g2, doublereal *g3,
565 doublereal *sig1, doublereal *sig2, doublereal *sig3)
566 {
567 /* System generated locals */
568 doublereal ret_val;
569
570 /* Builtin functions */
571 double sqrt(doublereal);
572
573 /* Local variables */
574 static doublereal f, g[3];
575 static integer i__;
576 static doublereal c1, c2, c3, q1[3], q2[3], q3[3], s1, s2, s3, u1[3], u2[
577 3], u3[3], ds, dv, u1n, u2n, u3n, sig, val, dum, sum;
578 extern doublereal hval_(doublereal *, doublereal *, doublereal *,
579 doublereal *, doublereal *, doublereal *);
580 extern /* Subroutine */ int arcint_(doublereal *, doublereal *,
581 doublereal *, doublereal *, doublereal *, doublereal *,
582 doublereal *, doublereal *, doublereal *, doublereal *,
583 doublereal *);
584
585
586 /* *********************************************************** */
587
588 /* From SSRFPACK */
589 /* Robert J. Renka */
590 /* Dept. of Computer Science */
591 /* Univ. of North Texas */
592 /* renka@cs.unt.edu */
593 /* 05/09/92 */
594
595 /* Given data values and gradients at the three vertices of */
596 /* a spherical triangle containing a point P, this routine */
597 /* computes the value of F at P where F interpolates the ver- */
598 /* tex data. Along the triangle sides, the interpolatory */
599 /* function F is the Hermite interpolatory tension spline */
600 /* defined by the values and tangential gradient components */
601 /* at the endpoints, and the gradient component normal to the */
602 /* triangle side varies linearly with respect to arc-length */
603 /* between the normal gradient components at the endpoints. */
604 /* A first-order C-1 blending method is used on the underly- */
605 /* ing planar triangle. Since values and gradients on an arc */
606 /* depend only on the vertex data, the method results in C-1 */
607 /* continuity when used to interpolate over a triangulation. */
608
609 /* The blending method consists of taking F(P) to be a */
610 /* weighted sum of the values at PP of three univariate Her- */
611 /* mite interpolatory tension splines defined on the line */
612 /* segments which join the vertices to the opposite sides and */
613 /* pass through PP: the central projection of P onto the */
614 /* underlying planar triangle. The tension factors for these */
615 /* splines are obtained by linear interpolation between the */
616 /* pair of tension factors associated with the triangle sides */
617 /* which join at the appropriate vertex. */
618
619 /* A tension factor SIGMA associated with a Hermite interp- */
620 /* olatory tension spline is a nonnegative parameter which */
621 /* determines the curviness of the spline. SIGMA = 0 results */
622 /* in a cubic spline, and the spline approaches the linear */
623 /* interpolant as SIGMA increases. */
624
625 /* On input: */
626
627 /* B1,B2,B3 = Barycentric coordinates of PP with re- */
628 /* spect to the (planar) underlying triangle */
629 /* (V1,V2,V3), where PP is the central */
630 /* projection of P onto this triangle. */
631
632 /* V1,V2,V3 = Cartesian coordinates of the vertices of */
633 /* a spherical triangle containing P. V3 */
634 /* Left V1->V2. */
635
636 /* F1,F2,F3 = Data values associated with the vertices. */
637
638 /* G1,G2,G3 = Gradients associated with the vertices. */
639 /* Gi is orthogonal to Vi for i = 1,2,3. */
640
641 /* SIG1,SIG2,SIG3 = Tension factors associated with the */
642 /* triangle sides opposite V1, V2, and */
643 /* V3, respectively. */
644
645 /* Input parameters are not altered by this function. */
646
647 /* On output: */
648
649 /* FVAL = Interpolated value at P. */
650
651 /* Each vector V above contains X, Y, and Z components in */
652 /* V(1), V(2), and V(3), respectively. */
653
654 /* SSRFPACK modules required by FVAL: ARCINT, ARCLEN, HVAL */
655
656 /* Intrinsic function called by FVAL: SQRT */
657
658 /* *********************************************************** */
659
660
661 /* Local parameters: */
662
663 /* C1,C2,C3 = Coefficients (weight functions) of partial */
664 /* interpolants. C1 = 1 on the edge opposite */
665 /* V1 and C1 = 0 on the other edges. Simi- */
666 /* larly for C2 and C3. C1+C2+C3 = 1. */
667 /* DS = Directional derivative (scaled by distnace) */
668 /* at U1, U2, or U3: DS = (G,U1-V1)/U1N = */
669 /* -(G,V1)/U1N on side opposite V1, where G/ */
670 /* U1N (plus an orthogonal component) is the */
671 /* projection of G onto the planar triangle */
672 /* DUM = Dummy variable for calls to ARCINT */
673 /* DV = Directional derivatives (scaled by distance) */
674 /* at a vertex: D1 = (G1,U1-V1) = (G1,U1) */
675 /* F,G = Value and gradient at Q1 Q2, or Q3 obtained */
676 /* by interpolation along one of the arcs of */
677 /* the spherical triangle */
678 /* I = DO-loop index */
679 /* Q1,Q2,Q3 = Central projections of U1, U2, and U3 onto */
680 /* the sphere and thus lying on an arc of the */
681 /* spherical triangle */
682 /* SIG = Tension factor for a side-vertex (partial) */
683 /* interpolant: obtained by linear interpo- */
684 /* lation applied to triangle side tensions */
685 /* SUM = Quantity used to normalize C1, C2, and C3 */
686 /* S1,S2,S3 = Sums of pairs of barycentric coordinates: */
687 /* used to compute U1, U2, U3, and SIG */
688 /* U1,U2,U3 = Points on the boundary of the planar trian- */
689 /* gle and lying on the lines containing PP */
690 /* and the vertices. U1 is opposite V1, etc. */
691 /* U1N,U2N,U3N = Quantities used to compute Q1, Q2, and Q3 */
692 /* (magnitudes of U1, U2, and U3) */
693 /* VAL = Local variable used to accumulate the con- */
694 /* tributions to FVAL */
695
696
697 /* Compute weight functions C1, C2, and C3. */
698
699 /* Parameter adjustments */
700 --g3;
701 --g2;
702 --g1;
703 --v3;
704 --v2;
705 --v1;
706
707 /* Function Body */
708 c1 = *b2 * *b3;
709 c2 = *b3 * *b1;
710 c3 = *b1 * *b2;
711 sum = c1 + c2 + c3;
712 if (sum <= 0.) {
713
714 /* P coincides with a vertex. */
715
716 ret_val = *b1 * *f1 + *b2 * *f2 + *b3 * *f3;
717 return ret_val;
718 }
719
720 /* Normalize C1, C2, and C3. */
721
722 c1 /= sum;
723 c2 /= sum;
724 c3 /= sum;
725
726 /* Compute (S1,S2,S3), (U1,U2,U3) and (U1N,U2N,U3N). */
727
728 s1 = *b2 + *b3;
729 s2 = *b3 + *b1;
730 s3 = *b1 + *b2;
731 u1n = 0.;
732 u2n = 0.;
733 u3n = 0.;
734 for (i__ = 1; i__ <= 3; ++i__) {
735 u1[i__ - 1] = (*b2 * v2[i__] + *b3 * v3[i__]) / s1;
736 u2[i__ - 1] = (*b3 * v3[i__] + *b1 * v1[i__]) / s2;
737 u3[i__ - 1] = (*b1 * v1[i__] + *b2 * v2[i__]) / s3;
738 u1n += u1[i__ - 1] * u1[i__ - 1];
739 u2n += u2[i__ - 1] * u2[i__ - 1];
740 u3n += u3[i__ - 1] * u3[i__ - 1];
741 /* L1: */
742 }
743
744 /* Compute Q1, Q2, and Q3. */
745
746 u1n = sqrt(u1n);
747 u2n = sqrt(u2n);
748 u3n = sqrt(u3n);
749 for (i__ = 1; i__ <= 3; ++i__) {
750 q1[i__ - 1] = u1[i__ - 1] / u1n;
751 q2[i__ - 1] = u2[i__ - 1] / u2n;
752 q3[i__ - 1] = u3[i__ - 1] / u3n;
753 /* L2: */
754 }
755
756 /* Compute interpolated value (VAL) at P by looping on */
757 /* triangle sides. */
758
759 val = 0.;
760
761 /* Contribution from side opposite V1: */
762
763 /* Compute value and gradient at Q1 by interpolating */
764 /* between V2 and V3. */
765
766 arcint_(q1, &v2[1], &v3[1], f2, f3, &g2[1], &g3[1], sig1, &f, g, &dum);
767
768 /* Add in the contribution. */
769
770 dv = g1[1] * u1[0] + g1[2] * u1[1] + g1[3] * u1[2];
771 ds = -(g[0] * v1[1] + g[1] * v1[2] + g[2] * v1[3]) / u1n;
772 sig = (*b2 * *sig3 + *b3 * *sig2) / s1;
773 val += c1 * hval_(b1, f1, &f, &dv, &ds, &sig);
774
775 /* Contribution from side opposite V2: */
776
777 /* Compute value and gradient at Q2 by interpolating */
778 /* between V3 and V1. */
779
780 arcint_(q2, &v3[1], &v1[1], f3, f1, &g3[1], &g1[1], sig2, &f, g, &dum);
781
782 /* Add in the contribution. */
783
784 dv = g2[1] * u2[0] + g2[2] * u2[1] + g2[3] * u2[2];
785 ds = -(g[0] * v2[1] + g[1] * v2[2] + g[2] * v2[3]) / u2n;
786 sig = (*b3 * *sig1 + *b1 * *sig3) / s2;
787 val += c2 * hval_(b2, f2, &f, &dv, &ds, &sig);
788
789 /* Contribution from side opposite V3: */
790
791 /* Compute interpolated value and gradient at Q3 */
792 /* by interpolating between V1 and V2. */
793
794 arcint_(q3, &v1[1], &v2[1], f1, f2, &g1[1], &g2[1], sig3, &f, g, &dum);
795
796 /* Add in the final contribution. */
797
798 dv = g3[1] * u3[0] + g3[2] * u3[1] + g3[3] * u3[2];
799 ds = -(g[0] * v3[1] + g[1] * v3[2] + g[2] * v3[3]) / u3n;
800 sig = (*b1 * *sig2 + *b2 * *sig1) / s3;
801 ret_val = val + c3 * hval_(b3, f3, &f, &dv, &ds, &sig);
802 return ret_val;
803 } /* fval_ */
804
getsig_(integer * n,doublereal * x,doublereal * y,doublereal * z__,doublereal * h__,integer * list,integer * lptr,integer * lend,doublereal * grad,doublereal * tol,doublereal * sigma,doublereal * dsmax,integer * ier)805 /* Subroutine */ int getsig_(integer *n, doublereal *x, doublereal *y,
806 doublereal *z__, doublereal *h__, integer *list, integer *lptr,
807 integer *lend, doublereal *grad, doublereal *tol, doublereal *sigma,
808 doublereal *dsmax, integer *ier)
809 {
810 /* Initialized data */
811
812 static doublereal sbig = 85.;
813
814 /* System generated locals */
815 integer i__1, i__2;
816 doublereal d__1, d__2;
817
818 /* Builtin functions */
819 double sqrt(doublereal), exp(doublereal), d_sign(doublereal *, doublereal
820 *);
821
822 /* Local variables */
823 static doublereal a, e, f, s, t, c1, c2, d0, d1, d2, f0;
824 static integer n1, n2;
825 static doublereal p1[3], p2[3], s1, s2, t0, t1, t2, al, fp, tm, un[3];
826 static integer nm1, lp1, lp2;
827 static doublereal tp1, d1d2, scm, dsm, ems, sig;
828 static integer lpl;
829 static doublereal sgn;
830 static integer nit;
831 static doublereal ssm, ems2, d1pd2, fneg, dsig, dmax__, fmax;
832 static integer icnt;
833 static doublereal ftol, rtol, stol, coshm, sigin, sinhm, ssinh;
834 extern doublereal store_(doublereal *);
835 static doublereal unorm;
836 extern doublereal arclen_(doublereal *, doublereal *);
837 static doublereal coshmm;
838 extern /* Subroutine */ int snhcsh_(doublereal *, doublereal *,
839 doublereal *, doublereal *);
840 extern integer lstptr_(integer *, integer *, integer *, integer *);
841
842
843 /* *********************************************************** */
844
845 /* From SSRFPACK */
846 /* Robert J. Renka */
847 /* Dept. of Computer Science */
848 /* Univ. of North Texas */
849 /* renka@cs.unt.edu */
850 /* 11/21/96 */
851
852 /* Given a triangulation of a set of nodes on the unit */
853 /* sphere, along with data values H and gradients GRAD at the */
854 /* nodes, this subroutine determines, for each triangulation */
855 /* arc, the smallest (nonnegative) tension factor SIGMA such */
856 /* that the Hermite interpolatory tension spline H(A), de- */
857 /* fined by SIGMA and the endpoint values and directional */
858 /* derivatives, preserves local shape properties of the data. */
859 /* In order to define the shape properties on an arc, it is */
860 /* convenient to map the arc to an interval (A1,A2). Then, */
861 /* denoting the endpoint data values by H1,H2 and the deriva- */
862 /* tives (tangential gradient components) by HP1,HP2, and */
863 /* letting S = (H2-H1)/(A2-A1), the data properties are */
864
865 /* Monotonicity: S, HP1, and HP2 are nonnegative or */
866 /* nonpositive, */
867 /* and */
868
869 /* Convexity: HP1 .LE. S .LE. HP2 or HP1 .GE. S */
870 /* .GE. HP2. */
871
872 /* The corresponding properties of H are constant sign of the */
873 /* first and second derivatives, respectively. Note that, */
874 /* unless HP1 = S = HP2, infinite tension is required (and H */
875 /* is linear on the interval) if S = 0 in the case of mono- */
876 /* tonicity, or if HP1 = S or HP2 = S in the case of */
877 /* convexity. */
878
879 /* Note that if gradients are to be computed by Subroutine */
880 /* GRADG or function values and gradients are computed by */
881 /* SMSURF, it may be desirable to alternate those computa- */
882 /* tions (which require tension factors) with calls to this */
883 /* subroutine. This iterative procedure should terminate */
884 /* with a call to GETSIG in order to ensure that the shape */
885 /* properties are preserved, and convergence can be achieved */
886 /* (at the cost of optimality) by allowing only increases in */
887 /* tension factors (refer to the parameter descriptions for */
888 /* SIGMA, DSMAX, and IER). */
889
890 /* Refer to functions SIG0, SIG1, and SIG2 for means of */
891 /* selecting minimum tension factors to preserve more general */
892 /* properties. */
893
894 /* On input: */
895
896 /* N = Number of nodes in the triangulation. N .GE. 3. */
897
898 /* X,Y,Z = Arrays of length N containing the Cartesian */
899 /* coordinates of the nodes. */
900
901 /* H = Array of length N containing data values at the */
902 /* nodes. H(I) is associated with (X(I),Y(I),Z(I)) */
903 /* for I = 1,...,N. */
904
905 /* LIST,LPTR,LEND = Data structure defining the tri- */
906 /* angulation. Refer to STRIPACK */
907 /* Subroutine TRMESH. */
908
909 /* GRAD = Array dimensioned 3 by N whose columns con- */
910 /* tain gradients at the nodes. GRAD( ,J) must */
911 /* be orthogonal to node J: GRAD(1,J)*X(J) + */
912 /* GRAD(2,J)*Y(J) + GRAD(3,J)*Z(J) = 0.. Refer */
913 /* to Subroutines GRADG, GRADL, and SMSURF. */
914
915 /* TOL = Tolerance whose magnitude determines how close */
916 /* each tension factor is to its optimal value */
917 /* when nonzero finite tension is necessary and */
918 /* sufficient to satisfy the constraint -- */
919 /* abs(TOL) is an upper bound on the magnitude */
920 /* of the smallest (nonnegative) or largest (non- */
921 /* positive) value of the first or second deriva- */
922 /* tive of H in the interval. Thus, the con- */
923 /* straint is satisfied, but possibly with more */
924 /* tension than necessary. */
925
926 /* The above parameters are not altered by this routine. */
927
928 /* SIGMA = Array of length 2*NA = 6*(N-1)-2*NB, where */
929 /* NA and NB are the numbers of arcs and boun- */
930 /* dary nodes, respectively, containing minimum */
931 /* values of the tension factors. The tension */
932 /* factors are associated with arcs in one-to- */
933 /* one correspondence with LIST entries. Note */
934 /* that each arc N1-N2 has two LIST entries and */
935 /* thus, the tension factor is stored in both */
936 /* SIGMA(I) and SIGMA(J) where LIST(I) = N2 (in */
937 /* the adjacency list for N1) and LIST(J) = N1 */
938 /* (in the list associated with N2). SIGMA */
939 /* should be set to all zeros if minimal ten- */
940 /* sion is desired, and should be unchanged */
941 /* from a previous call in order to ensure con- */
942 /* vergence of the iterative procedure describ- */
943 /* ed in the header comments. */
944
945 /* On output: */
946
947 /* SIGMA = Array containing tension factors for which */
948 /* H(A) preserves the local data properties on */
949 /* each triangulation arc, with the restriction */
950 /* that SIGMA(I) .LE. 85 for all I (unless the */
951 /* input value is larger). The factors are as */
952 /* small as possible (within the tolerance) but */
953 /* not less than their input values. If infin- */
954 /* ite tension is required on an arc, the cor- */
955 /* responding factor is SIGMA(I) = 85 (and H */
956 /* is an approximation to the linear inter- */
957 /* polant on the arc), and if neither property */
958 /* is satisfied by the data, then SIGMA(I) = 0 */
959 /* (assuming its input value is 0), and thus H */
960 /* is cubic on the arc. */
961
962 /* DSMAX = Maximum increase in a component of SIGMA */
963 /* from its input value. */
964
965 /* IER = Error indicator and information flag: */
966 /* IER = I if no errors were encountered and I */
967 /* components of SIGMA were altered from */
968 /* their input values for I .GE. 0. */
969 /* IER = -1 if N < 3. SIGMA is not altered in */
970 /* this case. */
971 /* IER = -2 if duplicate nodes were encountered. */
972
973 /* STRIPACK modules required by GETSIG: LSTPTR, STORE */
974
975 /* SSRFPACK modules required by GETSIG: ARCLEN, SNHCSH */
976
977 /* Intrinsic functions called by GETSIG: ABS, EXP, MAX, MIN, */
978 /* SIGN, SQRT */
979
980 /* *********************************************************** */
981
982
983 /* Parameter adjustments */
984 grad -= 4;
985 --lend;
986 --h__;
987 --z__;
988 --y;
989 --x;
990 --list;
991 --lptr;
992 --sigma;
993
994 /* Function Body */
995 nm1 = *n - 1;
996 if (nm1 < 2) {
997 goto L11;
998 }
999
1000 /* Compute an absolute tolerance FTOL = abs(TOL) and a */
1001 /* relative tolerance RTOL = 100*Macheps. */
1002
1003 ftol = fabs(*tol);
1004 rtol = 1.;
1005 L1:
1006 rtol /= 2.;
1007 d__1 = rtol + 1.;
1008 if (store_(&d__1) > 1.) {
1009 goto L1;
1010 }
1011 rtol *= 200.;
1012
1013 /* Print a heading. */
1014
1015 /* IF (LUN .GE. 0) WRITE (LUN,100) N, FTOL */
1016 /* 100 FORMAT ('1',13X,'GETSIG -- N =',I4,', TOL = ',E10.3//) */
1017 if (dbg_verbose) fprintf (stderr, "GETSIG -- N = %4d TOL = %g\n", *n, ftol);
1018
1019 /* Initialize change counter ICNT and maximum change DSM for */
1020 /* the loop on arcs. */
1021
1022 icnt = 0;
1023 dsm = 0.;
1024
1025 /* Loop on arcs N1-N2 for which N2 > N1. LPL points to the */
1026 /* last neighbor of N1. */
1027
1028 i__1 = nm1;
1029 for (n1 = 1; n1 <= i__1; ++n1) {
1030 lpl = lend[n1];
1031 lp1 = lpl;
1032
1033 /* Top of loop on neighbors N2 of N1. */
1034
1035 L2:
1036 lp1 = lptr[lp1];
1037 n2 = (i__2 = list[lp1], abs(i__2));
1038 if (n2 <= n1) {
1039 goto L9;
1040 }
1041
1042 /* Print a message and compute parameters for the arc: */
1043 /* nodal coordinates P1 and P2, arc-length AL, */
1044 /* UNORM = magnitude of P1 X P2, and */
1045 /* SIGIN = input SIGMA value. */
1046
1047 /* IF (LUN .GE. 0) WRITE (LUN,110) N1, N2 */
1048 /* 110 FORMAT (/1X,'ARC',I4,' -',I4) */
1049 if (dbg_verbose) fprintf (stderr, "ARC %d - %d\n", n1, n2);
1050 p1[0] = x[n1];
1051 p1[1] = y[n1];
1052 p1[2] = z__[n1];
1053 p2[0] = x[n2];
1054 p2[1] = y[n2];
1055 p2[2] = z__[n2];
1056 al = arclen_(p1, p2);
1057 un[0] = p1[1] * p2[2] - p1[2] * p2[1];
1058 un[1] = p1[2] * p2[0] - p1[0] * p2[2];
1059 un[2] = p1[0] * p2[1] - p1[1] * p2[0];
1060 unorm = sqrt(un[0] * un[0] + un[1] * un[1] + un[2] * un[2]);
1061 if (unorm == 0. || al == 0.) {
1062 goto L12;
1063 }
1064 sigin = sigma[lp1];
1065 if (sigin >= sbig) {
1066 goto L9;
1067 }
1068
1069 /* Compute scaled directional derivatives S1,S2 at the end- */
1070 /* points (for the direction N1->N2), first difference S, */
1071 /* and second differences D1,D2. */
1072
1073 s1 = al * (grad[n1 * 3 + 1] * p2[0] + grad[n1 * 3 + 2] * p2[1] + grad[
1074 n1 * 3 + 3] * p2[2]) / unorm;
1075 s2 = -al * (grad[n2 * 3 + 1] * p1[0] + grad[n2 * 3 + 2] * p1[1] +
1076 grad[n2 * 3 + 3] * p1[2]) / unorm;
1077 s = h__[n2] - h__[n1];
1078 d1 = s - s1;
1079 d2 = s2 - s;
1080 d1d2 = d1 * d2;
1081
1082 /* Test for infinite tension required to satisfy either */
1083 /* property. */
1084
1085 sig = sbig;
1086 if ((d1d2 == 0. && s1 != s2) || (s == 0. && s1 * s2 > 0.)) {
1087 goto L8;
1088 }
1089
1090 /* Test for SIGMA = 0 sufficient. The data satisfies convex- */
1091 /* ity iff D1D2 .GE. 0, and D1D2 = 0 implies S1 = S = S2. */
1092
1093 sig = 0.;
1094 if (d1d2 < 0.) {
1095 goto L4;
1096 }
1097 if (d1d2 == 0.) {
1098 goto L8;
1099 }
1100 /* Computing MAX */
1101 d__1 = d1 / d2, d__2 = d2 / d1;
1102 t = max(d__1,d__2);
1103 if (t <= 2.) {
1104 goto L8;
1105 }
1106 tp1 = t + 1.;
1107
1108 /* Convexity: find a zero of F(SIG) = SIG*coshm(SIG)/ */
1109 /* sinhm(SIG) - TP1. */
1110
1111 /* F(0) = 2-T < 0, F(TP1) .GE. 0, the derivative of F */
1112 /* vanishes at SIG = 0, and the second derivative of F is */
1113 /* .2 at SIG = 0. A quadratic approximation is used to */
1114 /* obtain a starting point for the Newton method. */
1115
1116 sig = sqrt(t * 10. - 20.);
1117 nit = 0;
1118
1119 /* Top of loop: */
1120
1121 L3:
1122 if (sig <= .5) {
1123 snhcsh_(&sig, &sinhm, &coshm, &coshmm);
1124 t1 = coshm / sinhm;
1125 fp = t1 + sig * (sig / sinhm - t1 * t1 + 1.);
1126 } else {
1127
1128 /* Scale sinhm and coshm by 2*exp(-SIG) in order to avoid */
1129 /* overflow with large SIG. */
1130
1131 ems = exp(-sig);
1132 ssm = 1. - ems * (ems + sig + sig);
1133 t1 = (1. - ems) * (1. - ems) / ssm;
1134 fp = t1 + sig * (sig * 2. * ems / ssm - t1 * t1 + 1.);
1135 }
1136
1137 f = sig * t1 - tp1;
1138 /* IF (LUN .GE. 0) WRITE (LUN,120) SIG, F, FP */
1139 /* 120 FORMAT (1X,'CONVEXITY -- SIG = ',E15.8, */
1140 /* . ', F(SIG) = ',E15.8/1X,35X,'FP(SIG) = ', */
1141 /* . E15.8) */
1142 if (dbg_verbose) fprintf (stderr, "CONVEXITY -- SIG = %g F(SIG) = %g FP(SIG) = %g\n", sig, f, fp);
1143 ++nit;
1144
1145 /* Test for convergence. */
1146
1147 if (fp <= 0.) {
1148 goto L8;
1149 }
1150 dsig = -f / fp;
1151 if (fabs(dsig) <= rtol * sig || (f >= 0. && f <= ftol) || fabs(f) <= rtol)
1152 {
1153 goto L8;
1154 }
1155
1156 /* Update SIG. */
1157
1158 sig += dsig;
1159 goto L3;
1160
1161 /* Convexity cannot be satisfied. Monotonicity can be satis- */
1162 /* fied iff S1*S .GE. 0 and S2*S .GE. 0 since S .NE. 0. */
1163
1164 L4:
1165 if (s1 * s < 0. || s2 * s < 0.) {
1166 goto L8;
1167 }
1168 t0 = s * 3. - s1 - s2;
1169 d0 = t0 * t0 - s1 * s2;
1170
1171 /* SIGMA = 0 is sufficient for monotonicity iff S*T0 .GE. 0 */
1172 /* or D0 .LE. 0. */
1173
1174 if (d0 <= 0. || s * t0 >= 0.) {
1175 goto L8;
1176 }
1177
1178 /* Monotonicity: find a zero of F(SIG) = sign(S)*HP(R), */
1179 /* where HPP(R) = 0 and HP, HPP denote derivatives of H. */
1180 /* F has a unique zero, F(0) < 0, and F approaches */
1181 /* abs(S) as SIG increases. */
1182
1183 /* Initialize parameters for the secant method. The method */
1184 /* uses three points: (SG0,F0), (SIG,F), and */
1185 /* (SNEG,FNEG), where SG0 and SNEG are defined implicitly */
1186 /* by DSIG = SIG - SG0 and DMAX = SIG - SNEG. */
1187
1188 sgn = d_sign(&c_b23, &s);
1189 sig = sbig;
1190 fmax = sgn * (sig * s - s1 - s2) / (sig - 2.);
1191 if (fmax <= 0.) {
1192 goto L8;
1193 }
1194 stol = rtol * sig;
1195 f = fmax;
1196 f0 = sgn * d0 / ((d1 - d2) * 3.);
1197 fneg = f0;
1198 dsig = sig;
1199 dmax__ = sig;
1200 d1pd2 = d1 + d2;
1201 nit = 0;
1202
1203 /* Top of loop: compute the change in SIG by linear */
1204 /* interpolation. */
1205
1206 L5:
1207 dsig = -f * dsig / (f - f0);
1208 /* IF (LUN .GE. 0) WRITE (LUN,130) DSIG */
1209 /* 130 FORMAT (1X,'MONOTONICITY -- DSIG = ',E15.8) */
1210 if (dbg_verbose) fprintf (stderr, "MONOTONICITY -- DSIG = %g\n", dsig);
1211 if (fabs(dsig) > fabs(dmax__) || dsig * dmax__ > 0.) {
1212 goto L7;
1213 }
1214
1215 /* Restrict the step-size such that abs(DSIG) .GE. STOL/2. */
1216 /* Note that DSIG and DMAX have opposite signs. */
1217
1218 if (fabs(dsig) < stol / 2.) {
1219 d__1 = stol / 2.;
1220 dsig = -d_sign(&d__1, &dmax__);
1221 }
1222
1223 /* Update SIG, F0, and F. */
1224
1225 sig += dsig;
1226 f0 = f;
1227 if (sig <= .5) {
1228
1229 /* Use approximations to the hyperbolic functions designed */
1230 /* to avoid cancellation error with small SIG. */
1231
1232 snhcsh_(&sig, &sinhm, &coshm, &coshmm);
1233 c1 = sig * coshm * d2 - sinhm * d1pd2;
1234 c2 = sig * (sinhm + sig) * d2 - coshm * d1pd2;
1235 a = c2 - c1;
1236 e = sig * sinhm - coshmm - coshmm;
1237 } else {
1238
1239 /* Scale sinhm and coshm by 2*exp(-SIG) in order to avoid */
1240 /* overflow with large SIG. */
1241
1242 ems = exp(-sig);
1243 ems2 = ems + ems;
1244 tm = 1. - ems;
1245 ssinh = tm * (ems + 1.);
1246 ssm = ssinh - sig * ems2;
1247 scm = tm * tm;
1248 c1 = sig * scm * d2 - ssm * d1pd2;
1249 c2 = sig * ssinh * d2 - scm * d1pd2;
1250
1251 /* R is in (0,1) and well-defined iff HPP(T1)*HPP(T2) < 0. */
1252
1253 f = fmax;
1254 if (c1 * (sig * scm * d1 - ssm * d1pd2) >= 0.) {
1255 goto L6;
1256 }
1257 a = ems2 * (sig * tm * d2 + (tm - sig) * d1pd2);
1258 if (a * (c2 + c1) < 0.) {
1259 goto L6;
1260 }
1261 e = sig * ssinh - scm - scm;
1262 }
1263
1264 f = (sgn * (e * s2 - c2) + sqrt(a * (c2 + c1))) / e;
1265
1266 /* Update the number of iterations NIT. */
1267
1268 L6:
1269 ++nit;
1270 /* IF (LUN .GE. 0) WRITE (LUN,140) NIT, SIG, F */
1271 /* 140 FORMAT (1X,11X,I2,' -- SIG = ',E15.8,', F = ', */
1272 /* . E15.8) */
1273 if (dbg_verbose) fprintf (stderr, "%d -- SIG = %g F = %g\n", nit, sig, f);
1274
1275 /* Test for convergence. */
1276
1277 stol = rtol * sig;
1278 if (fabs(dmax__) <= stol || (f >= 0. && f <= ftol) || fabs(f) <= rtol) {
1279 goto L8;
1280 }
1281 dmax__ += dsig;
1282 if (f0 * f > 0. && fabs(f) >= fabs(f0)) {
1283 goto L7;
1284 }
1285 if (f0 * f <= 0.) {
1286
1287 /* F and F0 have opposite signs. Update (SNEG,FNEG) to */
1288 /* (SG0,F0) so that F and FNEG always have opposite */
1289 /* signs. If SIG is closer to SNEG than SG0 and abs(F) */
1290 /* < abs(FNEG), then swap (SNEG,FNEG) with (SG0,F0). */
1291
1292 t1 = dmax__;
1293 t2 = fneg;
1294 dmax__ = dsig;
1295 fneg = f0;
1296 if (fabs(dsig) > fabs(t1) && fabs(f) < fabs(t2)) {
1297
1298 dsig = t1;
1299 f0 = t2;
1300 }
1301 }
1302 goto L5;
1303
1304 /* Bottom of loop: F0*F > 0 and the new estimate would */
1305 /* be outside of the bracketing interval of length */
1306 /* abs(DMAX). Reset (SG0,F0) to (SNEG,FNEG). */
1307
1308 L7:
1309 dsig = dmax__;
1310 f0 = fneg;
1311 goto L5;
1312
1313 /* Update SIGMA, ICNT, and DSM if necessary. */
1314
1315 L8:
1316 sig = min(sig,sbig);
1317 if (sig > sigin) {
1318 sigma[lp1] = sig;
1319 lp2 = lstptr_(&lend[n2], &n1, &list[1], &lptr[1]);
1320 sigma[lp2] = sig;
1321 ++icnt;
1322 /* Computing MAX */
1323 d__1 = dsm, d__2 = sig - sigin;
1324 dsm = max(d__1,d__2);
1325 }
1326
1327 /* Bottom of loop on neighbors N2 of N1. */
1328
1329 L9:
1330 if (lp1 != lpl) {
1331 goto L2;
1332 }
1333 /* L10: */
1334 }
1335
1336 /* No errors encountered. */
1337
1338 *dsmax = dsm;
1339 *ier = icnt;
1340 return 0;
1341
1342 /* N < 3. */
1343
1344 L11:
1345 *dsmax = 0.;
1346 *ier = -1;
1347 return 0;
1348
1349 /* Nodes N1 and N2 coincide. */
1350
1351 L12:
1352 *dsmax = dsm;
1353 *ier = -2;
1354 return 0;
1355 } /* getsig_ */
1356
givens_(doublereal * a,doublereal * b,doublereal * c__,doublereal * s)1357 /* Subroutine */ int givens_(doublereal *a, doublereal *b, doublereal *c__,
1358 doublereal *s)
1359 {
1360 /* Builtin functions */
1361 double sqrt(doublereal);
1362
1363 /* Local variables */
1364 static doublereal r__, u, v, aa, bb;
1365
1366
1367 /* *********************************************************** */
1368
1369 /* From SSRFPACK */
1370 /* Robert J. Renka */
1371 /* Dept. of Computer Science */
1372 /* Univ. of North Texas */
1373 /* renka@cs.unt.edu */
1374 /* 05/09/92 */
1375
1376 /* This subroutine constructs the Givens plane rotation, */
1377
1378 /* ( C S) */
1379 /* G = ( ) , where C*C + S*S = 1, */
1380 /* (-S C) */
1381
1382 /* which zeros the second component of the vector (A,B)**T */
1383 /* (transposed). Subroutine ROTATE may be called to apply */
1384 /* the transformation to a 2 by N matrix. */
1385
1386 /* This routine is identical to Subroutine SROTG from the */
1387 /* LINPACK BLAS (Basic Linear Algebra Subroutines). */
1388
1389 /* On input: */
1390
1391 /* A,B = Components of the vector defining the rota- */
1392 /* tion. These are overwritten by values R */
1393 /* and Z (described below) which define C and S. */
1394
1395 /* On output: */
1396
1397 /* A = Signed Euclidean norm R of the input vector: */
1398 /* R = +/-SQRT(A*A + B*B) */
1399
1400 /* B = Value Z such that: */
1401 /* C = SQRT(1-Z*Z) and S=Z if ABS(Z) .LE. 1, and */
1402 /* C = 1/Z and S = SQRT(1-C*C) if ABS(Z) > 1. */
1403
1404 /* C = +/-(A/R) or 1 if R = 0. */
1405
1406 /* S = +/-(B/R) or 0 if R = 0. */
1407
1408 /* Modules required by GIVENS: None */
1409
1410 /* Intrinsic functions called by GIVENS: ABS, SQRT */
1411
1412 /* *********************************************************** */
1413
1414
1415 /* Local parameters: */
1416
1417 /* AA,BB = Local copies of A and B */
1418 /* R = C*A + S*B = +/-SQRT(A*A+B*B) */
1419 /* U,V = Variables used to scale A and B for computing R */
1420
1421 aa = *a;
1422 bb = *b;
1423 if (fabs(aa) > fabs(bb)) {
1424
1425 /* ABS(A) > ABS(B). */
1426
1427 u = aa + aa;
1428 v = bb / u;
1429 r__ = sqrt(v * v + .25) * u;
1430 *c__ = aa / r__;
1431 *s = v * (*c__ + *c__);
1432
1433 /* Note that R has the sign of A, C > 0, and S has */
1434 /* SIGN(A)*SIGN(B). */
1435
1436 *b = *s;
1437 *a = r__;
1438 } else if (bb != 0.) {
1439
1440 /* ABS(A) .LE. ABS(B). */
1441
1442 u = bb + bb;
1443 v = aa / u;
1444
1445 /* Store R in A. */
1446
1447 *a = sqrt(v * v + .25) * u;
1448 *s = bb / *a;
1449 *c__ = v * (*s + *s);
1450
1451 /* Note that R has the sign of B, S > 0, and C has */
1452 /* SIGN(A)*SIGN(B). */
1453
1454 *b = 1.;
1455 if (*c__ != 0.) {
1456 *b = 1. / *c__;
1457 }
1458 } else {
1459
1460 /* A = B = 0. */
1461
1462 *c__ = 1.;
1463 *s = 0.;
1464 }
1465 return 0;
1466 } /* givens_ */
1467
gradg_(integer * n,doublereal * x,doublereal * y,doublereal * z__,doublereal * f,integer * list,integer * lptr,integer * lend,integer * iflgs,doublereal * sigma,integer * nit,doublereal * dgmax,doublereal * grad,integer * ier)1468 /* Subroutine */ int gradg_(integer *n, doublereal *x, doublereal *y,
1469 doublereal *z__, doublereal *f, integer *list, integer *lptr, integer
1470 *lend, integer *iflgs, doublereal *sigma, integer *nit, doublereal *
1471 dgmax, doublereal *grad, integer *ier)
1472 {
1473 /* System generated locals */
1474 integer i__1, i__2;
1475 doublereal d__1, d__2;
1476
1477 /* Builtin functions */
1478 double sqrt(doublereal), atan(doublereal);
1479
1480 /* Local variables */
1481 static doublereal d__;
1482 static integer j, k;
1483 static doublereal t, g1, g2, g3, r1, r2, a11, a12, a22, fk, sd, cx;
1484 static integer nn;
1485 static doublereal cy, xj, xk, yk, zk, yj, zj, sx, sy, xs, ys, dg1, dg2,
1486 dgk[3], den;
1487 static integer ifl;
1488 static doublereal det, sig;
1489 static integer lpj, lpl;
1490 static doublereal tol, alfa, dgmx;
1491 static integer iter;
1492 static doublereal sinal;
1493 static integer maxit;
1494 extern /* Subroutine */ int grcoef_(doublereal *, doublereal *,
1495 doublereal *), constr_(doublereal *, doublereal *, doublereal *,
1496 doublereal *, doublereal *, doublereal *, doublereal *), aplyrt_(
1497 doublereal *, doublereal *, doublereal *, doublereal *,
1498 doublereal *, doublereal *, doublereal *);
1499
1500
1501 /* *********************************************************** */
1502
1503 /* From SSRFPACK */
1504 /* Robert J. Renka */
1505 /* Dept. of Computer Science */
1506 /* Univ. of North Texas */
1507 /* renka@cs.unt.edu */
1508 /* 07/24/96 */
1509
1510 /* Given a triangulation of N nodes on the unit sphere with */
1511 /* data values F at the nodes and tension factors SIGMA asso- */
1512 /* ciated with the arcs, this routine uses a global method */
1513 /* to compute estimated gradients at the nodes. The method */
1514 /* consists of minimizing a quadratic functional Q(G) over */
1515 /* the N-vector G of gradients, where Q approximates the */
1516 /* linearized curvature of the restriction to arcs of the */
1517 /* interpolatory function F defined by Function FVAL. The */
1518 /* restriction of F to an arc of the triangulation is the */
1519 /* Hermite interpolatory tension spline defined by the data */
1520 /* values and tangential gradient components at the endpoints */
1521 /* of the arc. Letting D1F(A) and D2F(A) denote first and */
1522 /* second derivatives of F with respect to a parameter A var- */
1523 /* ying along a triangulation arc, Q is the sum of integrals */
1524 /* over the arcs of D2F(A)**2 + ((SIGMA/L)*(D1F(A)-S))**2, */
1525 /* where L denotes arc-length, SIGMA is the appropriate ten- */
1526 /* sion factor, and S is the slope of the linear function of */
1527 /* A which interpolates the values of F at the endpoints of */
1528 /* the arc. */
1529
1530 /* Since the gradient at node K lies in the plane tangent */
1531 /* to the sphere surface at K, it is effectively defined by */
1532 /* only two components -- its X and Y components in the coor- */
1533 /* dinate system obtained by rotating K to the north pole. */
1534 /* Thus, the minimization problem corresponds to an order-2N */
1535 /* symmetric positive-definite sparse linear system which is */
1536 /* solved by a block Gauss-Seidel method with 2 by 2 blocks. */
1537
1538 /* An alternative method, Subroutine GRADL, computes a */
1539 /* local approximation to the gradient at a single node and, */
1540 /* although less efficient when all gradients are needed, was */
1541 /* found to be generally more accurate (in the case of uni- */
1542 /* form zero tension) when the nodal distribution is very */
1543 /* dense, varies greatly, or does not cover the sphere. */
1544 /* GRADG, on the other hand, was found to be slightly more */
1545 /* accurate on a uniform distribution of 514 nodes. */
1546
1547 /* On input: */
1548
1549 /* N = Number of nodes in the triangulation. N .GE. 3. */
1550
1551 /* X,Y,Z = Arrays of length N containing Cartesian */
1552 /* coordinates of the nodes. X(I)**2 + Y(I)**2 */
1553 /* + Z(I)**2 = 1 for I = 1,...,N. */
1554
1555 /* F = Array of length N containing data values at the */
1556 /* nodes. F(I) is associated with (X(I),Y(I),Z(I)) */
1557 /* for I = 1,...,N. */
1558
1559 /* LIST,LPTR,LEND = Data structure defining the trian- */
1560 /* gulation. Refer to STRIPACK */
1561 /* Subroutine TRMESH. */
1562
1563 /* IFLGS = Tension factor option: */
1564 /* IFLGS .LE. 0 if a single uniform tension */
1565 /* factor is to be used. */
1566 /* IFLGS .GE. 1 if variable tension is desired. */
1567
1568 /* SIGMA = Uniform tension factor (IFLGS .LE. 0), or */
1569 /* array containing tension factors associated */
1570 /* with arcs in one-to-one correspondence with */
1571 /* LIST entries (IFLGS .GE. 1). Refer to Sub- */
1572 /* programs GETSIG, SIG0, SIG1, and SIG2. */
1573
1574 /* The above parameters are not altered by this routine. */
1575
1576 /* NIT = Maximum number of Gauss-Seidel iterations to */
1577 /* be applied. This maximum will likely be a- */
1578 /* chieved if DGMAX is smaller than the machine */
1579 /* precision. NIT .GE. 0. */
1580
1581 /* DGMAX = Nonnegative convergence criterion. The */
1582 /* method is terminated when the maximum change */
1583 /* in a gradient between iterations is at most */
1584 /* DGMAX. The change in a gradient is taken to */
1585 /* be the Euclidean norm of the difference (in */
1586 /* the rotated coordinate system) relative to 1 */
1587 /* plus the norm of the old gradient value. */
1588
1589 /* GRAD = 3 by N array whose columns contain initial */
1590 /* solution estimates (zero vectors are suffici- */
1591 /* ent). GRAD(I,J) contains component I of the */
1592 /* gradient at node J for I = 1,2,3 (X,Y,Z) and */
1593 /* J = 1,...,N. GRAD( ,J) must be orthogonal to */
1594 /* node J -- GRAD(1,J)*X(J) + GRAD(2,J)*Y(J) + */
1595 /* GRAD(3,J)*Z(J) = 0. */
1596
1597 /* On output: */
1598
1599 /* NIT = Number of Gauss-Seidel iterations employed. */
1600
1601 /* DGMAX = Maximum change in a gradient at the last */
1602 /* iteration. */
1603
1604 /* GRAD = Estimated gradients. See the description */
1605 /* under input parameters. GRAD is not changed */
1606 /* if IER = -1. */
1607
1608 /* IER = Error indicator: */
1609 /* IER = 0 if no errors were encountered and the */
1610 /* convergence criterion was achieved. */
1611 /* IER = 1 if no errors were encountered but con- */
1612 /* vergence was not achieved within NIT */
1613 /* iterations. */
1614 /* IER = -1 if N or DGMAX is outside its valid */
1615 /* range or NIT .LT. 0 on input. */
1616 /* IER = -2 if all nodes are collinear or the */
1617 /* triangulation is invalid. */
1618 /* IER = -3 if duplicate nodes were encountered. */
1619
1620 /* SSRFPACK modules required by GRADG: APLYRT, CONSTR, */
1621 /* GRCOEF, SNHCSH */
1622
1623 /* Intrinsic functions called by GRADG: ATAN, MAX, SQRT */
1624
1625 /* *********************************************************** */
1626
1627
1628 /* Local parameters: */
1629
1630 /* ALFA = Arc-length between nodes K and J */
1631 /* A11,A12,A22 = Matrix components of the 2 by 2 block A*DG */
1632 /* = R where A is symmetric, (DG1,DG2,0) is */
1633 /* the change in the gradient at K, and R is */
1634 /* the residual */
1635 /* CX,CY = Components of a rotation mapping K to the */
1636 /* north pole (0,0,1) */
1637 /* D = Function of SIG computed by GRCOEF -- factor */
1638 /* in the order-2 system */
1639 /* DEN = ALFA*SINAL**2 -- factor in the 2 by 2 system */
1640 /* DET = Determinant of the order-2 matrix */
1641 /* DGK = Change in GRAD( ,K) from the previous esti- */
1642 /* mate in the original coordinate system */
1643 /* DGMX = Maximum change in a gradient between itera- */
1644 /* tions */
1645 /* DG1,DG2 = Solution of the 2 by 2 system -- first 2 */
1646 /* components of DGK in the rotated coordi- */
1647 /* nate system */
1648 /* FK = Data value F(K) */
1649 /* G1,G2,G3 = Components of GRAD( ,K) */
1650 /* IFL = Local copy of IFLGS */
1651 /* ITER = Number of iterations used */
1652 /* J = Neighbor of K */
1653 /* K = DO-loop and node index */
1654 /* LPJ = LIST pointer of node J as a neighbor of K */
1655 /* LPL = Pointer to the last neighbor of K */
1656 /* MAXIT = Input value of NIT */
1657 /* NN = Local copy of N */
1658 /* R1,R2 = Components of the residual -- derivatives of */
1659 /* Q with respect to the components of the */
1660 /* gradient at node K */
1661 /* SD = Function of SIG computed by GRCOEF -- factor */
1662 /* in the order-2 system */
1663 /* SIG = Tension factor associated with ARC K-J */
1664 /* SINAL = SIN(ALFA) -- magnitude of the vector cross */
1665 /* product between nodes K and J */
1666 /* SX,SY = Components of a rotation mapping K to the */
1667 /* north pole (0,0,1) */
1668 /* T = Temporary storage for factors in the system */
1669 /* components */
1670 /* TOL = Local copy of DGMAX */
1671 /* XK,YK,ZK = Coordinates of node K -- X(K), Y(K), Z(K) */
1672 /* XJ,YJ,ZJ = Coordinates of node J in the rotated coor- */
1673 /* dinate system */
1674 /* XS,YS = XJ**2, YJ**2 */
1675
1676 /* Parameter adjustments */
1677 grad -= 4;
1678 --lend;
1679 --f;
1680 --z__;
1681 --y;
1682 --x;
1683 --list;
1684 --lptr;
1685 --sigma;
1686
1687 /* Function Body */
1688 nn = *n;
1689 ifl = *iflgs;
1690 maxit = *nit;
1691 tol = *dgmax;
1692
1693 /* Test for errors in input, and initialize iteration count, */
1694 /* tension factor, and output value of DGMAX. */
1695
1696 if (nn < 3 || maxit < 0 || tol < 0.) {
1697 goto L11;
1698 }
1699 iter = 0;
1700 sig = sigma[1];
1701 dgmx = 0.;
1702
1703 /* Top of iteration loop. */
1704
1705 L1:
1706 if (iter == maxit) {
1707 goto L4;
1708 }
1709 dgmx = 0.;
1710
1711 /* Loop on nodes. */
1712
1713 i__1 = nn;
1714 for (k = 1; k <= i__1; ++k) {
1715 xk = x[k];
1716 yk = y[k];
1717 zk = z__[k];
1718 fk = f[k];
1719 g1 = grad[k * 3 + 1];
1720 g2 = grad[k * 3 + 2];
1721 g3 = grad[k * 3 + 3];
1722
1723 /* Construct the rotation mapping node K to the north pole. */
1724
1725 constr_(&xk, &yk, &zk, &cx, &sx, &cy, &sy);
1726
1727 /* Initialize components of the 2 by 2 system for the */
1728 /* change (DG1,DG2,0) in the K-th solution components */
1729 /* (symmetric matrix in A and residual in R). */
1730
1731 a11 = 0.;
1732 a12 = 0.;
1733 a22 = 0.;
1734 r1 = 0.;
1735 r2 = 0.;
1736
1737 /* Loop on neighbors J of node K. */
1738
1739 lpl = lend[k];
1740 lpj = lpl;
1741 L2:
1742 lpj = lptr[lpj];
1743 j = (i__2 = list[lpj], abs(i__2));
1744
1745 /* Compute the coordinates of J in the rotated system. */
1746
1747 t = sx * y[j] + cx * z__[j];
1748 yj = cx * y[j] - sx * z__[j];
1749 zj = sy * x[j] + cy * t;
1750 xj = cy * x[j] - sy * t;
1751
1752 /* Compute arc-length ALFA between J and K, SINAL = */
1753 /* SIN(ALFA), and DEN = ALFA*SIN(ALFA)**2. */
1754
1755 alfa = atan(sqrt((1. - zj) / (zj + 1.))) * 2.;
1756 xs = xj * xj;
1757 ys = yj * yj;
1758 sinal = sqrt(xs + ys);
1759 den = alfa * (xs + ys);
1760
1761 /* Test for coincident nodes and compute functions of SIG: */
1762 /* D = SIG*(SIG*COSHM-SINHM)/E and SD = SIG*SINHM/E for */
1763 /* E = SIG*SINH-2*COSHM. */
1764
1765 if (den == 0.) {
1766 goto L13;
1767 }
1768 if (ifl >= 1) {
1769 sig = sigma[lpj];
1770 }
1771 grcoef_(&sig, &d__, &sd);
1772
1773 /* Update the system components for node J. */
1774
1775 t = d__ / den;
1776 a11 += t * xs;
1777 a12 += t * xj * yj;
1778 a22 += t * ys;
1779 t = (d__ + sd) * (fk - f[j]) / (alfa * alfa * sinal) + (d__ * (g1 * x[
1780 j] + g2 * y[j] + g3 * z__[j]) - sd * (grad[j * 3 + 1] * xk +
1781 grad[j * 3 + 2] * yk + grad[j * 3 + 3] * zk)) / den;
1782 r1 -= t * xj;
1783 r2 -= t * yj;
1784
1785 /* Bottom of loop on neighbors. */
1786
1787 if (lpj != lpl) {
1788 goto L2;
1789 }
1790
1791 /* Solve the 2 by 2 system and update DGMAX. */
1792
1793 det = a11 * a22 - a12 * a12;
1794 if (det == 0. || a11 == 0.) {
1795 goto L12;
1796 }
1797 dg2 = (a11 * r2 - a12 * r1) / det;
1798 dg1 = (r1 - a12 * dg2) / a11;
1799 /* Computing MAX */
1800 d__1 = dgmx, d__2 = sqrt(dg1 * dg1 + dg2 * dg2) / (sqrt(g1 * g1 + g2 *
1801 g2 + g3 * g3) + 1.);
1802 dgmx = max(d__1,d__2);
1803
1804 /* Rotate (DG1,DG2,0) back to the original coordinate */
1805 /* system and update GRAD( ,K). */
1806
1807 aplyrt_(&dg1, &dg2, &cx, &sx, &cy, &sy, dgk);
1808 grad[k * 3 + 1] = g1 + dgk[0];
1809 grad[k * 3 + 2] = g2 + dgk[1];
1810 grad[k * 3 + 3] = g3 + dgk[2];
1811 /* L3: */
1812 }
1813
1814 /* Increment ITER and test for convergence. */
1815
1816 ++iter;
1817 if (dgmx > tol) {
1818 goto L1;
1819 }
1820
1821 /* The method converged. */
1822
1823 *nit = iter;
1824 *dgmax = dgmx;
1825 *ier = 0;
1826 return 0;
1827
1828 /* The method failed to converge within NIT iterations. */
1829
1830 L4:
1831 *dgmax = dgmx;
1832 *ier = 1;
1833 return 0;
1834
1835 /* Invalid input parameter. */
1836
1837 L11:
1838 *nit = 0;
1839 *dgmax = 0.;
1840 *ier = -1;
1841 return 0;
1842
1843 /* Node K and its neighbors are collinear. */
1844
1845 L12:
1846 *nit = 0;
1847 *dgmax = dgmx;
1848 *ier = -2;
1849 return 0;
1850
1851 /* Nodes K and J coincide. */
1852
1853 L13:
1854 *nit = 0;
1855 *dgmax = dgmx;
1856 *ier = -3;
1857 return 0;
1858 } /* gradg_ */
1859
gradl_(integer * n,integer * k,doublereal * x,doublereal * y,doublereal * z__,doublereal * w,integer * list,integer * lptr,integer * lend,doublereal * g,integer * ier)1860 /* Subroutine */ int gradl_(integer *n, integer *k, doublereal *x, doublereal
1861 *y, doublereal *z__, doublereal *w, integer *list, integer *lptr,
1862 integer *lend, doublereal *g, integer *ier)
1863 {
1864 /* Initialized data */
1865
1866 static doublereal rtol = 1e-6;
1867 static doublereal dtol = .01;
1868 static doublereal sf = 1.;
1869
1870 /* System generated locals */
1871 integer i__1;
1872 doublereal d__1, d__2;
1873
1874 /* Builtin functions */
1875 double sqrt(doublereal);
1876
1877 /* Local variables */
1878 static doublereal a[36] /* was [6][6] */, c__;
1879 static integer i__, j, l;
1880 static doublereal s, df;
1881 static integer kk;
1882 static doublereal av, rf, cx;
1883 static integer nn;
1884 static doublereal cy;
1885 static integer np;
1886 static doublereal dx, dy, wk, xp, yp, zp, sx, sy, wt;
1887 static integer im1, ip1, jp1, lm1;
1888 static doublereal rin;
1889 static integer lnp;
1890 static doublereal sum, dmin__;
1891 static integer lmin, ierr, lmax;
1892 static doublereal avsq;
1893 static integer npts[30];
1894 extern /* Subroutine */ int sph_getnp_(doublereal *, doublereal *, doublereal
1895 *, integer *, integer *, integer *, integer *, integer *,
1896 doublereal *, integer *), aplyr_(doublereal *, doublereal *,
1897 doublereal *, doublereal *, doublereal *, doublereal *,
1898 doublereal *, doublereal *, doublereal *, doublereal *), setup_(
1899 doublereal *, doublereal *, doublereal *, doublereal *,
1900 doublereal *, doublereal *, doublereal *, doublereal *), givens_(
1901 doublereal *, doublereal *, doublereal *, doublereal *), rotate_(
1902 integer *, doublereal *, doublereal *, doublereal *, doublereal *)
1903 , constr_(doublereal *, doublereal *, doublereal *, doublereal *,
1904 doublereal *, doublereal *, doublereal *), aplyrt_(doublereal *,
1905 doublereal *, doublereal *, doublereal *, doublereal *,
1906 doublereal *, doublereal *);
1907
1908
1909 /* *********************************************************** */
1910
1911 /* From SSRFPACK */
1912 /* Robert J. Renka */
1913 /* Dept. of Computer Science */
1914 /* Univ. of North Texas */
1915 /* renka@cs.unt.edu */
1916 /* 07/24/96 */
1917
1918 /* Given a triangulation of a set of nodes on the unit */
1919 /* sphere with their associated data values W, this routine */
1920 /* estimates a gradient vector at node K as follows: the */
1921 /* coordinate system is rotated so that K becomes the north */
1922 /* pole, node K and a set of nearby nodes are projected */
1923 /* orthogonally onto the X-Y plane (in the new coordinate */
1924 /* system), a quadratic is fitted in a weighted least squares */
1925 /* sense to the data values at the projected nodes such that */
1926 /* the value (associated with K) at (0,0) is interpolated, X */
1927 /* and Y-partial derivative estimates DX and DY are computed */
1928 /* by differentiating the quadratic at (0,0), and the esti- */
1929 /* mated gradient G is obtained by rotating (DX,DY,0) back to */
1930 /* the original coordinate system. Note that G lies in the */
1931 /* plane tangent to the sphere at node K, i.e., G is orthogo- */
1932 /* nal to the unit vector represented by node K. A Marquardt */
1933 /* stabilization factor is used if necessary to ensure a */
1934 /* well-conditioned least squares system, and a unique solu- */
1935 /* tion exists unless the nodes are collinear. */
1936
1937 /* On input: */
1938
1939 /* N = Number of nodes in the triangulation. N .GE. 7. */
1940
1941 /* K = Node at which the gradient is sought. 1 .LE. K */
1942 /* .LE. N. */
1943
1944 /* X,Y,Z = Arrays containing the Cartesian coordinates */
1945 /* of the nodes. */
1946
1947 /* W = Array containing the data values at the nodes. */
1948 /* W(I) is associated with (X(I),Y(I),Z(I)) for */
1949 /* I = 1,...,N. */
1950
1951 /* LIST,LPTR,LEND = Data structure defining the trian- */
1952 /* gulation. Refer to STRIPACK */
1953 /* Subroutine TRMESH. */
1954
1955 /* Input parameters are not altered by this routine. */
1956
1957 /* On output: */
1958
1959 /* G = X, Y, and Z components (in that order) of the */
1960 /* estimated gradient at node K unless IER < 0. */
1961
1962 /* IER = Error indicator: */
1963 /* IER .GE. 6 if no errors were encountered. */
1964 /* IER contains the number of nodes */
1965 /* (including K) used in the least */
1966 /* squares fit. */
1967 /* IER = -1 if N or K is outside its valid range. */
1968 /* IER = -2 if the least squares system has no */
1969 /* unique solution due to duplicate or */
1970 /* collinear nodes. */
1971
1972 /* STRIPACK module required by GRADL: GETNP */
1973
1974 /* SSRFPACK modules required by GRADL: APLYR, APLYRT, */
1975 /* CONSTR, GIVENS, */
1976 /* ROTATE, SETUP */
1977
1978 /* Intrinsic functions called by GRADL: ABS, MIN, REAL, SQRT */
1979
1980 /* *********************************************************** */
1981
1982
1983 /* Parameter adjustments */
1984 --lend;
1985 --w;
1986 --z__;
1987 --y;
1988 --x;
1989 --list;
1990 --lptr;
1991 --g;
1992
1993 /* Function Body */
1994
1995 /* Local parameters: */
1996
1997 /* A = Transpose of the (upper triangle of the) aug- */
1998 /* mented regression matrix */
1999 /* AV = Root-mean-square distance (in the rotated */
2000 /* coordinate system) between the origin and */
2001 /* the nodes (other than K) in the least */
2002 /* squares fit. The first 3 columns of A**T */
2003 /* are scaled by 1/AVSQ, the next 2 by 1/AV. */
2004 /* AVSQ = AV*AV: accumulated in SUM */
2005 /* C,S = Components of the plane rotation used to */
2006 /* triangularize the regression matrix */
2007 /* CX,SX = Components of a plane rotation about the X- */
2008 /* axis which, together with CY and SY, define */
2009 /* a mapping from node K to the north pole */
2010 /* (0,0,1) */
2011 /* CY,SY = Components of a plane rotation about the Y- */
2012 /* axis */
2013 /* DF = Negative Z component (in the rotated coordi- */
2014 /* nate system) of an element NP of NPTS -- */
2015 /* increasing function of the angular distance */
2016 /* between K and NP. DF lies in the interval */
2017 /* (-1,1). */
2018 /* DMIN = Minimum of the magnitudes of the diagonal */
2019 /* elements of the triangularized regression */
2020 /* matrix */
2021 /* DTOL = Tolerance for detecting an ill-conditioned */
2022 /* system (DMIN is required to be at least */
2023 /* DTOL) */
2024 /* DX,DY = X and Y components of the estimated gradient */
2025 /* in the rotated coordinate system */
2026 /* I,J = Loop indexes */
2027 /* IERR = Error flag for calls to GETNP (not checked) */
2028 /* IM1,IP1 = I-1, I+1 */
2029 /* JP1 = J+1 */
2030 /* KK = Local copy of K */
2031 /* L = Number of columns of A**T to which a rotation */
2032 /* is applied */
2033 /* LM1 = LMIN-1 */
2034 /* LMIN,LMAX = Min(LMN,N), Min(LMX,N) */
2035 /* LMN,LMX = Minimum and maximum values of LNP for N */
2036 /* sufficiently large. In most cases LMN-1 */
2037 /* nodes are used in the fit. 7 .LE. LMN .LE. */
2038 /* LMX. */
2039 /* LNP = Length of NPTS or LMAX+1 */
2040 /* NN = Local copy of N */
2041 /* NP = Element of NPTS to be added to the system */
2042 /* NPTS = Array containing the indexes of a sequence of */
2043 /* nodes ordered by angular distance from K. */
2044 /* NPTS(1)=K and the first LNP-1 elements of */
2045 /* NPTS are used in the least squares fit. */
2046 /* unless LNP = LMAX+1, NPTS(LNP) determines R */
2047 /* (see RIN). */
2048 /* RF = Value of DF associated with NPTS(LNP) unless */
2049 /* LNP = LMAX+1 (see RIN) */
2050 /* RIN = Inverse of a radius of influence R which */
2051 /* enters into WT: R = 1+RF unless all ele- */
2052 /* ments of NPTS are used in the fit (LNP = */
2053 /* LMAX+1), in which case R is the distance */
2054 /* function associated with some point more */
2055 /* distant from K than NPTS(LMAX) */
2056 /* RTOL = Tolerance for determining LNP (and hence R): */
2057 /* if the increase in DF between two successive */
2058 /* elements of NPTS is less than RTOL, they are */
2059 /* treated as being the same distance from node */
2060 /* K and an additional node is added */
2061 /* SF = Marquardt stabilization factor used to damp */
2062 /* out the first 3 solution components (second */
2063 /* partials of the quadratic) when the system */
2064 /* is ill-conditioned. Increasing SF results */
2065 /* in more damping (a more nearly linear fit). */
2066 /* SUM = Sum of squared Euclidean distances (in the */
2067 /* rotated coordinate system) between the */
2068 /* origin and the nodes used in the least */
2069 /* squares fit */
2070 /* WK = W(K) -- data value at node K */
2071 /* WT = Weight for the equation coreesponding to NP: */
2072 /* WT = (R-D)/(R*D) = 1/D - RIN, where D = 1-ZP */
2073 /* is associated with NP */
2074 /* XP,YP,ZP = Coordinates of NP in the rotated coordinate */
2075 /* system unless ZP < 0, in which case */
2076 /* (XP,YP,0) lies on the equator */
2077
2078 nn = *n;
2079 kk = *k;
2080 wk = w[kk];
2081
2082 /* Check for errors and initialize LMIN, LMAX. */
2083
2084 if (nn < 7 || kk < 1 || kk > nn) {
2085 goto L13;
2086 }
2087 lmin = min(10,nn);
2088 lmax = min(30,nn);
2089
2090 /* Compute NPTS, LNP, AVSQ, AV, and R. */
2091 /* Set NPTS to the closest LMIN-1 nodes to K. DF contains */
2092 /* the negative Z component (in the rotated coordinate */
2093 /* system) of the new node on return from GETNP. */
2094
2095 sum = 0.;
2096 npts[0] = kk;
2097 lm1 = lmin - 1;
2098 i__1 = lm1;
2099 for (lnp = 2; lnp <= i__1; ++lnp) {
2100 sph_getnp_(&x[1], &y[1], &z__[1], &list[1], &lptr[1], &lend[1], &lnp,
2101 npts, &df, &ierr);
2102 sum = sum + 1. - df * df;
2103 /* L1: */
2104 }
2105
2106 /* Add additional nodes to NPTS until the increase in */
2107 /* R = 1+RF is at least RTOL. */
2108
2109 i__1 = lmax;
2110 for (lnp = lmin; lnp <= i__1; ++lnp) {
2111 sph_getnp_(&x[1], &y[1], &z__[1], &list[1], &lptr[1], &lend[1], &lnp,
2112 npts, &rf, &ierr);
2113 if (rf - df >= rtol) {
2114 goto L3;
2115 }
2116 sum = sum + 1. - rf * rf;
2117 /* L2: */
2118 }
2119
2120 /* Use all LMAX nodes in the least squares fit. R is */
2121 /* arbitrarily increased by 5 percent. */
2122
2123 rf = rf * 1.05 + .05;
2124 lnp = lmax + 1;
2125
2126 /* There are LNP-2 equations corresponding to nodes */
2127 /* NPTS(2),...,NPTS(LNP-1). */
2128
2129 L3:
2130 avsq = sum / (doublereal) (lnp - 2);
2131 av = sqrt(avsq);
2132 rin = 1. / (rf + 1.);
2133
2134 /* Construct the rotation. */
2135
2136 constr_(&x[kk], &y[kk], &z__[kk], &cx, &sx, &cy, &sy);
2137
2138 /* Set up the first 5 equations of the augmented regression */
2139 /* matrix (transposed) as the columns of A, and zero out */
2140 /* the lower triangle (upper triangle of A) with Givens */
2141 /* rotations. */
2142
2143 for (i__ = 1; i__ <= 5; ++i__) {
2144 np = npts[i__];
2145 aplyr_(&x[np], &y[np], &z__[np], &cx, &sx, &cy, &sy, &xp, &yp, &zp);
2146 wt = 1. / (1. - zp) - rin;
2147 setup_(&xp, &yp, &w[np], &wk, &av, &avsq, &wt, &a[i__ * 6 - 6]);
2148 if (i__ == 1) {
2149 goto L5;
2150 }
2151 im1 = i__ - 1;
2152 i__1 = im1;
2153 for (j = 1; j <= i__1; ++j) {
2154 jp1 = j + 1;
2155 l = 6 - j;
2156 givens_(&a[j + j * 6 - 7], &a[j + i__ * 6 - 7], &c__, &s);
2157 rotate_(&l, &c__, &s, &a[jp1 + j * 6 - 7], &a[jp1 + i__ * 6 - 7]);
2158 /* L4: */
2159 }
2160 L5:
2161 ;
2162 }
2163
2164 /* Add the additional equations to the system using */
2165 /* the last column of A. I .LE. LNP. */
2166
2167 i__ = 7;
2168 L6:
2169 if (i__ == lnp) {
2170 goto L8;
2171 }
2172 np = npts[i__ - 1];
2173 aplyr_(&x[np], &y[np], &z__[np], &cx, &sx, &cy, &sy, &xp, &yp, &zp);
2174 wt = 1. / (1. - zp) - rin;
2175 setup_(&xp, &yp, &w[np], &wk, &av, &avsq, &wt, &a[30]);
2176 for (j = 1; j <= 5; ++j) {
2177 jp1 = j + 1;
2178 l = 6 - j;
2179 givens_(&a[j + j * 6 - 7], &a[j + 29], &c__, &s);
2180 rotate_(&l, &c__, &s, &a[jp1 + j * 6 - 7], &a[jp1 + 29]);
2181 /* L7: */
2182 }
2183 ++i__;
2184 goto L6;
2185
2186 /* Test the system for ill-conditioning. */
2187
2188 L8:
2189 /* Computing MIN */
2190 d__1 = fabs(a[0]), d__2 = fabs(a[7]), d__1 = min(d__1,d__2), d__2 = fabs(a[
2191 14]), d__1 = min(d__1,d__2), d__2 = fabs(a[21]), d__1 = min(d__1,
2192 d__2), d__2 = fabs(a[28]);
2193 dmin__ = min(d__1,d__2);
2194 if (dmin__ >= dtol) {
2195 goto L12;
2196 }
2197 if (lnp <= lmax) {
2198
2199 /* Add another node to the system and increase R. */
2200 /* I = LNP. */
2201
2202 ++lnp;
2203 if (lnp <= lmax) {
2204 sph_getnp_(&x[1], &y[1], &z__[1], &list[1], &lptr[1], &lend[1], &lnp,
2205 npts, &rf, &ierr);
2206 }
2207 rin = 1. / ((rf + 1.) * 1.05);
2208 goto L6;
2209 }
2210
2211 /* Stabilize the system by damping second partials. Add */
2212 /* multiples of the first three unit vectors to the first */
2213 /* three equations. */
2214
2215 for (i__ = 1; i__ <= 3; ++i__) {
2216 a[i__ + 29] = sf;
2217 ip1 = i__ + 1;
2218 for (j = ip1; j <= 6; ++j) {
2219 a[j + 29] = 0.;
2220 /* L9: */
2221 }
2222 for (j = i__; j <= 5; ++j) {
2223 jp1 = j + 1;
2224 l = 6 - j;
2225 givens_(&a[j + j * 6 - 7], &a[j + 29], &c__, &s);
2226 rotate_(&l, &c__, &s, &a[jp1 + j * 6 - 7], &a[jp1 + 29]);
2227 /* L10: */
2228 }
2229 /* L11: */
2230 }
2231
2232 /* Test the linear portion of the stabilized system for */
2233 /* ill-conditioning. */
2234
2235 /* Computing MIN */
2236 d__1 = fabs(a[21]), d__2 = fabs(a[28]);
2237 dmin__ = min(d__1,d__2);
2238 if (dmin__ < dtol) {
2239 goto L14;
2240 }
2241
2242 /* Solve the 2 by 2 triangular system for the estimated */
2243 /* partial derivatives. */
2244
2245 L12:
2246 dy = a[29] / a[28];
2247 dx = (a[23] - a[22] * dy) / a[21] / av;
2248 dy /= av;
2249
2250 /* Rotate the gradient (DX,DY,0) back into the original */
2251 /* coordinate system. */
2252
2253 aplyrt_(&dx, &dy, &cx, &sx, &cy, &sy, &g[1]);
2254 *ier = lnp - 1;
2255 return 0;
2256
2257 /* N or K is outside its valid range. */
2258
2259 L13:
2260 *ier = -1;
2261 return 0;
2262
2263 /* No unique solution due to collinear nodes. */
2264
2265 L14:
2266 *ier = -2;
2267 return 0;
2268 } /* gradl_ */
2269
grcoef_(doublereal * sigma,doublereal * d__,doublereal * sd)2270 /* Subroutine */ int grcoef_(doublereal *sigma, doublereal *d__, doublereal *
2271 sd)
2272 {
2273 /* Builtin functions */
2274 double exp(doublereal);
2275
2276 /* Local variables */
2277 static doublereal e, scm, sig, ems, ssm, coshm, sinhm, ssinh, coshmm;
2278 extern /* Subroutine */ int snhcsh_(doublereal *, doublereal *,
2279 doublereal *, doublereal *);
2280
2281
2282 /* *********************************************************** */
2283
2284 /* From SSRFPACK */
2285 /* Robert J. Renka */
2286 /* Dept. of Computer Science */
2287 /* Univ. of North Texas */
2288 /* renka@cs.unt.edu */
2289 /* 11/21/96 */
2290
2291 /* This subroutine computes factors involved in the linear */
2292 /* systems solved by Subroutines GRADG and SMSGS. */
2293
2294 /* On input: */
2295
2296 /* SIGMA = Nonnegative tension factor associated with a */
2297 /* triangulation arc. */
2298
2299 /* SIGMA is not altered by this routine. */
2300
2301 /* On output: */
2302
2303 /* D = Diagonal factor. D = SIG*(SIG*Coshm(SIG) - */
2304 /* Sinhm(SIG))/E where E = SIG*Sinh(SIG) - 2* */
2305 /* Coshm(SIG). D > 0, and D = 4 at SIG = 0. */
2306
2307 /* SD = Off-diagonal factor. SD = SIG*Sinhm(SIG)/E. */
2308 /* SD > 0, and SD = 2 at SIG = 0. */
2309
2310 /* SSRFPACK module required by GRCOEF: SNHCSH */
2311
2312 /* Intrinsic function called by GRCOEF: EXP */
2313
2314 /* *********************************************************** */
2315
2316 sig = *sigma;
2317 if (sig < 1e-9) {
2318
2319 /* Cubic function: */
2320
2321 *d__ = 4.;
2322 *sd = 2.;
2323 } else if (sig <= .5) {
2324
2325 /* 0 < SIG .LE. .5. */
2326
2327 /* Use approximations designed to avoid cancellation error */
2328 /* in the hyperbolic functions when SIGMA is small. */
2329
2330 snhcsh_(&sig, &sinhm, &coshm, &coshmm);
2331 e = sig * sinhm - coshmm - coshmm;
2332 *d__ = sig * (sig * coshm - sinhm) / e;
2333 *sd = sig * sinhm / e;
2334 } else {
2335
2336 /* SIG > .5. */
2337
2338 /* Scale SINHM, COSHM, and E by 2*EXP(-SIG) in order to */
2339 /* avoid overflow when SIGMA is large. */
2340
2341 ems = exp(-sig);
2342 ssinh = 1. - ems * ems;
2343 ssm = ssinh - sig * 2. * ems;
2344 scm = (1. - ems) * (1. - ems);
2345 e = sig * ssinh - scm - scm;
2346 *d__ = sig * (sig * scm - ssm) / e;
2347 *sd = sig * ssm / e;
2348 }
2349 return 0;
2350 } /* grcoef_ */
2351
hval_(doublereal * b,doublereal * h1,doublereal * h2,doublereal * hp1,doublereal * hp2,doublereal * sigma)2352 doublereal hval_(doublereal *b, doublereal *h1, doublereal *h2, doublereal *
2353 hp1, doublereal *hp2, doublereal *sigma)
2354 {
2355 /* System generated locals */
2356 doublereal ret_val;
2357
2358 /* Builtin functions */
2359 double exp(doublereal);
2360
2361 /* Local variables */
2362 static doublereal e, s, b1, b2, d1, d2, e1, e2, cm, sm, tm, ts, cm2, sb1,
2363 sb2, sm2, tm1, tm2, cmm, sig, ems, dummy;
2364 extern /* Subroutine */ int snhcsh_(doublereal *, doublereal *,
2365 doublereal *, doublereal *);
2366
2367
2368 /* *********************************************************** */
2369
2370 /* From SSRFPACK */
2371 /* Robert J. Renka */
2372 /* Dept. of Computer Science */
2373 /* Univ. of North Texas */
2374 /* renka@cs.unt.edu */
2375 /* 11/21/96 */
2376
2377 /* Given a line segment P1-P2 containing a point P, along */
2378 /* with values and derivatives at the endpoints, this func- */
2379 /* tion returns the value H(P), where H is the Hermite inter- */
2380 /* polatory tension spline defined by the endpoint data. */
2381
2382 /* On input: */
2383
2384 /* B = Local coordinate of P with respect to P1-P2: */
2385 /* P = B*P1 + (1-B)*P2, and thus B = d(P,P2)/ */
2386 /* d(P1,P2), where d(P1,P2) is the distance between */
2387 /* P1 and P2. B < 0 or B > 1 results in extrapola- */
2388 /* tion. */
2389
2390 /* H1,H2 = Values interpolated at P1 and P2, respec- */
2391 /* tively. */
2392
2393 /* HP1,HP2 = Products of d(P1,P2) with first order der- */
2394 /* ivatives at P1 and P2, respectively. HP1 */
2395 /* may, for example, be the scalar product of */
2396 /* P2-P1 with a gradient at P1. */
2397
2398 /* SIGMA = Nonnegative tension factor associated with */
2399 /* the spline. SIGMA = 0 corresponds to a */
2400 /* cubic spline, and H approaches the linear */
2401 /* interpolant of H1 and H2 as SIGMA increases. */
2402
2403 /* Input parameters are not altered by this function. */
2404
2405 /* On output: */
2406
2407 /* HVAL = Interpolated value H(P). */
2408
2409 /* SSRFPACK module required by HVAL: SNHCSH */
2410
2411 /* Intrinsic functions called by HVAL: ABS, EXP */
2412
2413 /* *********************************************************** */
2414
2415 b1 = *b;
2416 b2 = 1. - b1;
2417
2418 /* Compute slope S and second differences D1 and D2 scaled */
2419 /* by the separation between P1 and P2. */
2420
2421 s = *h2 - *h1;
2422 d1 = s - *hp1;
2423 d2 = *hp2 - s;
2424
2425 /* Test the range of SIGMA. */
2426
2427 sig = fabs(*sigma);
2428 if (sig < 1e-9) {
2429
2430 /* Hermite cubic interpolation: */
2431
2432 ret_val = *h1 + b2 * (*hp1 + b2 * (d1 + b1 * (d1 - d2)));
2433 } else if (sig <= .5) {
2434
2435 /* 0 < SIG .LE. .5. Use approximations designed to avoid */
2436 /* cancellation error in the hyperbolic functions. */
2437
2438 sb2 = sig * b2;
2439 snhcsh_(&sig, &sm, &cm, &cmm);
2440 snhcsh_(&sb2, &sm2, &cm2, &dummy);
2441 e = sig * sm - cmm - cmm;
2442 ret_val = *h1 + b2 * *hp1 + ((cm * sm2 - sm * cm2) * (d1 + d2) + sig *
2443 (cm * cm2 - (sm + sig) * sm2) * d1) / (sig * e);
2444 } else {
2445
2446 /* SIG > .5. Use negative exponentials in order to avoid */
2447 /* overflow. Note that EMS = EXP(-SIG). */
2448
2449 sb1 = sig * b1;
2450 sb2 = sig - sb1;
2451 e1 = exp(-sb1);
2452 e2 = exp(-sb2);
2453 ems = e1 * e2;
2454 tm = 1. - ems;
2455 ts = tm * tm;
2456 tm1 = 1. - e1;
2457 tm2 = 1. - e2;
2458 e = tm * (sig * (ems + 1.) - tm - tm);
2459 ret_val = *h1 + b2 * s + (tm * tm1 * tm2 * (d1 + d2) + sig * ((e2 *
2460 tm1 * tm1 - b1 * ts) * d1 + (e1 * tm2 * tm2 - b2 * ts) * d2))
2461 / (sig * e);
2462 }
2463 return ret_val;
2464 } /* hval_ */
2465
intrc0_(integer * n,doublereal * plat,doublereal * plon,doublereal * x,doublereal * y,doublereal * z__,doublereal * w,integer * list,integer * lptr,integer * lend,integer * ist,doublereal * pw,integer * ier)2466 /* Subroutine */ int intrc0_(integer *n, doublereal *plat, doublereal *plon,
2467 doublereal *x, doublereal *y, doublereal *z__, doublereal *w, integer
2468 *list, integer *lptr, integer *lend, integer *ist, doublereal *pw,
2469 integer *ier)
2470 {
2471 /* Builtin functions */
2472 double cos(doublereal), sin(doublereal);
2473
2474 /* Local variables */
2475 static doublereal p[3], b1, b2, b3;
2476 static integer i1, i2, i3, n1, n2;
2477 static doublereal s12;
2478 static integer lp;
2479 static doublereal cos_plat;
2480 static doublereal sum, ptn1, ptn2;
2481 extern /* Subroutine */ int sph_trfind_(integer *, doublereal *, integer *,
2482 doublereal *, doublereal *, doublereal *, integer *, integer *,
2483 integer *, doublereal *, doublereal *, doublereal *, integer *,
2484 integer *, integer *);
2485
2486
2487 /* *********************************************************** */
2488
2489 /* From SSRFPACK */
2490 /* Robert J. Renka */
2491 /* Dept. of Computer Science */
2492 /* Univ. of North Texas */
2493 /* renka@cs.unt.edu */
2494 /* 07/24/96 */
2495
2496 /* Given a triangulation of a set of nodes on the unit */
2497 /* sphere, along with data values at the nodes, this sub- */
2498 /* routine computes the value at a point P of a continuous */
2499 /* function which interpolates the data values. The interp- */
2500 /* olatory function is linear on each underlying triangle */
2501 /* (planar triangle with the same vertices as a spherical */
2502 /* triangle). If P is not contained in a triangle, an ex- */
2503 /* trapolated value is taken to be the interpolated value at */
2504 /* the nearest point of the triangulation boundary. */
2505
2506 /* On input: */
2507
2508 /* N = Number of nodes in the triangulation. N .GE. 3. */
2509
2510 /* PLAT,PLON = Latitude and longitude of P in radians. */
2511
2512 /* X,Y,Z = Arrays containing Cartesian coordinates of */
2513 /* the nodes. */
2514
2515 /* W = Array containing data values at the nodes. W(I) */
2516 /* is associated with (X(I),Y(I),Z(I)) for I = */
2517 /* 1,...,N. */
2518
2519 /* LIST,LPTR,LEND = Data structure defining the trian- */
2520 /* gulation. Refer to STRIPACK */
2521 /* Subroutine TRMESH. */
2522
2523 /* IST = Index of the starting node in the search for a */
2524 /* triangle containing P. 1 .LE. IST .LE. N. */
2525 /* The output value of IST from a previous call */
2526 /* may be a good choice. */
2527
2528 /* Input parameters other than IST are not altered by this */
2529 /* routine. */
2530
2531 /* On output: */
2532
2533 /* IST = Index of one of the vertices of the triangle */
2534 /* containing P (or nearest P) unless IER = -1 */
2535 /* or IER = -2. */
2536
2537 /* PW = Value of the interpolatory function at P if */
2538 /* IER .GE. 0. */
2539
2540 /* IER = Error indicator: */
2541 /* IER = 0 if interpolation was performed */
2542 /* successfully. */
2543 /* IER = 1 if extrapolation was performed */
2544 /* successfully. */
2545 /* IER = -1 if N < 3 or IST is outside its valid */
2546 /* range. */
2547 /* IER = -2 if the nodes are collinear. */
2548 /* IER = -3 if P is not in a triangle and the */
2549 /* angle between P and the nearest boun- */
2550 /* dary point is at least 90 degrees. */
2551
2552 /* STRIPACK modules required by INTRC0: JRAND, LSTPTR, */
2553 /* STORE, TRFIND */
2554
2555 /* Intrinsic functions called by INTRC0: COS, SIN */
2556
2557 /* *********************************************************** */
2558
2559
2560 /* Local parameters: */
2561
2562 /* B1,B2,B3 = Barycentric coordinates of the central projec- */
2563 /* tion of P onto the underlying planar trian- */
2564 /* gle, or (B1 and B2) projection of Q onto the */
2565 /* underlying line segment N1-N2 when P is */
2566 /* exterior. Unnormalized coordinates are */
2567 /* computed by TRFIND when P is in a triangle. */
2568 /* I1,I2,I3 = Vertex indexes returned by TRFIND */
2569 /* LP = LIST pointer to N1 as a neighbor of N2 or N2 */
2570 /* as a neighbor of N1 */
2571 /* N1,N2 = Endpoints of a boundary arc which is visible */
2572 /* from P when P is not contained in a triangle */
2573 /* P = Cartesian coordinates of P */
2574 /* PTN1 = Scalar product (P,N1) */
2575 /* PTN2 = Scalar product (P,N2) */
2576 /* S12 = Scalar product (N1,N2) */
2577 /* SUM = Quantity used to normalize the barycentric */
2578 /* coordinates */
2579
2580 /* Parameter adjustments */
2581 --lend;
2582 --w;
2583 --z__;
2584 --y;
2585 --x;
2586 --list;
2587 --lptr;
2588
2589 /* Function Body */
2590 if (*n < 3 || *ist < 1 || *ist > *n) {
2591 goto L11;
2592 }
2593
2594 /* Transform (PLAT,PLON) to Cartesian coordinates. */
2595
2596 sincos (*plat, &p[2], &cos_plat); sincos (*plon, &p[1], &p[0]); p[0] *= cos_plat; p[1] *= cos_plat;
2597
2598 /* Find the vertex indexes of a triangle containing P. */
2599
2600 sph_trfind_(ist, p, n, &x[1], &y[1], &z__[1], &list[1], &lptr[1], &lend[1], &
2601 b1, &b2, &b3, &i1, &i2, &i3);
2602 if (i1 == 0) {
2603 goto L12;
2604 }
2605 *ist = i1;
2606 if (i3 != 0) {
2607
2608 /* P is contained in the triangle (I1,I2,I3). Normalize the */
2609 /* barycentric coordinates. */
2610
2611 sum = b1 + b2 + b3;
2612 b1 /= sum;
2613 b2 /= sum;
2614 b3 /= sum;
2615 *pw = b1 * w[i1] + b2 * w[i2] + b3 * w[i3];
2616 *ier = 0;
2617 return 0;
2618 }
2619
2620 /* P is exterior to the triangulation, and I1 and I2 are */
2621 /* boundary nodes which are visible from P. Set PW to the */
2622 /* interpolated value at Q, where Q is the closest boundary */
2623 /* point to P. */
2624
2625 /* Traverse the boundary starting from the rightmost visible */
2626 /* node I1. */
2627
2628 n1 = i1;
2629 ptn1 = p[0] * x[n1] + p[1] * y[n1] + p[2] * z__[n1];
2630 if (i1 != i2) {
2631 goto L2;
2632 }
2633
2634 /* All boundary nodes are visible from P. Find a boundary */
2635 /* arc N1->N2 such that P Left (N2 X N1)->N1. */
2636
2637 /* Counterclockwise boundary traversal: */
2638 /* Set N2 to the first neighbor of N1. */
2639
2640 L1:
2641 lp = lend[n1];
2642 lp = lptr[lp];
2643 n2 = list[lp];
2644
2645 /* Compute inner products (N1,N2) and (P,N2), and compute */
2646 /* B2 = DET(P,N1,N2 X N1). */
2647
2648 s12 = x[n1] * x[n2] + y[n1] * y[n2] + z__[n1] * z__[n2];
2649 ptn2 = p[0] * x[n2] + p[1] * y[n2] + p[2] * z__[n2];
2650 b2 = ptn2 - s12 * ptn1;
2651 if (b2 <= 0.) {
2652 goto L2;
2653 }
2654
2655 /* P Right (N2 X N1)->N1 -- Iterate. */
2656
2657 n1 = n2;
2658 i1 = n1;
2659 ptn1 = ptn2;
2660 goto L1;
2661
2662 /* P Left (N2 X N1)->N1, where N2 is the first neighbor of P1. */
2663 /* Clockwise boundary traversal: */
2664
2665 L2:
2666 n2 = n1;
2667 ptn2 = ptn1;
2668
2669 /* Set N1 to the last neighbor of N2 and test for */
2670 /* termination. */
2671
2672 lp = lend[n2];
2673 n1 = -list[lp];
2674 if (n1 == i1) {
2675 goto L13;
2676 }
2677
2678 /* Compute inner products (N1,N2) and (P,N1). */
2679
2680 s12 = x[n1] * x[n2] + y[n1] * y[n2] + z__[n1] * z__[n2];
2681 ptn1 = p[0] * x[n1] + p[1] * y[n1] + p[2] * z__[n1];
2682
2683 /* Compute B2 = DET(P,N1,N2 X N1) = DET(Q,N1,N2 X N1)*(P,Q). */
2684
2685 b2 = ptn2 - s12 * ptn1;
2686 if (b2 <= 0.) {
2687 goto L2;
2688 }
2689
2690 /* Compute B1 = DET(P,N2 X N1,N2) = DET(Q,N2 X N1,N2)*(P,Q). */
2691
2692 b1 = ptn1 - s12 * ptn2;
2693 if (b1 <= 0.) {
2694
2695 /* Q = N2. */
2696
2697 *pw = w[n2];
2698 } else {
2699
2700 /* P Strictly Left (N2 X N1)->N2 and P Strictly Left */
2701 /* N1->(N2 X N1). Thus Q lies on the interior of N1->N2. */
2702 /* Normalize the coordinates and compute PW. */
2703
2704 sum = b1 + b2;
2705 *pw = (b1 * w[n1] + b2 * w[n2]) / sum;
2706 }
2707 *ier = 1;
2708 return 0;
2709
2710 /* N or IST is outside its valid range. */
2711
2712 L11:
2713 *ier = -1;
2714 return 0;
2715
2716 /* Collinear nodes. */
2717
2718 L12:
2719 *ier = -2;
2720 return 0;
2721
2722 /* The angular distance between P and the closest boundary */
2723 /* point to P is at least 90 degrees. */
2724
2725 L13:
2726 *ier = -3;
2727 return 0;
2728 } /* intrc0_ */
2729
intrc1_(integer * n,doublereal * plat,doublereal * plon,doublereal * x,doublereal * y,doublereal * z__,doublereal * f,integer * list,integer * lptr,integer * lend,integer * iflgs,doublereal * sigma,integer * iflgg,doublereal * grad,integer * ist,doublereal * fp,integer * ier)2730 /* Subroutine */ int intrc1_(integer *n, doublereal *plat, doublereal *plon,
2731 doublereal *x, doublereal *y, doublereal *z__, doublereal *f, integer
2732 *list, integer *lptr, integer *lend, integer *iflgs, doublereal *
2733 sigma, integer *iflgg, doublereal *grad, integer *ist, doublereal *fp,
2734 integer *ier)
2735 {
2736 /* Builtin functions */
2737 double cos(doublereal), sin(doublereal), sqrt(doublereal);
2738
2739 /* Local variables */
2740 static doublereal a;
2741 static integer i__;
2742 static doublereal p[3], q[3], b1, b2, b3, g1[3], g2[3];
2743 static integer i1, i2, i3;
2744 static doublereal g3[3];
2745 static integer n1, n2;
2746 static doublereal p1[3], p2[3], p3[3], s1, s2, s3, fq, gq[3], s12;
2747 static integer lp, nn;
2748 static doublereal dum[3], gqn, sum, ptn1, ptn2;
2749 extern doublereal fval_(doublereal *, doublereal *, doublereal *,
2750 doublereal *, doublereal *, doublereal *, doublereal *,
2751 doublereal *, doublereal *, doublereal *, doublereal *,
2752 doublereal *, doublereal *, doublereal *, doublereal *);
2753 static integer ierr;
2754 static doublereal ptgq;
2755 extern /* Subroutine */ int gradl_(integer *, integer *, doublereal *,
2756 doublereal *, doublereal *, doublereal *, integer *, integer *,
2757 integer *, doublereal *, integer *);
2758 static doublereal qnorm;
2759 extern doublereal arclen_(doublereal *, doublereal *);
2760 extern /* Subroutine */ int arcint_(doublereal *, doublereal *,
2761 doublereal *, doublereal *, doublereal *, doublereal *,
2762 doublereal *, doublereal *, doublereal *, doublereal *,
2763 doublereal *), sph_trfind_(integer *, doublereal *, integer *,
2764 doublereal *, doublereal *, doublereal *, integer *, integer *,
2765 integer *, doublereal *, doublereal *, doublereal *, integer *,
2766 integer *, integer *);
2767 extern integer lstptr_(integer *, integer *, integer *, integer *);
2768
2769
2770 /* *********************************************************** */
2771
2772 /* From SSRFPACK */
2773 /* Robert J. Renka */
2774 /* Dept. of Computer Science */
2775 /* Univ. of North Texas */
2776 /* renka@cs.unt.edu */
2777 /* 07/25/96 */
2778
2779 /* Given a triangulation of a set of nodes on the unit */
2780 /* sphere, along with data values and gradients at the nodes, */
2781 /* this routine computes a value F(P), where F interpolates */
2782 /* the nodal data and is once-continuously differentiable */
2783 /* over the convex hull of the nodes. Refer to Function FVAL */
2784 /* for further details. If P is not contained in a triangle, */
2785 /* an extrapolated value is computed by extending F beyond */
2786 /* the boundary in a continuous fashion. */
2787
2788 /* On input: */
2789
2790 /* N = Number of nodes in the triangulation. N .GE. 3 */
2791 /* and N .GE. 7 if IFLGG .LE. 0. */
2792
2793 /* PLAT,PLON = Latitude and longitude in radians of the */
2794 /* point P at which F is to be evaluated. */
2795
2796 /* X,Y,Z = Arrays of length N containing Cartesian */
2797 /* coordinates of the nodes. */
2798
2799 /* F = Array of length N containing values of F at the */
2800 /* nodes: F(I) = F(X(I),Y(I),Z(I)). */
2801
2802 /* LIST,LPTR,LEND = Data structure defining the trian- */
2803 /* gulation. Refer to STRIPACK */
2804 /* Subroutine TRMESH. */
2805
2806 /* IFLGS = Tension factor option: */
2807 /* IFLGS .LE. 0 if a single uniform tension */
2808 /* factor is to be used. */
2809 /* IFLGS .GE. 1 if variable tension is desired. */
2810
2811 /* SIGMA = Uniform tension factor (IFLGS .LE. 0), or */
2812 /* array containing tension factors associated */
2813 /* with arcs in one-to-one correspondence with */
2814 /* LIST entries (IFLGS .GE. 1). Refer to Sub- */
2815 /* programs FVAL, GETSIG, SIG0, SIG1, and SIG2. */
2816
2817 /* IFLGG = Gradient option: */
2818 /* IFLGG .LE. 0 if INTRC1 is to provide grad- */
2819 /* ient estimates as needed (from */
2820 /* GRADL). */
2821 /* IFLGG .GE. 1 if gradients are user-provided */
2822 /* in GRAD. This is more effici- */
2823 /* ent if INTRC1 is to be called */
2824 /* several times. */
2825
2826 /* GRAD = 3 by N array whose I-th column contains */
2827 /* an estimated gradient at node I if IFLGG .GE. */
2828 /* 1, or unused dummy parameter if IFLGG .LE. 0. */
2829 /* Refer to Subroutines GRADL and GRADG. */
2830
2831 /* IST = Index of the starting node in the search for a */
2832 /* triangle containing P. The output value of */
2833 /* IST from a previous call may be a good choice. */
2834 /* 1 .LE. IST .LE. N. */
2835
2836 /* Input parameters other than IST are not altered by this */
2837 /* routine. */
2838
2839 /* On output: */
2840
2841 /* IST = Index of one of the vertices of the triangle */
2842 /* containing P (or a boundary node if P is not */
2843 /* contained in a triangle) unless IER = -1 or */
2844 /* IER = -2. */
2845
2846 /* FP = Value of F at P unless IER < 0, in which case */
2847 /* FP is not defined. */
2848
2849 /* IER = Error indicator and information flag: */
2850 /* IER = 0 if no errors were encountered and P is */
2851 /* contained in a triangle. */
2852 /* IER = 1 if no errors were encountered and */
2853 /* extrapolation was required. */
2854 /* IER = -1 if N or IST is outside its valid */
2855 /* range. */
2856 /* IER = -2 if the nodes are collinear. */
2857 /* IER = -3 if the angular distance between P and */
2858 /* the nearest point of the triangula- */
2859 /* tion is at least 90 degrees. */
2860
2861 /* STRIPACK modules required by INTRC1: JRAND, LSTPTR, STORE, */
2862 /* TRFIND */
2863 /* (and optionally) GETNP if IFLGG .LE. 0 */
2864
2865 /* SSRFPACK modules required by INTRC1: ARCINT, ARCLEN, */
2866 /* FVAL, HVAL, SNHCSH */
2867 /* (and if IFLGG .LE. 0) APLYR, APLYRT, CONSTR, */
2868 /* GIVENS, GRADL, */
2869 /* ROTATE, SETUP */
2870
2871 /* Intrinsic functions called by INTRC1: COS, SIN, SQRT */
2872
2873 /* *********************************************************** */
2874
2875
2876 /* Local parameters: */
2877
2878 /* A = Angular separation between P and Q */
2879 /* B1,B2,B3 = Barycentric coordinates of the central projec- */
2880 /* tion of P onto the underlying planar triangle, */
2881 /* or (B1 and B2) projection of Q onto the */
2882 /* underlying line segment N1-N2 when P is */
2883 /* exterior. Unnormalized coordinates are */
2884 /* computed by TRFIND when P is in a triangle. */
2885 /* DUM = Dummy parameter for ARCINT */
2886 /* FQ,GQ = Interpolated value and gradient at Q */
2887 /* GQN = Negative of the component of GQ in the direction */
2888 /* Q->P */
2889 /* G1,G2,G3 = Gradients at I1, I2, and I3, or (G1 and G2) at */
2890 /* N1 and N2 */
2891 /* I = DO-loop index */
2892 /* IERR = Error flag for calls to GRADL */
2893 /* I1,I2,I3 = Vertex indexes returned by TRFIND */
2894 /* LP = LIST pointer */
2895 /* N1,N2 = Indexes of the endpoints of a boundary arc when */
2896 /* P is exterior (not contained in a triangle) */
2897 /* NN = Local copy of N */
2898 /* P = Cartesian coordinates of P */
2899 /* P1,P2,P3 = Cartesian coordinates of the vertices I1, I2, */
2900 /* and I3, or (P1 and P2) coordinates of N1 and */
2901 /* N2 if P is exterior */
2902 /* PTGQ = Scalar product (P,GQ) -- factor of the component */
2903 /* of GQ in the direction Q->P */
2904 /* PTN1 = Scalar product (P,N1) -- factor of B1 and B2 */
2905 /* PTN2 = Scalar product (P,N2) -- factor of B1 and B2 */
2906 /* Q = Closest boundary point to P when P is exterior */
2907 /* QNORM = Factor used to normalize Q */
2908 /* S1,S2,S3 = Tension factors associated with the triangle */
2909 /* sides opposite I1, I2, and I3, or (S1) the */
2910 /* boundary arc N1-N2 */
2911 /* S12 = Scalar product (N1,N2) -- factor of B1 and B2 */
2912 /* SUM = Quantity used to normalize the barycentric */
2913 /* coordinates */
2914
2915 /* Parameter adjustments */
2916 grad -= 4;
2917 --lend;
2918 --f;
2919 --z__;
2920 --y;
2921 --x;
2922 --list;
2923 --lptr;
2924 --sigma;
2925
2926 /* Function Body */
2927 nn = *n;
2928 if (nn < 3 || (*iflgg <= 0 && nn < 7) || *ist < 1 || *ist > nn) {
2929 goto L11;
2930 }
2931
2932 /* Transform (PLAT,PLON) to Cartesian coordinates. */
2933
2934 p[0] = cos(*plat) * cos(*plon);
2935 p[1] = cos(*plat) * sin(*plon);
2936 p[2] = sin(*plat);
2937
2938 /* Locate P with respect to the triangulation. */
2939
2940 sph_trfind_(ist, p, &nn, &x[1], &y[1], &z__[1], &list[1], &lptr[1], &lend[1],
2941 &b1, &b2, &b3, &i1, &i2, &i3);
2942 if (i1 == 0) {
2943 goto L12;
2944 }
2945 *ist = i1;
2946 if (i3 != 0) {
2947
2948 /* P is contained in the triangle (I1,I2,I3). Store the */
2949 /* vertex coordinates, gradients, and tension factors in */
2950 /* local variables. */
2951
2952 p1[0] = x[i1];
2953 p1[1] = y[i1];
2954 p1[2] = z__[i1];
2955 p2[0] = x[i2];
2956 p2[1] = y[i2];
2957 p2[2] = z__[i2];
2958 p3[0] = x[i3];
2959 p3[1] = y[i3];
2960 p3[2] = z__[i3];
2961 if (*iflgg > 0) {
2962
2963 /* Gradients are user-provided. */
2964
2965 for (i__ = 1; i__ <= 3; ++i__) {
2966 g1[i__ - 1] = grad[i__ + i1 * 3];
2967 g2[i__ - 1] = grad[i__ + i2 * 3];
2968 g3[i__ - 1] = grad[i__ + i3 * 3];
2969 /* L1: */
2970 }
2971 } else {
2972
2973 /* Compute gradient estimates at the vertices. */
2974
2975 gradl_(&nn, &i1, &x[1], &y[1], &z__[1], &f[1], &list[1], &lptr[1],
2976 &lend[1], g1, &ierr);
2977 if (ierr < 0) {
2978 goto L12;
2979 }
2980 gradl_(&nn, &i2, &x[1], &y[1], &z__[1], &f[1], &list[1], &lptr[1],
2981 &lend[1], g2, &ierr);
2982 if (ierr < 0) {
2983 goto L12;
2984 }
2985 gradl_(&nn, &i3, &x[1], &y[1], &z__[1], &f[1], &list[1], &lptr[1],
2986 &lend[1], g3, &ierr);
2987 if (ierr < 0) {
2988 goto L12;
2989 }
2990 }
2991
2992 if (*iflgs > 0) {
2993
2994 /* Variable tension: */
2995
2996 lp = lstptr_(&lend[i2], &i3, &list[1], &lptr[1]);
2997 s1 = sigma[lp];
2998 lp = lstptr_(&lend[i3], &i1, &list[1], &lptr[1]);
2999 s2 = sigma[lp];
3000 lp = lstptr_(&lend[i1], &i2, &list[1], &lptr[1]);
3001 s3 = sigma[lp];
3002 } else {
3003
3004 /* Uniform tension: */
3005
3006 s1 = sigma[1];
3007 s2 = s1;
3008 s3 = s1;
3009 }
3010
3011 /* Normalize the coordinates. */
3012
3013 sum = b1 + b2 + b3;
3014 b1 /= sum;
3015 b2 /= sum;
3016 b3 /= sum;
3017 *fp = fval_(&b1, &b2, &b3, p1, p2, p3, &f[i1], &f[i2], &f[i3], g1, g2,
3018 g3, &s1, &s2, &s3);
3019 *ier = 0;
3020 return 0;
3021 }
3022
3023 /* P is exterior to the triangulation, and I1 and I2 are */
3024 /* boundary nodes which are visible from P. Extrapolate to */
3025 /* P by linear (with respect to arc-length) interpolation */
3026 /* of the value and directional derivative (gradient comp- */
3027 /* onent in the direction Q->P) of the interpolatory */
3028 /* surface at Q where Q is the closest boundary point to P. */
3029
3030 /* Determine Q by traversing the boundary starting from I1. */
3031
3032 n1 = i1;
3033 ptn1 = p[0] * x[n1] + p[1] * y[n1] + p[2] * z__[n1];
3034 if (i1 != i2) {
3035 goto L3;
3036 }
3037
3038 /* All boundary nodes are visible from P. Find a boundary */
3039 /* arc N1->N2 such that P Left (N2 X N1)->N1. */
3040
3041 /* Counterclockwise boundary traversal: */
3042 /* Set N2 to the first neighbor of N1. */
3043
3044 L2:
3045 lp = lend[n1];
3046 lp = lptr[lp];
3047 n2 = list[lp];
3048
3049 /* Compute inner products (N1,N2) and (P,N2), and compute */
3050 /* B2 = Det(P,N1,N2 X N1). */
3051
3052 s12 = x[n1] * x[n2] + y[n1] * y[n2] + z__[n1] * z__[n2];
3053 ptn2 = p[0] * x[n2] + p[1] * y[n2] + p[2] * z__[n2];
3054 b2 = ptn2 - s12 * ptn1;
3055 if (b2 <= 0.) {
3056 goto L3;
3057 }
3058
3059 /* P Right (N2 X N1)->N1: iterate. */
3060
3061 n1 = n2;
3062 i1 = n1;
3063 ptn1 = ptn2;
3064 goto L2;
3065
3066 /* P Left (N2 X N1)->N1 where N2 is the first neighbor of N1. */
3067 /* Clockwise boundary traversal: */
3068
3069 L3:
3070 n2 = n1;
3071 ptn2 = ptn1;
3072
3073 /* Set N1 to the last neighbor of N2 and test for */
3074 /* termination. */
3075
3076 lp = lend[n2];
3077 n1 = -list[lp];
3078 if (n1 == i1) {
3079 goto L13;
3080 }
3081
3082 /* Compute inner products (N1,N2) and (P,N1). */
3083
3084 s12 = x[n1] * x[n2] + y[n1] * y[n2] + z__[n1] * z__[n2];
3085 ptn1 = p[0] * x[n1] + p[1] * y[n1] + p[2] * z__[n1];
3086
3087 /* Compute B2 = Det(P,N1,N2 X N1) = Det(Q,N1,N2 X N1)*(P,Q). */
3088
3089 b2 = ptn2 - s12 * ptn1;
3090 if (b2 <= 0.) {
3091 goto L3;
3092 }
3093
3094 /* Compute B1 = Det(P,N2 X N1,N2) = Det(Q,N2 X N1,N2)*(P,Q). */
3095
3096 b1 = ptn1 - s12 * ptn2;
3097 if (b1 <= 0.) {
3098
3099 /* Q = N2. Store value, coordinates, and gradient at Q. */
3100
3101 fq = f[n2];
3102 q[0] = x[n2];
3103 q[1] = y[n2];
3104 q[2] = z__[n2];
3105 if (*iflgg > 0) {
3106 for (i__ = 1; i__ <= 3; ++i__) {
3107 gq[i__ - 1] = grad[i__ + n2 * 3];
3108 /* L4: */
3109 }
3110 } else {
3111 gradl_(&nn, &n2, &x[1], &y[1], &z__[1], &f[1], &list[1], &lptr[1],
3112 &lend[1], gq, &ierr);
3113 if (ierr < 0) {
3114 goto L12;
3115 }
3116 }
3117
3118 /* Extrapolate to P: FP = FQ + A*(GQ,Q X (PXQ)/SIN(A)), */
3119 /* where A is the angular separation between Q and P, */
3120 /* and Sin(A) is the magnitude of P X Q. */
3121
3122 a = arclen_(q, p);
3123 ptgq = p[0] * gq[0] + p[1] * gq[1] + p[2] * gq[2];
3124 *fp = fq;
3125 if (a != 0.) {
3126 *fp += ptgq * a / sin(a);
3127 }
3128 *ier = 1;
3129 return 0;
3130 }
3131
3132 /* P Strictly Left (N2 X N1)->N2 and P Strictly Left */
3133 /* N1->(N2 X N1). Thus Q lies on the interior of N1->N2. */
3134 /* Store coordinates of N1 and N2 in local variables. */
3135
3136 p1[0] = x[n1];
3137 p1[1] = y[n1];
3138 p1[2] = z__[n1];
3139 p2[0] = x[n2];
3140 p2[1] = y[n2];
3141 p2[2] = z__[n2];
3142
3143 /* Compute the central projection of Q onto P2-P1 and */
3144 /* normalize to obtain Q. */
3145
3146 qnorm = 0.;
3147 for (i__ = 1; i__ <= 3; ++i__) {
3148 q[i__ - 1] = b1 * p1[i__ - 1] + b2 * p2[i__ - 1];
3149 qnorm += q[i__ - 1] * q[i__ - 1];
3150 /* L5: */
3151 }
3152 qnorm = sqrt(qnorm);
3153 for (i__ = 1; i__ <= 3; ++i__) {
3154 q[i__ - 1] /= qnorm;
3155 /* L6: */
3156 }
3157
3158 /* Store gradients at N1 and N2 and tension factor S1. */
3159
3160 if (*iflgg > 0) {
3161 for (i__ = 1; i__ <= 3; ++i__) {
3162 g1[i__ - 1] = grad[i__ + n1 * 3];
3163 g2[i__ - 1] = grad[i__ + n2 * 3];
3164 /* L7: */
3165 }
3166 } else {
3167 gradl_(&nn, &n1, &x[1], &y[1], &z__[1], &f[1], &list[1], &lptr[1], &
3168 lend[1], g1, &ierr);
3169 if (ierr < 0) {
3170 goto L12;
3171 }
3172 gradl_(&nn, &n2, &x[1], &y[1], &z__[1], &f[1], &list[1], &lptr[1], &
3173 lend[1], g2, &ierr);
3174 if (ierr < 0) {
3175 goto L12;
3176 }
3177 }
3178
3179 if (*iflgs <= 0) {
3180 s1 = sigma[1];
3181 }
3182 if (*iflgs >= 1) {
3183 s1 = sigma[lp];
3184 }
3185
3186 /* Compute an interpolated value and normal gradient */
3187 /* component at Q. */
3188
3189 arcint_(q, p1, p2, &f[n1], &f[n2], g1, g2, &s1, &fq, dum, &gqn);
3190
3191 /* Extrapolate to P: the normal gradient component GQN is */
3192 /* the negative of the component in the direction Q->P. */
3193
3194 *fp = fq - gqn * arclen_(q, p);
3195 *ier = 1;
3196 return 0;
3197
3198 /* N or IST is outside its valid range. */
3199
3200 L11:
3201 *ier = -1;
3202 return 0;
3203
3204 /* Collinear nodes encountered. */
3205
3206 L12:
3207 *ier = -2;
3208 return 0;
3209
3210 /* The distance between P and the closest boundary point */
3211 /* is at least 90 degrees. */
3212
3213 L13:
3214 *ier = -3;
3215 return 0;
3216 } /* intrc1_ */
3217
rotate_(integer * n,doublereal * c__,doublereal * s,doublereal * x,doublereal * y)3218 /* Subroutine */ int rotate_(integer *n, doublereal *c__, doublereal *s,
3219 doublereal *x, doublereal *y)
3220 {
3221 /* System generated locals */
3222 integer i__1;
3223
3224 /* Local variables */
3225 static integer i__;
3226 static doublereal xi, yi;
3227
3228
3229 /* *********************************************************** */
3230
3231 /* From SSRFPACK */
3232 /* Robert J. Renka */
3233 /* Dept. of Computer Science */
3234 /* Univ. of North Texas */
3235 /* renka@cs.unt.edu */
3236 /* 09/01/88 */
3237
3238 /* ( C S) */
3239 /* This subroutine applies the Givens rotation ( ) to */
3240 /* (-S C) */
3241 /* (X(1) ... X(N)) */
3242 /* the 2 by N matrix ( ) . */
3243 /* (Y(1) ... Y(N)) */
3244
3245 /* This routine is identical to Subroutine SROT from the */
3246 /* LINPACK BLAS (Basic Linear Algebra Subroutines). */
3247
3248 /* On input: */
3249
3250 /* N = Number of columns to be rotated. */
3251
3252 /* C,S = Elements of the Givens rotation. Refer to */
3253 /* Subroutine GIVENS. */
3254
3255 /* The above parameters are not altered by this routine. */
3256
3257 /* X,Y = Arrays of length .GE. N containing the compo- */
3258 /* nents of the vectors to be rotated. */
3259
3260 /* On output: */
3261
3262 /* X,Y = Arrays containing the rotated vectors (not */
3263 /* altered if N < 1). */
3264
3265 /* Modules required by ROTATE: None */
3266
3267 /* *********************************************************** */
3268
3269
3270 /* Parameter adjustments */
3271 --y;
3272 --x;
3273
3274 /* Function Body */
3275 i__1 = *n;
3276 for (i__ = 1; i__ <= i__1; ++i__) {
3277 xi = x[i__];
3278 yi = y[i__];
3279 x[i__] = *c__ * xi + *s * yi;
3280 y[i__] = -(*s) * xi + *c__ * yi;
3281 /* L1: */
3282 }
3283 return 0;
3284 } /* rotate_ */
3285
setup_(doublereal * xi,doublereal * yi,doublereal * wi,doublereal * wk,doublereal * s1,doublereal * s2,doublereal * wt,doublereal * row)3286 /* Subroutine */ int setup_(doublereal *xi, doublereal *yi, doublereal *wi,
3287 doublereal *wk, doublereal *s1, doublereal *s2, doublereal *wt,
3288 doublereal *row)
3289 {
3290 static doublereal w1, w2;
3291
3292
3293 /* *********************************************************** */
3294
3295 /* From SSRFPACK */
3296 /* Robert J. Renka */
3297 /* Dept. of Computer Science */
3298 /* Univ. of North Texas */
3299 /* renka@cs.unt.edu */
3300 /* 05/09/92 */
3301
3302 /* This subroutine sets up the I-th row of an augmented */
3303 /* regression matrix for a weighted least squares fit of a */
3304 /* quadratic function Q(X,Y) to a set of data values Wi, */
3305 /* where Q(0,0) = Wk. The first 3 columns (quadratic terms) */
3306 /* are scaled by 1/S2 and the fourth and fifth columns (lin- */
3307 /* ear terms) are scaled by 1/S1. */
3308
3309 /* On input: */
3310
3311 /* XI,YI = Coordinates of node I. */
3312
3313 /* WI = Data value at node I. */
3314
3315 /* WK = Data value interpolated by Q at the origin. */
3316
3317 /* S1,S2 = Inverse scale factors. */
3318
3319 /* WT = Weight factor corresponding to the I-th */
3320 /* equation. */
3321
3322 /* ROW = Array of length 6. */
3323
3324 /* Input parameters are not altered by this routine. */
3325
3326 /* On output: */
3327
3328 /* ROW = Array containing a row of the augmented re- */
3329 /* gression matrix. */
3330
3331 /* Modules required by SETUP: None */
3332
3333 /* *********************************************************** */
3334
3335
3336 /* Local parameters: */
3337
3338 /* W1 = Weighted scale factor for the linear terms */
3339 /* W2 = Weighted scale factor for the quadratic terms */
3340
3341 /* Parameter adjustments */
3342 --row;
3343
3344 /* Function Body */
3345 w1 = *wt / *s1;
3346 w2 = *wt / *s2;
3347 row[1] = *xi * *xi * w2;
3348 row[2] = *xi * *yi * w2;
3349 row[3] = *yi * *yi * w2;
3350 row[4] = *xi * w1;
3351 row[5] = *yi * w1;
3352 row[6] = (*wi - *wk) * *wt;
3353 return 0;
3354 } /* setup_ */
3355
sig0_(integer * n1,integer * n2,integer * n,doublereal * x,doublereal * y,doublereal * z__,doublereal * h__,integer * list,integer * lptr,integer * lend,doublereal * grad,integer * iflgb,doublereal * hbnd,doublereal * tol,integer * iflgs,doublereal * sigma,integer * ier)3356 doublereal sig0_(integer *n1, integer *n2, integer *n, doublereal *x,
3357 doublereal *y, doublereal *z__, doublereal *h__, integer *list,
3358 integer *lptr, integer *lend, doublereal *grad, integer *iflgb,
3359 doublereal *hbnd, doublereal *tol, integer *iflgs, doublereal *sigma,
3360 integer *ier)
3361 {
3362 /* Initialized data */
3363
3364 static doublereal sbig = 85.;
3365
3366 /* System generated locals */
3367 integer i__1;
3368 doublereal ret_val, d__1, d__2, d__3, d__4, d__5, d__6;
3369
3370 /* Builtin functions */
3371 double sqrt(doublereal), d_sign(doublereal *, doublereal *), exp(
3372 doublereal), log(doublereal);
3373
3374 /* Local variables */
3375 static doublereal a, b, c__, d__, e, f, r__, s, t, a0, b0, c1, c2, d0, d2,
3376 f0, h1, h2, p1[3], p2[3], s1, s2, t0, t1, t2, aa, al, rf, tm, un[
3377 3];
3378 static integer lp1, lp2;
3379 static doublereal bnd, scm, sig, ems;
3380 static integer lpl, nit;
3381 static doublereal ssm, d1pd2, fneg, dsig, dmax__, fmax, sneg, ftol, rsig,
3382 rtol, stol, coshm, sinhm, ssinh;
3383 extern doublereal store_(doublereal *);
3384 static doublereal unorm;
3385 extern doublereal arclen_(doublereal *, doublereal *);
3386 static doublereal coshmm;
3387 extern /* Subroutine */ int snhcsh_(doublereal *, doublereal *,
3388 doublereal *, doublereal *);
3389
3390
3391 /* *********************************************************** */
3392
3393 /* From SSRFPACK */
3394 /* Robert J. Renka */
3395 /* Dept. of Computer Science */
3396 /* Univ. of North Texas */
3397 /* renka@cs.unt.edu */
3398 /* 11/21/96 */
3399
3400 /* Given a triangulation of a set of nodes on the unit */
3401 /* sphere, along with data values H and gradients GRAD at the */
3402 /* nodes, this function determines the smallest tension fac- */
3403 /* tor SIG0 such that the Hermite interpolatory tension */
3404 /* spline H(A), defined by SIG0 and the endpoint values and */
3405 /* directional derivatives associated with an arc N1-N2, is */
3406 /* bounded (either above or below) by HBND for all A in */
3407 /* (A1,A2), where (A1,A2) denotes an interval corresponding */
3408 /* to the arc and A is the arc-length. */
3409
3410 /* On input: */
3411
3412 /* N1,N2 = Nodal indexes of the endpoints of an arc for */
3413 /* which the tension factor is to be computed. */
3414 /* The indexes must be distinct and lie in the */
3415 /* range 1 to N, and if IFLGS .GE. 1, they must */
3416 /* correspond to adjacent nodes in the triangu- */
3417 /* lation. */
3418
3419 /* N = Number of nodes in the triangulation. N .GE. 3. */
3420
3421 /* X,Y,Z = Arrays of length N containing coordinates of */
3422 /* the nodes. X(I)**2 + Y(I)**2 + Z(I)**2 = 1. */
3423
3424 /* H = Array of length N containing data values at the */
3425 /* nodes. H(I) is associated with (X(I),Y(I),Z(I)) */
3426 /* for I = 1 to N. */
3427
3428 /* LIST,LPTR,LEND = Data structure defining the trian- */
3429 /* gulation. Refer to STRIPACK */
3430 /* Subroutine TRMESH. */
3431
3432 /* GRAD = Array dimensioned 3 by N whose columns con- */
3433 /* tain gradients at the nodes. GRAD( ,J) must */
3434 /* be orthogonal to node J: GRAD(1,J)*X(J) + */
3435 /* GRAD(2,J)*Y(J) + GRAD(3,J)*Z(J) = 0. Refer */
3436 /* to Subroutines GRADG, GRADL, and SMSURF. */
3437
3438 /* IFLGB = Bound option indicator: */
3439 /* IFLGB = -1 if HBND is a lower bound on H. */
3440 /* IFLGB = 1 if HBND is an upper bound on H. */
3441
3442 /* HBND = Bound on H. HBND .LE. min(H1,H2) if IFLGB = */
3443 /* -1 and HBND .GE. max(H1,H2) if IFLGB = 1, */
3444 /* where H1 and H2 are the data values at the */
3445 /* endpoints of the arc N1-N2. */
3446
3447 /* TOL = Tolerance whose magnitude determines how close */
3448 /* SIG0 is to its optimal value when nonzero */
3449 /* finite tension is necessary and sufficient to */
3450 /* satisfy the constraint. For a lower bound, */
3451 /* SIG0 is chosen so that HBND .LE. HMIN .LE. */
3452 /* HBND + abs(TOL), where HMIN is the minimum */
3453 /* value of H on the arc, and for an upper bound, */
3454 /* the maximum of H satisfies HBND - abs(TOL) */
3455 /* .LE. HMAX .LE. HBND. Thus, the constraint is */
3456 /* satisfied but possibly with more tension than */
3457 /* necessary. */
3458
3459 /* IFLGS = Tension array option indicator: */
3460 /* IFLGS .LE. 0 if SIGMA is not to be used. */
3461 /* IFLGS .GE. 1 if SIGMA is to be updated by */
3462 /* storing SIG0 in the appropriate */
3463 /* locations. */
3464
3465 /* The above parameters are not altered by this function. */
3466
3467 /* SIGMA = Dummy parameter (IFLGS .LE. 0) or array con- */
3468 /* taining tension factors associated with arcs */
3469 /* in one-to-one correspondence with LIST */
3470 /* entries (IFLGS .GE. 1). Refer to Subroutine */
3471 /* GETSIG. */
3472
3473 /* On output: */
3474
3475 /* SIGMA = Tension factor array updated with the new */
3476 /* value if and only if IFLGS .GE. 1 and IER */
3477 /* .GE. 0. */
3478
3479 /* IER = Error indicator: */
3480 /* IER = 0 if no errors were encountered and the */
3481 /* constraint can be satisfied with fin- */
3482 /* ite tension. */
3483 /* IER = 1 if no errors were encountered but in- */
3484 /* finite tension is required to satisfy */
3485 /* the constraint (e.g., IFLGB = -1, HBND */
3486 /* = H(A1), and the directional deriva- */
3487 /* tive of H at A1 is negative). */
3488 /* IER = -1 if N1, N2, N, or IFLGB is outside its */
3489 /* valid range. */
3490 /* IER = -2 if nodes N1 and N2 coincide or IFLGS */
3491 /* .GE. 1 and the nodes are not adja- */
3492 /* cent. */
3493 /* IER = -3 if HBND is outside its valid range. */
3494
3495 /* SIG0 = Minimum tension factor defined above unless */
3496 /* IER < 0, in which case SIG0 = -1. If IER */
3497 /* = 1, SIG0 is set to 85, resulting in an */
3498 /* approximation to the linear interpolant of */
3499 /* the endpoint values. */
3500
3501 /* STRIPACK module required by SIG0: STORE */
3502
3503 /* SSRFPACK modules required by SIG0: ARCLEN, SNHCSH */
3504
3505 /* Intrinsic functions called by SIG0: ABS, EXP, LOG, MAX, */
3506 /* MIN, REAL, SIGN, */
3507 /* SQRT */
3508
3509 /* *********************************************************** */
3510
3511
3512 /* Parameter adjustments */
3513 grad -= 4;
3514 --lend;
3515 --h__;
3516 --z__;
3517 --y;
3518 --x;
3519 --list;
3520 --lptr;
3521 --sigma;
3522
3523 /* Function Body */
3524 rf = (doublereal) (*iflgb);
3525 bnd = *hbnd;
3526
3527 /* Print a heading. */
3528
3529 /* IF (LUN .GE. 0 .AND. RF .LT. 0.) WRITE (LUN,100) N1, */
3530 /* . N2, BND */
3531 /* IF (LUN .GE. 0 .AND. RF .GT. 0.) WRITE (LUN,110) N1, */
3532 /* . N2, BND */
3533 /* 100 FORMAT (//1X,'SIG0 -- N1 =',I4,', N2 =',I4, */
3534 /* . ', LOWER BOUND = ',E15.8) */
3535 /* 110 FORMAT (//1X,'SIG0 -- N1 =',I4,', N2 =',I4, */
3536 /* . ', UPPER BOUND = ',E15.8) */
3537 if (dbg_verbose) if (rf < 0.0) fprintf (stderr, "SIG0 -- N1 = %d N2 = %d LOWER BOUND = %g\n", *n1, *n2, bnd);
3538 if (dbg_verbose) if (rf > 0.0) fprintf (stderr, "SIG0 -- N1 = %d N2 = %d UPPER BOUND = %g\n", *n1, *n2, bnd);
3539
3540 /* Test for errors and store local parameters. */
3541
3542 *ier = -1;
3543 /* Computing MAX */
3544 i__1 = max(*n1,*n2);
3545 if (min(*n1,*n2) < 1 || *n1 == *n2 || max(i__1,3) > *n || fabs(rf) != 1.) {
3546 goto L11;
3547 }
3548 *ier = -2;
3549 if (*iflgs > 0) {
3550
3551 /* Set LP1 and LP2 to the pointers to N2 as a neighbor of */
3552 /* N1 and N1 as a neighbor of N2, respectively. */
3553
3554 lpl = lend[*n1];
3555 lp1 = lptr[lpl];
3556 L1:
3557 if (list[lp1] == *n2) {
3558 goto L2;
3559 }
3560 lp1 = lptr[lp1];
3561 if (lp1 != lpl) {
3562 goto L1;
3563 }
3564 if ((i__1 = list[lp1], abs(i__1)) != *n2) {
3565 goto L11;
3566 }
3567
3568 L2:
3569 lpl = lend[*n2];
3570 lp2 = lptr[lpl];
3571 L3:
3572 if (list[lp2] == *n1) {
3573 goto L4;
3574 }
3575 lp2 = lptr[lp2];
3576 if (lp2 != lpl) {
3577 goto L3;
3578 }
3579 if ((i__1 = list[lp2], abs(i__1)) != *n1) {
3580 goto L11;
3581 }
3582 }
3583
3584 /* Store nodal coordinates P1 and P2, compute arc-length AL */
3585 /* and unit normal UN = (P1 X P2)/UNORM, and test for */
3586 /* coincident nodes. */
3587
3588 L4:
3589 p1[0] = x[*n1];
3590 p1[1] = y[*n1];
3591 p1[2] = z__[*n1];
3592 p2[0] = x[*n2];
3593 p2[1] = y[*n2];
3594 p2[2] = z__[*n2];
3595 al = arclen_(p1, p2);
3596 un[0] = p1[1] * p2[2] - p1[2] * p2[1];
3597 un[1] = p1[2] * p2[0] - p1[0] * p2[2];
3598 un[2] = p1[0] * p2[1] - p1[1] * p2[0];
3599 unorm = sqrt(un[0] * un[0] + un[1] * un[1] + un[2] * un[2]);
3600 if (unorm == 0. || al == 0.) {
3601 goto L11;
3602 }
3603
3604 /* Store endpoint data values and test for valid constraint. */
3605
3606 h1 = h__[*n1];
3607 h2 = h__[*n2];
3608 *ier = -3;
3609 if ((rf < 0. && min(h1,h2) < bnd) || (rf > 0. && bnd < max(h1,h2))) {
3610 goto L11;
3611 }
3612
3613 /* Compute scaled directional derivatives S1,S2 at the end- */
3614 /* points (for the direction N1->N2) and test for infinite */
3615 /* tension required. */
3616
3617 s1 = al * (grad[*n1 * 3 + 1] * p2[0] + grad[*n1 * 3 + 2] * p2[1] + grad[*
3618 n1 * 3 + 3] * p2[2]) / unorm;
3619 s2 = -al * (grad[*n2 * 3 + 1] * p1[0] + grad[*n2 * 3 + 2] * p1[1] + grad[*
3620 n2 * 3 + 3] * p1[2]) / unorm;
3621 *ier = 1;
3622 sig = sbig;
3623 if ((h1 == bnd && rf * s1 > 0.) || (h2 == bnd && rf * s2 < 0.)) {
3624 goto L10;
3625 }
3626
3627 /* Test for SIG = 0 sufficient. */
3628
3629 *ier = 0;
3630 sig = 0.;
3631 if (rf * s1 <= 0. && rf * s2 >= 0.) {
3632 goto L10;
3633 }
3634
3635 /* Compute first difference S and coefficients A0 and B0 */
3636 /* of the Hermite cubic interpolant H0(A) = H2 - (S2*R + */
3637 /* B0*R**2 + (A0/3)*R**3), where R(A) = (A2-A)/AL. */
3638
3639 s = h2 - h1;
3640 t0 = s * 3. - s1 - s2;
3641 a0 = (s - t0) * 3.;
3642 b0 = t0 - s2;
3643 d0 = t0 * t0 - s1 * s2;
3644
3645 /* H0 has local extrema in (A1,A2) iff S1*S2 < 0 or */
3646 /* (T0*(S1+S2) < 0 and D0 .GE. 0). */
3647
3648 if (s1 * s2 >= 0. && (t0 * (s1 + s2) >= 0. || d0 < 0.)) {
3649 goto L10;
3650 }
3651 if (a0 == 0.) {
3652
3653 /* H0 is quadratic and has an extremum at R = -S2/(2*B0). */
3654 /* H0(R) = H2 + S2**2/(4*B0). Note that A0 = 0 implies */
3655 /* 2*B0 = S1-S2, and S1*S2 < 0 implies B0 .NE. 0. */
3656 /* Also, the extremum is a min iff HBND is a lower bound. */
3657
3658 f0 = (bnd - h2 - s2 * s2 / (b0 * 4.)) * rf;
3659 } else {
3660
3661 /* A0 .NE. 0 and H0 has extrema at R = (-B0 +/- SQRT(D0))/ */
3662 /* A0 = S2/(-B0 -/+ SQRT(D0)), where the negative root */
3663 /* corresponds to a min. The expression for R is chosen */
3664 /* to avoid cancellation error. H0(R) = H2 + (S2*B0 + */
3665 /* 2*D0*R)/(3*A0). */
3666
3667 d__1 = sqrt(d0);
3668 t = -b0 - d_sign(&d__1, &b0);
3669 r__ = t / a0;
3670 if (rf * b0 > 0.) {
3671 r__ = s2 / t;
3672 }
3673 f0 = (bnd - h2 - (s2 * b0 + d0 * 2. * r__) / (a0 * 3.)) * rf;
3674 }
3675
3676 /* F0 .GE. 0 iff SIG = 0 is sufficient to satisfy the */
3677 /* constraint. */
3678
3679 if (f0 >= 0.) {
3680 goto L10;
3681 }
3682
3683 /* Find a zero of F(SIG) = (BND-H(R))*RF where the derivative */
3684 /* of H, HP, vanishes at R. F is a nondecreasing function, */
3685 /* F(0) < 0, and F = FMAX for SIG sufficiently large. */
3686
3687 /* Initialize parameters for the secant method. The method */
3688 /* uses three points: (SG0,F0), (SIG,F), and (SNEG,FNEG), */
3689 /* where SG0 and SNEG are defined implicitly by DSIG = SIG */
3690 /* - SG0 and DMAX = SIG - SNEG. SG0 is initially zero and */
3691 /* SNEG is initialized to a sufficiently large value that */
3692 /* FNEG > 0. This value is used only if the initial value */
3693 /* of F is negative. */
3694
3695 /* Computing MAX */
3696 /* Computing MIN */
3697 d__5 = (d__1 = h1 - bnd, fabs(d__1)), d__6 = (d__2 = h2 - bnd, fabs(d__2));
3698 d__3 = .001, d__4 = min(d__5,d__6);
3699 fmax = max(d__3,d__4);
3700 /* Computing MAX */
3701 d__3 = (d__1 = h1 - bnd, fabs(d__1)), d__4 = (d__2 = h2 - bnd, fabs(d__2));
3702 t = max(d__3,d__4);
3703 /* Computing MAX */
3704 d__1 = fabs(s1), d__2 = fabs(s2);
3705 sig = max(d__1,d__2) / t;
3706 dmax__ = sig * (1. - t / fmax);
3707 sneg = sig - dmax__;
3708 /* IF (LUN .GE. 0) WRITE (LUN,120) SIG, SNEG, F0, FMAX */
3709 /* 120 FORMAT (1X,8X,'SIG = ',E15.8,', SNEG = ',E15.8/ */
3710 /* . 1X,9X,'F0 = ',E15.8,', FMAX = ',E15.8/) */
3711 if (dbg_verbose) fprintf (stderr, "SIG = %g SNEG = %g F0 = %g FMAX = %g\n", sig, sneg, f0, fmax);
3712 dsig = sig;
3713 fneg = fmax;
3714 d2 = s2 - s;
3715 d1pd2 = s2 - s1;
3716 nit = 0;
3717
3718 /* Compute an absolute tolerance FTOL = abs(TOL) and a */
3719 /* relative tolerance RTOL = 100*Macheps. */
3720
3721 ftol = fabs(*tol);
3722 rtol = 1.;
3723 L5:
3724 rtol /= 2.;
3725 d__1 = rtol + 1.;
3726 if (store_(&d__1) > 1.) {
3727 goto L5;
3728 }
3729 rtol *= 200.;
3730
3731 /* Top of loop: compute F. */
3732
3733 L6:
3734 ems = exp(-sig);
3735 if (sig <= .5) {
3736
3737 /* Use approximations designed to avoid cancellation error */
3738 /* (associated with small SIG) in the modified hyperbolic */
3739 /* functions. */
3740
3741 snhcsh_(&sig, &sinhm, &coshm, &coshmm);
3742 c1 = sig * coshm * d2 - sinhm * d1pd2;
3743 c2 = sig * (sinhm + sig) * d2 - coshm * d1pd2;
3744 a = c2 - c1;
3745 aa = a / ems;
3746 e = sig * sinhm - coshmm - coshmm;
3747 } else {
3748
3749 /* Scale SINHM and COSHM by 2*exp(-SIG) in order to avoid */
3750 /* overflow. */
3751
3752 tm = 1. - ems;
3753 ssinh = tm * (ems + 1.);
3754 ssm = ssinh - sig * 2. * ems;
3755 scm = tm * tm;
3756 c1 = sig * scm * d2 - ssm * d1pd2;
3757 c2 = sig * ssinh * d2 - scm * d1pd2;
3758 aa = (sig * tm * d2 + (tm - sig) * d1pd2) * 2.;
3759 a = ems * aa;
3760 e = sig * ssinh - scm - scm;
3761 }
3762
3763 /* HP(R) = (S2 - (C1*sinh(SIG*R) - C2*coshm(SIG*R))/E)/DT */
3764 /* = 0 for ESR = (-B +/- sqrt(D))/A = C/(-B -/+ sqrt(D)) */
3765 /* where ESR = exp(SIG*R), A = C2-C1, D = B**2 - A*C, and */
3766 /* B and C are defined below. */
3767
3768 b = e * s2 - c2;
3769 c__ = c2 + c1;
3770 d__ = b * b - a * c__;
3771 f = 0.;
3772 if (aa * c__ == 0. && b == 0.) {
3773 goto L7;
3774 }
3775 f = fmax;
3776 if (d__ < 0.) {
3777 goto L7;
3778 }
3779 t1 = sqrt(d__);
3780 t = -b - d_sign(&t1, &b);
3781 rsig = 0.;
3782 if (rf * b < 0. && aa != 0.) {
3783 if (t / aa > 0.) {
3784 rsig = sig + log(t / aa);
3785 }
3786 }
3787 if ((rf * b > 0. || aa == 0.) && c__ / t > 0.) {
3788 rsig = log(c__ / t);
3789 }
3790 if ((rsig <= 0. || rsig >= sig) && b != 0.) {
3791 goto L7;
3792 }
3793
3794 /* H(R) = H2 - (B*SIG*R + C1 + RF*sqrt(D))/(SIG*E). */
3795
3796 f = (bnd - h2 + (b * rsig + c1 + rf * t1) / (sig * e)) * rf;
3797
3798 /* Update the number of iterations NIT. */
3799
3800 L7:
3801 ++nit;
3802 /* IF (LUN .GE. 0) WRITE (LUN,130) NIT, SIG, F */
3803 /* 130 FORMAT (1X,3X,I2,' -- SIG = ',E15.8,', F = ', */
3804 /* . E15.8) */
3805 if (dbg_verbose) fprintf (stderr, "%d -- SIG = %g F = %g\n", nit, sig, f);
3806 if (f0 * f < 0.) {
3807
3808 /* F0*F < 0. Update (SNEG,FNEG) to (SG0,F0) so that F and */
3809 /* FNEG always have opposite signs. If SIG is closer to */
3810 /* SNEG than SG0, then swap (SNEG,FNEG) with (SG0,F0). */
3811
3812 t1 = dmax__;
3813 t2 = fneg;
3814 dmax__ = dsig;
3815 fneg = f0;
3816 if (fabs(dsig) > fabs(t1)) {
3817
3818 dsig = t1;
3819 f0 = t2;
3820 }
3821 }
3822
3823 /* Test for convergence. */
3824
3825 stol = rtol * sig;
3826 if (fabs(dmax__) <= stol || (f >= 0. && f <= ftol) || fabs(f) <= rtol) {
3827 goto L10;
3828 }
3829
3830 /* Test for F0 = F = FMAX or F < 0 on the first iteration. */
3831
3832 if (f0 != f && (nit > 1 || f > 0.)) {
3833 goto L9;
3834 }
3835
3836 /* F*F0 > 0 and either the new estimate would be outside */
3837 /* of the bracketing interval of length abs(DMAX) or */
3838 /* F < 0 on the first iteration. Reset (SG0,F0) to */
3839 /* (SNEG,FNEG). */
3840
3841 L8:
3842 dsig = dmax__;
3843 f0 = fneg;
3844
3845 /* Compute the change in SIG by linear interpolation */
3846 /* between (SG0,F0) and (SIG,F). */
3847
3848 L9:
3849 dsig = -f * dsig / (f - f0);
3850 /* IF (LUN .GE. 0) WRITE (LUN,140) DSIG */
3851 /* 140 FORMAT (1X,8X,'DSIG = ',E15.8) */
3852 if (dbg_verbose) fprintf (stderr, "DSIG = %g\n", dsig);
3853 if (fabs(dsig) > fabs(dmax__) || dsig * dmax__ > 0.) {
3854 goto L8;
3855 }
3856
3857 /* Restrict the step-size such that abs(DSIG) .GE. STOL/2. */
3858 /* Note that DSIG and DMAX have opposite signs. */
3859
3860 if (fabs(dsig) < stol / 2.) {
3861 d__1 = stol / 2.;
3862 dsig = -d_sign(&d__1, &dmax__);
3863 }
3864
3865 /* Bottom of loop: Update SIG, DMAX, and F0. */
3866
3867 sig += dsig;
3868 dmax__ += dsig;
3869 f0 = f;
3870 goto L6;
3871
3872 /* No errors encountered. */
3873
3874 L10:
3875 ret_val = sig;
3876 if (*iflgs <= 0) {
3877 return ret_val;
3878 }
3879 sigma[lp1] = sig;
3880 sigma[lp2] = sig;
3881 return ret_val;
3882
3883 /* Error termination. */
3884
3885 L11:
3886 ret_val = -1.;
3887 return ret_val;
3888 } /* sig0_ */
3889
sig1_(integer * n1,integer * n2,integer * n,doublereal * x,doublereal * y,doublereal * z__,doublereal * h__,integer * list,integer * lptr,integer * lend,doublereal * grad,integer * iflgb,doublereal * hpbnd,doublereal * tol,integer * iflgs,doublereal * sigma,integer * ier)3890 doublereal sig1_(integer *n1, integer *n2, integer *n, doublereal *x,
3891 doublereal *y, doublereal *z__, doublereal *h__, integer *list,
3892 integer *lptr, integer *lend, doublereal *grad, integer *iflgb,
3893 doublereal *hpbnd, doublereal *tol, integer *iflgs, doublereal *sigma,
3894 integer *ier)
3895 {
3896 /* Initialized data */
3897
3898 static doublereal sbig = 85.;
3899
3900 /* System generated locals */
3901 integer i__1;
3902 doublereal ret_val, d__1, d__2;
3903
3904 /* Builtin functions */
3905 double sqrt(doublereal), exp(doublereal), d_sign(doublereal *, doublereal
3906 *);
3907
3908 /* Local variables */
3909 static doublereal a, e, f, s, a0, b0, c0, c1, c2, d0, d1, d2, f0, p1[3],
3910 p2[3], s1, s2, t0, t1, t2, al, rf, tm, un[3];
3911 static integer lp1, lp2;
3912 static doublereal bnd, sig, ems;
3913 static integer lpl, nit;
3914 static doublereal ems2, d1pd2, fneg, dsig, dmax__, fmax, sinh__, ftol,
3915 rtol, stol, coshm, sinhm;
3916 extern doublereal store_(doublereal *);
3917 static doublereal unorm;
3918 extern doublereal arclen_(doublereal *, doublereal *);
3919 static doublereal coshmm;
3920 extern /* Subroutine */ int snhcsh_(doublereal *, doublereal *,
3921 doublereal *, doublereal *);
3922
3923
3924 /* *********************************************************** */
3925
3926 /* From SSRFPACK */
3927 /* Robert J. Renka */
3928 /* Dept. of Computer Science */
3929 /* Univ. of North Texas */
3930 /* renka@cs.unt.edu */
3931 /* 11/21/96 */
3932
3933 /* Given a triangulation of a set of nodes on the unit */
3934 /* sphere, along with data values H and gradients GRAD at the */
3935 /* nodes, this function determines the smallest tension fac- */
3936 /* tor SIG1 such that the first derivative HP(A) of the */
3937 /* Hermite interpolatory tension spline H(A), defined by SIG1 */
3938 /* and the endpoint values and directional derivatives asso- */
3939 /* ciated with an arc N1-N2, is bounded (either above or */
3940 /* below) by HPBND for all A in (A1,A2), where (A1,A2) de- */
3941 /* notes an interval corresponding to the arc and A denotes */
3942 /* arc-length. */
3943
3944 /* On input: */
3945
3946 /* N1,N2 = Nodal indexes of the endpoints of an arc for */
3947 /* which the tension factor is to be computed. */
3948 /* The indexes must be distinct and lie in the */
3949 /* range 1 to N, and if IFLGS .GE. 1, they must */
3950 /* correspond to adjacent nodes in the triangu- */
3951 /* lation. */
3952
3953 /* N = Number of nodes in the triangulation. N .GE. 3. */
3954
3955 /* X,Y,Z = Arrays of length N containing coordinates of */
3956 /* the nodes. X(I)**2 + Y(I)**2 + Z(I)**2 = 1. */
3957
3958 /* H = Array of length N containing data values at the */
3959 /* nodes. H(I) is associated with (X(I),Y(I),Z(I)) */
3960 /* for I = 1 to N. */
3961
3962 /* LIST,LPTR,LEND = Data structure defining the trian- */
3963 /* gulation. Refer to STRIPACK */
3964 /* Subroutine TRMESH. */
3965
3966 /* GRAD = Array dimensioned 3 by N whose columns con- */
3967 /* gradients at the nodes. GRAD( ,J) must be */
3968 /* orthogonal to node J: GRAD(1,J)*X(J) + */
3969 /* GRAD(2,J)*Y(J) + GRAD(3,J)*Z(J) = 0. Refer */
3970 /* to Subroutines GRADG, GRADL, and SMSURF. */
3971
3972 /* IFLGB = Bound option indicator: */
3973 /* IFLGB = -1 if HPBND is a lower bound on HP. */
3974 /* IFLGB = 1 if HPBND is an upper bound on HP. */
3975
3976 /* HPBND = Bound on HP. HPBND .LE. min(HP1,HP2,S) if */
3977 /* IFLGB = -1 and HPBND .GE. max(HP1,HP2,S) if */
3978 /* IFLGB = 1, where HP1 and HP2 are the direc- */
3979 /* tional derivatives at the endpoints of the */
3980 /* arc N1-N2, and S is the slope of the linear */
3981 /* interpolant of the endpoint data values. */
3982
3983 /* TOL = Tolerance whose magnitude determines how close */
3984 /* SIG1 is to its optimal value when nonzero */
3985 /* finite tension is necessary and sufficient to */
3986 /* satisfy the constraint. For a lower bound, */
3987 /* SIG1 is chosen so that HPBND .LE. HPMIN .LE. */
3988 /* HPBND + abs(TOL), where HPMIN is the minimum */
3989 /* value of HP on the arc. For an upper bound, */
3990 /* the maximum of HP satisfies HPBND - abs(TOL) */
3991 /* .LE. HPMAX .LE. HPBND. Thus, the constraint */
3992 /* is satisfied but possibly with more tension */
3993 /* than necessary. */
3994
3995 /* IFLGS = Tension array option indicator: */
3996 /* IFLGS .LE. 0 if SIGMA is not to be used. */
3997 /* IFLGS .GE. 1 if SIGMA is to be updated by */
3998 /* storing SIG1 in the appropriate */
3999 /* locations. */
4000
4001 /* The above parameters are not altered by this function. */
4002
4003 /* SIGMA = Dummy parameter (IFLGS .LE. 0) or array */
4004 /* containing tension factors associated with */
4005 /* arcs in one-to-one correspondence with LIST */
4006 /* entries (IFLGS .GE. 1). Refer to Subroutine */
4007 /* GETSIG. */
4008
4009 /* On output: */
4010
4011 /* SIGMA = Tension factor array updated with the new */
4012 /* value if and only if IFLGS .GE. 1 and IER */
4013 /* .GE. 0. */
4014
4015 /* IER = Error indicator: */
4016 /* IER = 0 if no errors were encountered and the */
4017 /* constraint can be satisfied with fin- */
4018 /* ite tension. */
4019 /* IER = 1 if no errors were encountered but in- */
4020 /* finite tension is required to satisfy */
4021 /* the constraint (e.g., IFLGB = -1, */
4022 /* HPBND = S, and HP1 > S). */
4023 /* IER = -1 if N1, N2, N, or IFLGB is outside its */
4024 /* valid range. */
4025 /* IER = -2 if nodes N1 and N2 coincide or IFLGS */
4026 /* .GE. 1 and the nodes are not adja- */
4027 /* cent. */
4028 /* IER = -3 if HPBND is outside its valid range. */
4029
4030 /* SIG1 = Minimum tension factor defined above unless */
4031 /* IER < 0, in which case SIG1 = -1. If IER */
4032 /* = 1, SIG1 is set to 85, resulting in an */
4033 /* approximation to the linear interpolant of */
4034 /* the endpoint values. */
4035
4036 /* STRIPACK module required by SIG1: STORE */
4037
4038 /* SSRFPACK modules required by SIG1: ARCLEN, SNHCSH */
4039
4040 /* Intrinsic functions called by SIG1: ABS, EXP, MAX, MIN, */
4041 /* REAL, SIGN, SQRT */
4042
4043 /* *********************************************************** */
4044
4045
4046 /* Parameter adjustments */
4047 grad -= 4;
4048 --lend;
4049 --h__;
4050 --z__;
4051 --y;
4052 --x;
4053 --list;
4054 --lptr;
4055 --sigma;
4056
4057 /* Function Body */
4058 rf = (doublereal) (*iflgb);
4059 bnd = *hpbnd;
4060
4061 /* Print a heading. */
4062
4063 /* IF (LUN .GE. 0 .AND. RF .LT. 0.) WRITE (LUN,100) N1, */
4064 /* . N2, BND */
4065 /* IF (LUN .GE. 0 .AND. RF .GT. 0.) WRITE (LUN,110) N1, */
4066 /* . N2, BND */
4067 /* 100 FORMAT (//1X,'SIG1 -- N1 =',I4,', N2 =',I4, */
4068 /* . ', LOWER BOUND = ',E15.8) */
4069 /* 110 FORMAT (//1X,'SIG1 -- N1 =',I4,', N2 =',I4, */
4070 /* . ', UPPER BOUND = ',E15.8) */
4071 if (dbg_verbose) if (rf < 0.0) fprintf (stderr, "SIG1 -- N1 = %d N2 = %d LOWER BOUND = %g\n", *n1, *n2, bnd);
4072 if (dbg_verbose) if (rf > 0.0) fprintf (stderr, "SIG1 -- N1 = %d N2 = %d UPPER BOUND = %g\n", *n1, *n2, bnd);
4073
4074 /* Test for errors and store local parameters. */
4075
4076 *ier = -1;
4077 /* Computing MAX */
4078 i__1 = max(*n1,*n2);
4079 if (min(*n1,*n2) < 1 || *n1 == *n2 || max(i__1,3) > *n || fabs(rf) != 1.) {
4080 goto L11;
4081 }
4082 *ier = -2;
4083 if (*iflgs > 0) {
4084
4085 /* Set LP1 and LP2 to the pointers to N2 as a neighbor of */
4086 /* N1 and N1 as a neighbor of N2, respectively. */
4087
4088 lpl = lend[*n1];
4089 lp1 = lptr[lpl];
4090 L1:
4091 if (list[lp1] == *n2) {
4092 goto L2;
4093 }
4094 lp1 = lptr[lp1];
4095 if (lp1 != lpl) {
4096 goto L1;
4097 }
4098 if ((i__1 = list[lp1], abs(i__1)) != *n2) {
4099 goto L11;
4100 }
4101
4102 L2:
4103 lpl = lend[*n2];
4104 lp2 = lptr[lpl];
4105 L3:
4106 if (list[lp2] == *n1) {
4107 goto L4;
4108 }
4109 lp2 = lptr[lp2];
4110 if (lp2 != lpl) {
4111 goto L3;
4112 }
4113 if ((i__1 = list[lp2], abs(i__1)) != *n1) {
4114 goto L11;
4115 }
4116 }
4117
4118 /* Store nodal coordinates P1 and P2, compute arc-length AL */
4119 /* and unit normal UN = (P1 X P2)/UNORM, and test for */
4120 /* coincident nodes. */
4121
4122 L4:
4123 p1[0] = x[*n1];
4124 p1[1] = y[*n1];
4125 p1[2] = z__[*n1];
4126 p2[0] = x[*n2];
4127 p2[1] = y[*n2];
4128 p2[2] = z__[*n2];
4129 al = arclen_(p1, p2);
4130 un[0] = p1[1] * p2[2] - p1[2] * p2[1];
4131 un[1] = p1[2] * p2[0] - p1[0] * p2[2];
4132 un[2] = p1[0] * p2[1] - p1[1] * p2[0];
4133 unorm = sqrt(un[0] * un[0] + un[1] * un[1] + un[2] * un[2]);
4134 if (unorm == 0. || al == 0.) {
4135 goto L11;
4136 }
4137
4138 /* Compute first difference S and scaled directional deriva- */
4139 /* tives S1,S2 at the endpoints (for the direction N1->N2). */
4140
4141 s = h__[*n2] - h__[*n1];
4142 s1 = al * (grad[*n1 * 3 + 1] * p2[0] + grad[*n1 * 3 + 2] * p2[1] + grad[*
4143 n1 * 3 + 3] * p2[2]) / unorm;
4144 s2 = -al * (grad[*n2 * 3 + 1] * p1[0] + grad[*n2 * 3 + 2] * p1[1] + grad[*
4145 n2 * 3 + 3] * p1[2]) / unorm;
4146
4147 /* Test for a valid constraint. */
4148
4149 *ier = -3;
4150 /* Computing MIN */
4151 d__1 = min(s1,s2);
4152 /* Computing MAX */
4153 d__2 = max(s1,s2);
4154 if ((rf < 0. && min(d__1,s) < bnd) || (rf > 0. && bnd < max(d__2,s))) {
4155 goto L11;
4156 }
4157
4158 /* Test for infinite tension required. */
4159
4160 *ier = 1;
4161 sig = sbig;
4162 if (s == bnd && (s1 != s || s2 != s)) {
4163 goto L10;
4164 }
4165
4166 /* Test for SIG = 0 sufficient. The Hermite cubic interpo- */
4167 /* lant H0 has derivative HP0(T) = (S2 + 2*B0*R + A0*R**2)/ */
4168 /* AL, where R = (T2-T)/AL. */
4169
4170 *ier = 0;
4171 sig = 0.;
4172 t0 = s * 3. - s1 - s2;
4173 b0 = t0 - s2;
4174 c0 = t0 - s1;
4175 a0 = -b0 - c0;
4176
4177 /* HP0(R) has an extremum (at R = -B0/A0) in (0,1) iff */
4178 /* B0*C0 > 0 and the third derivative of H0 has the */
4179 /* sign of A0. */
4180
4181 if (b0 * c0 <= 0. || a0 * rf > 0.) {
4182 goto L10;
4183 }
4184
4185 /* A0*RF < 0 and HP0(R) = -D0/(DT*A0) at R = -B0/A0. */
4186
4187 d0 = t0 * t0 - s1 * s2;
4188 f0 = (bnd + d0 / (a0 * al)) * rf;
4189 if (f0 >= 0.) {
4190 goto L10;
4191 }
4192
4193 /* Find a zero of F(SIG) = (BND-HP(R))*RF, where HP has an */
4194 /* extremum at R. F has a unique zero, F(0) = F0 < 0, and */
4195 /* F = (BND-S)*RF > 0 for SIG sufficiently large. */
4196
4197 /* Initialize parameters for the secant method. The method */
4198 /* uses three points: (SG0,F0), (SIG,F), and (SNEG,FNEG), */
4199 /* where SG0 and SNEG are defined implicitly by DSIG = SIG */
4200 /* - SG0 and DMAX = SIG - SNEG. SG0 is initially zero and */
4201 /* SIG is initialized to the zero of (BND - (SIG*S-S1-S2)/ */
4202 /* (AL*(SIG-2.)))*RF -- a value for which F(SIG) .GE. 0 and */
4203 /* F(SIG) = 0 for SIG sufficiently large that 2*SIG is in- */
4204 /* significant relative to exp(SIG). */
4205
4206 fmax = (bnd - s / al) * rf;
4207 sig = 2. - a0 / ((al * bnd - s) * 3.);
4208 /* IF (LUN .GE. 0) WRITE (LUN,120) F0, FMAX, SIG */
4209 /* 120 FORMAT (1X,9X,'F0 = ',E15.8,', FMAX = ',E15.8/ */
4210 /* . 1X,8X,'SIG = ',E15.8/) */
4211 if (dbg_verbose) fprintf (stderr, "F0 = %g FMAX = %g SIG = %g\n", f0, fmax, sig);
4212 d__1 = sig * exp(-sig) + .5;
4213 if (store_(&d__1) == .5) {
4214 goto L10;
4215 }
4216 dsig = sig;
4217 dmax__ = sig * -2.;
4218 fneg = fmax;
4219 d1 = s - s1;
4220 d2 = s2 - s;
4221 d1pd2 = d1 + d2;
4222 nit = 0;
4223
4224 /* Compute an absolute tolerance FTOL = abs(TOL), and a */
4225 /* relative tolerance RTOL = 100*Macheps. */
4226
4227 ftol = fabs(*tol);
4228 rtol = 1.;
4229 L5:
4230 rtol /= 2.;
4231 d__1 = rtol + 1.;
4232 if (store_(&d__1) > 1.) {
4233 goto L5;
4234 }
4235 rtol *= 200.;
4236
4237 /* Top of loop: compute F. */
4238
4239 L6:
4240 if (sig <= .5) {
4241
4242 /* Use approximations designed to avoid cancellation */
4243 /* error (associated with small SIG) in the modified */
4244 /* hyperbolic functions. */
4245
4246 snhcsh_(&sig, &sinhm, &coshm, &coshmm);
4247 c1 = sig * coshm * d2 - sinhm * d1pd2;
4248 c2 = sig * (sinhm + sig) * d2 - coshm * d1pd2;
4249 a = c2 - c1;
4250 e = sig * sinhm - coshmm - coshmm;
4251 } else {
4252
4253 /* Scale SINHM and COSHM by 2*exp(-SIG) in order to avoid */
4254 /* overflow. */
4255
4256 ems = exp(-sig);
4257 ems2 = ems + ems;
4258 tm = 1. - ems;
4259 sinh__ = tm * (ems + 1.);
4260 sinhm = sinh__ - sig * ems2;
4261 coshm = tm * tm;
4262 c1 = sig * coshm * d2 - sinhm * d1pd2;
4263 c2 = sig * sinh__ * d2 - coshm * d1pd2;
4264 a = ems2 * (sig * tm * d2 + (tm - sig) * d1pd2);
4265 e = sig * sinh__ - coshm - coshm;
4266 }
4267
4268 /* The second derivative HPP of H(R) has a zero at exp(SIG* */
4269 /* R) = SQRT((C2+C1)/A) and R is in (0,1) and well- */
4270 /* defined iff HPP(T1)*HPP(T2) < 0. */
4271
4272 f = fmax;
4273 t1 = a * (c2 + c1);
4274 if (t1 >= 0.) {
4275 if (c1 * (sig * coshm * d1 - sinhm * d1pd2) < 0.) {
4276
4277 /* HP(R) = (B+SIGN(A)*SQRT(A*C))/(AL*E) at the critical */
4278 /* value of R, where A = C2-C1, B = E*S2-C2, and C = C2 + */
4279 /* C1. Note that RF*A < 0. */
4280
4281 f = (bnd - (e * s2 - c2 - rf * sqrt(t1)) / (al * e)) * rf;
4282 }
4283 }
4284
4285 /* Update the number of iterations NIT. */
4286
4287 ++nit;
4288 /* IF (LUN .GE. 0) WRITE (LUN,130) NIT, SIG, F */
4289 /* 130 FORMAT (1X,3X,I2,' -- SIG = ',E15.8,', F = ', */
4290 /* . E15.8) */
4291 if (dbg_verbose) fprintf (stderr, "%d -- SIG = %g F = %g\n", nit, sig, f);
4292 if (f0 * f < 0.) {
4293
4294 /* F0*F < 0. Update (SNEG,FNEG) to (SG0,F0) so that F */
4295 /* and FNEG always have opposite signs. If SIG is closer */
4296 /* to SNEG than SG0 and abs(F) < abs(FNEG), then swap */
4297 /* (SNEG,FNEG) with (SG0,F0). */
4298
4299 t1 = dmax__;
4300 t2 = fneg;
4301 dmax__ = dsig;
4302 fneg = f0;
4303 if (fabs(dsig) > fabs(t1) && fabs(f) < fabs(t2)) {
4304
4305 dsig = t1;
4306 f0 = t2;
4307 }
4308 }
4309
4310 /* Test for convergence. */
4311
4312 stol = rtol * sig;
4313 if (fabs(dmax__) <= stol || (f >= 0. && f <= ftol) || fabs(f) <= rtol) {
4314 goto L10;
4315 }
4316 if (f0 * f < 0. || fabs(f) < fabs(f0)) {
4317 goto L8;
4318 }
4319
4320 /* F*F0 > 0 and the new estimate would be outside of the */
4321 /* bracketing interval of length abs(DMAX). Reset */
4322 /* (SG0,F0) to (SNEG,FNEG). */
4323
4324 L7:
4325 dsig = dmax__;
4326 f0 = fneg;
4327
4328 /* Compute the change in SIG by linear interpolation */
4329 /* between (SG0,F0) and (SIG,F). */
4330
4331 L8:
4332 dsig = -f * dsig / (f - f0);
4333 /* IF (LUN .GE. 0) WRITE (LUN,140) DSIG */
4334 /* 140 FORMAT (1X,8X,'DSIG = ',E15.8) */
4335 if (dbg_verbose) fprintf (stderr, "DSIG = %g\n", dsig);
4336 if (fabs(dsig) > fabs(dmax__) || dsig * dmax__ > 0.) {
4337 goto L7;
4338 }
4339
4340 /* Restrict the step-size such that abs(DSIG) .GE. STOL/2. */
4341 /* Note that DSIG and DMAX have opposite signs. */
4342
4343 if (fabs(dsig) < stol / 2.) {
4344 d__1 = stol / 2.;
4345 dsig = -d_sign(&d__1, &dmax__);
4346 }
4347
4348 /* Bottom of loop: update SIG, DMAX, and F0. */
4349
4350 sig += dsig;
4351 dmax__ += dsig;
4352 f0 = f;
4353 goto L6;
4354
4355 /* No errors encountered. */
4356
4357 L10:
4358 ret_val = sig;
4359 if (*iflgs <= 0) {
4360 return ret_val;
4361 }
4362 sigma[lp1] = sig;
4363 sigma[lp2] = sig;
4364 return ret_val;
4365
4366 /* Error termination. */
4367
4368 L11:
4369 ret_val = -1.;
4370 return ret_val;
4371 } /* sig1_ */
4372
sig2_(integer * n1,integer * n2,integer * n,doublereal * x,doublereal * y,doublereal * z__,doublereal * h__,integer * list,integer * lptr,integer * lend,doublereal * grad,doublereal * tol,integer * iflgs,doublereal * sigma,integer * ier)4373 doublereal sig2_(integer *n1, integer *n2, integer *n, doublereal *x,
4374 doublereal *y, doublereal *z__, doublereal *h__, integer *list,
4375 integer *lptr, integer *lend, doublereal *grad, doublereal *tol,
4376 integer *iflgs, doublereal *sigma, integer *ier)
4377 {
4378 /* Initialized data */
4379
4380 static doublereal sbig = 85.;
4381
4382 /* System generated locals */
4383 integer i__1;
4384 doublereal ret_val, d__1, d__2;
4385
4386 /* Builtin functions */
4387 double sqrt(doublereal), exp(doublereal);
4388
4389 /* Local variables */
4390 static doublereal f, s, t, d1, d2, p1[3], p2[3], t1, al, fp, un[3];
4391 static integer lp1, lp2;
4392 static doublereal tp1, d1d2, sig, ems;
4393 static integer lpl, nit;
4394 static doublereal ssm, dsig, ftol, rtol, coshm, sinhm, dummy;
4395 extern doublereal store_(doublereal *);
4396 static doublereal unorm;
4397 extern doublereal arclen_(doublereal *, doublereal *);
4398 extern /* Subroutine */ int snhcsh_(doublereal *, doublereal *,
4399 doublereal *, doublereal *);
4400
4401
4402 /* *********************************************************** */
4403
4404 /* From SSRFPACK */
4405 /* Robert J. Renka */
4406 /* Dept. of Computer Science */
4407 /* Univ. of North Texas */
4408 /* renka@cs.unt.edu */
4409 /* 07/25/96 */
4410
4411 /* Given a triangulation of a set of nodes on the unit */
4412 /* sphere, along with data values H and gradients GRAD at the */
4413 /* nodes, this function determines the smallest tension fac- */
4414 /* tor SIG2 such that the Hermite interpolatory tension */
4415 /* spline H(A), defined by SIG2 and the endpoint values and */
4416 /* directional derivatives associated with an arc N1-N2, */
4417 /* preserves convexity (or concavity) of the data: */
4418
4419 /* HP1 .LE. S .LE. HP2 implies HPP(A) .GE. 0, and */
4420 /* HP1 .GE. S .GE. HP2 implies HPP(A) .LE. 0 */
4421
4422 /* for all A in the open interval (A1,A2) corresponding to */
4423 /* the arc, where HP1 and HP2 are the derivative values of H */
4424 /* at the endpoints, S is the slope of the linear interpolant */
4425 /* of the endpoint data values, HPP denotes the second deriv- */
4426 /* ative of H, and A is arc-length. Note, however, that */
4427 /* infinite tension is required if HP1 = S or HP2 = S (unless */
4428 /* HP1 = HP2 = S). */
4429
4430 /* On input: */
4431
4432 /* N1,N2 = Nodal indexes of the endpoints of an arc for */
4433 /* which the tension factor is to be computed. */
4434 /* The indexes must be distinct and lie in the */
4435 /* range 1 to N, and if IFLGS .GE. 1, they must */
4436 /* correspond to adjacent nodes in the triangu- */
4437 /* lation. */
4438
4439 /* N = Number of nodes in the triangulation. N .GE. 3. */
4440
4441 /* X,Y,Z = Arrays of length N containing coordinates of */
4442 /* the nodes. X(I)**2 + Y(I)**2 + Z(I)**2 = 1. */
4443
4444 /* H = Array of length N containing data values at the */
4445 /* nodes. H(I) is associated with (X(I),Y(I),Z(I)) */
4446 /* for I = 1 to N. */
4447
4448 /* LIST,LPTR,LEND = Data structure defining the trian- */
4449 /* gulation. Refer to STRIPACK */
4450 /* Subroutine TRMESH. */
4451
4452 /* GRAD = Array dimensioned 3 by N whose columns con- */
4453 /* gradients at the nodes. GRAD( ,J) must be */
4454 /* orthogonal to node J: GRAD(1,J)*X(J) + */
4455 /* GRAD(2,J)*Y(J) + GRAD(3,J)*Z(J) = 0. Refer */
4456 /* to Subroutines GRADG, GRADL, and SMSURF. */
4457
4458 /* TOL = Tolerance whose magnitude determines how close */
4459 /* SIG2 is to its optimal value when nonzero */
4460 /* finite tension is necessary and sufficient to */
4461 /* satisfy convexity or concavity. In the case */
4462 /* convexity, SIG2 is chosen so that 0 .LE. */
4463 /* HPPMIN .LE. abs(TOL), where HPPMIN is the */
4464 /* minimum value of HPP on the arc. In the case */
4465 /* of concavity, the maximum value of HPP satis- */
4466 /* fies -abs(TOL) .LE. HPPMAX .LE. 0. Thus, the */
4467 /* constraint is satisfied but possibly with more */
4468 /* tension than necessary. */
4469
4470 /* IFLGS = Tension array option indicator: */
4471 /* IFLGS .LE. 0 if SIGMA is not to be used. */
4472 /* IFLGS .GE. 1 if SIGMA is to be updated by */
4473 /* storing SIG2 in the appropriate */
4474 /* locations. */
4475
4476 /* The above parameters are not altered by this function. */
4477
4478 /* SIGMA = Dummy parameter (IFLGS .LE. 0) or array */
4479 /* containing tension factors associated with */
4480 /* arcs in one-to-one correspondence with LIST */
4481 /* entries (IFLGS .GE. 1). Refer to Subroutine */
4482 /* GETSIG. */
4483
4484 /* On output: */
4485
4486 /* SIGMA = Tension factor array updated with the new */
4487 /* value if and only if IFLGS .GE. 1 and IER */
4488 /* .GE. 0. */
4489
4490 /* IER = Error indicator: */
4491 /* IER = 0 if no errors were encountered and fin- */
4492 /* ite tension is sufficient to satisfy */
4493 /* convexity (or concavity). */
4494 /* IER = 1 if no errors were encountered but in- */
4495 /* finite tension is required to satisfy */
4496 /* convexity. */
4497 /* IER = 2 if the data does not satisfy convexity */
4498 /* or concavity. */
4499 /* IER = -1 if N1, N2, or N is outside its valid */
4500 /* range. */
4501 /* IER = -2 if nodes N1 and N2 coincide or IFLGS */
4502 /* .GE. 1 and the nodes are not adja- */
4503 /* cent. */
4504
4505 /* SIG2 = Minimum tension factor defined above unless */
4506 /* IER < 0, in which case SIG2 = -1. If IER */
4507 /* = 1, SIG2 is set to 85, resulting in an */
4508 /* approximation to the linear interpolant of */
4509 /* the endpoint values. If IER = 2, SIG2 = 0, */
4510 /* resulting in the Hermite cubic interpolant. */
4511
4512 /* STRIPACK module required by SIG2: STORE */
4513
4514 /* SSRFPACK modules required by SIG2: ARCLEN, SNHCSH */
4515
4516 /* Intrinsic functions called by SIG2: ABS, EXP, MAX, MIN, */
4517 /* SQRT */
4518
4519 /* *********************************************************** */
4520
4521
4522 /* Parameter adjustments */
4523 grad -= 4;
4524 --lend;
4525 --h__;
4526 --z__;
4527 --y;
4528 --x;
4529 --list;
4530 --lptr;
4531 --sigma;
4532
4533 /* Function Body */
4534
4535 /* Print a heading. */
4536
4537 /* IF (LUN .GE. 0) WRITE (LUN,100) N1, N2 */
4538 /* 100 FORMAT (//1X,'SIG2 -- N1 =',I4,', N2 =',I4) */
4539 if (dbg_verbose) fprintf (stderr, "SIG2 -- N1 = %d N2 = %d\n", *n1, *n2);
4540
4541 /* Test for errors and set local parameters. */
4542
4543 *ier = -1;
4544 /* Computing MAX */
4545 i__1 = max(*n1,*n2);
4546 if (min(*n1,*n2) < 1 || *n1 == *n2 || max(i__1,3) > *n) {
4547 goto L11;
4548 }
4549 *ier = -2;
4550 if (*iflgs > 0) {
4551
4552 /* Set LP1 and LP2 to the pointers to N2 as a neighbor of */
4553 /* N1 and N1 as a neighbor of N2, respectively. */
4554
4555 lpl = lend[*n1];
4556 lp1 = lptr[lpl];
4557 L1:
4558 if (list[lp1] == *n2) {
4559 goto L2;
4560 }
4561 lp1 = lptr[lp1];
4562 if (lp1 != lpl) {
4563 goto L1;
4564 }
4565 if ((i__1 = list[lp1], abs(i__1)) != *n2) {
4566 goto L11;
4567 }
4568
4569 L2:
4570 lpl = lend[*n2];
4571 lp2 = lptr[lpl];
4572 L3:
4573 if (list[lp2] == *n1) {
4574 goto L4;
4575 }
4576 lp2 = lptr[lp2];
4577 if (lp2 != lpl) {
4578 goto L3;
4579 }
4580 if ((i__1 = list[lp2], abs(i__1)) != *n1) {
4581 goto L11;
4582 }
4583 }
4584
4585 /* Store nodal coordinates P1 and P2, compute arc-length AL */
4586 /* and unit normal UN = (P1 X P2)/UNORM, and test for */
4587 /* coincident nodes. */
4588
4589 L4:
4590 p1[0] = x[*n1];
4591 p1[1] = y[*n1];
4592 p1[2] = z__[*n1];
4593 p2[0] = x[*n2];
4594 p2[1] = y[*n2];
4595 p2[2] = z__[*n2];
4596 al = arclen_(p1, p2);
4597 un[0] = p1[1] * p2[2] - p1[2] * p2[1];
4598 un[1] = p1[2] * p2[0] - p1[0] * p2[2];
4599 un[2] = p1[0] * p2[1] - p1[1] * p2[0];
4600 unorm = sqrt(un[0] * un[0] + un[1] * un[1] + un[2] * un[2]);
4601 if (unorm == 0. || al == 0.) {
4602 goto L11;
4603 }
4604
4605 /* Compute first and second differences and test for infinite */
4606 /* tension required. */
4607
4608 s = h__[*n2] - h__[*n1];
4609 d1 = s - al * (grad[*n1 * 3 + 1] * p2[0] + grad[*n1 * 3 + 2] * p2[1] +
4610 grad[*n1 * 3 + 3] * p2[2]) / unorm;
4611 d2 = -al * (grad[*n2 * 3 + 1] * p1[0] + grad[*n2 * 3 + 2] * p1[1] + grad[*
4612 n2 * 3 + 3] * p1[2]) / unorm - s;
4613 d1d2 = d1 * d2;
4614 *ier = 1;
4615 sig = sbig;
4616 if (d1d2 == 0. && d1 != d2) {
4617 goto L10;
4618 }
4619
4620 /* Test for a valid constraint. */
4621
4622 *ier = 2;
4623 sig = 0.;
4624 if (d1d2 < 0.) {
4625 goto L10;
4626 }
4627
4628 /* Test for SIG = 0 sufficient. */
4629
4630 *ier = 0;
4631 if (d1d2 == 0.) {
4632 goto L10;
4633 }
4634 /* Computing MAX */
4635 d__1 = d1 / d2, d__2 = d2 / d1;
4636 t = max(d__1,d__2);
4637 if (t <= 2.) {
4638 goto L10;
4639 }
4640
4641 /* Find a zero of F(SIG) = SIG*COSHM(SIG)/SINHM(SIG) - (T+1). */
4642 /* Since the derivative of F vanishes at the origin, a */
4643 /* quadratic approximation is used to obtain an initial */
4644 /* estimate for the Newton method. */
4645
4646 tp1 = t + 1.;
4647 sig = sqrt(t * 10. - 20.);
4648 nit = 0;
4649
4650 /* Compute an absolute tolerance FTOL = abs(TOL) and a */
4651 /* relative tolerance RTOL = 100*Macheps. */
4652
4653 ftol = fabs(*tol);
4654 rtol = 1.;
4655 L5:
4656 rtol /= 2.;
4657 d__1 = rtol + 1.;
4658 if (store_(&d__1) > 1.) {
4659 goto L5;
4660 }
4661 rtol *= 200.;
4662
4663 /* Top of loop: evaluate F and its derivative FP. */
4664
4665 L6:
4666 if (sig <= .5) {
4667
4668 /* Use approximations designed to avoid cancellation error */
4669 /* in the hyperbolic functions. */
4670
4671 snhcsh_(&sig, &sinhm, &coshm, &dummy);
4672 t1 = coshm / sinhm;
4673 fp = t1 + sig * (sig / sinhm - t1 * t1 + 1.);
4674 } else {
4675
4676 /* Scale SINHM and COSHM by 2*exp(-SIG) in order to avoid */
4677 /* overflow. */
4678
4679 ems = exp(-sig);
4680 ssm = 1. - ems * (ems + sig + sig);
4681 t1 = (1. - ems) * (1. - ems) / ssm;
4682 fp = t1 + sig * (sig * 2. * ems / ssm - t1 * t1 + 1.);
4683 }
4684
4685 f = sig * t1 - tp1;
4686
4687 /* Update the number of iterations NIT. */
4688
4689 ++nit;
4690 /* IF (LUN .GE. 0) WRITE (LUN,110) NIT, SIG, F, FP */
4691 /* 110 FORMAT (1X,3X,I2,' -- SIG = ',E15.8,', F = ', */
4692 /* . E15.8/1X,31X,'FP = ',E15.8) */
4693 if (dbg_verbose) fprintf (stderr, "%d -- SIG = %g F = %g FP = %g\n", nit, sig, f, fp);
4694
4695 /* Test for convergence. */
4696
4697 if (fp <= 0.) {
4698 goto L10;
4699 }
4700 dsig = -f / fp;
4701 if (fabs(dsig) <= rtol * sig || (f >= 0. && f <= ftol) || fabs(f) <= rtol) {
4702 goto L10;
4703 }
4704
4705 /* Bottom of loop: update SIG. */
4706
4707 sig += dsig;
4708 goto L6;
4709
4710 /* No errors encountered. */
4711
4712 L10:
4713 ret_val = sig;
4714 if (*iflgs <= 0) {
4715 return ret_val;
4716 }
4717 sigma[lp1] = sig;
4718 sigma[lp2] = sig;
4719 return ret_val;
4720
4721 /* Error termination. */
4722
4723 L11:
4724 ret_val = -1.;
4725 return ret_val;
4726 } /* sig2_ */
4727
smsgs_(integer * n,doublereal * x,doublereal * y,doublereal * z__,doublereal * u,integer * list,integer * lptr,integer * lend,integer * iflgs,doublereal * sigma,doublereal * w,doublereal * p,integer * nit,doublereal * dfmax,doublereal * f,doublereal * grad,integer * ier)4728 /* Subroutine */ int smsgs_(integer *n, doublereal *x, doublereal *y,
4729 doublereal *z__, doublereal *u, integer *list, integer *lptr, integer
4730 *lend, integer *iflgs, doublereal *sigma, doublereal *w, doublereal *
4731 p, integer *nit, doublereal *dfmax, doublereal *f, doublereal *grad,
4732 integer *ier)
4733 {
4734 /* System generated locals */
4735 integer i__1, i__2;
4736 doublereal d__1, d__2;
4737
4738 /* Builtin functions */
4739 double sqrt(doublereal), atan(doublereal);
4740
4741 /* Local variables */
4742 static integer j, k;
4743 static doublereal t, g1, g2, g3, r1, r2, r3, t1, t2, t3, t4, t5, t6, c11,
4744 c12, c13, c22, c23, c33, df, fk, cx;
4745 static integer nn;
4746 static doublereal cy, pp, xj, xk, yj, yk, zj, zk, sx, sy, xs, ys, rr2,
4747 rr3, cc22, cc23, cc33, dgk[3];
4748 static integer ifl;
4749 static doublereal gjk, det, gkj, dgx, dgy, sig;
4750 static integer lpj, lpl;
4751 static doublereal tol, den1, den2, alfa, dfmx;
4752 static integer iter;
4753 static doublereal alfsq, sinal;
4754 static integer itmax;
4755 extern /* Subroutine */ int grcoef_(doublereal *, doublereal *,
4756 doublereal *), constr_(doublereal *, doublereal *, doublereal *,
4757 doublereal *, doublereal *, doublereal *, doublereal *), aplyrt_(
4758 doublereal *, doublereal *, doublereal *, doublereal *,
4759 doublereal *, doublereal *, doublereal *);
4760
4761
4762 /* *********************************************************** */
4763
4764 /* From SSRFPACK */
4765 /* Robert J. Renka */
4766 /* Dept. of Computer Science */
4767 /* Univ. of North Texas */
4768 /* renka@cs.unt.edu */
4769 /* 07/25/96 */
4770
4771 /* This subroutine solves the symmetric positive definite */
4772 /* linear system associated with minimizing the quadratic */
4773 /* functional Q(F,FX,FY,FZ) described in Subroutine SMSURF. */
4774 /* Since the gradient at node K lies in the plane tangent to */
4775 /* the sphere surface at K, it is effectively defined by only */
4776 /* two components -- its X and Y components in the coordinate */
4777 /* system obtained by rotating K to the north pole. Thus, */
4778 /* the minimization problem corresponds to an order-3N system */
4779 /* which is solved by the block Gauss-Seidel method with 3 by */
4780 /* 3 blocks. */
4781
4782 /* On input: */
4783
4784 /* N,X,Y,Z,U,LIST,LPTR,LEND,IFLGS,SIGMA,W = Parameters */
4785 /* as described in Subroutine SMSURF. */
4786
4787 /* P = Positive smoothing parameter defining Q. */
4788
4789 /* The above parameters are not altered by this routine. */
4790
4791 /* NIT = Maximum number of iterations to be used. This */
4792 /* maximum will likely be achieved if DFMAX is */
4793 /* smaller than the machine precision. NIT .GE. */
4794 /* 0. */
4795
4796 /* DFMAX = Nonnegative convergence criterion. The */
4797 /* method is terminated when the maximum */
4798 /* change in a solution F-component between */
4799 /* iterations is at most DFMAX. The change in */
4800 /* a component is taken to be the absolute */
4801 /* difference relative to 1 plus the old value. */
4802
4803 /* F = Initial estimate of the first N solution compo- */
4804 /* nents. */
4805
4806 /* GRAD = 3 by N array containing initial estimates of */
4807 /* the last 3N solution components (the gradi- */
4808 /* ent with FX, FY, and FZ in rows 1, 2, and 3, */
4809 /* respectively). */
4810
4811 /* On output: */
4812
4813 /* NIT = Number of Gauss-Seidel iterations employed. */
4814
4815 /* DFMAX = Maximum relative change in a solution F- */
4816 /* component at the last iteration. */
4817
4818 /* F = First N solution components -- function values */
4819 /* at the nodes. */
4820
4821 /* GRAD = Last 3N solution components -- gradients at */
4822 /* the nodes. */
4823
4824 /* IER = Error indicator: */
4825 /* IER = 0 if no errors were encountered and the */
4826 /* convergence criterion was achieved. */
4827 /* IER = 1 if no errors were encountered but con- */
4828 /* vergence was not achieved within NIT */
4829 /* iterations. */
4830 /* IER = -1 if N, P, NIT, or DFMAX is outside its */
4831 /* valid range on input. F and GRAD are */
4832 /* not altered in this case. */
4833 /* IER = -2 if all nodes are collinear or the */
4834 /* triangulation is invalid. */
4835 /* IER = -3 if duplicate nodes were encountered. */
4836
4837 /* SSRFPACK modules required by SMSGS: APLYRT, CONSTR, */
4838 /* GRCOEF, SNHCSH */
4839
4840 /* Intrinsic functions called by SMSGS: ABS, ATAN, MAX, SQRT */
4841
4842 /* *********************************************************** */
4843
4844
4845 /* Parameter adjustments */
4846 grad -= 4;
4847 --f;
4848 --w;
4849 --lend;
4850 --u;
4851 --z__;
4852 --y;
4853 --x;
4854 --list;
4855 --lptr;
4856 --sigma;
4857
4858 /* Function Body */
4859 nn = *n;
4860 ifl = *iflgs;
4861 pp = *p;
4862 itmax = *nit;
4863 tol = *dfmax;
4864
4865 /* Test for errors in input and initialize iteration count, */
4866 /* tension factor, and output value of DFMAX. */
4867
4868 if (nn < 3 || pp <= 0. || itmax < 0 || tol < 0.) {
4869 goto L5;
4870 }
4871 iter = 0;
4872 sig = sigma[1];
4873 dfmx = 0.;
4874
4875 /* Top of iteration loop. */
4876
4877 L1:
4878 if (iter == itmax) {
4879 goto L4;
4880 }
4881 dfmx = 0.;
4882
4883 /* Loop on nodes. */
4884
4885 i__1 = nn;
4886 for (k = 1; k <= i__1; ++k) {
4887 xk = x[k];
4888 yk = y[k];
4889 zk = z__[k];
4890 fk = f[k];
4891 g1 = grad[k * 3 + 1];
4892 g2 = grad[k * 3 + 2];
4893 g3 = grad[k * 3 + 3];
4894
4895 /* Construct the rotation mapping node K to the north pole. */
4896
4897 constr_(&xk, &yk, &zk, &cx, &sx, &cy, &sy);
4898
4899 /* Initialize components of the order-3 system for the */
4900 /* change (DF,DGX,DGY) in the K-th solution components. */
4901
4902 c11 = pp * w[k];
4903 c12 = 0.;
4904 c13 = 0.;
4905 c22 = 0.;
4906 c23 = 0.;
4907 c33 = 0.;
4908 r1 = c11 * (u[k] - fk);
4909 r2 = 0.;
4910 r3 = 0.;
4911
4912 /* Loop on neighbors J of node K. */
4913
4914 lpl = lend[k];
4915 lpj = lpl;
4916 L2:
4917 lpj = lptr[lpj];
4918 j = (i__2 = list[lpj], abs(i__2));
4919
4920 /* Compute the coordinates of J in the rotated system. */
4921
4922 t = sx * y[j] + cx * z__[j];
4923 yj = cx * y[j] - sx * z__[j];
4924 zj = sy * x[j] + cy * t;
4925 xj = cy * x[j] - sy * t;
4926
4927 /* Compute arc-length ALFA between K and J, ALFSQ = ALFA* */
4928 /* ALFA, SINAL = SIN(ALFA), DEN1 = ALFA*SIN(ALFA)**2, and */
4929 /* DEN2 = ALFSQ*SINAL. */
4930
4931 alfa = atan(sqrt((1. - zj) / (zj + 1.))) * 2.;
4932 alfsq = alfa * alfa;
4933 xs = xj * xj;
4934 ys = yj * yj;
4935 sinal = sqrt(xs + ys);
4936 den1 = alfa * (xs + ys);
4937 den2 = alfsq * sinal;
4938
4939 /* Test for coincident nodes and compute functions of SIG: */
4940 /* T1 = SIG*SIG*COSHM/E, T2 = SIG*SINHM/E, and T3 = SIG* */
4941 /* (SIG*COSHM-SINHM)/E for E = SIG*SINH - 2*COSHM. */
4942
4943 if (den1 == 0.) {
4944 goto L7;
4945 }
4946 if (ifl >= 1) {
4947 sig = sigma[lpj];
4948 }
4949 grcoef_(&sig, &t3, &t2);
4950 t1 = t2 + t3;
4951
4952 /* Update system components for node J. */
4953
4954 t4 = t1 * 2. / (alfa * alfsq);
4955 t5 = t1 / den2;
4956 t6 = t3 / den1;
4957 c11 += t4;
4958 c12 += t5 * xj;
4959 c13 += t5 * yj;
4960 c22 += t6 * xs;
4961 c23 += t6 * xj * yj;
4962 c33 += t6 * ys;
4963 gkj = g1 * x[j] + g2 * y[j] + g3 * z__[j];
4964 gjk = grad[j * 3 + 1] * xk + grad[j * 3 + 2] * yk + grad[j * 3 + 3] *
4965 zk;
4966 r1 = r1 + t4 * (f[j] - fk) + t5 * (gjk - gkj);
4967 t = t5 * (f[j] - fk) - t6 * gkj + t2 * gjk / den1;
4968 r2 += t * xj;
4969 r3 += t * yj;
4970
4971 /* Bottom of loop on neighbors. */
4972
4973 if (lpj != lpl) {
4974 goto L2;
4975 }
4976
4977 /* Solve the system associated with the K-th block. */
4978
4979 cc22 = c11 * c22 - c12 * c12;
4980 cc23 = c11 * c23 - c12 * c13;
4981 cc33 = c11 * c33 - c13 * c13;
4982 rr2 = c11 * r2 - c12 * r1;
4983 rr3 = c11 * r3 - c13 * r1;
4984 det = cc22 * cc33 - cc23 * cc23;
4985 if (det == 0. || cc22 == 0. || c11 == 0.) {
4986 goto L6;
4987 }
4988 dgy = (cc22 * rr3 - cc23 * rr2) / det;
4989 dgx = (rr2 - cc23 * dgy) / cc22;
4990 df = (r1 - c12 * dgx - c13 * dgy) / c11;
4991
4992 /* Rotate (DGX,DGY,0) back to the original coordinate */
4993 /* system, and update GRAD( ,K), F(K), and DFMX. */
4994
4995 aplyrt_(&dgx, &dgy, &cx, &sx, &cy, &sy, dgk);
4996 grad[k * 3 + 1] = g1 + dgk[0];
4997 grad[k * 3 + 2] = g2 + dgk[1];
4998 grad[k * 3 + 3] = g3 + dgk[2];
4999 f[k] = fk + df;
5000 /* Computing MAX */
5001 d__1 = dfmx, d__2 = fabs(df) / (fabs(fk) + 1.);
5002 dfmx = max(d__1,d__2);
5003 /* L3: */
5004 }
5005
5006 /* Increment ITER and test for convergence. */
5007
5008 ++iter;
5009 if (dfmx > tol) {
5010 goto L1;
5011 }
5012
5013 /* The method converged. */
5014
5015 *nit = iter;
5016 *dfmax = dfmx;
5017 *ier = 0;
5018 return 0;
5019
5020 /* The method failed to converge within NIT iterations. */
5021
5022 L4:
5023 *dfmax = dfmx;
5024 *ier = 1;
5025 return 0;
5026
5027 /* Invalid input parameter. */
5028
5029 L5:
5030 *nit = 0;
5031 *dfmax = 0.;
5032 *ier = -1;
5033 return 0;
5034
5035 /* Node K and its neighbors are collinear. */
5036
5037 L6:
5038 *nit = 0;
5039 *dfmax = dfmx;
5040 *ier = -2;
5041 return 0;
5042
5043 /* Nodes J and K coincide. */
5044
5045 L7:
5046 *nit = 0;
5047 *dfmax = dfmx;
5048 *ier = -3;
5049 return 0;
5050 } /* smsgs_ */
5051
smsurf_(integer * n,doublereal * x,doublereal * y,doublereal * z__,doublereal * u,integer * list,integer * lptr,integer * lend,integer * iflgs,doublereal * sigma,doublereal * w,doublereal * sm,doublereal * smtol,doublereal * gstol,integer * lprnt,doublereal * f,doublereal * grad,integer * ier)5052 /* Subroutine */ int smsurf_(integer *n, doublereal *x, doublereal *y,
5053 doublereal *z__, doublereal *u, integer *list, integer *lptr, integer
5054 *lend, integer *iflgs, doublereal *sigma, doublereal *w, doublereal *
5055 sm, doublereal *smtol, doublereal *gstol, integer *lprnt, doublereal *
5056 f, doublereal *grad, integer *ier)
5057 {
5058 /* Initialized data */
5059
5060 static integer itmax = 50;
5061 static integer nitmax = 40;
5062
5063 /* System generated locals */
5064 integer i__1;
5065 doublereal d__1;
5066
5067 /* Builtin functions */
5068 double sqrt(doublereal);
5069
5070 /* Local variables */
5071 static doublereal c__, g;
5072 static integer i__;
5073 static doublereal p, s, g0, q2, dp;
5074 static integer nn;
5075 static doublereal wi;
5076 static integer nit, lun;
5077 static doublereal tol, gneg, dmax__;
5078 static integer ierr, iter;
5079 static doublereal sumw, q2min, q2max, dfmax;
5080 extern /* Subroutine */ int smsgs_(integer *, doublereal *, doublereal *,
5081 doublereal *, doublereal *, integer *, integer *, integer *,
5082 integer *, doublereal *, doublereal *, doublereal *, integer *,
5083 doublereal *, doublereal *, doublereal *, integer *);
5084
5085
5086 /* *********************************************************** */
5087
5088 /* From SSRFPACK */
5089 /* Robert J. Renka */
5090 /* Dept. of Computer Science */
5091 /* Univ. of North Texas */
5092 /* renka@cs.unt.edu */
5093 /* 07/21/98 */
5094
5095 /* Given a triangulation of N nodes on the unit sphere with */
5096 /* data values U at the nodes and tension factors SIGMA */
5097 /* associated with the arcs, this routine determines a set of */
5098 /* nodal function values F and gradients GRAD = (FX,FY,FZ) */
5099 /* such that a quadratic functional Q1(F,GRAD) is minimized */
5100 /* subject to the constraint Q2(F) .LE. SM for Q2(F) = */
5101 /* (U-F)**T*W*(U-F), where W is a diagonal matrix of positive */
5102 /* weights. The functional Q1 is an approximation to the */
5103 /* linearized curvature over the triangulation of a C-1 fun- */
5104 /* ction F(V), V a unit vector, which interpolates the nodal */
5105 /* values and gradients. Subroutines INTRC1 and UNIF may be */
5106 /* called to evaluate F at arbitrary points. */
5107
5108 /* The smoothing procedure is an extension of the method */
5109 /* for cubic spline smoothing due to C. Reinsch -- Numer. */
5110 /* Math., 10 (1967) and 16 (1971). Refer to Function FVAL */
5111 /* for a further description of the interpolant F. Letting */
5112 /* D1F(T) and D2F(T) denote first and second derivatives of F */
5113 /* with respect to a parameter T varying along a triangula- */
5114 /* tion arc, Q1 is the sum of integrals over the arcs of */
5115 /* D2F(T)**2 + ((SIGMA/L)*(D1F(T)-S))**2 where L denotes arc- */
5116 /* length, SIGMA is the appropriate tension factor, and S is */
5117 /* the slope of the linear function of T which interpolates */
5118 /* the values of F at the endpoints of the arc. Introducing */
5119 /* a smoothing parameter P, and assuming the constraint is */
5120 /* active, the problem is equivalent to minimizing Q(P,F, */
5121 /* GRAD) = Q1(F,GRAD) + P*(Q2(F)-SM). The secant method is */
5122 /* used to find a zero of G(P) = 1/SQRT(Q2) - 1/SQRT(SM) */
5123 /* where F(P) satisfies the order-3N symmetric positive def- */
5124 /* inite linear system obtained by setting the gradient of Q */
5125 /* (treated as a function of F and GRAD with GRAD tangent to */
5126 /* the sphere surface) to zero. The linear system is solved */
5127 /* by the block Gauss-Seidel method (refer to SMSGS). */
5128
5129 /* Note that the method can also be used to select grad- */
5130 /* ients for the interpolation problem (F = U, SM = 0, and P */
5131 /* infinite). This is achieved by a call to Subroutine */
5132 /* GRADG. */
5133
5134 /* On input: */
5135
5136 /* N = Number of nodes in the triangulation. N .GE. 3. */
5137
5138 /* X,Y,Z = Arrays of length N containing Cartesian */
5139 /* coordinates of the nodes. */
5140
5141 /* U = Array of length N containing data values at the */
5142 /* nodes. */
5143
5144 /* LIST,LPTR,LEND = Data structure defining the trian- */
5145 /* gulation. Refer to STRIPACK */
5146 /* Subroutine TRMESH. */
5147
5148 /* IFLGS = Tension factor option: */
5149 /* IFLGS .LE. 0 if a single uniform tension */
5150 /* factor is to be used. */
5151 /* IFLGS .GE. 1 if variable tension is desired. */
5152
5153 /* SIGMA = Uniform tension factor (IFLGS .LE. 0), or */
5154 /* array containing tension factors associated */
5155 /* with arcs in one-to-one correspondence with */
5156 /* LIST entries (IFLGS .GE. 1). Refer to Sub- */
5157 /* programs GETSIG, SIG0, SIG1, and SIG2. */
5158
5159 /* W = Array of length N containing positive weights */
5160 /* associated with the data values. The recommend- */
5161 /* ed value of W(I) is 1/DU**2 where DU is the */
5162 /* standard deviation associated with U(I). DU**2 */
5163 /* is the expected value of the squared error in */
5164 /* the measurement of U(I). (The mean error is */
5165 /* assumed to be zero.) */
5166
5167 /* SM = Positive parameter specifying an upper bound on */
5168 /* Q2(F). Note that F is constant (and Q2(F) */
5169 /* is minimized) if SM is sufficiently large that */
5170 /* the constraint is not active. It is recommend- */
5171 /* ed that SM satisfy N-SQRT(2N) .LE. SM .LE. N+ */
5172 /* SQRT(2N). */
5173
5174 /* SMTOL = Parameter in the open interval (0,1) speci- */
5175 /* fying the relative error allowed in satisfy- */
5176 /* ing the constraint -- the constraint is */
5177 /* assumed to be satisfied if SM*(1-SMTOL) .LE. */
5178 /* Q2 .LE. SM*(1+SMTOL). A reasonable value */
5179 /* for SMTOL is SQRT(2/N). */
5180
5181 /* GSTOL = Nonnegative tolerance defining the conver- */
5182 /* gence criterion for the Gauss-Seidel method. */
5183 /* Refer to parameter DFMAX in Subroutine */
5184 /* SMSGS. A recommended value is .05*DU**2, */
5185 /* where DU is an average standard deviation */
5186 /* in the data values. */
5187
5188 /* LPRNT = Logical unit on which diagnostic messages */
5189 /* are printed, or negative integer specifying */
5190 /* no diagnostics. For each secant iteration, */
5191 /* the following values are printed: P, G(P), */
5192 /* NIT, DFMAX, and DP, where NIT denotes the */
5193 /* number of Gauss-Seidel iterations used in */
5194 /* the computation of G, DFMAX denotes the max- */
5195 /* imum relative change in a solution component */
5196 /* in the last Gauss-Seidel iteration, and DP */
5197 /* is the change in P computed by linear inter- */
5198 /* polation between the current point (P,G) and */
5199 /* a previous point. */
5200
5201 /* Input parameters are not altered by this routine. */
5202
5203 /* On output: */
5204
5205 /* F = Array of length N containing nodal function val- */
5206 /* ues unless IER < 0. */
5207
5208 /* GRAD = 3 by N array whose columns contain gradients */
5209 /* of F at the nodes unless IER < 0. */
5210
5211 /* IER = Error indicator and information flag: */
5212 /* IER = 0 if no errors were encountered and the */
5213 /* constraint is active -- Q2(F) is ap- */
5214 /* proximately equal to SM. */
5215 /* IER = 1 if no errors were encountered but the */
5216 /* constraint is not active -- F and GRAD */
5217 /* are the values and gradients of a con- */
5218 /* stant function which minimizes Q2(F), */
5219 /* and Q1 = 0. */
5220 /* IER = 2 if the constraint could not be satis- */
5221 /* fied to within SMTOL due to */
5222 /* ill-conditioned linear systems. */
5223 /* IER = -1 if N, W, SM, SMTOL, or GSTOL is out- */
5224 /* side its valid range on input. */
5225 /* IER = -2 if all nodes are collinear or the */
5226 /* triangulation is invalid. */
5227 /* IER = -3 if duplicate nodes were encountered. */
5228
5229 /* SSRFPACK modules required by SMSURF: APLYRT, CONSTR, */
5230 /* GRCOEF, SMSGS, */
5231 /* SNHCSH */
5232
5233 /* Intrinsic functions called by SMSURF: ABS, SQRT */
5234
5235 /* *********************************************************** */
5236
5237
5238 /* Local parameters: */
5239
5240 /* ITMAX = Maximum number of secant iterations. */
5241 /* LUN = Local copy of LPRNT. */
5242 /* NITMAX = Maximum number of Gauss-Seidel iterations for */
5243 /* each secant iteration. */
5244 /* NN = Local copy of N. */
5245 /* TOL = Local copy of GSTOL. */
5246
5247 /* Parameter adjustments */
5248 grad -= 4;
5249 --f;
5250 --w;
5251 --lend;
5252 --u;
5253 --z__;
5254 --y;
5255 --x;
5256 --list;
5257 --lptr;
5258 --sigma;
5259
5260 /* Function Body */
5261
5262 nn = *n;
5263 tol = *gstol;
5264 lun = *lprnt;
5265 if (lun > 99) {
5266 lun = -1;
5267 }
5268
5269 /* Test for errors and initialize F to the weighted least */
5270 /* squares fit of a constant function to the data. */
5271
5272 *ier = -1;
5273 if (nn < 3 || *sm <= 0. || *smtol <= 0. || *smtol >= 1. || tol <= 0.) {
5274 return 0;
5275 }
5276 c__ = 0.;
5277 sumw = 0.;
5278 i__1 = nn;
5279 for (i__ = 1; i__ <= i__1; ++i__) {
5280 wi = w[i__];
5281 if (wi <= 0.) {
5282 return 0;
5283 }
5284 c__ += wi * u[i__];
5285 sumw += wi;
5286 /* L1: */
5287 }
5288 c__ /= sumw;
5289
5290 /* Compute nodal values and gradients, and accumulate Q2 = */
5291 /* (U-F)**T*W*(U-F). */
5292
5293 q2 = 0.;
5294 i__1 = nn;
5295 for (i__ = 1; i__ <= i__1; ++i__) {
5296 f[i__] = c__;
5297 grad[i__ * 3 + 1] = 0.;
5298 grad[i__ * 3 + 2] = 0.;
5299 grad[i__ * 3 + 3] = 0.;
5300 /* Computing 2nd power */
5301 d__1 = u[i__] - f[i__];
5302 q2 += w[i__] * (d__1 * d__1);
5303 /* L2: */
5304 }
5305
5306 /* Compute bounds on Q2 defined by SMTOL, and test for the */
5307 /* constraint satisfied by the constant fit. */
5308
5309 q2min = *sm * (1. - *smtol);
5310 q2max = *sm * (*smtol + 1.);
5311 if (q2 <= q2max) {
5312
5313 /* The constraint is satisfied by a constant function. */
5314
5315 *ier = 1;
5316 /* IF (LUN .GE. 0) WRITE (LUN,100) */
5317 /* 100 FORMAT (///1X,'SMSURF -- THE CONSTRAINT IS NOT ', */
5318 /* . 'ACTIVE AND THE FITTING FCN IS CONSTANT.') */
5319 if (dbg_verbose) fprintf (stderr, "SMSURF -- THE CONSTRAINT IS NOT ACTIVE AND THE FITTING FCN IS CONSTANT\n");
5320 return 0;
5321 }
5322
5323 /* Compute G0 = G(0) and print a heading. */
5324
5325 *ier = 0;
5326 s = 1. / sqrt(*sm);
5327 g0 = 1. / sqrt(q2) - s;
5328 /* IF (LUN .GE. 0) WRITE (LUN,110) SM, TOL, NITMAX, G0 */
5329 /* 110 FORMAT (///1X,'SMSURF -- SM = ',E10.4,', GSTOL = ', */
5330 /* . E7.1,', NITMAX = ',I2,', G(0) = ',E15.8) */
5331 if (dbg_verbose) fprintf (stderr, "SMSURF -- SM = %g GSTOL = %g NITMAX = %d G(0) = %g\n", *sm, tol, nitmax, g0);
5332
5333 /* G(P) is strictly increasing and concave, and G(0) .LT. 0. */
5334 /* Initialize parameters for the secant method. The method */
5335 /* uses three points -- (P0,G0), (P,G), and (PNEG,GNEG) */
5336 /* where P0 and PNEG are defined implicitly by DP = P - P0 */
5337 /* and DMAX = P - PNEG. */
5338
5339 p = *sm * 10.;
5340 dp = p;
5341 dmax__ = 0.;
5342 iter = 0;
5343
5344 /* Top of loop -- compute G. */
5345
5346 L3:
5347 nit = nitmax;
5348 dfmax = tol;
5349 smsgs_(&nn, &x[1], &y[1], &z__[1], &u[1], &list[1], &lptr[1], &lend[1],
5350 iflgs, &sigma[1], &w[1], &p, &nit, &dfmax, &f[1], &grad[4], &ierr)
5351 ;
5352 if (ierr < 0) {
5353 *ier = ierr;
5354 }
5355
5356 /* IERR = -1 in SMSGS could be caused by P = 0 as a result */
5357 /* of inaccurate solutions to ill-conditioned systems. */
5358
5359 if (ierr == -1) {
5360 *ier = 2;
5361 }
5362 if (ierr < 0) {
5363 return 0;
5364 }
5365 q2 = 0.;
5366 i__1 = nn;
5367 for (i__ = 1; i__ <= i__1; ++i__) {
5368 /* Computing 2nd power */
5369 d__1 = u[i__] - f[i__];
5370 q2 += w[i__] * (d__1 * d__1);
5371 /* L4: */
5372 }
5373 g = 1. / sqrt(q2) - s;
5374 ++iter;
5375 /* IF (LUN .GE. 0) WRITE (LUN,120) ITER, P, G, NIT, DFMAX */
5376 /* 120 FORMAT (/1X,I2,' -- P = ',E15.8,', G = ',E15.8, */
5377 /* . ', NIT = ',I2,', DFMAX = ',E12.6) */
5378 if (dbg_verbose) fprintf (stderr, " %d -- P = %g G = %g NIT = %d DFMAX = %g\n", iter, p, g, nit, dfmax);
5379
5380 /* Test for convergence. */
5381
5382 if (q2min <= q2 && q2 <= q2max) {
5383 return 0;
5384 }
5385 if (iter >= itmax) {
5386 *ier = 2;
5387 return 0;
5388 }
5389 if (dmax__ == 0. && g <= 0.) {
5390
5391 /* Increase P until G(P) > 0. */
5392
5393 p *= 10.;
5394 dp = p;
5395 goto L3;
5396 }
5397
5398 /* A bracketing interval [P0,P] has been found. */
5399
5400 if (g0 * g <= 0.) {
5401
5402 /* G0*G < 0. Update (PNEG,GNEG) to (P0,G0) so that G */
5403 /* and GNEG always have opposite signs. */
5404
5405 dmax__ = dp;
5406 gneg = g0;
5407 }
5408
5409 /* Compute the change in P by linear interpolation between */
5410 /* (P0,G0) and (P,G). */
5411
5412 L5:
5413 dp = -g * dp / (g - g0);
5414 /* IF (LUN .GE. 0) WRITE (LUN,130) DP */
5415 /* 130 FORMAT (1X,5X,'DP = ',E15.8) */
5416 if (dbg_verbose) fprintf (stderr, "DP = %g\n", dp);
5417 if (fabs(dp) > fabs(dmax__)) {
5418
5419 /* G0*G .GT. 0 and the new estimate would be outside of the */
5420 /* bracketing interval of length ABS(DMAX). Reset */
5421 /* (P0,G0) to (PNEG,GNEG). */
5422
5423 dp = dmax__;
5424 g0 = gneg;
5425 goto L5;
5426 }
5427
5428 /* Bottom of loop -- update P, DMAX, and G0. */
5429
5430 p += dp;
5431 dmax__ += dp;
5432 g0 = g;
5433 goto L3;
5434 } /* smsurf_ */
5435
snhcsh_(doublereal * x,doublereal * sinhm,doublereal * coshm,doublereal * coshmm)5436 /* Subroutine */ int snhcsh_(doublereal *x, doublereal *sinhm, doublereal *
5437 coshm, doublereal *coshmm)
5438 {
5439 /* Initialized data */
5440
5441 static doublereal c1 = .1666666666659;
5442 static doublereal c2 = .008333333431546;
5443 static doublereal c3 = 1.984107350948e-4;
5444 static doublereal c4 = 2.768286868175e-6;
5445
5446 /* Builtin functions */
5447 double exp(doublereal);
5448
5449 /* Local variables */
5450 static doublereal f, ax, xc, xs, xsd2, xsd4, expx;
5451
5452
5453 /* *********************************************************** */
5454
5455 /* From SSRFPACK */
5456 /* Robert J. Renka */
5457 /* Dept. of Computer Science */
5458 /* Univ. of North Texas */
5459 /* renka@cs.unt.edu */
5460 /* 03/18/90 */
5461
5462 /* This subroutine computes approximations to the modified */
5463 /* hyperbolic functions defined below with relative error */
5464 /* bounded by 4.7E-12 for a floating point number system with */
5465 /* sufficient precision. For IEEE standard single precision, */
5466 /* the relative error is less than 1.E-5 for all x. */
5467
5468 /* Note that the 13-digit constants in the data statements */
5469 /* below may not be acceptable to all compilers. */
5470
5471 /* On input: */
5472
5473 /* X = Point at which the functions are to be */
5474 /* evaluated. */
5475
5476 /* X is not altered by this routine. */
5477
5478 /* On output: */
5479
5480 /* SINHM = sinh(X) - X. */
5481
5482 /* COSHM = cosh(X) - 1. */
5483
5484 /* COSHMM = cosh(X) - 1 - X*X/2. */
5485
5486 /* Modules required by SNHCSH: None */
5487
5488 /* Intrinsic functions called by SNHCSH: ABS, EXP */
5489
5490 /* *********************************************************** */
5491
5492
5493 ax = fabs(*x);
5494 xs = ax * ax;
5495 if (ax <= .5) {
5496
5497 /* Approximations for small X: */
5498
5499 xc = *x * xs;
5500 *sinhm = xc * (((c4 * xs + c3) * xs + c2) * xs + c1);
5501 xsd4 = xs * .25;
5502 xsd2 = xsd4 + xsd4;
5503 f = (((c4 * xsd4 + c3) * xsd4 + c2) * xsd4 + c1) * xsd4;
5504 *coshmm = xsd2 * f * (f + 2.);
5505 *coshm = *coshmm + xsd2;
5506 } else {
5507
5508 /* Approximations for large X: */
5509
5510 expx = exp(ax);
5511 *sinhm = -(1. / expx + ax + ax - expx) / 2.;
5512 if (*x < 0.) {
5513 *sinhm = -(*sinhm);
5514 }
5515 *coshm = (1. / expx - 2. + expx) / 2.;
5516 *coshmm = *coshm - xs / 2.;
5517 }
5518 return 0;
5519 } /* snhcsh_ */
5520
unif_(integer * n,doublereal * x,doublereal * y,doublereal * z__,doublereal * f,integer * list,integer * lptr,integer * lend,integer * iflgs,doublereal * sigma,integer * nrow,integer * ni,integer * nj,doublereal * plat,doublereal * plon,integer * iflgg,doublereal * grad,doublereal * ff,integer * ier)5521 /* Subroutine */ int unif_(integer *n, doublereal *x, doublereal *y,
5522 doublereal *z__, doublereal *f, integer *list, integer *lptr, integer
5523 *lend, integer *iflgs, doublereal *sigma, integer *nrow, integer *ni,
5524 integer *nj, doublereal *plat, doublereal *plon, integer *iflgg,
5525 doublereal *grad, doublereal *ff, integer *ier)
5526 {
5527 /* Initialized data */
5528
5529 static integer nst = 1;
5530
5531 /* System generated locals */
5532 integer ff_dim1, ff_offset, i__1, i__2;
5533
5534 /* Local variables */
5535 static integer i__, j, nn, nx, ny, ifl, nex, ist, ierr;
5536 extern /* Subroutine */ int gradl_(integer *, integer *, doublereal *,
5537 doublereal *, doublereal *, doublereal *, integer *, integer *,
5538 integer *, doublereal *, integer *), intrc1_(integer *,
5539 doublereal *, doublereal *, doublereal *, doublereal *,
5540 doublereal *, doublereal *, integer *, integer *, integer *,
5541 integer *, doublereal *, integer *, doublereal *, integer *,
5542 doublereal *, integer *);
5543
5544
5545 /* *********************************************************** */
5546
5547 /* From SSRFPACK */
5548 /* Robert J. Renka */
5549 /* Dept. of Computer Science */
5550 /* Univ. of North Texas */
5551 /* renka@cs.unt.edu */
5552 /* 07/25/96 */
5553
5554 /* Given a Delaunay triangulation of a set of nodes on the */
5555 /* unit sphere, along with data values and tension factors */
5556 /* associated with the triangulation arcs, this routine */
5557 /* interpolates the data values to a uniform grid for such */
5558 /* applications as contouring. The interpolant is once con- */
5559 /* tinuously differentiable. Extrapolation is performed at */
5560 /* grid points exterior to the triangulation when the nodes */
5561 /* do not cover the entire sphere. */
5562
5563 /* On input: */
5564
5565 /* N = Number of nodes. N .GE. 3 and N .GE. 7 if */
5566 /* IFLAG .NE. 1. */
5567
5568 /* X,Y,Z = Arrays containing Cartesian coordinates of */
5569 /* the nodes. X(I)**2 + Y(I)**2 + Z(I)**2 = 1 */
5570 /* for I = 1 to N. */
5571
5572 /* F = Array containing data values. F(I) is associ- */
5573 /* ated with (X(I),Y(I),Z(I)). */
5574
5575 /* LIST,LPTR,LEND = Data structure defining the trian- */
5576 /* gulation. Refer to STRIPACK */
5577 /* Subroutine TRMESH. */
5578
5579 /* IFLGS = Tension factor option: */
5580 /* IFLGS .LE. 0 if a single uniform tension */
5581 /* factor is to be used. */
5582 /* IFLGS .GE. 1 if variable tension is desired. */
5583
5584 /* SIGMA = Uniform tension factor (IFLGS .LE. 0), or */
5585 /* array containing tension factors associated */
5586 /* with arcs in one-to-one correspondence with */
5587 /* LIST entries (IFLGS .GE. 1). Refer to Sub- */
5588 /* programs GETSIG, SIG0, SIG1, and SIG2. */
5589
5590 /* NROW = Number of rows in the dimension statement of */
5591 /* FF. */
5592
5593 /* NI,NJ = Number of rows and columns in the uniform */
5594 /* grid. 1 .LE. NI .LE. NROW and 1 .LE. NJ. */
5595
5596 /* PLAT,PLON = Arrays of length NI and NJ, respective- */
5597 /* ly, containing the latitudes and */
5598 /* longitudes of the grid lines. */
5599
5600 /* IFLGG = Option indicator: */
5601 /* IFLGG = 0 if gradient estimates at the ver- */
5602 /* tices of a triangle are to be */
5603 /* recomputed for each grid point in */
5604 /* the triangle and not saved. */
5605 /* IFLGG = 1 if gradient estimates are input in */
5606 /* GRAD. */
5607 /* IFLGG = 2 if gradient estimates are to be */
5608 /* computed once for each node (by */
5609 /* GRADL) and saved in GRAD. */
5610
5611 /* The above parameters are not altered by this routine. */
5612
5613 /* GRAD = 3 by N array whose columns contain the X, Y, */
5614 /* and Z components (in that order) of the grad- */
5615 /* ients at the nodes if IFLGG = 1, array of */
5616 /* sufficient size if IFLGG = 2, or dummy para- */
5617 /* meter if IFLGG = 0. */
5618
5619 /* Gradient estimates may be computed by Subroutines GRADL or */
5620 /* GRADG if IFLGG = 1. */
5621
5622 /* FF = NROW by NCOL array with NROW .GE. NI and NCOL */
5623 /* .GE. NJ. */
5624
5625 /* On output: */
5626
5627 /* GRAD = Array containing estimated gradients as de- */
5628 /* fined above if IFLGG = 2 and IER .GE. 0. */
5629 /* GRAD is not altered if IFLGG .NE. 2. */
5630
5631 /* FF = Interpolated values at the grid points if IER */
5632 /* .GE. 0. FF(I,J) = F(PLAT(I),PLON(J)) for I = */
5633 /* 1,...,NI and J = 1,...,NJ. */
5634
5635 /* IER = Error indicator: */
5636 /* IER = K if no errors were encountered and K */
5637 /* grid points required extrapolation for */
5638 /* K .GE. 0. */
5639 /* IER = -1 if N, NI, NJ, or IFLGG is outside its */
5640 /* valid range. */
5641 /* IER = -2 if the nodes are collinear. */
5642 /* IER = -3 if extrapolation failed due to the */
5643 /* uniform grid extending too far beyond */
5644 /* the triangulation boundary. */
5645
5646 /* STRIPACK modules required by UNIF: GETNP, JRAND, LSTPTR, */
5647 /* STORE, TRFIND */
5648
5649 /* SSRFPACK modules required by UNIF: APLYR, APLYRT, ARCINT, */
5650 /* ARCLEN, CONSTR, */
5651 /* FVAL, GIVENS, GRADL, */
5652 /* HVAL, INTRC1, */
5653 /* ROTATE, SETUP, */
5654 /* SNHCSH */
5655
5656 /* *********************************************************** */
5657
5658 /* Parameter adjustments */
5659 grad -= 4;
5660 --lend;
5661 --f;
5662 --z__;
5663 --y;
5664 --x;
5665 --list;
5666 --lptr;
5667 --sigma;
5668 --plat;
5669 ff_dim1 = *nrow;
5670 ff_offset = 1 + ff_dim1;
5671 ff -= ff_offset;
5672 --plon;
5673
5674 /* Function Body */
5675
5676 /* Local parameters: */
5677
5678 /* I,J = DO-loop indexes */
5679 /* IERR = Error flag for calls to GRADL and INTRC1 */
5680 /* IFL = Local copy of IFLGG */
5681 /* IST = Parameter for INTRC1 */
5682 /* NEX = Number of grid points exterior to the triangula- */
5683 /* tion boundary (number of extrapolated values) */
5684 /* NN = Local copy of N */
5685 /* NST = Initial value for IST */
5686 /* NX,NY = Local copies of NI and NJ */
5687
5688 nn = *n;
5689 nx = *ni;
5690 ny = *nj;
5691 ifl = *iflgg;
5692 if (nx < 1 || nx > *nrow || ny < 1 || ifl < 0 || ifl > 2) {
5693 goto L4;
5694 }
5695 ist = nst;
5696 if (ifl == 2) {
5697
5698 /* Compute gradient estimates at the nodes. */
5699
5700 i__1 = nn;
5701 for (i__ = 1; i__ <= i__1; ++i__) {
5702 gradl_(&nn, &i__, &x[1], &y[1], &z__[1], &f[1], &list[1], &lptr[1]
5703 , &lend[1], &grad[i__ * 3 + 1], &ierr);
5704 if (ierr < 0) {
5705 goto L5;
5706 }
5707 /* L1: */
5708 }
5709 ifl = 1;
5710 }
5711
5712 /* Compute uniform grid points and interpolated values. */
5713
5714 nex = 0;
5715 i__1 = ny;
5716 for (j = 1; j <= i__1; ++j) {
5717 i__2 = nx;
5718 for (i__ = 1; i__ <= i__2; ++i__) {
5719 intrc1_(&nn, &plat[i__], &plon[j], &x[1], &y[1], &z__[1], &f[1], &
5720 list[1], &lptr[1], &lend[1], iflgs, &sigma[1], &ifl, &
5721 grad[4], &ist, &ff[i__ + j * ff_dim1], &ierr);
5722 if (ierr < 0) {
5723 goto L5;
5724 }
5725 nex += ierr;
5726 /* L2: */
5727 }
5728 /* L3: */
5729 }
5730 *ier = nex;
5731 return 0;
5732
5733 /* NI, NJ, or IFLGG is outside its valid range. */
5734
5735 L4:
5736 *ier = -1;
5737 return 0;
5738
5739 /* Error in GRADL or INTRC1. */
5740
5741 L5:
5742 *ier = ierr;
5743 return 0;
5744 } /* unif_ */
5745
5746 }
5747 #ifdef __cplusplus
5748 }
5749 #endif
5750