roll projCtx through various other low level functions

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

73 files changed, 344 insertions, 330 deletions

diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index 98ec0a6b..7729a13e 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -153,10 +153,10 @@ setup(PJ *P) {
 	return P;
 }
 ENTRY1(aea,en)
-	P->phi1 = pj_param(P->params, "rlat_1").f;
-	P->phi2 = pj_param(P->params, "rlat_2").f;
+	P->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
+	P->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
 ENDENTRY(setup(P))
 ENTRY1(leac,en)
-	P->phi2 = pj_param(P->params, "rlat_1").f;
-	P->phi1 = pj_param(P->params, "bsouth").i ? - HALFPI: HALFPI;
+	P->phi2 = pj_param(P->ctx, P->params, "rlat_1").f;
+	P->phi1 = pj_param(P->ctx, P->params, "bsouth").i ? - HALFPI: HALFPI;
 ENDENTRY(setup(P))
diff --git a/src/PJ_aeqd.c b/src/PJ_aeqd.c
index 53942c82..fa7c3f63 100644
--- a/src/PJ_aeqd.c
+++ b/src/PJ_aeqd.c
@@ -87,7 +87,7 @@ FORWARD(e_forward); /* elliptical */
 		ct = cos(t); st = sin(t);
 		Az = atan2(sin(lp.lam) * ct, P->cosph0 * st - P->sinph0 * coslam * ct);
 		cA = cos(Az); sA = sin(Az);
-		s = aasin( fabs(sA) < TOL ?
+		s = aasin( P->ctx, fabs(sA) < TOL ?
 			(P->cosph0 * st - P->sinph0 * coslam * ct) / cA :
 			sin(lp.lam) * ct / sA );
 		H = P->He * cA;
@@ -147,7 +147,7 @@ INVERSE(e_guam_inv); /* Guam elliptical */
 	lp.phi = P->phi0;
 	for (i = 0; i < 3; ++i) {
 		t = P->e * sin(lp.phi);
-		lp.phi = pj_inv_mlfn(P->M1 + xy.y -
+		lp.phi = pj_inv_mlfn(P->ctx, P->M1 + xy.y -
 			x2 * tan(lp.phi) * (t = sqrt(1. - t * t)), P->es, P->en);
 	}
 	lp.lam = xy.x * t / cos(lp.phi);
@@ -170,8 +170,8 @@ INVERSE(e_inverse); /* elliptical */
 		D = c / P->N1;
 		E = D * (1. - D * D * (A * (1. + A) / 6. + B * (1. + 3.*A) * D / 24.));
 		F = 1. - E * E * (A / 2. + B * E / 6.);
-		psi = aasin(P->sinph0 * cos(E) + t * sin(E));
-		lp.lam = aasin(sin(Az) * sin(E) / cos(psi));
+		psi = aasin(P->ctx, P->sinph0 * cos(E) + t * sin(E));
+		lp.lam = aasin(P->ctx, sin(Az) * sin(E) / cos(psi));
 		if ((t = fabs(psi)) < EPS10)
 			lp.phi = 0.;
 		else if (fabs(t - HALFPI) < 0.)
@@ -180,7 +180,7 @@ INVERSE(e_inverse); /* elliptical */
 			lp.phi = atan((1. - P->es * F * P->sinph0 / sin(psi)) * tan(psi) /
 				P->one_es);
 	} else { /* Polar */
-		lp.phi = pj_inv_mlfn(P->mode == N_POLE ? P->Mp - c : P->Mp + c,
+		lp.phi = pj_inv_mlfn(P->ctx, P->mode == N_POLE ? P->Mp - c : P->Mp + c,
 			P->es, P->en);
 		lp.lam = atan2(xy.x, P->mode == N_POLE ? -xy.y : xy.y);
 	}
@@ -201,11 +201,11 @@ INVERSE(s_inverse); /* spherical */
 		sinc = sin(c_rh);
 		cosc = cos(c_rh);
 		if (P->mode == EQUIT) {
-			lp.phi = aasin(xy.y * sinc / c_rh);
+                        lp.phi = aasin(P->ctx, xy.y * sinc / c_rh);
 			xy.x *= sinc;
 			xy.y = cosc * c_rh;
 		} else {
-			lp.phi = aasin(cosc * P->sinph0 + xy.y * sinc * P->cosph0 /
+			lp.phi = aasin(P->ctx,cosc * P->sinph0 + xy.y * sinc * P->cosph0 /
 				c_rh);
 			xy.y = (cosc - P->sinph0 * sin(lp.phi)) * c_rh;
 			xy.x *= sinc * P->cosph0;
@@ -228,7 +228,7 @@ FREEUP;
 	}
 }
 ENTRY1(aeqd, en)
-	P->phi0 = pj_param(P->params, "rlat_0").f;
+	P->phi0 = pj_param(P->ctx, P->params, "rlat_0").f;
 	if (fabs(fabs(P->phi0) - HALFPI) < EPS10) {
 		P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
 		P->sinph0 = P->phi0 < 0. ? -1. : 1.;
@@ -246,7 +246,7 @@ ENTRY1(aeqd, en)
 		P->inv = s_inverse; P->fwd = s_forward;
 	} else {
 		if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
-		if (pj_param(P->params, "bguam").i) {
+		if (pj_param(P->ctx, P->params, "bguam").i) {
 			P->M1 = pj_mlfn(P->phi0, P->sinph0, P->cosph0, P->en);
 			P->inv = e_guam_inv; P->fwd = e_guam_fwd;
 		} else {
diff --git a/src/PJ_airy.c b/src/PJ_airy.c
index cc764175..1675df1e 100644
--- a/src/PJ_airy.c
+++ b/src/PJ_airy.c
@@ -94,8 +94,8 @@ FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(airy)
 	double beta;
 
-	P->no_cut = pj_param(P->params, "bno_cut").i;
-	beta = 0.5 * (HALFPI - pj_param(P->params, "rlat_b").f);
+	P->no_cut = pj_param(P->ctx, P->params, "bno_cut").i;
+	beta = 0.5 * (HALFPI - pj_param(P->ctx, P->params, "rlat_b").f);
 	if (fabs(beta) < EPS)
 		P->Cb = -0.5;
 	else {
diff --git a/src/PJ_aitoff.c b/src/PJ_aitoff.c
index c0d958ba..b01578e0 100644
--- a/src/PJ_aitoff.c
+++ b/src/PJ_aitoff.c
@@ -66,9 +66,9 @@ ENTRY0(aitoff)
 ENDENTRY(setup(P))
 ENTRY0(wintri)
 	P->mode = 1;
-	if (pj_param(P->params, "tlat_1").i)
+	if (pj_param(P->ctx, P->params, "tlat_1").i)
         {
-		if ((P->cosphi1 = cos(pj_param(P->params, "rlat_1").f)) == 0.)
+		if ((P->cosphi1 = cos(pj_param(P->ctx, P->params, "rlat_1").f)) == 0.)
 			E_ERROR(-22)
         }
 	else /* 50d28' or acos(2/pi) */
diff --git a/src/PJ_bipc.c b/src/PJ_bipc.c
index 5cb74256..471e7266 100644
--- a/src/PJ_bipc.c
+++ b/src/PJ_bipc.c
@@ -125,7 +125,7 @@ INVERSE(s_inverse); /* spheroid */
 }
 FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(bipc)
-	P->noskew = pj_param(P->params, "bns").i;
+	P->noskew = pj_param(P->ctx, P->params, "bns").i;
 	P->inv = s_inverse;
 	P->fwd = s_forward;
 	P->es = 0.;
diff --git a/src/PJ_bonne.c b/src/PJ_bonne.c
index 81a5bf98..4b1e7866 100644
--- a/src/PJ_bonne.c
+++ b/src/PJ_bonne.c
@@ -45,7 +45,7 @@ INVERSE(e_inverse); /* ellipsoid */
 	double s, rh;
 
 	rh = hypot(xy.x, xy.y = P->am1 - xy.y);
-	lp.phi = pj_inv_mlfn(P->am1 + P->m1 - rh, P->es, P->en);
+	lp.phi = pj_inv_mlfn(P->ctx, P->am1 + P->m1 - rh, P->es, P->en);
 	if ((s = fabs(lp.phi)) < HALFPI) {
 		s = sin(lp.phi);
 		lp.lam = rh * atan2(xy.x, xy.y) *
@@ -65,7 +65,7 @@ FREEUP;
 ENTRY1(bonne, en)
 	double c;
 
-	P->phi1 = pj_param(P->params, "rlat_1").f;
+	P->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
 	if (fabs(P->phi1) < EPS10) E_ERROR(-23);
 	if (P->es) {
 		P->en = pj_enfn(P->es);
diff --git a/src/PJ_cass.c b/src/PJ_cass.c
index 4539c66b..38fa9db5 100644
--- a/src/PJ_cass.c
+++ b/src/PJ_cass.c
@@ -40,7 +40,7 @@ FORWARD(s_forward); /* spheroid */
 INVERSE(e_inverse); /* ellipsoid */
 	double ph1;
 
-	ph1 = pj_inv_mlfn(P->m0 + xy.y, P->es, P->en);
+	ph1 = pj_inv_mlfn(P->ctx, P->m0 + xy.y, P->es, P->en);
 	P->tn = tan(ph1); P->t = P->tn * P->tn;
 	P->n = sin(ph1);
 	P->r = 1. / (1. - P->es * P->n * P->n);
diff --git a/src/PJ_cea.c b/src/PJ_cea.c
index 97e4f3aa..786eadcb 100644
--- a/src/PJ_cea.c
+++ b/src/PJ_cea.c
@@ -42,8 +42,8 @@ FREEUP;
 ENTRY1(cea, apa)
 	double t;
 
-	if (pj_param(P->params, "tlat_ts").i &&
-		(P->k0 = cos(t = pj_param(P->params, "rlat_ts").f)) < 0.) E_ERROR(-24)
+	if (pj_param(P->ctx, P->params, "tlat_ts").i &&
+		(P->k0 = cos(t = pj_param(P->ctx, P->params, "rlat_ts").f)) < 0.) E_ERROR(-24)
 	else
 		t = 0.;
 	if (P->es) {
diff --git a/src/PJ_chamb.c b/src/PJ_chamb.c
index c17a7133..65f21129 100644
--- a/src/PJ_chamb.c
+++ b/src/PJ_chamb.c
@@ -17,17 +17,17 @@ PROJ_HEAD(chamb, "Chamberlin Trimetric") "\n\tMisc Sph, no inv."
 #define THIRD 0.333333333333333333
 #define TOL 1e-9
 	static VECT /* distance and azimuth from point 1 to point 2 */
-vect(double dphi, double c1, double s1, double c2, double s2, double dlam) {
+vect(projCtx ctx, double dphi, double c1, double s1, double c2, double s2, double dlam) {
 	VECT v;
 	double cdl, dp, dl;
 
 	cdl = cos(dlam);
 	if (fabs(dphi) > 1. || fabs(dlam) > 1.)
-		v.r = aacos(s1 * s2 + c1 * c2 * cdl);
+		v.r = aacos(ctx, s1 * s2 + c1 * c2 * cdl);
 	else { /* more accurate for smaller distances */
 		dp = sin(.5 * dphi);
 		dl = sin(.5 * dlam);
-		v.r = 2. * aasin(sqrt(dp * dp + c1 * c2 * dl * dl));
+		v.r = 2. * aasin(ctx,sqrt(dp * dp + c1 * c2 * dl * dl));
 	}
 	if (fabs(v.r) > TOL)
 		v.Az = atan2(c2 * sin(dlam), c1 * s2 - s1 * c2 * cdl);
@@ -36,8 +36,8 @@ vect(double dphi, double c1, double s1, double c2, double s2, double dlam) {
 	return v;
 }
 	static double /* law of cosines */
-lc(double b,double c,double a) {
-	return aacos(.5 * (b * b + c * c - a * a) / (b * c));
+lc(projCtx ctx, double b,double c,double a) {
+	return aacos(ctx, .5 * (b * b + c * c - a * a) / (b * c));
 }
 FORWARD(s_forward); /* spheroid */
 	double sinphi, cosphi, a;
@@ -47,7 +47,7 @@ FORWARD(s_forward); /* spheroid */
 	sinphi = sin(lp.phi);
 	cosphi = cos(lp.phi);
 	for (i = 0; i < 3; ++i) { /* dist/azimiths from control */
-		v[i] = vect(lp.phi - P->c[i].phi, P->c[i].cosphi, P->c[i].sinphi,
+		v[i] = vect(P->ctx, lp.phi - P->c[i].phi, P->c[i].cosphi, P->c[i].sinphi,
 			cosphi, sinphi, lp.lam - P->c[i].lam);
 		if ( ! v[i].r)
 			break;
@@ -59,7 +59,7 @@ FORWARD(s_forward); /* spheroid */
 		xy = P->p;
 		for (i = 0; i < 3; ++i) {
 			j = i == 2 ? 0 : i + 1;
-			a = lc(P->c[i].v.r, v[i].r, v[j].r);
+			a = lc(P->ctx,P->c[i].v.r, v[i].r, v[j].r);
 			if (v[i].Az < 0.)
 				a = -a;
 			if (! i) { /* coord comp unique to each arc */
@@ -87,22 +87,22 @@ ENTRY0(chamb)
 
 	for (i = 0; i < 3; ++i) { /* get control point locations */
 		(void)sprintf(line, "rlat_%d", i+1);
-		P->c[i].phi = pj_param(P->params, line).f;
+		P->c[i].phi = pj_param(P->ctx, P->params, line).f;
 		(void)sprintf(line, "rlon_%d", i+1);
-		P->c[i].lam = pj_param(P->params, line).f;
+		P->c[i].lam = pj_param(P->ctx, P->params, line).f;
 		P->c[i].lam = adjlon(P->c[i].lam - P->lam0);
 		P->c[i].cosphi = cos(P->c[i].phi);
 		P->c[i].sinphi = sin(P->c[i].phi);
 	}
 	for (i = 0; i < 3; ++i) { /* inter ctl pt. distances and azimuths */
 		j = i == 2 ? 0 : i + 1;
-		P->c[i].v = vect(P->c[j].phi - P->c[i].phi, P->c[i].cosphi, P->c[i].sinphi,
+		P->c[i].v = vect(P->ctx,P->c[j].phi - P->c[i].phi, P->c[i].cosphi, P->c[i].sinphi,
 			P->c[j].cosphi, P->c[j].sinphi, P->c[j].lam - P->c[i].lam);
 		if (! P->c[i].v.r) E_ERROR(-25);
 		/* co-linearity problem ignored for now */
 	}
-	P->beta_0 = lc(P->c[0].v.r, P->c[2].v.r, P->c[1].v.r);
-	P->beta_1 = lc(P->c[0].v.r, P->c[1].v.r, P->c[2].v.r);
+	P->beta_0 = lc(P->ctx,P->c[0].v.r, P->c[2].v.r, P->c[1].v.r);
+	P->beta_1 = lc(P->ctx,P->c[0].v.r, P->c[1].v.r, P->c[2].v.r);
 	P->beta_2 = PI - P->beta_0;
 	P->p.y = 2. * (P->c[0].p.y = P->c[1].p.y = P->c[2].v.r * sin(P->beta_0));
 	P->c[2].p.y = 0.;
diff --git a/src/PJ_eck4.c b/src/PJ_eck4.c
index 2b5c8101..f90267dd 100644
--- a/src/PJ_eck4.c
+++ b/src/PJ_eck4.c
@@ -35,9 +35,9 @@ FORWARD(s_forward); /* spheroid */
 INVERSE(s_inverse); /* spheroid */
 	double c;
 
-	lp.phi = aasin(xy.y / C_y);
+	lp.phi = aasin(P->ctx,xy.y / C_y);
 	lp.lam = xy.x / (C_x * (1. + (c = cos(lp.phi))));
-	lp.phi = aasin((lp.phi + sin(lp.phi) * (c + 2.)) / C_p);
+	lp.phi = aasin(P->ctx,(lp.phi + sin(lp.phi) * (c + 2.)) / C_p);
 	return (lp);
 }
 FREEUP; if (P) pj_dalloc(P); }
diff --git a/src/PJ_eqc.c b/src/PJ_eqc.c
index 533fb63e..c69a6baa 100644
--- a/src/PJ_eqc.c
+++ b/src/PJ_eqc.c
@@ -16,7 +16,7 @@ INVERSE(s_inverse); /* spheroid */
 }
 FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(eqc)
-	if ((P->rc = cos(pj_param(P->params, "rlat_ts").f)) <= 0.) E_ERROR(-24);
+	if ((P->rc = cos(pj_param(P->ctx, P->params, "rlat_ts").f)) <= 0.) E_ERROR(-24);
 	P->inv = s_inverse;
 	P->fwd = s_forward;
 	P->es = 0.;
diff --git a/src/PJ_eqdc.c b/src/PJ_eqdc.c
index eda95b0d..01a0f03f 100644
--- a/src/PJ_eqdc.c
+++ b/src/PJ_eqdc.c
@@ -28,7 +28,7 @@ INVERSE(e_inverse); /* sphere & ellipsoid */
 		}
 		lp.phi = P->c - P->rho;
 		if (P->ellips)
-			lp.phi = pj_inv_mlfn(lp.phi, P->es, P->en);
+			lp.phi = pj_inv_mlfn(P->ctx, lp.phi, P->es, P->en);
 		lp.lam = atan2(xy.x, xy.y) / P->n;
 	} else {
 		lp.lam = 0.;
@@ -51,8 +51,8 @@ ENTRY1(eqdc, en)
 	double cosphi, sinphi;
 	int secant;
 
-	P->phi1 = pj_param(P->params, "rlat_1").f;
-	P->phi2 = pj_param(P->params, "rlat_2").f;
+	P->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
+	P->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
 	if (fabs(P->phi1 + P->phi2) < EPS10) E_ERROR(-21);
 	if (!(P->en = pj_enfn(P->es)))
 		E_ERROR_0;
diff --git a/src/PJ_fouc_s.c b/src/PJ_fouc_s.c
index dee05a96..b84b3f82 100644
--- a/src/PJ_fouc_s.c
+++ b/src/PJ_fouc_s.c
@@ -28,14 +28,14 @@ INVERSE(s_inverse); /* spheroid */
 		if (!i)
 			lp.phi = xy.y < 0. ? -HALFPI : HALFPI;
 	} else
-		lp.phi = aasin(xy.y);
+		lp.phi = aasin(P->ctx,xy.y);
 	V = cos(lp.phi);
 	lp.lam = xy.x * (P->n + P->n1 * V) / V;
 	return (lp);
 }
 FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(fouc_s)
-	P->n = pj_param(P->params, "dn").f;
+	P->n = pj_param(P->ctx, P->params, "dn").f;
 	if (P->n < 0. || P->n > 1.)
 		E_ERROR(-99)
 	P->n1 = 1. - P->n;
diff --git a/src/PJ_geos.c b/src/PJ_geos.c
index ff99d7b3..9b9d36f4 100644
--- a/src/PJ_geos.c
+++ b/src/PJ_geos.c
@@ -124,7 +124,7 @@ INVERSE(e_inverse); /* ellipsoid */
 }
 FREEUP; if (P) free(P); }
 ENTRY0(geos)
-	if ((P->h = pj_param(P->params, "dh").f) <= 0.) E_ERROR(-30);
+	if ((P->h = pj_param(P->ctx, P->params, "dh").f) <= 0.) E_ERROR(-30);
 	if (P->phi0) E_ERROR(-46);
 	P->radius_g = 1. + (P->radius_g_1 = P->h / P->a);
 	P->C  = P->radius_g * P->radius_g - 1.0;
diff --git a/src/PJ_gn_sinu.c b/src/PJ_gn_sinu.c
index e08ba83c..bfd8bc2d 100644
--- a/src/PJ_gn_sinu.c
+++ b/src/PJ_gn_sinu.c
@@ -21,7 +21,7 @@ FORWARD(e_forward); /* ellipsoid */
 INVERSE(e_inverse); /* ellipsoid */
 	double s;
 
-	if ((s = fabs(lp.phi = pj_inv_mlfn(xy.y, P->es, P->en))) < HALFPI) {
+	if ((s = fabs(lp.phi = pj_inv_mlfn(P->ctx, xy.y, P->es, P->en))) < HALFPI) {
 		s = sin(lp.phi);
 		lp.lam = xy.x * sqrt(1. - P->es * s * s) / cos(lp.phi);
 	} else if ((s - EPS10) < HALFPI)
@@ -32,7 +32,7 @@ INVERSE(e_inverse); /* ellipsoid */
 /* General spherical sinusoidals */
 FORWARD(s_forward); /* sphere */
 	if (!P->m)
-		lp.phi = P->n != 1. ? aasin(P->n * sin(lp.phi)): lp.phi;
+		lp.phi = P->n != 1. ? aasin(P->ctx,P->n * sin(lp.phi)): lp.phi;
 	else {
 		double k, V;
 		int i;
@@ -53,8 +53,8 @@ FORWARD(s_forward); /* sphere */
 }
 INVERSE(s_inverse); /* sphere */
 	xy.y /= P->C_y;
-	lp.phi = P->m ? aasin((P->m * xy.y + sin(xy.y)) / P->n) :
-		( P->n != 1. ? aasin(sin(xy.y) / P->n) : xy.y );
+	lp.phi = P->m ? aasin(P->ctx,(P->m * xy.y + sin(xy.y)) / P->n) :
+		( P->n != 1. ? aasin(P->ctx,sin(xy.y) / P->n) : xy.y );
 	lp.lam = xy.x / (P->C_x * (P->m + cos(xy.y)));
 	return (lp);
 }
@@ -89,9 +89,9 @@ ENTRY1(mbtfps, en)
 	setup(P);
 ENDENTRY(P)
 ENTRY1(gn_sinu, en)
-	if (pj_param(P->params, "tn").i && pj_param(P->params, "tm").i) {
-		P->n = pj_param(P->params, "dn").f;
-		P->m = pj_param(P->params, "dm").f;
+	if (pj_param(P->ctx, P->params, "tn").i && pj_param(P->ctx, P->params, "tm").i) {
+		P->n = pj_param(P->ctx, P->params, "dn").f;
+		P->m = pj_param(P->ctx, P->params, "dm").f;
 	} else
 		E_ERROR(-99)
 	setup(P);
diff --git a/src/PJ_gstmerc.c b/src/PJ_gstmerc.c
index 23b86ee8..bffe0b26 100644
--- a/src/PJ_gstmerc.c
+++ b/src/PJ_gstmerc.c
@@ -28,7 +28,7 @@ INVERSE(s_inverse); /* spheroid */
     sinC= sin((xy.y*P->a - P->YS)/P->n2)/cosh((xy.x*P->a - P->XS)/P->n2);
     LC= log(pj_tsfn(-1.0*asin(sinC),0.0,0.0));
     lp.lam= L/P->n1;
-    lp.phi= -1.0*pj_phi2(exp((LC-P->c)/P->n1),P->e);
+    lp.phi= -1.0*pj_phi2(P->ctx, exp((LC-P->c)/P->n1),P->e);
     /*fprintf(stderr,"inv:\nL      =%16.13f\nsinC   =%16.13f\nLC     =%16.13f\nXY(%16.4f,%16.4f)=LP(%16.13f,%16.13f)\n",L,sinC,LC,((xy.x/P->ra)+P->x0)/P->to_meter,((xy.y/P->ra)+P->y0)/P->to_meter,lp.lam+P->lam0,lp.phi);*/
 	return (lp);
 }
diff --git a/src/PJ_hammer.c b/src/PJ_hammer.c
index a2ba02aa..719a2ce4 100644
--- a/src/PJ_hammer.c
+++ b/src/PJ_hammer.c
@@ -15,12 +15,12 @@ FORWARD(s_forward); /* spheroid */
 }
 FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(hammer)
-	if (pj_param(P->params, "tW").i) {
-		if ((P->w = fabs(pj_param(P->params, "dW").f)) <= 0.) E_ERROR(-27);
+	if (pj_param(P->ctx, P->params, "tW").i) {
+		if ((P->w = fabs(pj_param(P->ctx, P->params, "dW").f)) <= 0.) E_ERROR(-27);
 	} else
 		P->w = .5;
-	if (pj_param(P->params, "tM").i) {
-		if ((P->m = fabs(pj_param(P->params, "dM").f)) <= 0.) E_ERROR(-27);
+	if (pj_param(P->ctx, P->params, "tM").i) {
+		if ((P->m = fabs(pj_param(P->ctx, P->params, "dM").f)) <= 0.) E_ERROR(-27);
 	} else
 		P->m = 1.;
 	P->rm = 1. / P->m;
diff --git a/src/PJ_imw_p.c b/src/PJ_imw_p.c
index 135b3dab..ae411116 100644
--- a/src/PJ_imw_p.c
+++ b/src/PJ_imw_p.c
@@ -13,12 +13,12 @@ PROJ_HEAD(imw_p, "International Map of the World Polyconic")
 phi12(PJ *P, double *del, double *sig) {
 	int err = 0;
 
-	if (!pj_param(P->params, "tlat_1").i ||
-		!pj_param(P->params, "tlat_2").i) {
+	if (!pj_param(P->ctx, P->params, "tlat_1").i ||
+		!pj_param(P->ctx, P->params, "tlat_2").i) {
 		err = -41;
 	} else {
-		P->phi_1 = pj_param(P->params, "rlat_1").f;
-		P->phi_2 = pj_param(P->params, "rlat_2").f;
+		P->phi_1 = pj_param(P->ctx, P->params, "rlat_1").f;
+		P->phi_2 = pj_param(P->ctx, P->params, "rlat_2").f;
 		*del = 0.5 * (P->phi_2 - P->phi_1);
 		*sig = 0.5 * (P->phi_2 + P->phi_1);
 		err = (fabs(*del) < EPS || fabs(*sig) < EPS) ? -42 : 0;
@@ -113,8 +113,8 @@ ENTRY1(imw_p, en)
 		P->phi_1 = P->phi_2;
 		P->phi_2 = del;
 	}
-	if (pj_param(P->params, "tlon_1").i)
-		P->lam_1 = pj_param(P->params, "rlon_1").f;
+	if (pj_param(P->ctx, P->params, "tlon_1").i)
+		P->lam_1 = pj_param(P->ctx, P->params, "rlon_1").f;
 	else { /* use predefined based upon latitude */
 		sig = fabs(sig * RAD_TO_DEG);
 		if (sig <= 60)		sig = 2.;
diff --git a/src/PJ_krovak.c b/src/PJ_krovak.c
index 0b46b8db..bd00aea6 100644
--- a/src/PJ_krovak.c
+++ b/src/PJ_krovak.c
@@ -121,7 +121,7 @@ FORWARD(e_forward); /* ellipsoid */
 	xy.y = ro * cos(eps) / a;
 	xy.x = ro * sin(eps) / a;
 
-        if( !pj_param(P -> params, "tczech").i )
+        if( !pj_param(P->ctx, P->params, "tczech").i )
 	  {
 	    xy.y *= -1.0;
 	    xy.x *= -1.0;
@@ -176,7 +176,7 @@ INVERSE(e_inverse); /* ellipsoid */
 	xy.x=xy.y;
 	xy.y=xy0;
 
-        if( !pj_param(P -> params, "tczech").i )
+        if( !pj_param(P->ctx, P->params, "tczech").i )
 	  {
 	    xy.x *= -1.0;
 	    xy.y *= -1.0;
@@ -221,7 +221,7 @@ ENTRY0(krovak)
 	/* read some Parameters,
 	 * here Latitude Truescale */
 
-	ts = pj_param(P->params, "rlat_ts").f;
+	ts = pj_param(P->ctx, P->params, "rlat_ts").f;
 	P->C_x = ts;
 	
 	/* we want Bessel as fixed ellipsoid */
@@ -229,17 +229,17 @@ ENTRY0(krovak)
 	P->e = sqrt(P->es = 0.006674372230614);
 
         /* if latitude of projection center is not set, use 49d30'N */
-	if (!pj_param(P->params, "tlat_0").i)
+	if (!pj_param(P->ctx, P->params, "tlat_0").i)
             P->phi0 = 0.863937979737193; 
 
         /* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */
         /* that will correspond to using longitudes relative to greenwich    */
         /* as input and output, instead of lat/long relative to Ferro */
-	if (!pj_param(P->params, "tlon_0").i)
+	if (!pj_param(P->ctx, P->params, "tlon_0").i)
             P->lam0 = 0.7417649320975901 - 0.308341501185665;
 
         /* if scale not set default to 0.9999 */
-	if (!pj_param(P->params, "tk").i)
+	if (!pj_param(P->ctx, P->params, "tk").i)
             P->k0 = 0.9999;
 
 	/* always the same */
diff --git a/src/PJ_labrd.c b/src/PJ_labrd.c
index 70cd2f42..4cb39ec8 100644
--- a/src/PJ_labrd.c
+++ b/src/PJ_labrd.c
@@ -86,8 +86,8 @@ FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(labrd)
 	double Az, sinp, R, N, t;
 
-	P->rot	= pj_param(P->params, "bno_rot").i == 0;
-	Az = pj_param(P->params, "razi").f;
+	P->rot	= pj_param(P->ctx, P->params, "bno_rot").i == 0;
+	Az = pj_param(P->ctx, P->params, "razi").f;
 	sinp = sin(P->phi0);
 	t = 1. - P->es * sinp * sinp;
 	N = 1. / sqrt(t);
diff --git a/src/PJ_lagrng.c b/src/PJ_lagrng.c
index d0e123a9..8a13b3da 100644
--- a/src/PJ_lagrng.c
+++ b/src/PJ_lagrng.c
@@ -26,9 +26,9 @@ FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(lagrng)
 	double phi1;
 
-	if ((P->rw = pj_param(P->params, "dW").f) <= 0) E_ERROR(-27);
+	if ((P->rw = pj_param(P->ctx, P->params, "dW").f) <= 0) E_ERROR(-27);
 	P->hrw = 0.5 * (P->rw = 1. / P->rw);
-	phi1 = pj_param(P->params, "rlat_1").f;
+	phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
 	if (fabs(fabs(phi1 = sin(phi1)) - 1.) < TOL) E_ERROR(-22);
 	P->a1 = pow((1. - phi1)/(1. + phi1), P->hrw);
 	P->es = 0.; P->fwd = s_forward;
diff --git a/src/PJ_lcc.c b/src/PJ_lcc.c
index 80a85d0b..9d3494bf 100644
--- a/src/PJ_lcc.c
+++ b/src/PJ_lcc.c
@@ -34,7 +34,7 @@ INVERSE(e_inverse); /* ellipsoid & spheroid */
 			xy.y = -xy.y;
 		}
 		if (P->ellips) {
-			if ((lp.phi = pj_phi2(pow(rho / P->c, 1./P->n), P->e))
+			if ((lp.phi = pj_phi2(P->ctx, pow(rho / P->c, 1./P->n), P->e))
 				== HUGE_VAL)
 				I_ERROR;
 		} else
@@ -64,12 +64,12 @@ ENTRY0(lcc)
 	double cosphi, sinphi;
 	int secant;
 
-	P->phi1 = pj_param(P->params, "rlat_1").f;
-	if (pj_param(P->params, "tlat_2").i)
-		P->phi2 = pj_param(P->params, "rlat_2").f;
+	P->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
+	if (pj_param(P->ctx, P->params, "tlat_2").i)
+		P->phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
 	else {
 		P->phi2 = P->phi1;
-		if (!pj_param(P->params, "tlat_0").i)
+		if (!pj_param(P->ctx, P->params, "tlat_0").i)
 			P->phi0 = P->phi1;
 	}
 	if (fabs(P->phi1 + P->phi2) < EPS10) E_ERROR(-21);
diff --git a/src/PJ_lcca.c b/src/PJ_lcca.c
index 9d8fe90a..dc27a282 100644
--- a/src/PJ_lcca.c
+++ b/src/PJ_lcca.c
@@ -46,7 +46,7 @@ INVERSE(e_inverse); /* ellipsoid & spheroid */
 		if (fabs(dif) < DEL_TOL) break;
 	}
 	if (!i) I_ERROR
-	lp.phi = pj_inv_mlfn(S + P->M0, P->es, P->en);
+	lp.phi = pj_inv_mlfn(P->ctx, S + P->M0, P->es, P->en);
 	return (lp);
 }
 FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
@@ -54,7 +54,7 @@ ENTRY0(lcca)
 	double s2p0, N0, R0, tan0, tan20;
 
 	if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
-	if (!pj_param(P->params, "tlat_0").i) E_ERROR(50);
+	if (!pj_param(P->ctx, P->params, "tlat_0").i) E_ERROR(50);
 	if (P->phi0 == 0.) E_ERROR(51);
 	P->l = sin(P->phi0);
 	P->M0 = pj_mlfn(P->phi0, P->l, cos(P->phi0), P->en);
diff --git a/src/PJ_loxim.c b/src/PJ_loxim.c
index 80c15318..595423f7 100644
--- a/src/PJ_loxim.c
+++ b/src/PJ_loxim.c
@@ -33,7 +33,7 @@ INVERSE(s_inverse); /* spheroid */
 }
 FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(loxim);
-	P->phi1 = pj_param(P->params, "rlat_1").f;
+	P->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
 	if ((P->cosphi1 = cos(P->phi1)) < EPS) E_ERROR(-22);
 	P->tanphi1 = tan(FORTPI + 0.5 * P->phi1);
 	P->inv = s_inverse; P->fwd = s_forward;
diff --git a/src/PJ_lsat.c b/src/PJ_lsat.c
index e8edb542..d11a5c14 100644
--- a/src/PJ_lsat.c
+++ b/src/PJ_lsat.c
@@ -66,7 +66,7 @@ FORWARD(e_forward); /* ellipsoid */
 	}
 	if (l) {
 		sp = sin(lp.phi);
-		phidp = aasin((P->one_es * P->ca * sp - P->sa * cos(lp.phi) * 
+		phidp = aasin(P->ctx,(P->one_es * P->ca * sp - P->sa * cos(lp.phi) * 
 			sin(lamt)) / sqrt(1. - P->es * sp * sp));
 		tanph = log(tan(FORTPI + .5 * phidp));
 		sd = sin(lamdp);
@@ -116,7 +116,7 @@ INVERSE(e_inverse); /* ellipsoid */
 	lamt -= HALFPI * (1. - scl) * sl;
 	lp.lam = lamt - P->p22 * lamdp;
 	if (fabs(P->sa) < TOL)
-	    lp.phi = aasin(spp / sqrt(P->one_es * P->one_es + P->es * sppsq));
+	    lp.phi = aasin(P->ctx,spp / sqrt(P->one_es * P->one_es + P->es * sppsq));
 	else
 		lp.phi = atan((tan(lamdp) * cos(lamt) - P->ca * sin(lamt)) /
 			(P->one_es * P->sa));
@@ -127,9 +127,9 @@ ENTRY0(lsat)
     int land, path;
     double lam, alf, esc, ess;
 
-	land = pj_param(P->params, "ilsat").i;
+	land = pj_param(P->ctx, P->params, "ilsat").i;
 	if (land <= 0 || land > 5) E_ERROR(-28);
-	path = pj_param(P->params, "ipath").i;
+	path = pj_param(P->ctx, P->params, "ipath").i;
 	if (path <= 0 || path > (land <= 3 ? 251 : 233)) E_ERROR(-29);
 	if (land <= 3) {
 		P->lam0 = DEG_TO_RAD * 128.87 - TWOPI / 251. * path;
diff --git a/src/PJ_mbt_fps.c b/src/PJ_mbt_fps.c
index c41f2e37..57a96e72 100644
--- a/src/PJ_mbt_fps.c
+++ b/src/PJ_mbt_fps.c
@@ -29,9 +29,9 @@ FORWARD(s_forward); /* spheroid */
 INVERSE(s_inverse); /* spheroid */
 	double t;
 
-	lp.phi = C2 * (t = aasin(xy.y / C_y));
+	lp.phi = C2 * (t = aasin(P->ctx,xy.y / C_y));
 	lp.lam = xy.x / (C_x * (1. + 3. * cos(lp.phi)/cos(t)));
-	lp.phi = aasin((C1 * sin(t) + sin(lp.phi)) / C3);
+	lp.phi = aasin(P->ctx,(C1 * sin(t) + sin(lp.phi)) / C3);
 	return (lp);
 }
 FREEUP; if (P) pj_dalloc(P); }
diff --git a/src/PJ_merc.c b/src/PJ_merc.c
index f8e0426a..4b991c1c 100644
--- a/src/PJ_merc.c
+++ b/src/PJ_merc.c
@@ -15,7 +15,7 @@ FORWARD(s_forward); /* spheroid */
 	return (xy);
 }
 INVERSE(e_inverse); /* ellipsoid */
-	if ((lp.phi = pj_phi2(exp(- xy.y / P->k0), P->e)) == HUGE_VAL) I_ERROR;
+	if ((lp.phi = pj_phi2(P->ctx, exp(- xy.y / P->k0), P->e)) == HUGE_VAL) I_ERROR;
 	lp.lam = xy.x / P->k0;
 	return (lp);
 }
@@ -29,8 +29,8 @@ ENTRY0(merc)
 	double phits=0.0;
 	int is_phits;
 
-	if( (is_phits = pj_param(P->params, "tlat_ts").i) ) {
-		phits = fabs(pj_param(P->params, "rlat_ts").f);
+	if( (is_phits = pj_param(P->ctx, P->params, "tlat_ts").i) ) {
+		phits = fabs(pj_param(P->ctx, P->params, "rlat_ts").f);
 		if (phits >= HALFPI) E_ERROR(-24);
 	}
 	if (P->es) { /* ellipsoid */
diff --git a/src/PJ_mod_ster.c b/src/PJ_mod_ster.c
index 854f35e1..185e95e3 100644
--- a/src/PJ_mod_ster.c
+++ b/src/PJ_mod_ster.c
@@ -60,7 +60,7 @@ INVERSE(e_inverse); /* ellipsoid */
 			lp.phi = P->phi0;
 			return lp;
 		}
-		chi = aasin(cosz * P->schio + p.i * sinz * P->cchio / rh);
+		chi = aasin(P->ctx, cosz * P->schio + p.i * sinz * P->cchio / rh);
 		phi = chi;
 		for (nn = 20; nn ;--nn) {
 			esphi = P->e * sin(phi);
diff --git a/src/PJ_moll.c b/src/PJ_moll.c
index f63efdbd..98f0457d 100644
--- a/src/PJ_moll.c
+++ b/src/PJ_moll.c
@@ -27,10 +27,10 @@ FORWARD(s_forward); /* spheroid */
 	return (xy);
 }
 INVERSE(s_inverse); /* spheroid */
-	lp.phi = aasin(xy.y / P->C_y);
+	lp.phi = aasin(P->ctx, xy.y / P->C_y);
 	lp.lam = xy.x / (P->C_x * cos(lp.phi));
 	lp.phi += lp.phi;
-	lp.phi = aasin((lp.phi + sin(lp.phi)) / P->C_p);
+	lp.phi = aasin(P->ctx, (lp.phi + sin(lp.phi)) / P->C_p);
 	return (lp);
 }
 FREEUP; if (P) pj_dalloc(P); }
diff --git a/src/PJ_nell.c b/src/PJ_nell.c
index ed275860..dd41b468 100644
--- a/src/PJ_nell.c
+++ b/src/PJ_nell.c
@@ -22,7 +22,7 @@ FORWARD(s_forward); /* spheroid */
 }
 INVERSE(s_inverse); /* spheroid */
 	lp.lam = 2. * xy.x / (1. + cos(xy.y));
-	lp.phi = aasin(0.5 * (xy.y + sin(xy.y)));
+	lp.phi = aasin(P->ctx,0.5 * (xy.y + sin(xy.y)));
 	return (lp);
 }
 FREEUP; if (P) pj_dalloc(P); }
diff --git a/src/PJ_nsper.c b/src/PJ_nsper.c
index a3ee264e..0f355b93 100644
--- a/src/PJ_nsper.c
+++ b/src/PJ_nsper.c
@@ -115,7 +115,7 @@ INVERSE(s_inverse); /* spheroid */
 FREEUP; if (P) pj_dalloc(P); }
 	static PJ *
 setup(PJ *P) {
-	if ((P->height = pj_param(P->params, "dh").f) <= 0.) E_ERROR(-30);
+	if ((P->height = pj_param(P->ctx, P->params, "dh").f) <= 0.) E_ERROR(-30);
 	if (fabs(fabs(P->phi0) - HALFPI) < EPS10)
 		P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
 	else if (fabs(P->phi0) < EPS10)
@@ -141,8 +141,8 @@ ENDENTRY(setup(P))
 ENTRY0(tpers)
 	double omega, gamma;
 
-	omega = pj_param(P->params, "dtilt").f * DEG_TO_RAD;
-	gamma = pj_param(P->params, "dazi").f * DEG_TO_RAD;
+	omega = pj_param(P->ctx, P->params, "dtilt").f * DEG_TO_RAD;
+	gamma = pj_param(P->ctx, P->params, "dazi").f * DEG_TO_RAD;
 	P->tilt = 1;
 	P->cg = cos(gamma); P->sg = sin(gamma);
 	P->cw = cos(omega); P->sw = sin(omega);
diff --git a/src/PJ_ob_tran.c b/src/PJ_ob_tran.c
index 460a600f..34e0f6b8 100644
--- a/src/PJ_ob_tran.c
+++ b/src/PJ_ob_tran.c
@@ -21,7 +21,7 @@ FORWARD(o_forward); /* spheroid */
 	cosphi = cos(lp.phi);
 	lp.lam = adjlon(aatan2(cosphi * sin(lp.lam), P->sphip * cosphi * coslam +
 		P->cphip * sinphi) + P->lamp);
-	lp.phi = aasin(P->sphip * sinphi - P->cphip * cosphi * coslam);
+	lp.phi = aasin(P->ctx,P->sphip * sinphi - P->cphip * cosphi * coslam);
 	return (P->link->fwd(lp, P->link));
 }
 FORWARD(t_forward); /* spheroid */
@@ -32,7 +32,7 @@ FORWARD(t_forward); /* spheroid */
 	cosphi = cos(lp.phi);
 	coslam = cos(lp.lam);
 	lp.lam = adjlon(aatan2(cosphi * sin(lp.lam), sin(lp.phi)) + P->lamp);
-	lp.phi = aasin(- cosphi * coslam);
+	lp.phi = aasin(P->ctx, - cosphi * coslam);
 	return (P->link->fwd(lp, P->link));
 }
 INVERSE(o_inverse); /* spheroid */
@@ -43,7 +43,7 @@ INVERSE(o_inverse); /* spheroid */
 		coslam = cos(lp.lam -= P->lamp);
 		sinphi = sin(lp.phi);
 		cosphi = cos(lp.phi);
-		lp.phi = aasin(P->sphip * sinphi + P->cphip * cosphi * coslam);
+		lp.phi = aasin(P->ctx,P->sphip * sinphi + P->cphip * cosphi * coslam);
 		lp.lam = aatan2(cosphi * sin(lp.lam), P->sphip * cosphi * coslam -
 			P->cphip * sinphi);
 	}
@@ -57,7 +57,7 @@ INVERSE(t_inverse); /* spheroid */
 		cosphi = cos(lp.phi);
 		t = lp.lam - P->lamp;
 		lp.lam = aatan2(cosphi * sin(t), - sin(lp.phi));
-		lp.phi = aasin(cosphi * cos(t));
+		lp.phi = aasin(P->ctx,cosphi * cos(t));
 	}
 	return (lp);
 }
@@ -74,7 +74,7 @@ ENTRY1(ob_tran, link)
 	char *name, *s;
 
 	/* get name of projection to be translated */
-	if (!(name = pj_param(P->params, "so_proj").s)) E_ERROR(-26);
+	if (!(name = pj_param(P->ctx, P->params, "so_proj").s)) E_ERROR(-26);
 	for (i = 0; (s = pj_list[i].id) && strcmp(name, s) ; ++i) ;
 	if (!s || !(P->link = (*pj_list[i].proj)(0))) E_ERROR(-37);
 	/* copy existing header into new */
@@ -97,12 +97,12 @@ ENTRY1(ob_tran, link)
 		freeup(P);
 		return 0;
 	}
-	if (pj_param(P->params, "to_alpha").i) {
+	if (pj_param(P->ctx, P->params, "to_alpha").i) {
 		double lamc, phic, alpha;
 
-		lamc	= pj_param(P->params, "ro_lon_c").f;
-		phic	= pj_param(P->params, "ro_lat_c").f;
-		alpha	= pj_param(P->params, "ro_alpha").f;
+		lamc	= pj_param(P->ctx, P->params, "ro_lon_c").f;
+		phic	= pj_param(P->ctx, P->params, "ro_lat_c").f;
+		alpha	= pj_param(P->ctx, P->params, "ro_alpha").f;
 /*
 		if (fabs(phic) <= TOL ||
 			fabs(fabs(phic) - HALFPI) <= TOL ||
@@ -111,17 +111,17 @@ ENTRY1(ob_tran, link)
 		if (fabs(fabs(phic) - HALFPI) <= TOL)
 			E_ERROR(-32);
 		P->lamp = lamc + aatan2(-cos(alpha), -sin(alpha) * sin(phic));
-		phip = aasin(cos(phic) * sin(alpha));
-	} else if (pj_param(P->params, "to_lat_p").i) { /* specified new pole */
-		P->lamp = pj_param(P->params, "ro_lon_p").f;
-		phip = pj_param(P->params, "ro_lat_p").f;
+		phip = aasin(P->ctx,cos(phic) * sin(alpha));
+	} else if (pj_param(P->ctx, P->params, "to_lat_p").i) { /* specified new pole */
+		P->lamp = pj_param(P->ctx, P->params, "ro_lon_p").f;
+		phip = pj_param(P->ctx, P->params, "ro_lat_p").f;
 	} else { /* specified new "equator" points */
 		double lam1, lam2, phi1, phi2, con;
 
-		lam1 = pj_param(P->params, "ro_lon_1").f;
-		phi1 = pj_param(P->params, "ro_lat_1").f;
-		lam2 = pj_param(P->params, "ro_lon_2").f;
-		phi2 = pj_param(P->params, "ro_lat_2").f;
+		lam1 = pj_param(P->ctx, P->params, "ro_lon_1").f;
+		phi1 = pj_param(P->ctx, P->params, "ro_lat_1").f;
+		lam2 = pj_param(P->ctx, P->params, "ro_lon_2").f;
+		phi2 = pj_param(P->ctx, P->params, "ro_lat_2").f;
 		if (fabs(phi1 - phi2) <= TOL ||
 			(con = fabs(phi1)) <= TOL ||
 			fabs(con - HALFPI) <= TOL ||
diff --git a/src/PJ_ocea.c b/src/PJ_ocea.c
index cdffd6e4..8c7bd161 100644
--- a/src/PJ_ocea.c
+++ b/src/PJ_ocea.c
@@ -41,16 +41,16 @@ ENTRY0(ocea)
 
 	P->rok = P->a / P->k0;
 	P->rtk = P->a * P->k0;
-	if ( pj_param(P->params, "talpha").i) {
-		alpha	= pj_param(P->params, "ralpha").f;
-		lonz = pj_param(P->params, "rlonc").f;
+	if ( pj_param(P->ctx, P->params, "talpha").i) {
+		alpha	= pj_param(P->ctx, P->params, "ralpha").f;
+		lonz = pj_param(P->ctx, P->params, "rlonc").f;
 		P->singam = atan(-cos(alpha)/(-sin(phi_0) * sin(alpha))) + lonz;
 		P->sinphi = asin(cos(phi_0) * sin(alpha));
 	} else {
-		phi_1 = pj_param(P->params, "rlat_1").f;
-		phi_2 = pj_param(P->params, "rlat_2").f;
-		lam_1 = pj_param(P->params, "rlon_1").f;
-		lam_2 = pj_param(P->params, "rlon_2").f;
+		phi_1 = pj_param(P->ctx, P->params, "rlat_1").f;
+		phi_2 = pj_param(P->ctx, P->params, "rlat_2").f;
+		lam_1 = pj_param(P->ctx, P->params, "rlon_1").f;
+		lam_2 = pj_param(P->ctx, P->params, "rlon_2").f;
 		P->singam = atan2(cos(phi_1) * sin(phi_2) * cos(lam_1) -
 			sin(phi_1) * cos(phi_2) * cos(lam_2),
 			sin(phi_1) * cos(phi_2) * sin(lam_2) -
diff --git a/src/PJ_oea.c b/src/PJ_oea.c
index b5e42a1f..b84a7ea6 100644
--- a/src/PJ_oea.c
+++ b/src/PJ_oea.c
@@ -13,9 +13,9 @@ FORWARD(s_forward); /* sphere */
 	sp = sin(lp.phi);
 	cl = cos(lp.lam);
 	Az = aatan2(cp * sin(lp.lam), P->cp0 * sp - P->sp0 * cp * cl) + P->theta;
-	shz = sin(0.5 * aacos(P->sp0 * sp + P->cp0 * cp * cl));
-	M = aasin(shz * sin(Az));
-	N = aasin(shz * cos(Az) * cos(M) / cos(M * P->two_r_m));
+	shz = sin(0.5 * aacos(P->ctx, P->sp0 * sp + P->cp0 * cp * cl));
+	M = aasin(P->ctx, shz * sin(Az));
+	N = aasin(P->ctx, shz * cos(Az) * cos(M) / cos(M * P->two_r_m));
 	xy.y = P->n * sin(N * P->two_r_n);
 	xy.x = P->m * sin(M * P->two_r_m) * cos(N) / cos(N * P->two_r_n);
 	return (xy);
@@ -23,26 +23,26 @@ FORWARD(s_forward); /* sphere */
 INVERSE(s_inverse); /* sphere */
 	double N, M, xp, yp, z, Az, cz, sz, cAz;
 
-	N = P->hn * aasin(xy.y * P->rn);
-	M = P->hm * aasin(xy.x * P->rm * cos(N * P->two_r_n) / cos(N));
+	N = P->hn * aasin(P->ctx,xy.y * P->rn);
+	M = P->hm * aasin(P->ctx,xy.x * P->rm * cos(N * P->two_r_n) / cos(N));
 	xp = 2. * sin(M);
 	yp = 2. * sin(N) * cos(M * P->two_r_m) / cos(M);
 	cAz = cos(Az = aatan2(xp, yp) - P->theta);
-	z = 2. * aasin(0.5 * hypot(xp, yp));
+	z = 2. * aasin(P->ctx, 0.5 * hypot(xp, yp));
 	sz = sin(z);
 	cz = cos(z);
-	lp.phi = aasin(P->sp0 * cz + P->cp0 * sz * cAz);
+	lp.phi = aasin(P->ctx, P->sp0 * cz + P->cp0 * sz * cAz);
 	lp.lam = aatan2(sz * sin(Az),
 		P->cp0 * cz - P->sp0 * sz * cAz);
 	return (lp);
 }
 FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(oea)
-	if (((P->n = pj_param(P->params, "dn").f) <= 0.) ||
-		((P->m = pj_param(P->params, "dm").f) <= 0.))
+	if (((P->n = pj_param(P->ctx, P->params, "dn").f) <= 0.) ||
+		((P->m = pj_param(P->ctx, P->params, "dm").f) <= 0.))
 		E_ERROR(-39)
 	else {
-		P->theta = pj_param(P->params, "rtheta").f;
+		P->theta = pj_param(P->ctx, P->params, "rtheta").f;
 		P->sp0 = sin(P->phi0);
 		P->cp0 = cos(P->phi0);
 		P->rn = 1./ P->n;
diff --git a/src/PJ_omerc.c b/src/PJ_omerc.c
index c8b61a71..51d398d2 100644
--- a/src/PJ_omerc.c
+++ b/src/PJ_omerc.c
@@ -84,7 +84,7 @@ INVERSE(e_inverse); /* ellipsoid */
 		lp.phi = Up < 0. ? -HALFPI : HALFPI;
 	} else {
 		lp.phi = P->E / sqrt((1. + Up) / (1. - Up));
-		if ((lp.phi = pj_phi2(pow(lp.phi, 1. / P->B), P->e)) == HUGE_VAL)
+		if ((lp.phi = pj_phi2(P->ctx, pow(lp.phi, 1. / P->B), P->e)) == HUGE_VAL)
 			I_ERROR;
 		lp.lam = - P->rB * atan2((Sp * P->cosgam -
 			Vp * P->singam), cos(P->BrA * u));
@@ -97,19 +97,19 @@ ENTRY0(omerc)
 		gamma0, lamc, lam1, lam2, phi1, phi2, alpha_c;
 	int alp, gam, no_off = 0;
 
-	P->no_rot = pj_param(P->params, "tno_rot").i;
-        if ((alp = pj_param(P->params, "talpha").i) != 0)
-		alpha_c = pj_param(P->params, "ralpha").f;
-        if ((gam = pj_param(P->params, "tgamma").i) != 0)
-		gamma = pj_param(P->params, "rgamma").f;
+	P->no_rot = pj_param(P->ctx, P->params, "tno_rot").i;
+        if ((alp = pj_param(P->ctx, P->params, "talpha").i) != 0)
+		alpha_c = pj_param(P->ctx, P->params, "ralpha").f;
+        if ((gam = pj_param(P->ctx, P->params, "tgamma").i) != 0)
+		gamma = pj_param(P->ctx, P->params, "rgamma").f;
 	if (alp || gam) {
-		lamc	= pj_param(P->params, "rlonc").f;
-		no_off = pj_param(P->params, "tno_off").i;
+		lamc	= pj_param(P->ctx, P->params, "rlonc").f;
+		no_off = pj_param(P->ctx, P->params, "tno_off").i;
 	} else {
-		lam1 = pj_param(P->params, "rlon_1").f;
-		phi1 = pj_param(P->params, "rlat_1").f;
-		lam2 = pj_param(P->params, "rlon_2").f;
-		phi2 = pj_param(P->params, "rlat_2").f;
+		lam1 = pj_param(P->ctx, P->params, "rlon_1").f;
+		phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
+		lam2 = pj_param(P->ctx, P->params, "rlon_2").f;
+		phi2 = pj_param(P->ctx, P->params, "rlat_2").f;
 		if (fabs(phi1 - phi2) <= TOL ||
 			(con = fabs(phi1)) <= TOL ||
 			fabs(con - HALFPI) <= TOL ||
diff --git a/src/PJ_putp2.c b/src/PJ_putp2.c
index e846f13b..33bc67ae 100644
--- a/src/PJ_putp2.c
+++ b/src/PJ_putp2.c
@@ -31,9 +31,9 @@ FORWARD(s_forward); /* spheroid */
 INVERSE(s_inverse); /* spheroid */
 	double c;
 
-	lp.phi = aasin(xy.y / C_y);
+	lp.phi = aasin(P->ctx,xy.y / C_y);
 	lp.lam = xy.x / (C_x * ((c = cos(lp.phi)) - 0.5));
-	lp.phi = aasin((lp.phi + sin(lp.phi) * (c - 1.)) / C_p);
+	lp.phi = aasin(P->ctx,(lp.phi + sin(lp.phi) * (c - 1.)) / C_p);
 	return (lp);
 }
 FREEUP; if (P) pj_dalloc(P); }
diff --git a/src/PJ_putp4p.c b/src/PJ_putp4p.c
index 35066c2a..32036bbe 100644
--- a/src/PJ_putp4p.c
+++ b/src/PJ_putp4p.c
@@ -5,18 +5,18 @@
 PROJ_HEAD(putp4p, "Putnins P4'") "\n\tPCyl., Sph.";
 PROJ_HEAD(weren, "Werenskiold I") "\n\tPCyl., Sph.";
 FORWARD(s_forward); /* spheroid */
-	lp.phi = aasin(0.883883476 * sin(lp.phi));
+	lp.phi = aasin(P->ctx,0.883883476 * sin(lp.phi));
 	xy.x = P->C_x * lp.lam * cos(lp.phi);
 	xy.x /= cos(lp.phi *= 0.333333333333333);
 	xy.y = P->C_y * sin(lp.phi);
 	return (xy);
 }
 INVERSE(s_inverse); /* spheroid */
-	lp.phi = aasin(xy.y / P->C_y);
+	lp.phi = aasin(P->ctx,xy.y / P->C_y);
 	lp.lam = xy.x * cos(lp.phi) / P->C_x;
 	lp.phi *= 3.;
 	lp.lam /= cos(lp.phi);
-	lp.phi = aasin(1.13137085 * sin(lp.phi));
+	lp.phi = aasin(P->ctx,1.13137085 * sin(lp.phi));
 	return (lp);
 }
 FREEUP; if (P) pj_dalloc(P); }
diff --git a/src/PJ_putp6.c b/src/PJ_putp6.c
index 28e82cfd..5909dd4c 100644
--- a/src/PJ_putp6.c
+++ b/src/PJ_putp6.c
@@ -32,7 +32,7 @@ INVERSE(s_inverse); /* spheroid */
 	lp.phi = xy.y / P->C_y;
 	r = sqrt(1. + lp.phi * lp.phi);
 	lp.lam = xy.x / (P->C_x * (P->D - r));
-	lp.phi = aasin( ( (P->A - r) * lp.phi - log(lp.phi + r) ) / P->B);
+	lp.phi = aasin( P->ctx, ( (P->A - r) * lp.phi - log(lp.phi + r) ) / P->B);
 	return (lp);
 }
 FREEUP; if (P) pj_dalloc(P); }
diff --git a/src/PJ_rpoly.c b/src/PJ_rpoly.c
index 1440327d..7a7670ed 100644
--- a/src/PJ_rpoly.c
+++ b/src/PJ_rpoly.c
@@ -27,7 +27,7 @@ FORWARD(s_forward); /* spheroid */
 }
 FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(rpoly)
-	if ((P->mode = (P->phi1 = fabs(pj_param(P->params, "rlat_ts").f)) > EPS)) {
+	if ((P->mode = (P->phi1 = fabs(pj_param(P->ctx, P->params, "rlat_ts").f)) > EPS)) {
 		P->fxb = 0.5 * sin(P->phi1);
 		P->fxa = 0.5 / P->fxb;
 	}
diff --git a/src/PJ_sconics.c b/src/PJ_sconics.c
index c65635d8..b999f95c 100644
--- a/src/PJ_sconics.c
+++ b/src/PJ_sconics.c
@@ -37,12 +37,12 @@ phi12(PJ *P, double *del) {
 	double p1, p2;
 	int err = 0;
 
-	if (!pj_param(P->params, "tlat_1").i ||
-		!pj_param(P->params, "tlat_2").i) {
+	if (!pj_param(P->ctx, P->params, "tlat_1").i ||
+		!pj_param(P->ctx, P->params, "tlat_2").i) {
 		err = -41;
 	} else {
-		p1 = pj_param(P->params, "rlat_1").f;
-		p2 = pj_param(P->params, "rlat_2").f;
+		p1 = pj_param(P->ctx, P->params, "rlat_1").f;
+		p2 = pj_param(P->ctx, P->params, "rlat_2").f;
 		*del = 0.5 * (p2 - p1);
 		P->sig = 0.5 * (p2 + p1);
 		err = (fabs(*del) < EPS || fabs(P->sig) < EPS) ? -42 : 0;
diff --git a/src/PJ_somerc.c b/src/PJ_somerc.c
index d419e653..f7d1f119 100644
--- a/src/PJ_somerc.c
+++ b/src/PJ_somerc.c
@@ -14,8 +14,8 @@ FORWARD(e_forward);
 		+ P->K)) - HALFPI;
 	lamp = P->c * lp.lam;
 	cp = cos(phip);
-	phipp = aasin(P->cosp0 * sin(phip) - P->sinp0 * cp * cos(lamp));
-	lampp = aasin(cp * sin(lamp) / cos(phipp));
+	phipp = aasin(P->ctx,P->cosp0 * sin(phip) - P->sinp0 * cp * cos(lamp));
+	lampp = aasin(P->ctx,cp * sin(lamp) / cos(phipp));
 	xy.x = P->kR * lampp;
 	xy.y = P->kR * log(tan(FORTPI + 0.5 * phipp));
 	return (xy);
@@ -27,8 +27,8 @@ INVERSE(e_inverse); /* ellipsoid & spheroid */
 	phipp = 2. * (atan(exp(xy.y / P->kR)) - FORTPI);
 	lampp = xy.x / P->kR;
 	cp = cos(phipp);
-	phip = aasin(P->cosp0 * sin(phipp) + P->sinp0 * cp * cos(lampp));
-	lamp = aasin(cp * sin(lampp) / cos(phip));
+	phip = aasin(P->ctx,P->cosp0 * sin(phipp) + P->sinp0 * cp * cos(lampp));
+	lamp = aasin(P->ctx,cp * sin(lampp) / cos(phip));
 	con = (P->K - log(tan(FORTPI + 0.5 * phip)))/P->c;
 	for (i = NITER; i ; --i) {
 		esp = P->e * sin(phip);
@@ -55,7 +55,7 @@ ENTRY0(somerc)
 	cp *= cp;
 	P->c = sqrt(1 + P->es * cp * cp * P->rone_es);
 	sp = sin(P->phi0);
-	P->cosp0 = cos( phip0 = aasin(P->sinp0 = sp / P->c) );
+	P->cosp0 = cos( phip0 = aasin(P->ctx, P->sinp0 = sp / P->c) );
 	sp *= P->e;
 	P->K = log(tan(FORTPI + 0.5 * phip0)) - P->c * (
 		log(tan(FORTPI + 0.5 * P->phi0)) - P->hlf_e *
diff --git a/src/PJ_stere.c b/src/PJ_stere.c
index a31bf97b..2ee63f53 100644
--- a/src/PJ_stere.c
+++ b/src/PJ_stere.c
@@ -226,12 +226,12 @@ setup(PJ *P) { /* general initialization */
 	return P;
 }
 ENTRY0(stere)
-	P->phits = pj_param(P->params, "tlat_ts").i ?
-		P->phits = pj_param(P->params, "rlat_ts").f : HALFPI;
+	P->phits = pj_param(P->ctx, P->params, "tlat_ts").i ?
+		P->phits = pj_param(P->ctx, P->params, "rlat_ts").f : HALFPI;
 ENDENTRY(setup(P))
 ENTRY0(ups)
 	/* International Ellipsoid */
-	P->phi0 = pj_param(P->params, "bsouth").i ? - HALFPI: HALFPI;
+	P->phi0 = pj_param(P->ctx, P->params, "bsouth").i ? - HALFPI: HALFPI;
 	if (!P->es) E_ERROR(-34);
 	P->k0 = .994;
 	P->x0 = 2000000.;
diff --git a/src/PJ_sterea.c b/src/PJ_sterea.c
index 49ab4b28..083fe271 100644
--- a/src/PJ_sterea.c
+++ b/src/PJ_sterea.c
@@ -42,7 +42,7 @@ PROJ_HEAD(sterea, "Oblique Stereographic Alternative")
 FORWARD(e_forward); /* ellipsoid */
 	double cosc, sinc, cosl, k;
 
-	lp = pj_gauss(lp, P->en);
+	lp = pj_gauss(P->ctx, lp, P->en);
 	sinc = sin(lp.phi);
 	cosc = cos(lp.phi);
 	cosl = cos(lp.lam);
@@ -67,7 +67,7 @@ INVERSE(e_inverse); /* ellipsoid */
 		lp.phi = P->phic0;
 		lp.lam = 0.;
 	}
-	return(pj_inv_gauss(lp, P->en));
+	return(pj_inv_gauss(P->ctx, lp, P->en));
 }
 FREEUP; if (P) { if (P->en) free(P->en); free(P); } }
 ENTRY0(sterea)
diff --git a/src/PJ_sts.c b/src/PJ_sts.c
index f1557feb..c35d5cfc 100644
--- a/src/PJ_sts.c
+++ b/src/PJ_sts.c
@@ -27,7 +27,7 @@ INVERSE(s_inverse); /* spheroid */
 	double c;
 	
 	xy.y /= P->C_y;
-	c = cos(lp.phi = P->tan_mode ? atan(xy.y) : aasin(xy.y));
+	c = cos(lp.phi = P->tan_mode ? atan(xy.y) : aasin(P->ctx,xy.y));
 	lp.phi /= P->C_p;
 	lp.lam = xy.x / (P->C_x * cos(lp.phi));
 	if (P->tan_mode)
diff --git a/src/PJ_tmerc.c b/src/PJ_tmerc.c
index 5104608a..3b1e0fb5 100644
--- a/src/PJ_tmerc.c
+++ b/src/PJ_tmerc.c
@@ -89,7 +89,7 @@ FORWARD(s_forward); /* sphere */
 INVERSE(e_inverse); /* ellipsoid */
 	double n, con, cosphi, d, ds, sinphi, t;
 
-	lp.phi = pj_inv_mlfn(P->ml0 + xy.y / P->k0, P->es, P->en);
+	lp.phi = pj_inv_mlfn(P->ctx, P->ml0 + xy.y / P->k0, P->es, P->en);
 	if (fabs(lp.phi) >= HALFPI) {
 		lp.phi = xy.y < 0. ? -HALFPI : HALFPI;
 		lp.lam = 0.;
@@ -157,10 +157,10 @@ ENTRY1(utm, en)
 	int zone;
 
 	if (!P->es) E_ERROR(-34);
-	P->y0 = pj_param(P->params, "bsouth").i ? 10000000. : 0.;
+	P->y0 = pj_param(P->ctx, P->params, "bsouth").i ? 10000000. : 0.;
 	P->x0 = 500000.;
-	if (pj_param(P->params, "tzone").i) /* zone input ? */
-		if ((zone = pj_param(P->params, "izone").i) > 0 && zone <= 60)
+	if (pj_param(P->ctx, P->params, "tzone").i) /* zone input ? */
+		if ((zone = pj_param(P->ctx, P->params, "izone").i) > 0 && zone <= 60)
 			--zone;
 		else
 			E_ERROR(-35)
diff --git a/src/PJ_tpeqd.c b/src/PJ_tpeqd.c
index b4e76aa0..4ab5cf4e 100644
--- a/src/PJ_tpeqd.c
+++ b/src/PJ_tpeqd.c
@@ -10,8 +10,8 @@ FORWARD(s_forward); /* sphere */
 
 	sp = sin(lp.phi);
 	cp = cos(lp.phi);
-	z1 = aacos(P->sp1 * sp + P->cp1 * cp * cos(dl1 = lp.lam + P->dlam2));
-	z2 = aacos(P->sp2 * sp + P->cp2 * cp * cos(dl2 = lp.lam - P->dlam2));
+	z1 = aacos(P->ctx,P->sp1 * sp + P->cp1 * cp * cos(dl1 = lp.lam + P->dlam2));
+	z2 = aacos(P->ctx,P->sp2 * sp + P->cp2 * cp * cos(dl2 = lp.lam - P->dlam2));
 	z1 *= z1;
 	z2 *= z2;
 	xy.x = P->r2z0 * (t = z1 - z2);
@@ -29,13 +29,13 @@ INVERSE(s_inverse); /* sphere */
 	s = cz1 + cz2;
 	d = cz1 - cz2;
 	lp.lam = - atan2(d, (s * P->thz0));
-	lp.phi = aacos(hypot(P->thz0 * s, d) * P->rhshz0);
+	lp.phi = aacos(P->ctx,hypot(P->thz0 * s, d) * P->rhshz0);
 	if ( xy.y < 0. )
 		lp.phi = - lp.phi;
 	/* lam--phi now in system relative to P1--P2 base equator */
 	sp = sin(lp.phi);
 	cp = cos(lp.phi);
-	lp.phi = aasin(P->sa * sp + P->ca * cp * (s = cos(lp.lam -= P->lp)));
+	lp.phi = aasin(P->ctx,P->sa * sp + P->ca * cp * (s = cos(lp.lam -= P->lp)));
 	lp.lam = atan2(cp * sin(lp.lam), P->sa * cp * s - P->ca * sp) + P->lamc;
 	return lp;
 }
@@ -44,10 +44,10 @@ ENTRY0(tpeqd)
 	double lam_1, lam_2, phi_1, phi_2, A12, pp;
 
 	/* get control point locations */
-	phi_1 = pj_param(P->params, "rlat_1").f;
-	lam_1 = pj_param(P->params, "rlon_1").f;
-	phi_2 = pj_param(P->params, "rlat_2").f;
-	lam_2 = pj_param(P->params, "rlon_2").f;
+	phi_1 = pj_param(P->ctx, P->params, "rlat_1").f;
+	lam_1 = pj_param(P->ctx, P->params, "rlon_1").f;
+	phi_2 = pj_param(P->ctx, P->params, "rlat_2").f;
+	lam_2 = pj_param(P->ctx, P->params, "rlon_2").f;
 	if (phi_1 == phi_2 && lam_1 == lam_2) E_ERROR(-25);
 	P->lam0 = adjlon(0.5 * (lam_1 + lam_2));
 	P->dlam2 = adjlon(lam_2 - lam_1);
@@ -58,11 +58,11 @@ ENTRY0(tpeqd)
 	P->cs = P->cp1 * P->sp2;
 	P->sc = P->sp1 * P->cp2;
 	P->ccs = P->cp1 * P->cp2 * sin(P->dlam2);
-	P->z02 = aacos(P->sp1 * P->sp2 + P->cp1 * P->cp2 * cos(P->dlam2));
+	P->z02 = aacos(P->ctx,P->sp1 * P->sp2 + P->cp1 * P->cp2 * cos(P->dlam2));
 	P->hz0 = .5 * P->z02;
 	A12 = atan2(P->cp2 * sin(P->dlam2),
 		P->cp1 * P->sp2 - P->sp1 * P->cp2 * cos(P->dlam2));
-	P->ca = cos(pp = aasin(P->cp1 * sin(A12)));
+	P->ca = cos(pp = aasin(P->ctx,P->cp1 * sin(A12)));
 	P->sa = sin(pp);
 	P->lp = adjlon(atan2(P->cp1 * cos(A12), P->sp1) - P->hz0);
 	P->dlam2 *= .5;
diff --git a/src/PJ_urm5.c b/src/PJ_urm5.c
index 70a5a1d6..6928bb34 100644
--- a/src/PJ_urm5.c
+++ b/src/PJ_urm5.c
@@ -6,7 +6,7 @@ PROJ_HEAD(urm5, "Urmaev V") "\n\tPCyl., Sph.\n\tn= q= alphi=";
 FORWARD(s_forward); /* spheroid */
 	double t;
 
-	t = lp.phi = aasin(P->n * sin(lp.phi));
+	t = lp.phi = aasin(P->ctx,P->n * sin(lp.phi));
 	xy.x = P->m * lp.lam * cos(lp.phi);
 	t *= t;
 	xy.y = lp.phi * (1. + t * P->q3) * P->rmn;
@@ -16,9 +16,9 @@ FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(urm5)
 	double alpha, t;
 
-	P->n = pj_param(P->params, "dn").f;
-	P->q3 = pj_param(P->params, "dq").f / 3.;
-	alpha = pj_param(P->params, "ralpha").f;
+	P->n = pj_param(P->ctx, P->params, "dn").f;
+	P->q3 = pj_param(P->ctx, P->params, "dq").f / 3.;
+	alpha = pj_param(P->ctx, P->params, "ralpha").f;
 	t = P->n * sin(alpha);
 	P->m = cos(alpha) / sqrt(1. - t * t);
 	P->rmn = 1. / (P->m * P->n);
diff --git a/src/PJ_urmfps.c b/src/PJ_urmfps.c
index 4a4f31a9..5c5918ae 100644
--- a/src/PJ_urmfps.c
+++ b/src/PJ_urmfps.c
@@ -7,14 +7,14 @@ PROJ_HEAD(wag1, "Wagner I (Kavraisky VI)") "\n\tPCyl, Sph.";
 #define C_x 0.8773826753
 #define Cy 1.139753528477
 FORWARD(s_forward); /* sphere */
-	lp.phi = aasin(P->n * sin(lp.phi));
+	lp.phi = aasin(P->ctx,P->n * sin(lp.phi));
 	xy.x = C_x * lp.lam * cos(lp.phi);
 	xy.y = P->C_y * lp.phi;
 	return (xy);
 }
 INVERSE(s_inverse); /* sphere */
 	xy.y /= P->C_y;
-	lp.phi = aasin(sin(xy.y) / P->n);
+	lp.phi = aasin(P->ctx,sin(xy.y) / P->n);
 	lp.lam = xy.x / (C_x * cos(xy.y));
 	return (lp);
 }
@@ -28,8 +28,8 @@ setup(PJ *P) {
 	return P;
 }
 ENTRY0(urmfps)
-	if (pj_param(P->params, "tn").i) {
-		P->n = pj_param(P->params, "dn").f;
+	if (pj_param(P->ctx, P->params, "tn").i) {
+		P->n = pj_param(P->ctx, P->params, "dn").f;
 		if (P->n <= 0. || P->n > 1.)
 			E_ERROR(-40)
 	} else
diff --git a/src/PJ_wag2.c b/src/PJ_wag2.c
index 96fd5fcb..f4942bb2 100644
--- a/src/PJ_wag2.c
+++ b/src/PJ_wag2.c
@@ -6,7 +6,7 @@ PROJ_HEAD(wag2, "Wagner II") "\n\tPCyl., Sph.";
 #define C_p1 0.88022
 #define C_p2 0.88550
 FORWARD(s_forward); /* spheroid */
-	lp.phi = aasin(C_p1 * sin(C_p2 * lp.phi));
+	lp.phi = aasin(P->ctx,C_p1 * sin(C_p2 * lp.phi));
 	xy.x = C_x * lp.lam * cos(lp.phi);
 	xy.y = C_y * lp.phi;
 	return (xy);
@@ -14,7 +14,7 @@ FORWARD(s_forward); /* spheroid */
 INVERSE(s_inverse); /* spheroid */
 	lp.phi = xy.y / C_y;
 	lp.lam = xy.x / (C_x * cos(lp.phi));
-	lp.phi = aasin(sin(lp.phi) / C_p1) / C_p2;
+	lp.phi = aasin(P->ctx,sin(lp.phi) / C_p1) / C_p2;
 	return (lp);
 }
 FREEUP; if (P) pj_dalloc(P); }
diff --git a/src/PJ_wag3.c b/src/PJ_wag3.c
index fbab4c49..482e389c 100644
--- a/src/PJ_wag3.c
+++ b/src/PJ_wag3.c
@@ -18,7 +18,7 @@ FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(wag3)
 	double ts;
 
-	ts = pj_param(P->params, "rlat_ts").f;
+	ts = pj_param(P->ctx, P->params, "rlat_ts").f;
 	P->C_x = cos(ts) / cos(2.*ts/3.);
 	P->es = 0.; P->inv = s_inverse; P->fwd = s_forward;
 ENDENTRY(P)
diff --git a/src/PJ_wink1.c b/src/PJ_wink1.c
index 214e6f4f..dfb455e0 100644
--- a/src/PJ_wink1.c
+++ b/src/PJ_wink1.c
@@ -15,6 +15,6 @@ INVERSE(s_inverse); /* spheroid */
 }
 FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(wink1)
-	P->cosphi1 = cos(pj_param(P->params, "rlat_ts").f);
+	P->cosphi1 = cos(pj_param(P->ctx, P->params, "rlat_ts").f);
 	P->es = 0.; P->inv = s_inverse; P->fwd = s_forward;
 ENDENTRY(P)
diff --git a/src/PJ_wink2.c b/src/PJ_wink2.c
index 8eace009..388c790a 100644
--- a/src/PJ_wink2.c
+++ b/src/PJ_wink2.c
@@ -29,6 +29,6 @@ FORWARD(s_forward); /* spheroid */
 }
 FREEUP; if (P) pj_dalloc(P); }
 ENTRY0(wink2)
-	P->cosphi1 = cos(pj_param(P->params, "rlat_1").f);
+	P->cosphi1 = cos(pj_param(P->ctx, P->params, "rlat_1").f);
 	P->es = 0.; P->inv = 0; P->fwd = s_forward;
 ENDENTRY(P)
diff --git a/src/aasincos.c b/src/aasincos.c
index 460590e4..6f39478e 100644
--- a/src/aasincos.c
+++ b/src/aasincos.c
@@ -3,24 +3,26 @@
 #define ONE_TOL	 1.00000000000001
 #define TOL	0.000000001
 #define ATOL 1e-50
+
 	double
-aasin(double v) {
+aasin(projCtx ctx,double v) {
 	double av;
-
+        
 	if ((av = fabs(v)) >= 1.) {
 		if (av > ONE_TOL)
-			pj_errno = -19;
+                        pj_ctx_set_errno( ctx, -19 );
 		return (v < 0. ? -HALFPI : HALFPI);
 	}
 	return asin(v);
 }
+
 	double
-aacos(double v) {
+aacos(projCtx ctx, double v) {
 	double av;
 
 	if ((av = fabs(v)) >= 1.) {
 		if (av > ONE_TOL)
-			pj_errno = -19;
+                        pj_ctx_set_errno( ctx, -19 );
 		return (v < 0. ? PI : 0.);
 	}
 	return acos(v);
diff --git a/src/dmstor.c b/src/dmstor.c
index 1cc99f62..f1da5379 100644
--- a/src/dmstor.c
+++ b/src/dmstor.c
@@ -17,6 +17,11 @@ vm[] = {
 };
 	double
 dmstor(const char *is, char **rs) {
+	return dmstor_ctx( pj_get_default_ctx(), is, rs );
+}
+
+	double
+dmstor_ctx(projCtx ctx, const char *is, char **rs) {
 	int sign, n, nl;
 	char *p, *s, work[MAX_WORK];
 	double v, tv;
@@ -46,7 +51,7 @@ dmstor(const char *is, char **rs) {
 			n = 2; break;
 		case 'r': case 'R':
 			if (nl) {
-				pj_errno = -16;
+				pj_ctx_set_errno( ctx, -16 );
 				return HUGE_VAL;
 			}
 			++s;
@@ -58,7 +63,7 @@ dmstor(const char *is, char **rs) {
 			continue;
 		}
 		if (n < nl) {
-			pj_errno = -16;
+			pj_ctx_set_errno( ctx, -16 );
 			return HUGE_VAL;
 		}
 		v += tv * vm[n];
diff --git a/src/geod_set.c b/src/geod_set.c
index ba2b415e..e43c26c6 100644
--- a/src/geod_set.c
+++ b/src/geod_set.c
@@ -23,7 +23,7 @@ geod_set(int argc, char **argv) {
 	/* set elliptical parameters */
 	if (pj_ell_set(pj_get_default_ctx(),start, &geod_a, &es)) emess(1,"ellipse setup failure");
 	/* set units */
-	if (name = pj_param(start, "sunits").s) {
+	if (name = pj_param(NULL,start, "sunits").s) {
 		char *s;
                 struct PJ_UNITS *unit_list = pj_get_units_ref();
 		for (i = 0; (s = unit_list[i].id) && strcmp(name, s) ; ++i) ;
@@ -43,27 +43,27 @@ geod_set(int argc, char **argv) {
 		geod_f = f2 = f4 = f64 = 0.;
 	}
 	/* check if line or arc mode */
-	if (pj_param(start, "tlat_1").i) {
+	if (pj_param(NULL,start, "tlat_1").i) {
 		double del_S;
 #undef f
-		phi1 = pj_param(start, "rlat_1").f;
-		lam1 = pj_param(start, "rlon_1").f;
-		if (pj_param(start, "tlat_2").i) {
-			phi2 = pj_param(start, "rlat_2").f;
-			lam2 = pj_param(start, "rlon_2").f;
+		phi1 = pj_param(NULL,start, "rlat_1").f;
+		lam1 = pj_param(NULL,start, "rlon_1").f;
+		if (pj_param(NULL,start, "tlat_2").i) {
+			phi2 = pj_param(NULL,start, "rlat_2").f;
+			lam2 = pj_param(NULL,start, "rlon_2").f;
 			geod_inv();
 			geod_pre();
-		} else if (geod_S = pj_param(start, "dS").f) {
-			al12 = pj_param(start, "rA").f;
+		} else if (geod_S = pj_param(NULL,start, "dS").f) {
+			al12 = pj_param(NULL,start, "rA").f;
 			geod_pre();
 			geod_for();
 		} else emess(1,"incomplete geodesic/arc info");
-		if ((n_alpha = pj_param(start, "in_A").i) > 0) {
-			if (!(del_alpha = pj_param(start, "rdel_A").f))
+		if ((n_alpha = pj_param(NULL,start, "in_A").i) > 0) {
+			if (!(del_alpha = pj_param(NULL,start, "rdel_A").f))
 				emess(1,"del azimuth == 0");
-		} else if (del_S = fabs(pj_param(start, "ddel_S").f)) {
+		} else if (del_S = fabs(pj_param(NULL,start, "ddel_S").f)) {
 			n_S = geod_S / del_S + .5;
-		} else if ((n_S = pj_param(start, "in_S").i) <= 0)
+		} else if ((n_S = pj_param(NULL,start, "in_S").i) <= 0)
 			emess(1,"no interval divisor selected");
 	}
 	/* free up linked list */
diff --git a/src/pj_apply_gridshift.c b/src/pj_apply_gridshift.c
index 9ce9c02c..16f59767 100644
--- a/src/pj_apply_gridshift.c
+++ b/src/pj_apply_gridshift.c
@@ -89,7 +89,7 @@ int pj_apply_gridshift_2( PJ *defn, int inverse,
     {
         defn->gridlist = 
             pj_gridlist_from_nadgrids( pj_get_ctx( defn ),
-                                       pj_param(defn->params,"snadgrids").s,
+                                       pj_param(defn->ctx, defn->params,"snadgrids").s,
                                        &(defn->gridlist_count) );
 
         if( defn->gridlist == NULL || defn->gridlist_count == 0 )
diff --git a/src/pj_apply_vgridshift.c b/src/pj_apply_vgridshift.c
index 2d514587..fb35308c 100644
--- a/src/pj_apply_vgridshift.c
+++ b/src/pj_apply_vgridshift.c
@@ -60,7 +60,7 @@ int pj_apply_vgridshift( PJ *defn, const char *listname,
     {
         *gridlist_p = 
             pj_gridlist_from_nadgrids( pj_get_ctx(defn), 
-                                       pj_param(defn->params,listname).s,
+                                       pj_param(defn->ctx,defn->params,listname).s,
                                        gridlist_count_p );
 
         if( *gridlist_p == NULL || *gridlist_count_p == 0 )
diff --git a/src/pj_datum_set.c b/src/pj_datum_set.c
index 37df66b2..8bfe796d 100644
--- a/src/pj_datum_set.c
+++ b/src/pj_datum_set.c
@@ -54,7 +54,7 @@ int pj_datum_set(projCtx ctx, paralist *pl, PJ *projdef)
 /*      definition will last into the pj_ell_set() function called      */
 /*      after this one.                                                 */
 /* -------------------------------------------------------------------- */
-    if( (name = pj_param(pl,"sdatum").s) != NULL )
+    if( (name = pj_param(ctx, pl,"sdatum").s) != NULL )
     {
         paralist *curr;
         const char *s;
@@ -84,7 +84,7 @@ int pj_datum_set(projCtx ctx, paralist *pl, PJ *projdef)
 /* -------------------------------------------------------------------- */
 /*      Check for nadgrids parameter.                                   */
 /* -------------------------------------------------------------------- */
-    if( (nadgrids = pj_param(pl,"snadgrids").s) != NULL )
+    if( (nadgrids = pj_param(ctx, pl,"snadgrids").s) != NULL )
     {
         /* We don't actually save the value separately.  It will continue
            to exist int he param list for use in pj_apply_gridshift.c */
@@ -95,7 +95,7 @@ int pj_datum_set(projCtx ctx, paralist *pl, PJ *projdef)
 /* -------------------------------------------------------------------- */
 /*      Check for towgs84 parameter.                                    */
 /* -------------------------------------------------------------------- */
-    else if( (towgs84 = pj_param(pl,"stowgs84").s) != NULL )
+    else if( (towgs84 = pj_param(ctx, pl,"stowgs84").s) != NULL )
     {
         int    parm_count = 0;
         const char *s;
diff --git a/src/pj_ell_set.c b/src/pj_ell_set.c
index c0d139d7..6ef537e1 100644
--- a/src/pj_ell_set.c
+++ b/src/pj_ell_set.c
@@ -16,12 +16,12 @@ pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
 		/* check for varying forms of ellipsoid input */
 	*a = *es = 0.;
 	/* R takes precedence */
-	if (pj_param(pl, "tR").i)
-		*a = pj_param(pl, "dR").f;
+	if (pj_param(NULL, pl, "tR").i)
+		*a = pj_param(NULL,pl, "dR").f;
 	else { /* probable elliptical figure */
 
 		/* check if ellps present and temporarily append its values to pl */
-                if ((name = pj_param(pl, "sellps").s) != NULL) {
+                if ((name = pj_param(NULL,pl, "sellps").s) != NULL) {
 			char *s;
 
 			for (start = pl; start && start->next ; start = start->next) ;
@@ -31,50 +31,50 @@ pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
 			curr = curr->next = pj_mkparam(pj_ellps[i].major);
 			curr = curr->next = pj_mkparam(pj_ellps[i].ell);
 		}
-		*a = pj_param(pl, "da").f;
-		if (pj_param(pl, "tes").i) /* eccentricity squared */
-			*es = pj_param(pl, "des").f;
-		else if (pj_param(pl, "te").i) { /* eccentricity */
-			e = pj_param(pl, "de").f;
+		*a = pj_param(NULL,pl, "da").f;
+		if (pj_param(NULL,pl, "tes").i) /* eccentricity squared */
+			*es = pj_param(NULL,pl, "des").f;
+		else if (pj_param(NULL,pl, "te").i) { /* eccentricity */
+			e = pj_param(NULL,pl, "de").f;
 			*es = e * e;
-		} else if (pj_param(pl, "trf").i) { /* recip flattening */
-			*es = pj_param(pl, "drf").f;
+		} else if (pj_param(NULL,pl, "trf").i) { /* recip flattening */
+			*es = pj_param(NULL,pl, "drf").f;
 			if (!*es) {
 				pj_ctx_set_errno( ctx, -10);
 				goto bomb;
 			}
 			*es = 1./ *es;
 			*es = *es * (2. - *es);
-		} else if (pj_param(pl, "tf").i) { /* flattening */
-			*es = pj_param(pl, "df").f;
+		} else if (pj_param(NULL,pl, "tf").i) { /* flattening */
+			*es = pj_param(NULL,pl, "df").f;
 			*es = *es * (2. - *es);
-		} else if (pj_param(pl, "tb").i) { /* minor axis */
-			b = pj_param(pl, "db").f;
+		} else if (pj_param(NULL,pl, "tb").i) { /* minor axis */
+			b = pj_param(NULL,pl, "db").f;
 			*es = 1. - (b * b) / (*a * *a);
 		}     /* else *es == 0. and sphere of radius *a */
 		if (!b)
 			b = *a * sqrt(1. - *es);
 		/* following options turn ellipsoid into equivalent sphere */
-		if (pj_param(pl, "bR_A").i) { /* sphere--area of ellipsoid */
+		if (pj_param(NULL,pl, "bR_A").i) { /* sphere--area of ellipsoid */
 			*a *= 1. - *es * (SIXTH + *es * (RA4 + *es * RA6));
 			*es = 0.;
-		} else if (pj_param(pl, "bR_V").i) { /* sphere--vol. of ellipsoid */
+		} else if (pj_param(NULL,pl, "bR_V").i) { /* sphere--vol. of ellipsoid */
 			*a *= 1. - *es * (SIXTH + *es * (RV4 + *es * RV6));
 			*es = 0.;
-		} else if (pj_param(pl, "bR_a").i) { /* sphere--arithmetic mean */
+		} else if (pj_param(NULL,pl, "bR_a").i) { /* sphere--arithmetic mean */
 			*a = .5 * (*a + b);
 			*es = 0.;
-		} else if (pj_param(pl, "bR_g").i) { /* sphere--geometric mean */
+		} else if (pj_param(NULL,pl, "bR_g").i) { /* sphere--geometric mean */
 			*a = sqrt(*a * b);
 			*es = 0.;
-		} else if (pj_param(pl, "bR_h").i) { /* sphere--harmonic mean */
+		} else if (pj_param(NULL,pl, "bR_h").i) { /* sphere--harmonic mean */
 			*a = 2. * *a * b / (*a + b);
 			*es = 0.;
-		} else if ((i = pj_param(pl, "tR_lat_a").i) || /* sphere--arith. */
-			pj_param(pl, "tR_lat_g").i) { /* or geom. mean at latitude */
+		} else if ((i = pj_param(NULL,pl, "tR_lat_a").i) || /* sphere--arith. */
+			pj_param(NULL,pl, "tR_lat_g").i) { /* or geom. mean at latitude */
 			double tmp;
 
-			tmp = sin(pj_param(pl, i ? "rR_lat_a" : "rR_lat_g").f);
+			tmp = sin(pj_param(NULL,pl, i ? "rR_lat_a" : "rR_lat_g").f);
 			if (fabs(tmp) > HALFPI) {
                                 pj_ctx_set_errno(ctx,-11);
 				goto bomb;
diff --git a/src/pj_factors.c b/src/pj_factors.c
index 2c478a81..6a4db92e 100644
--- a/src/pj_factors.c
+++ b/src/pj_factors.c
@@ -13,7 +13,7 @@ pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) {
 
 	/* check for forward and latitude or longitude overange */
 	if ((t = fabs(lp.phi)-HALFPI) > EPS || fabs(lp.lam) > 10.) {
-                pj_ctx_set_errno( P->ctx->last_errno, -14);
+                pj_ctx_set_errno( P->ctx, -14);
 		return 1;
 	} else { /* proceed */
 		errno = pj_errno = 0;
@@ -72,7 +72,7 @@ pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) {
 		fac->s = (fac->der.y_p * fac->der.x_l - fac->der.x_p * fac->der.y_l) *
 			r / cosphi;
 		/* meridian-parallel angle theta prime */
-		fac->thetap = aasin(fac->s / (fac->h * fac->k));
+		fac->thetap = aasin(P->ctx,fac->s / (fac->h * fac->k));
 		/* Tissot ellips axis */
 		t = fac->k * fac->k + fac->h * fac->h;
 		fac->a = sqrt(t + 2. * fac->s);
@@ -80,7 +80,7 @@ pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) {
 		fac->b = 0.5 * (fac->a - t);
 		fac->a = 0.5 * (fac->a + t);
 		/* omega */
-		fac->omega = 2. * aasin((fac->a - fac->b)/(fac->a + fac->b));
+		fac->omega = 2. * aasin(P->ctx,(fac->a - fac->b)/(fac->a + fac->b));
 	}
 	return 0;
 }
diff --git a/src/pj_gauss.c b/src/pj_gauss.c
index bc25f444..c852c097 100644
--- a/src/pj_gauss.c
+++ b/src/pj_gauss.c
@@ -64,7 +64,7 @@ pj_gauss_ini(double e, double phi0, double *chi, double *rc) {
 	return ((void *)en);
 }
 	LP
-pj_gauss(LP elp, const void *en) {
+pj_gauss(projCtx ctx, LP elp, const void *en) {
 	LP slp;
 
 	slp.phi = 2. * atan( EN->K *
@@ -74,7 +74,7 @@ pj_gauss(LP elp, const void *en) {
 	return(slp);
 }
 	LP
-pj_inv_gauss(LP slp, const void *en) {
+pj_inv_gauss(projCtx ctx, LP slp, const void *en) {
 	LP elp;
 	double num;
 	int i;
@@ -89,6 +89,6 @@ pj_inv_gauss(LP slp, const void *en) {
 	}	
 	/* convergence failed */
 	if (!i)
-		pj_errno = -17;
+		pj_ctx_set_errno( ctx, -17 );
 	return (elp);
 }
diff --git a/src/pj_gridinfo.c b/src/pj_gridinfo.c
index e543d028..3bc11ed6 100644
--- a/src/pj_gridinfo.c
+++ b/src/pj_gridinfo.c
@@ -760,7 +760,6 @@ PJ_GRIDINFO *pj_gridinfo_init( projCtx ctx, const char *gridname )
 /* -------------------------------------------------------------------- */
     strcpy(fname, gridname);
     if (!(fp = pj_open_lib(ctx, fname, "rb"))) {
-        pj_errno = errno;
         return gilist;
     }
 
diff --git a/src/pj_gridlist.c b/src/pj_gridlist.c
index 3277bb5c..9f48146f 100644
--- a/src/pj_gridlist.c
+++ b/src/pj_gridlist.c
@@ -194,7 +194,7 @@ PJ_GRIDINFO **pj_gridlist_from_nadgrids( projCtx ctx, const char *nadgrids,
 
         if( end_char > sizeof(name) )
         {
-            pj_errno = -38;
+            pj_ctx_set_errno( ctx, -38 );
             pj_release_lock();
             return NULL;
         }
@@ -210,7 +210,7 @@ PJ_GRIDINFO **pj_gridlist_from_nadgrids( projCtx ctx, const char *nadgrids,
                                          &grid_max) 
             && required )
         {
-            pj_errno = -38;
+            pj_ctx_set_errno( ctx, -38 );
             pj_release_lock();
             return NULL;
         }
diff --git a/src/pj_init.c b/src/pj_init.c
index 073807d6..84f818f1 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -42,7 +42,7 @@ PJ_CVSID("$Id$");
 /*                              get_opt()                               */
 /************************************************************************/
 static paralist *
-get_opt(paralist **start, FILE *fid, char *name, paralist *next) {
+get_opt(projCtx ctx, paralist **start, FILE *fid, char *name, paralist *next) {
     char sword[302], *word = sword+1;
     int first = 1, len, c;
 
@@ -59,16 +59,16 @@ get_opt(paralist **start, FILE *fid, char *name, paralist *next) {
                 while((c = fgetc(fid)) != EOF && c != '\n') ;
                 break;
             }
-        } else if (!first && !pj_param(*start, sword).i) {
+        } else if (!first && !pj_param(ctx, *start, sword).i) {
             /* don't default ellipse if datum, ellps or any earth model
                information is set. */
             if( strncmp(word,"ellps=",6) != 0 
-                || (!pj_param(*start, "tdatum").i 
-                    && !pj_param(*start, "tellps").i 
-                    && !pj_param(*start, "ta").i 
-                    && !pj_param(*start, "tb").i 
-                    && !pj_param(*start, "trf").i 
-                    && !pj_param(*start, "tf").i) )
+                || (!pj_param(ctx, *start, "tdatum").i 
+                    && !pj_param(ctx, *start, "tellps").i 
+                    && !pj_param(ctx, *start, "ta").i 
+                    && !pj_param(ctx, *start, "tb").i 
+                    && !pj_param(ctx, *start, "trf").i 
+                    && !pj_param(ctx, *start, "tf").i) )
             {
                 next = next->next = pj_mkparam(word);
             }
@@ -88,9 +88,9 @@ get_defaults(projCtx ctx, paralist **start, paralist *next, char *name) {
     FILE *fid;
 
     if ( (fid = pj_open_lib(ctx,"proj_def.dat", "rt")) != NULL) {
-        next = get_opt(start, fid, "general", next);
+        next = get_opt(ctx, start, fid, "general", next);
         rewind(fid);
-        next = get_opt(start, fid, name, next);
+        next = get_opt(ctx, start, fid, name, next);
         (void)fclose(fid);
     }
     if (errno)
@@ -131,7 +131,7 @@ get_init(projCtx ctx, paralist **start, paralist *next, char *name) {
     else { pj_ctx_set_errno(ctx,-3); return NULL; }
 
     if ( (fid = pj_open_lib(ctx,fname, "rt")) != NULL)
-        next = get_opt(start, fid, opt, next);
+        next = get_opt(ctx, start, fid, opt, next);
     else
         return NULL;
     (void)fclose(fid);
@@ -254,22 +254,22 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
     if (ctx->last_errno) goto bum_call;
 
     /* check if +init present */
-    if (pj_param(start, "tinit").i) {
+    if (pj_param(ctx, start, "tinit").i) {
         paralist *last = curr;
 
-        if (!(curr = get_init(ctx,&start, curr, pj_param(start, "sinit").s)))
+        if (!(curr = get_init(ctx,&start, curr, pj_param(ctx, start, "sinit").s)))
             goto bum_call;
         if (curr == last) { pj_ctx_set_errno( ctx, -2); goto bum_call; }
     }
 
     /* find projection selection */
-    if (!(name = pj_param(start, "sproj").s))
+    if (!(name = pj_param(ctx, start, "sproj").s))
     { pj_ctx_set_errno( ctx, -4 ); goto bum_call; }
     for (i = 0; (s = pj_list[i].id) && strcmp(name, s) ; ++i) ;
     if (!s) { pj_ctx_set_errno( ctx, -5 ); goto bum_call; }
 
     /* set defaults, unless inhibited */
-    if (!pj_param(start, "bno_defs").i)
+    if (!pj_param(ctx, start, "bno_defs").i)
         curr = get_defaults(ctx,&start, curr, name);
     proj = (PJ *(*)(PJ *)) pj_list[i].proj;
 
@@ -316,24 +316,24 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
     }
         
     /* set PIN->geoc coordinate system */
-    PIN->geoc = (PIN->es && pj_param(start, "bgeoc").i);
+    PIN->geoc = (PIN->es && pj_param(ctx, start, "bgeoc").i);
 
     /* over-ranging flag */
-    PIN->over = pj_param(start, "bover").i;
+    PIN->over = pj_param(ctx, start, "bover").i;
 
     /* longitude center for wrapping */
-    PIN->has_geoid_vgrids = pj_param(start, "tgeoidgrids").i;
+    PIN->has_geoid_vgrids = pj_param(ctx, start, "tgeoidgrids").i;
 
     /* longitude center for wrapping */
-    PIN->is_long_wrap_set = pj_param(start, "tlon_wrap").i;
+    PIN->is_long_wrap_set = pj_param(ctx, start, "tlon_wrap").i;
     if (PIN->is_long_wrap_set)
-        PIN->long_wrap_center = pj_param(start, "rlon_wrap").f;
+        PIN->long_wrap_center = pj_param(ctx, start, "rlon_wrap").f;
 
     /* axis orientation */
-    if( (pj_param(start,"saxis").s) != NULL )
+    if( (pj_param(ctx, start,"saxis").s) != NULL )
     {
         static const char *axis_legal = "ewnsud";
-        const char *axis_arg = pj_param(start,"saxis").s;
+        const char *axis_arg = pj_param(ctx, start,"saxis").s;
         if( strlen(axis_arg) != 3 )
         {
             pj_ctx_set_errno( ctx, PJD_ERR_AXIS );
@@ -352,25 +352,25 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
         strcpy( PIN->axis, axis_arg );
     }
 
-    PIN->is_long_wrap_set = pj_param(start, "tlon_wrap").i;
+    PIN->is_long_wrap_set = pj_param(ctx, start, "tlon_wrap").i;
     if (PIN->is_long_wrap_set)
-        PIN->long_wrap_center = pj_param(start, "rlon_wrap").f;
+        PIN->long_wrap_center = pj_param(ctx, start, "rlon_wrap").f;
 
     /* central meridian */
-    PIN->lam0=pj_param(start, "rlon_0").f;
+    PIN->lam0=pj_param(ctx, start, "rlon_0").f;
 
     /* central latitude */
-    PIN->phi0 = pj_param(start, "rlat_0").f;
+    PIN->phi0 = pj_param(ctx, start, "rlat_0").f;
 
     /* false easting and northing */
-    PIN->x0 = pj_param(start, "dx_0").f;
-    PIN->y0 = pj_param(start, "dy_0").f;
+    PIN->x0 = pj_param(ctx, start, "dx_0").f;
+    PIN->y0 = pj_param(ctx, start, "dy_0").f;
 
     /* general scaling factor */
-    if (pj_param(start, "tk_0").i)
-        PIN->k0 = pj_param(start, "dk_0").f;
-    else if (pj_param(start, "tk").i)
-        PIN->k0 = pj_param(start, "dk").f;
+    if (pj_param(ctx, start, "tk_0").i)
+        PIN->k0 = pj_param(ctx, start, "dk_0").f;
+    else if (pj_param(ctx, start, "tk").i)
+        PIN->k0 = pj_param(ctx, start, "dk").f;
     else
         PIN->k0 = 1.;
     if (PIN->k0 <= 0.) {
@@ -380,12 +380,12 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
 
     /* set units */
     s = 0;
-    if ((name = pj_param(start, "sunits").s) != NULL) { 
+    if ((name = pj_param(ctx, start, "sunits").s) != NULL) { 
         for (i = 0; (s = pj_units[i].id) && strcmp(name, s) ; ++i) ;
         if (!s) { pj_ctx_set_errno( ctx, -7 ); goto bum_call; }
         s = pj_units[i].to_meter;
     }
-    if (s || (s = pj_param(start, "sto_meter").s)) {
+    if (s || (s = pj_param(ctx, start, "sto_meter").s)) {
         PIN->to_meter = strtod(s, &s);
         if (*s == '/') /* ratio number */
             PIN->to_meter /= strtod(++s, 0);
@@ -395,7 +395,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
 
     /* prime meridian */
     s = 0;
-    if ((name = pj_param(start, "spm").s) != NULL) { 
+    if ((name = pj_param(ctx, start, "spm").s) != NULL) { 
         const char *value = NULL;
         char *next_str = NULL;
 
@@ -409,12 +409,12 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
         }
             
         if( value == NULL 
-            && (dmstor(name,&next_str) != 0.0  || *name == '0')
+            && (dmstor_ctx(ctx,name,&next_str) != 0.0  || *name == '0')
             && *next_str == '\0' )
             value = name;
 
         if (!value) { pj_ctx_set_errno( ctx, -46 ); goto bum_call; }
-        PIN->from_greenwich = dmstor(value,NULL);
+        PIN->from_greenwich = dmstor_ctx(ctx,value,NULL);
     }
     else
         PIN->from_greenwich = 0.0;
diff --git a/src/pj_mlfn.c b/src/pj_mlfn.c
index d62bd485..aeab6117 100644
--- a/src/pj_mlfn.c
+++ b/src/pj_mlfn.c
@@ -40,7 +40,7 @@ pj_mlfn(double phi, double sphi, double cphi, double *en) {
 		+ sphi*(en[3] + sphi*en[4]))));
 }
 	double
-pj_inv_mlfn(double arg, double es, double *en) {
+pj_inv_mlfn(projCtx ctx, double arg, double es, double *en) {
 	double s, t, phi, k = 1./(1.-es);
 	int i;
 
@@ -52,6 +52,6 @@ pj_inv_mlfn(double arg, double es, double *en) {
 		if (fabs(t) < EPS)
 			return phi;
 	}
-	pj_errno = -17;
+	pj_ctx_set_errno( ctx, -17 );
 	return phi;
 }
diff --git a/src/pj_open_lib.c b/src/pj_open_lib.c
index b15cbdf6..22c180a8 100644
--- a/src/pj_open_lib.c
+++ b/src/pj_open_lib.c
@@ -166,6 +166,9 @@ pj_open_lib(projCtx ctx, char *name, char *mode) {
             name, sysname,
             fid == NULL ? "failed" : "succeeded" );
 
+    if( ctx->last_errno == 0 && errno != 0 )
+        pj_ctx_set_errno( ctx, errno );
+
     return(fid);
 #else
     return NULL;
diff --git a/src/pj_param.c b/src/pj_param.c
index c18554a0..119006e1 100644
--- a/src/pj_param.c
+++ b/src/pj_param.c
@@ -34,11 +34,15 @@ pj_mkparam(char *str) {
 /************************************************************************/
 
 	PVALUE /* test for presence or get parameter value */
-pj_param(paralist *pl, const char *opt) {
+pj_param(projCtx ctx, paralist *pl, const char *opt) {
+
 	int type;
 	unsigned l;
 	PVALUE value;
 
+	if( ctx == NULL )
+		ctx = pj_get_default_ctx();
+
 	type = *opt++;
 	/* simple linear lookup */
 	l = strlen(opt);
@@ -60,7 +64,7 @@ pj_param(paralist *pl, const char *opt) {
 			value.f = atof(opt);
 			break;
 		case 'r':	/* degrees input */
-			value.f = dmstor(opt, 0);
+			value.f = dmstor_ctx(ctx, opt, 0);
 			break;
 		case 's':	/* char string */
                         value.s = (char *) opt;
@@ -74,7 +78,7 @@ pj_param(paralist *pl, const char *opt) {
 				value.i = 1;
 				break;
 			default:
-				pj_errno = -8;
+				pj_ctx_set_errno(ctx, -8);
 				value.i = 0;
 				break;
 			}
diff --git a/src/pj_phi2.c b/src/pj_phi2.c
index 0900350e..5c81da48 100644
--- a/src/pj_phi2.c
+++ b/src/pj_phi2.c
@@ -6,7 +6,7 @@
 #define N_ITER 15
 
 	double
-pj_phi2(double ts, double e) {
+pj_phi2(projCtx ctx, double ts, double e) {
 	double eccnth, Phi, con, dphi;
 	int i;
 
@@ -20,6 +20,6 @@ pj_phi2(double ts, double e) {
 		Phi += dphi;
 	} while ( fabs(dphi) > TOL && --i);
 	if (i <= 0)
-		pj_errno = -18;
+		pj_ctx_set_errno( ctx, -18 );
 	return Phi;
 }
diff --git a/src/pj_transform.c b/src/pj_transform.c
index 60f6dd72..4f0585ea 100644
--- a/src/pj_transform.c
+++ b/src/pj_transform.c
@@ -470,8 +470,8 @@ int pj_compare_datums( PJ *srcdefn, PJ *dstdefn )
     }
     else if( srcdefn->datum_type == PJD_GRIDSHIFT )
     {
-        return strcmp( pj_param(srcdefn->params,"snadgrids").s,
-                       pj_param(dstdefn->params,"snadgrids").s ) == 0;
+        return strcmp( pj_param(srcdefn->ctx, srcdefn->params,"snadgrids").s,
+                       pj_param(dstdefn->ctx, dstdefn->params,"snadgrids").s ) == 0;
     }
     else
         return 1;
diff --git a/src/pj_utils.c b/src/pj_utils.c
index 89f5bda2..15adcd93 100644
--- a/src/pj_utils.c
+++ b/src/pj_utils.c
@@ -77,31 +77,31 @@ PJ *pj_latlong_from_proj( PJ *pj_in )
     pj_errno = 0;
     strcpy( defn, "+proj=latlong" );
 
-    if( pj_param(pj_in->params, "tdatum").i )
+    if( pj_param(pj_in->ctx, pj_in->params, "tdatum").i )
     {
         got_datum = TRUE;
         sprintf( defn+strlen(defn), " +datum=%s", 
-                 pj_param(pj_in->params,"sdatum").s );
+                 pj_param(pj_in->ctx, pj_in->params,"sdatum").s );
     }
-    else if( pj_param(pj_in->params, "tellps").i )
+    else if( pj_param(pj_in->ctx, pj_in->params, "tellps").i )
     {
         sprintf( defn+strlen(defn), " +ellps=%s", 
-                 pj_param(pj_in->params,"sellps").s );
+                 pj_param(pj_in->ctx, pj_in->params,"sellps").s );
     }
-    else if( pj_param(pj_in->params, "ta").i )
+    else if( pj_param(pj_in->ctx,pj_in->params, "ta").i )
     {
         sprintf( defn+strlen(defn), " +a=%s", 
-                 pj_param(pj_in->params,"sa").s );
+                 pj_param(pj_in->ctx,pj_in->params,"sa").s );
             
-        if( pj_param(pj_in->params, "tb").i )
+        if( pj_param(pj_in->ctx,pj_in->params, "tb").i )
             sprintf( defn+strlen(defn), " +b=%s", 
-                     pj_param(pj_in->params,"sb").s );
-        else if( pj_param(pj_in->params, "tes").i )
+                     pj_param(pj_in->ctx,pj_in->params,"sb").s );
+        else if( pj_param(pj_in->ctx,pj_in->params, "tes").i )
             sprintf( defn+strlen(defn), " +es=%s", 
-                     pj_param(pj_in->params,"ses").s );
-        else if( pj_param(pj_in->params, "tf").i )
+                     pj_param(pj_in->ctx,pj_in->params,"ses").s );
+        else if( pj_param(pj_in->ctx,pj_in->params, "tf").i )
             sprintf( defn+strlen(defn), " +f=%s", 
-                     pj_param(pj_in->params,"sf").s );
+                     pj_param(pj_in->ctx,pj_in->params,"sf").s );
         else
             sprintf( defn+strlen(defn), " +es=%.16g", 
                      pj_in->es );
@@ -115,41 +115,41 @@ PJ *pj_latlong_from_proj( PJ *pj_in )
 
     if( !got_datum )
     {
-        if( pj_param(pj_in->params, "ttowgs84").i )
+        if( pj_param(pj_in->ctx,pj_in->params, "ttowgs84").i )
             sprintf( defn+strlen(defn), " +towgs84=%s", 
-                     pj_param(pj_in->params,"stowgs84").s );
+                     pj_param(pj_in->ctx,pj_in->params,"stowgs84").s );
 
-        if( pj_param(pj_in->params, "tnadgrids").i )
+        if( pj_param(pj_in->ctx,pj_in->params, "tnadgrids").i )
             sprintf( defn+strlen(defn), " +nadgrids=%s", 
-                     pj_param(pj_in->params,"snadgrids").s );
+                     pj_param(pj_in->ctx,pj_in->params,"snadgrids").s );
     }
 
     /* copy over some other information related to ellipsoid */
-    if( pj_param(pj_in->params, "tR").i )
+    if( pj_param(pj_in->ctx,pj_in->params, "tR").i )
         sprintf( defn+strlen(defn), " +R=%s", 
-                 pj_param(pj_in->params,"sR").s );
+                 pj_param(pj_in->ctx,pj_in->params,"sR").s );
 
-    if( pj_param(pj_in->params, "tR_A").i )
+    if( pj_param(pj_in->ctx,pj_in->params, "tR_A").i )
         sprintf( defn+strlen(defn), " +R_A" );
 
-    if( pj_param(pj_in->params, "tR_V").i )
+    if( pj_param(pj_in->ctx,pj_in->params, "tR_V").i )
         sprintf( defn+strlen(defn), " +R_V" );
 
-    if( pj_param(pj_in->params, "tR_a").i )
+    if( pj_param(pj_in->ctx,pj_in->params, "tR_a").i )
         sprintf( defn+strlen(defn), " +R_a" );
 
-    if( pj_param(pj_in->params, "tR_lat_a").i )
+    if( pj_param(pj_in->ctx,pj_in->params, "tR_lat_a").i )
         sprintf( defn+strlen(defn), " +R_lat_a=%s", 
-                 pj_param(pj_in->params,"sR_lat_a").s );
+                 pj_param(pj_in->ctx,pj_in->params,"sR_lat_a").s );
 
-    if( pj_param(pj_in->params, "tR_lat_g").i )
+    if( pj_param(pj_in->ctx,pj_in->params, "tR_lat_g").i )
         sprintf( defn+strlen(defn), " +R_lat_g=%s", 
-                 pj_param(pj_in->params,"sR_lat_g").s );
+                 pj_param(pj_in->ctx,pj_in->params,"sR_lat_g").s );
 
     /* copy over prime meridian */
-    if( pj_param(pj_in->params, "tpm").i )
+    if( pj_param(pj_in->ctx,pj_in->params, "tpm").i )
         sprintf( defn+strlen(defn), " +pm=%s", 
-                 pj_param(pj_in->params,"spm").s );
+                 pj_param(pj_in->ctx,pj_in->params,"spm").s );
 
     return pj_init_plus_ctx( pj_in->ctx, defn );
 }
diff --git a/src/proj_mdist.c b/src/proj_mdist.c
index 9e14a37e..9c588a27 100644
--- a/src/proj_mdist.c
+++ b/src/proj_mdist.c
@@ -102,7 +102,7 @@ proj_mdist(double phi, double sphi, double cphi, const void *b) {
 	return(D + sc * sum);
 }
 	double
-proj_inv_mdist(double dist, const void *b) {
+proj_inv_mdist(projCtx ctx, double dist, const void *b) {
 	double s, t, phi, k;
 	int i;
 
@@ -118,6 +118,6 @@ proj_inv_mdist(double dist, const void *b) {
 			return phi;
 	}
 		/* convergence failed */
-	pj_errno = -17;
+	pj_ctx_set_errno(ctx, -17);
 	return phi;
 }
diff --git a/src/proj_rouss.c b/src/proj_rouss.c
index 637fe36f..64075abe 100644
--- a/src/proj_rouss.c
+++ b/src/proj_rouss.c
@@ -61,7 +61,7 @@ INVERSE(e_inverse); /* ellipsoid */
 	s = P->s0 + y*(1.+y2*(-P->D2+P->D8*y2))+
 		x2*(-P->D1+y*(-P->D3+y*(-P->D5+y*(-P->D7+y*P->D11)))+
 		x2*(P->D4+y*(P->D6+y*P->D10)-x2*P->D9));
-	lp.phi=proj_inv_mdist(s, P->en);
+	lp.phi=proj_inv_mdist(P->ctx, s, P->en);
 	s = sin(lp.phi);
 	lp.lam=al * sqrt(1. - P->es * s * s)/cos(lp.phi);
 	return (lp);
diff --git a/src/projects.h b/src/projects.h
index b7d69bfb..9535ffe7 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -352,11 +352,12 @@ typedef struct _pj_gi {
 
 /* procedure prototypes */
 double dmstor(const char *, char **);
+double dmstor_ctx(projCtx ctx, const char *, char **);
 void set_rtodms(int, int);
 char *rtodms(char *, double, int, int);
 double adjlon(double);
-double aacos(double), aasin(double), asqrt(double), aatan2(double, double);
-PVALUE pj_param(paralist *, const char *);
+double aacos(projCtx,double), aasin(projCtx,double), asqrt(double), aatan2(double, double);
+PVALUE pj_param(projCtx ctx, paralist *, const char *);
 paralist *pj_mkparam(char *);
 int pj_ell_set(projCtx ctx, paralist *, double *, double *);
 int pj_datum_set(projCtx,paralist *, PJ *);
@@ -370,11 +371,11 @@ void pj_insert_initcache( const char *filekey, const paralist *list);
 
 double *pj_enfn(double);
 double pj_mlfn(double, double, double, double *);
-double pj_inv_mlfn(double, double, double *);
+double pj_inv_mlfn(projCtx, double, double, double *);
 double pj_qsfn(double, double, double);
 double pj_tsfn(double, double, double);
 double pj_msfn(double, double, double);
-double pj_phi2(double, double);
+double pj_phi2(projCtx, double, double);
 double pj_qsfn_(double, PJ *);
 double *pj_authset(double);
 double pj_authlat(double, double *);
@@ -434,10 +435,10 @@ void pj_gridinfo_free( projCtx, PJ_GRIDINFO * );
 
 void *proj_mdist_ini(double);
 double proj_mdist(double, double, double, const void *);
-double proj_inv_mdist(double, const void *);
+double proj_inv_mdist(projCtx ctx, double, const void *);
 void *pj_gauss_ini(double, double, double *,double *);
-LP pj_gauss(LP, const void *);
-LP pj_inv_gauss(LP, const void *);
+LP pj_gauss(projCtx, LP, const void *);
+LP pj_inv_gauss(projCtx, LP, const void *);
 
 extern char const pj_release[];