/Proj4/biveval.c

http://github.com/route-me/route-me · C · 93 lines · 86 code · 5 blank · 2 comment · 16 complexity · e3c1974751411554bfa8af4794380368 MD5 · raw file

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