/Proj4/biveval.c

 1/* procedures for evaluating Tseries */
 2#ifndef lint
 3static const char SCCSID[]="@(#)biveval.c	4.4	93/06/12	GIE	REL";
 4#endif
 5# include "projects.h"
 6# define NEAR_ONE	1.00001
 7	static projUV
 8w2, w;
 9static double ceval(struct PW_COEF *C, int n) {
10	double d=0, dd=0, vd, vdd, tmp, *c;
11	int j;
12
13	for (C += n ; n-- ; --C ) {
14		j = C->m;
15		if (j) {
16			vd = vdd = 0.;
17			for (c = C->c + --j; j ; --j ) {
18				vd = w2.v * (tmp = vd) - vdd + *c--;
19				vdd = tmp;
20			}
21			d = w2.u * (tmp = d) - dd + w.v * vd - vdd + 0.5 * *c;
22		} else
23			d = w2.u * (tmp = d) - dd;
24		dd = tmp;
25	}
26	j = C->m;
27	if (j) {
28		vd = vdd = 0.;
29		for (c = C->c + --j; j ; --j ) {
30			vd = w2.v * (tmp = vd) - vdd + *c--;
31			vdd = tmp;
32		}
33		return (w.u * d - dd + 0.5 * ( w.v * vd - vdd + 0.5 * *c ));
34	} else
35		return (w.u * d - dd);
36}
37	projUV /* bivariate Chebyshev polynomial entry point */
38bcheval(projUV in, Tseries *T) {
39	projUV out;
40		/* scale to +-1 */
41 	w.u = ( in.u + in.u - T->a.u ) * T->b.u;
42 	w.v = ( in.v + in.v - T->a.v ) * T->b.v;
43	if (fabs(w.u) > NEAR_ONE || fabs(w.v) > NEAR_ONE) {
44		out.u = out.v = HUGE_VAL;
45		pj_errno = -36;
46	} else { /* double evaluation */
47		w2.u = w.u + w.u;
48		w2.v = w.v + w.v;
49		out.u = ceval(T->cu, T->mu);
50		out.v = ceval(T->cv, T->mv);
51	}
52	return out;
53}
54	projUV /* bivariate power polynomial entry point */
55bpseval(projUV in, Tseries *T) {
56	projUV out;
57	double *c, row;
58	int i, m;
59
60	out.u = out.v = 0.;
61	for (i = T->mu; i >= 0; --i) {
62		row = 0.;
63		m = T->cu[i].m;
64		if (m) {
65			c = T->cu[i].c + m;
66			while (m--)
67				row = *--c + in.v * row;
68		}
69		out.u = row + in.u * out.u;
70	}
71	for (i = T->mv; i >= 0; --i) {
72		row = 0.;
73		m = T->cv[i].m;
74		if (m) {
75			c = T->cv[i].c + m;
76			while (m--)
77				row = *--c + in.v * row;
78		}
79		out.v = row + in.u * out.v;
80	}
81	return out;
82}
83
84projUV /* general entry point selecting evaluation mode */
85biveval(projUV in, Tseries *T) {
86
87    if (T->power) {
88        return bpseval(in, T);
89    } else {
90        return bcheval(in, T);
91    }
92}
93