/Proj4/biveval.c
C | 93 lines | 86 code | 5 blank | 2 comment | 16 complexity | e3c1974751411554bfa8af4794380368 MD5 | raw file
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