restore the functions from biveval.c as they are used by some folks (#24)

git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@1514 4e78687f-474d-0410-85f9-8d5e500ac6b2

4 files changed, 98 insertions, 6 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 9b4a5603..9a4095cf 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -48,7 +48,7 @@ libproj_la_SOURCES = \
 	PJ_sts.c PJ_urm5.c PJ_urmfps.c PJ_wag2.c \
 	PJ_wag3.c PJ_wink1.c PJ_wink2.c pj_latlong.c pj_geocent.c \
 	aasincos.c adjlon.c bch2bps.c bchgen.c \
-	dmstor.c mk_cheby.c pj_auth.c \
+	biveval.c dmstor.c mk_cheby.c pj_auth.c \
 	pj_deriv.c pj_ell_set.c pj_ellps.c pj_errno.c \
 	pj_factors.c pj_fwd.c pj_init.c pj_inv.c \
 	pj_list.c pj_malloc.c pj_mlfn.c pj_msfn.c proj_mdist.c \
diff --git a/src/Makefile.in b/src/Makefile.in
index a9a18b42..3eff10e7 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -77,10 +77,10 @@ am_libproj_la_OBJECTS = PJ_aeqd.lo PJ_gnom.lo PJ_laea.lo \
 	PJ_putp6.lo PJ_robin.lo PJ_sts.lo PJ_urm5.lo PJ_urmfps.lo \
 	PJ_wag2.lo PJ_wag3.lo PJ_wink1.lo PJ_wink2.lo pj_latlong.lo \
 	pj_geocent.lo aasincos.lo adjlon.lo bch2bps.lo bchgen.lo \
-	dmstor.lo mk_cheby.lo pj_auth.lo pj_deriv.lo pj_ell_set.lo \
-	pj_ellps.lo pj_errno.lo pj_factors.lo pj_fwd.lo pj_init.lo \
-	pj_inv.lo pj_list.lo pj_malloc.lo pj_mlfn.lo pj_msfn.lo \
-	proj_mdist.lo pj_open_lib.lo pj_param.lo pj_phi2.lo \
+	biveval.lo dmstor.lo mk_cheby.lo pj_auth.lo pj_deriv.lo \
+	pj_ell_set.lo pj_ellps.lo pj_errno.lo pj_factors.lo pj_fwd.lo \
+	pj_init.lo pj_inv.lo pj_list.lo pj_malloc.lo pj_mlfn.lo \
+	pj_msfn.lo proj_mdist.lo pj_open_lib.lo pj_param.lo pj_phi2.lo \
 	pj_pr_list.lo pj_qsfn.lo pj_strerrno.lo pj_tsfn.lo pj_units.lo \
 	pj_zpoly1.lo rtodms.lo vector1.lo pj_release.lo pj_gauss.lo \
 	nad_cvt.lo nad_init.lo nad_intr.lo emess.lo \
@@ -280,7 +280,7 @@ libproj_la_SOURCES = \
 	PJ_sts.c PJ_urm5.c PJ_urmfps.c PJ_wag2.c \
 	PJ_wag3.c PJ_wink1.c PJ_wink2.c pj_latlong.c pj_geocent.c \
 	aasincos.c adjlon.c bch2bps.c bchgen.c \
-	dmstor.c mk_cheby.c pj_auth.c \
+	biveval.c dmstor.c mk_cheby.c pj_auth.c \
 	pj_deriv.c pj_ell_set.c pj_ellps.c pj_errno.c \
 	pj_factors.c pj_fwd.c pj_init.c pj_inv.c \
 	pj_list.c pj_malloc.c pj_mlfn.c pj_msfn.c proj_mdist.c \
@@ -516,6 +516,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/adjlon.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/bch2bps.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/bchgen.Plo@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/biveval.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cs2cs.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/dmstor.Plo@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/emess.Plo@am__quote@
diff --git a/src/biveval.c b/src/biveval.c
new file mode 100644
index 00000000..d0f0e743
--- /dev/null
+++ b/src/biveval.c
@@ -0,0 +1,88 @@
+/* procedures for evaluating Tseries */
+#ifndef lint
+static const char SCCSID[]="@(#)biveval.c	4.4	93/06/12	GIE	REL";
+#endif
+# include <projects.h>
+# define NEAR_ONE	1.00001
+static double ceval(struct PW_COEF *C, int n, projUV w, projUV w2) {
+	double d=0, dd=0, vd, vdd, tmp, *c;
+	int j;
+
+	for (C += n ; n-- ; --C ) {
+		if (j = C->m) {
+			vd = vdd = 0.;
+			for (c = C->c + --j; j ; --j ) {
+				vd = w2.v * (tmp = vd) - vdd + *c--;
+				vdd = tmp;
+			}
+			d = w2.u * (tmp = d) - dd + w.v * vd - vdd + 0.5 * *c;
+		} else
+			d = w2.u * (tmp = d) - dd;
+		dd = tmp;
+	}
+	if (j = C->m) {
+		vd = vdd = 0.;
+		for (c = C->c + --j; j ; --j ) {
+			vd = w2.v * (tmp = vd) - vdd + *c--;
+			vdd = tmp;
+		}
+		return (w.u * d - dd + 0.5 * ( w.v * vd - vdd + 0.5 * *c ));
+	} else
+		return (w.u * d - dd);
+}
+	projUV /* bivariate Chebyshev polynomial entry point */
+bcheval(projUV in, Tseries *T) {
+        projUV w2, w;
+	projUV out;
+		/* scale to +-1 */
+ 	w.u = ( in.u + in.u - T->a.u ) * T->b.u;
+ 	w.v = ( in.v + in.v - T->a.v ) * T->b.v;
+	if (fabs(w.u) > NEAR_ONE || fabs(w.v) > NEAR_ONE) {
+		out.u = out.v = HUGE_VAL;
+		pj_errno = -36;
+	} else { /* double evaluation */
+		w2.u = w.u + w.u;
+		w2.v = w.v + w.v;
+                out.u = ceval(T->cu, T->mu, w, w2);
+                out.v = ceval(T->cv, T->mv, w, w2);
+	}
+	return out;
+}
+	projUV /* bivariate power polynomial entry point */
+bpseval(projUV in, Tseries *T) {
+	projUV out;
+	double *c, row;
+	int i, m;
+
+	out.u = out.v = 0.;
+	for (i = T->mu; i >= 0; --i) {
+		row = 0.;
+		if (m = T->cu[i].m) {
+			c = T->cu[i].c + m;
+			while (m--)
+				row = *--c + in.v * row;
+		}
+		out.u = row + in.u * out.u;
+	}
+	for (i = T->mv; i >= 0; --i) {
+		row = 0.;
+		if (m = T->cv[i].m) {
+			c = T->cv[i].c + m;
+			while (m--)
+				row = *--c + in.v * row;
+		}
+		out.v = row + in.u * out.v;
+	}
+	return out;
+}
+
+projUV /* general entry point selecting evaluation mode */
+biveval(projUV in, Tseries *T) {
+
+    if (T->power) {
+        return bpseval(in, T);
+    } else {
+        return bcheval(in, T);
+    }
+}
+
diff --git a/src/projects.h b/src/projects.h
index 7cdbe2a9..6ed75648 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -373,6 +373,9 @@ typedef struct {	/* Chebyshev or Power series structure */
 	int power;		/* != 0 if power series, else Chebyshev */
 } Tseries;
 Tseries *mk_cheby(projUV, projUV, double, projUV *, projUV (*)(projUV), int, int, int);
+projUV bpseval(projUV, Tseries *);
+projUV bcheval(projUV, Tseries *);
+projUV biveval(projUV, Tseries *);
 void *vector1(int, int);
 void **vector2(int, int, int);
 void freev2(void **v, int nrows);