/Proj4/mk_cheby.c

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