/Proj4/mk_cheby.c

http://github.com/route-me/route-me · C · 169 lines · 161 code · 4 blank · 4 comment · 51 complexity · 748c507d3fee6170cae5727ab0183175 MD5 · raw file

  1. #ifndef lint
  2. static const char SCCSID[]="@(#)mk_cheby.c 4.5 94/03/22 GIE REL";
  3. #endif
  4. #include "projects.h"
  5. static void /* sum coefficients less than res */
  6. eval(projUV **w, int nu, int nv, double res, projUV *resid) {
  7. int i, j;
  8. double ab;
  9. projUV *s;
  10. resid->u = resid->v = 0.;
  11. for (i = 0; i < nu; ++i)
  12. for (s = w[i], j = 0; j < nv; ++j, ++s) {
  13. if ((ab = fabs(s->u)) < res)
  14. resid->u += ab;
  15. if ((ab = fabs(s->v)) < res)
  16. resid->v += ab;
  17. }
  18. }
  19. static Tseries * /* create power series structure */
  20. makeT(int nru, int nrv) {
  21. Tseries *Ts;
  22. int i;
  23. if ((Ts = (Tseries *)pj_malloc(sizeof(Tseries))) &&
  24. (Ts->cu = (struct PW_COEF *)pj_malloc(
  25. sizeof(struct PW_COEF) * nru)) &&
  26. (Ts->cv = (struct PW_COEF *)pj_malloc(
  27. sizeof(struct PW_COEF) * nrv))) {
  28. for (i = 0; i < nru; ++i)
  29. Ts->cu[i].c = 0;
  30. for (i = 0; i < nrv; ++i)
  31. Ts->cv[i].c = 0;
  32. return Ts;
  33. } else
  34. return 0;
  35. }
  36. Tseries *
  37. mk_cheby(projUV a, projUV b, double res, projUV *resid, projUV (*func)(projUV),
  38. int nu, int nv, int power) {
  39. int j, i, nru, nrv, *ncu, *ncv;
  40. Tseries *Ts = 0;
  41. projUV **w;
  42. double cutres;
  43. if (!(w = (projUV **)vector2(nu, nv, sizeof(projUV))) ||
  44. !(ncu = (int *)vector1(nu + nv, sizeof(int))))
  45. return 0;
  46. ncv = ncu + nu;
  47. if (!bchgen(a, b, nu, nv, w, func)) {
  48. projUV *s;
  49. double *p;
  50. /* analyse coefficients and adjust until residual OK */
  51. cutres = res;
  52. for (i = 4; i ; --i) {
  53. eval(w, nu, nv, cutres, resid);
  54. if (resid->u < res && resid->v < res)
  55. break;
  56. cutres *= 0.5;
  57. }
  58. if (i <= 0) /* warn of too many tries */
  59. resid->u = - resid->u;
  60. /* apply cut resolution and set pointers */
  61. nru = nrv = 0;
  62. for (j = 0; j < nu; ++j) {
  63. ncu[j] = ncv[j] = 0; /* clear column maxes */
  64. for (s = w[j], i = 0; i < nv; ++i, ++s) {
  65. if (fabs(s->u) < cutres) /* < resolution ? */
  66. s->u = 0.; /* clear coefficient */
  67. else
  68. ncu[j] = i + 1; /* update column max */
  69. if (fabs(s->v) < cutres) /* same for v coef's */
  70. s->v = 0.;
  71. else
  72. ncv[j] = i + 1;
  73. }
  74. if (ncu[j]) nru = j + 1; /* update row max */
  75. if (ncv[j]) nrv = j + 1;
  76. }
  77. if (power) { /* convert to bivariate power series */
  78. if (!bch2bps(a, b, w, nu, nv))
  79. goto error;
  80. /* possible change in some row counts, so readjust */
  81. nru = nrv = 0;
  82. for (j = 0; j < nu; ++j) {
  83. ncu[j] = ncv[j] = 0; /* clear column maxes */
  84. for (s = w[j], i = 0; i < nv; ++i, ++s) {
  85. if (s->u)
  86. ncu[j] = i + 1; /* update column max */
  87. if (s->v)
  88. ncv[j] = i + 1;
  89. }
  90. if (ncu[j]) nru = j + 1; /* update row max */
  91. if (ncv[j]) nrv = j + 1;
  92. }
  93. Ts = makeT(nru, nrv);
  94. if (Ts) {
  95. Ts->a = a;
  96. Ts->b = b;
  97. Ts->mu = nru - 1;
  98. Ts->mv = nrv - 1;
  99. Ts->power = 1;
  100. for (i = 0; i < nru; ++i) /* store coefficient rows for u */
  101. Ts->cu[i].m = ncu[i];
  102. if (Ts->cu[i].m)
  103. if ((p = Ts->cu[i].c =
  104. (double *)pj_malloc(sizeof(double) * ncu[i])))
  105. for (j = 0; j < ncu[i]; ++j)
  106. *p++ = (w[i] + j)->u;
  107. else
  108. goto error;
  109. for (i = 0; i < nrv; ++i) /* same for v */
  110. Ts->cv[i].m = ncv[i];
  111. if (Ts->cv[i].m)
  112. if ((p = Ts->cv[i].c =
  113. (double *)pj_malloc(sizeof(double) * ncv[i])))
  114. for (j = 0; j < ncv[i]; ++j)
  115. *p++ = (w[i] + j)->v;
  116. else
  117. goto error;
  118. }
  119. } else if ((Ts = makeT(nru, nrv))) {
  120. /* else make returned Chebyshev coefficient structure */
  121. Ts->mu = nru - 1; /* save row degree */
  122. Ts->mv = nrv - 1;
  123. Ts->a.u = a.u + b.u; /* set argument scaling */
  124. Ts->a.v = a.v + b.v;
  125. Ts->b.u = 1. / (b.u - a.u);
  126. Ts->b.v = 1. / (b.v - a.v);
  127. Ts->power = 0;
  128. for (i = 0; i < nru; ++i) /* store coefficient rows for u */
  129. Ts->cu[i].m = ncu[i];
  130. if (Ts->cu[i].m)
  131. if ((p = Ts->cu[i].c =
  132. (double *)pj_malloc(sizeof(double) * ncu[i])))
  133. for (j = 0; j < ncu[i]; ++j)
  134. *p++ = (w[i] + j)->u;
  135. else
  136. goto error;
  137. for (i = 0; i < nrv; ++i) /* same for v */
  138. Ts->cv[i].m = ncv[i];
  139. if (Ts->cv[i].m)
  140. if ((p = Ts->cv[i].c =
  141. (double *)pj_malloc(sizeof(double) * ncv[i])))
  142. for (j = 0; j < ncv[i]; ++j)
  143. *p++ = (w[i] + j)->v;
  144. else
  145. goto error;
  146. } else
  147. goto error;
  148. }
  149. goto gohome;
  150. error:
  151. if (Ts) { /* pj_dalloc up possible allocations */
  152. for (i = 0; i <= Ts->mu; ++i)
  153. if (Ts->cu[i].c)
  154. pj_dalloc(Ts->cu[i].c);
  155. for (i = 0; i <= Ts->mv; ++i)
  156. if (Ts->cv[i].c)
  157. pj_dalloc(Ts->cv[i].c);
  158. pj_dalloc(Ts);
  159. }
  160. Ts = 0;
  161. gohome:
  162. freev2((void **) w, nu);
  163. pj_dalloc(ncu);
  164. return Ts;
  165. }