--- ray/src/common/bsdf_m.c	2011/02/24 20:14:26	3.8
+++ ray/src/common/bsdf_m.c	2014/03/21 17:49:53	3.28
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdf_m.c,v 3.8 2011/02/24 20:14:26 greg Exp $";
+static const char RCSid[] = "$Id: bsdf_m.c,v 3.28 2014/03/21 17:49:53 greg Exp $";
 #endif
 /*
  *  bsdf_m.c
@@ -11,6 +11,7 @@ static const char RCSid[] = "$Id: bsdf_m.c,v 3.8 2011/
  *
  */
 
+#define	_USE_MATH_DEFINES
 #include "rtio.h"
 #include <stdlib.h>
 #include <math.h>
@@ -28,24 +29,10 @@ static const char RCSid[] = "$Id: bsdf_m.c,v 3.8 2011/
 #define RC_INTERR	(-4)
 #define RC_MEMERR	(-5)
 
-#define MAXLATS		46		/* maximum number of latitudes */
-
-/* BSDF angle specification */
-typedef struct {
-	char	name[64];		/* basis name */
-	int	nangles;		/* total number of directions */
-	struct {
-		float	tmin;			/* starting theta */
-		int	nphis;			/* number of phis (0 term) */
-	} lat[MAXLATS+1];		/* latitudes */
-} ANGLE_BASIS;
-
-#define	MAXABASES	7		/* limit on defined bases */
-
-static ANGLE_BASIS	abase_list[MAXABASES] = {
+ANGLE_BASIS	abase_list[MAXABASES] = {
 	{
 		"LBNL/Klems Full", 145,
-		{ {-5., 1},
+		{ {0., 1},
 		{5., 8},
 		{15., 16},
 		{25., 20},
@@ -57,7 +44,7 @@ static ANGLE_BASIS	abase_list[MAXABASES] = {
 		{90., 0} }
 	}, {
 		"LBNL/Klems Half", 73,
-		{ {-6.5, 1},
+		{ {0., 1},
 		{6.5, 8},
 		{19.5, 12},
 		{32.5, 16},
@@ -67,7 +54,7 @@ static ANGLE_BASIS	abase_list[MAXABASES] = {
 		{90., 0} }
 	}, {
 		"LBNL/Klems Quarter", 41,
-		{ {-9., 1},
+		{ {0., 1},
 		{9., 8},
 		{27., 12},
 		{46., 12},
@@ -76,29 +63,17 @@ static ANGLE_BASIS	abase_list[MAXABASES] = {
 	}
 };
 
-static int	nabases = 3;	/* current number of defined bases */
+int		nabases = 3;		/* current number of defined bases */
 
 static int
 fequal(double a, double b)
 {
-	if (b != .0)
+	if (b != 0)
 		a = a/b - 1.;
 	return (a <= 1e-6) & (a >= -1e-6);
 }
 
-/* returns the name of the given tag */
-#ifdef ezxml_name
-#undef ezxml_name
-static char *
-ezxml_name(ezxml_t xml)
-{
-	if (xml == NULL)
-		return NULL;
-	return xml->name;
-}
-#endif
-
-/* returns the given tag's character content or empty string if none */
+/* Returns the given tag's character content or empty string if none */
 #ifdef ezxml_txt
 #undef ezxml_txt
 static char *
@@ -157,16 +132,18 @@ SDnewMatrix(int ni, int no)
 /* Free a BSDF matrix */
 #define	SDfreeMatrix		free
 
-/* get vector for this angle basis index (front exiting) */
-static int
-fo_getvec(FVECT v, int ndx, double randX, void *p)
+/* Get vector for this angle basis index (front exiting) */
+int
+fo_getvec(FVECT v, double ndxr, void *p)
 {
-	ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
+	ANGLE_BASIS	*ab = (ANGLE_BASIS *)p;
+	int		ndx = (int)ndxr;
+	double		randX = ndxr - ndx;
 	double		rx[2];
 	int		li;
 	double		pol, azi, d;
 	
-	if ((ndx < 0) | (ndx >= ab->nangles))
+	if ((ndxr < 0) | (ndx >= ab->nangles))
 		return RC_FAIL;
 	for (li = 0; ndx >= ab->lat[li].nphis; li++)
 		ndx -= ab->lat[li].nphis;
@@ -181,19 +158,19 @@ fo_getvec(FVECT v, int ndx, double randX, void *p)
 	return RC_GOOD;
 }
 
-/* get index corresponding to the given vector (front exiting) */
-static int
+/* Get index corresponding to the given vector (front exiting) */
+int
 fo_getndx(const FVECT v, void *p)
 {
-	ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
+	ANGLE_BASIS	*ab = (ANGLE_BASIS *)p;
 	int	li, ndx;
-	double	pol, azi, d;
+	double	pol, azi;
 
 	if (v == NULL)
 		return -1;
-	if ((v[2] < .0) | (v[2] > 1.0))
+	if ((v[2] < 0) | (v[2] > 1.))
 		return -1;
-	pol = 180.0/M_PI*acos(v[2]);
+	pol = 180.0/M_PI*Acos(v[2]);
 	azi = 180.0/M_PI*atan2(v[1], v[0]);
 	if (azi < 0.0) azi += 360.0;
 	for (li = 1; ab->lat[li].tmin <= pol; li++)
@@ -210,8 +187,8 @@ fo_getndx(const FVECT v, void *p)
 /* compute square of real value */
 static double sq(double x) { return x*x; }
 
-/* get projected solid angle for this angle basis index (universal) */
-static double
+/* Get projected solid angle for this angle basis index (universal) */
+double
 io_getohm(int ndx, void *p)
 {
 	static int	last_li = -1;
@@ -227,19 +204,17 @@ io_getohm(int ndx, void *p)
 	if (li == last_li)			/* cached latitude? */
 		return last_ohm;
 	last_li = li;
-	theta1 = M_PI/180. * ab->lat[li+1].tmin;
-	if (ab->lat[li].nphis == 1)		/* special case */
-		return last_ohm = M_PI*(1. - sq(cos(theta1)));
 	theta = M_PI/180. * ab->lat[li].tmin;
+	theta1 = M_PI/180. * ab->lat[li+1].tmin;
 	return last_ohm = M_PI*(sq(cos(theta)) - sq(cos(theta1))) /
 				(double)ab->lat[li].nphis;
 }
 
-/* get vector for this angle basis index (back incident) */
-static int
-bi_getvec(FVECT v, int ndx, double randX, void *p)
+/* Get vector for this angle basis index (back incident) */
+int
+bi_getvec(FVECT v, double ndxr, void *p)
 {
-	if (!fo_getvec(v, ndx, randX, p))
+	if (!fo_getvec(v, ndxr, p))
 		return RC_FAIL;
 
 	v[0] = -v[0];
@@ -249,8 +224,8 @@ bi_getvec(FVECT v, int ndx, double randX, void *p)
 	return RC_GOOD;
 }
 
-/* get index corresponding to the vector (back incident) */
-static int
+/* Get index corresponding to the vector (back incident) */
+int
 bi_getndx(const FVECT v, void *p)
 {
 	FVECT  v2;
@@ -262,11 +237,11 @@ bi_getndx(const FVECT v, void *p)
 	return fo_getndx(v2, p);
 }
 
-/* get vector for this angle basis index (back exiting) */
-static int
-bo_getvec(FVECT v, int ndx, double randX, void *p)
+/* Get vector for this angle basis index (back exiting) */
+int
+bo_getvec(FVECT v, double ndxr, void *p)
 {
-	if (!fo_getvec(v, ndx, randX, p))
+	if (!fo_getvec(v, ndxr, p))
 		return RC_FAIL;
 
 	v[2] = -v[2];
@@ -274,8 +249,8 @@ bo_getvec(FVECT v, int ndx, double randX, void *p)
 	return RC_GOOD;
 }
 
-/* get index corresponding to the vector (back exiting) */
-static int
+/* Get index corresponding to the vector (back exiting) */
+int
 bo_getndx(const FVECT v, void *p)
 {
 	FVECT  v2;
@@ -287,11 +262,11 @@ bo_getndx(const FVECT v, void *p)
 	return fo_getndx(v2, p);
 }
 
-/* get vector for this angle basis index (front incident) */
-static int
-fi_getvec(FVECT v, int ndx, double randX, void *p)
+/* Get vector for this angle basis index (front incident) */
+int
+fi_getvec(FVECT v, double ndxr, void *p)
 {
-	if (!fo_getvec(v, ndx, randX, p))
+	if (!fo_getvec(v, ndxr, p))
 		return RC_FAIL;
 
 	v[0] = -v[0];
@@ -300,8 +275,8 @@ fi_getvec(FVECT v, int ndx, double randX, void *p)
 	return RC_GOOD;
 }
 
-/* get index corresponding to the vector (front incident) */
-static int
+/* Get index corresponding to the vector (front incident) */
+int
 fi_getndx(const FVECT v, void *p)
 {
 	FVECT  v2;
@@ -344,13 +319,12 @@ load_angle_basis(ezxml_t wab)
 				ezxml_child(ezxml_child(wbb,
 					"ThetaBounds"), "UpperTheta")));
 		if (!i)
-			abase_list[nabases].lat[i].tmin =
-					-abase_list[nabases].lat[i+1].tmin;
+			abase_list[nabases].lat[0].tmin = 0;
 		else if (!fequal(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
 					"ThetaBounds"), "LowerTheta"))),
 				abase_list[nabases].lat[i].tmin)) {
 			sprintf(SDerrorDetail, "Theta values disagree in '%s'",
-				abname);
+					abname);
 			return RC_DATERR;
 		}
 		abase_list[nabases].nangles +=
@@ -360,7 +334,7 @@ load_angle_basis(ezxml_t wab)
 				(abase_list[nabases].lat[i].nphis == 1 &&
 				abase_list[nabases].lat[i].tmin > FTINY)) {
 			sprintf(SDerrorDetail, "Illegal phi count in '%s'",
-				abname);
+					abname);
 			return RC_DATERR;
 		}
 	}
@@ -411,30 +385,35 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
 	SDSpectralDF	*df;
 	SDMat		*dp;
 	char		*sdata;
-	int		tfront;
 	int		inbi, outbi;
 	int		i;
 					/* allocate BSDF component */
 	sdata = ezxml_txt(ezxml_child(wdb, "WavelengthDataDirection"));
+	if (!sdata)
+		return RC_FAIL;
 	/*
 	 * Remember that front and back are reversed from WINDOW 6 orientations
 	 */
-	if ((tfront = !strcasecmp(sdata, "Transmission Back")) ||
-			(sd->tf == NULL &&
-				!strcasecmp(sdata, "Transmission Front"))) {
+	if (!strcasecmp(sdata, "Transmission Front")) {
+		if (sd->tb != NULL)
+			SDfreeSpectralDF(sd->tb);
+		if ((sd->tb = SDnewSpectralDF(1)) == NULL)
+			return RC_MEMERR;
+		df = sd->tb;
+	} else if (!strcasecmp(sdata, "Transmission Back")) {
 		if (sd->tf != NULL)
 			SDfreeSpectralDF(sd->tf);
 		if ((sd->tf = SDnewSpectralDF(1)) == NULL)
 			return RC_MEMERR;
 		df = sd->tf;
 	} else if (!strcasecmp(sdata, "Reflection Front")) {
-		if (sd->rb != NULL)	/* note back-front reversal */
+		if (sd->rb != NULL)
 			SDfreeSpectralDF(sd->rb);
 		if ((sd->rb = SDnewSpectralDF(1)) == NULL)
 			return RC_MEMERR;
 		df = sd->rb;
 	} else if (!strcasecmp(sdata, "Reflection Back")) {
-		if (sd->rf != NULL)	/* note front-back reversal */
+		if (sd->rf != NULL)
 			SDfreeSpectralDF(sd->rf);
 		if ((sd->rf = SDnewSpectralDF(1)) == NULL)
 			return RC_MEMERR;
@@ -476,17 +455,15 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
 	dp->ib_priv = &abase_list[inbi];
 	dp->ob_priv = &abase_list[outbi];
 	if (df == sd->tf) {
-		if (tfront) {
-			dp->ib_vec = &fi_getvec;
-			dp->ib_ndx = &fi_getndx;
-			dp->ob_vec = &bo_getvec;
-			dp->ob_ndx = &bo_getndx;
-		} else {
-			dp->ib_vec = &bi_getvec;
-			dp->ib_ndx = &bi_getndx;
-			dp->ob_vec = &fo_getvec;
-			dp->ob_ndx = &fo_getndx;
-		}
+		dp->ib_vec = &fi_getvec;
+		dp->ib_ndx = &fi_getndx;
+		dp->ob_vec = &bo_getvec;
+		dp->ob_ndx = &bo_getndx;
+	} else if (df == sd->tb) {
+		dp->ib_vec = &bi_getvec;
+		dp->ib_ndx = &bi_getndx;
+		dp->ob_vec = &fo_getvec;
+		dp->ob_ndx = &fo_getndx;
 	} else if (df == sd->rf) {
 		dp->ib_vec = &fi_getvec;
 		dp->ib_ndx = &fi_getndx;
@@ -504,29 +481,32 @@ load_bsdf_data(SDData *sd, ezxml_t wdb, int rowinc)
 	df->comp[0].dist = dp;
 	df->comp[0].func = &SDhandleMtx;
 					/* read BSDF data */
-	sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
+	sdata = ezxml_txt(ezxml_child(wdb, "ScatteringData"));
 	if (!sdata || !*sdata) {
 		sprintf(SDerrorDetail, "Missing BSDF ScatteringData in '%s'",
 				sd->name);
 		return RC_FORMERR;
 	}
 	for (i = 0; i < dp->ninc*dp->nout; i++) {
-		char  *sdnext = fskip(sdata);
+		char	*sdnext = fskip(sdata);
+		double	val;
 		if (sdnext == NULL) {
 			sprintf(SDerrorDetail,
 				"Bad/missing BSDF ScatteringData in '%s'",
 					sd->name);
 			return RC_FORMERR;
 		}
-		while (*sdnext && isspace(*sdnext))
+		while (isspace(*sdnext))
 			sdnext++;
 		if (*sdnext == ',') sdnext++;
+		if ((val = atof(sdata)) < 0)
+			val = 0;	/* don't allow negative values */
 		if (rowinc) {
 			int	r = i/dp->nout;
-			int	c = i - c*dp->nout;
-			mBSDF_value(dp,r,c) = atof(sdata);
+			int	c = i - r*dp->nout;
+			mBSDF_value(dp,r,c) = val;
 		} else
-			dp->bsdf[i] = atof(sdata);
+			dp->bsdf[i] = val;
 		sdata = sdnext;
 	}
 	return get_extrema(df);
@@ -543,6 +523,10 @@ subtract_min(SDMat *sm)
 	for (i = n; --i; )
 		if (sm->bsdf[i] < minv)
 			minv = sm->bsdf[i];
+	
+	if (minv <= FTINY)
+		return .0;
+
 	for (i = n; i--; )
 		sm->bsdf[i] -= minv;
 
@@ -568,7 +552,7 @@ extract_diffuse(SDValue *dv, SDSpectralDF *df)
 		c_cmix(&dv->spec, dv->cieY, &dv->spec, ymin, &df->comp[n].cspec[0]);
 		dv->cieY += ymin;
 	}
-	df->maxHemi -= dv->cieY;	/* adjust minimum hemispherical */
+	df->maxHemi -= dv->cieY;	/* adjust maximum hemispherical */
 					/* make sure everything is set */
 	c_ccvt(&dv->spec, C_CSXY+C_CSSPEC);
 }
@@ -577,12 +561,11 @@ extract_diffuse(SDValue *dv, SDSpectralDF *df)
 SDError
 SDloadMtx(SDData *sd, ezxml_t wtl)
 {
-	ezxml_t			wld, wdb;
-	int			rowIn;
-	struct BSDF_data	*dp;
-	char			*txt;
-	int			rval;
-	
+	ezxml_t		wld, wdb;
+	int		rowIn;
+	char		*txt;
+	int		rval;
+					/* basic checks and data ordering */
 	txt = ezxml_txt(ezxml_child(ezxml_child(wtl,
 			"DataDefinition"), "IncidentDataStructure"));
 	if (txt == NULL || !*txt) {
@@ -601,12 +584,14 @@ SDloadMtx(SDData *sd, ezxml_t wtl)
 				sd->name);
 		return SDEsupport;
 	}
-				/* get angle basis */
-	rval = load_angle_basis(ezxml_child(ezxml_child(wtl,
-				"DataDefinition"), "AngleBasis"));
-	if (rval < 0)
-		return convert_errcode(rval);
-				/* load BSDF components */
+					/* get angle bases */
+	for (wld = ezxml_child(ezxml_child(wtl, "DataDefinition"), "AngleBasis");
+				wld != NULL; wld = wld->next) {
+		rval = load_angle_basis(wld);
+		if (rval < 0)
+			return convert_errcode(rval);
+	}
+					/* load BSDF components */
 	for (wld = ezxml_child(wtl, "WavelengthData");
 				wld != NULL; wld = wld->next) {
 		if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")),
@@ -617,21 +602,28 @@ SDloadMtx(SDData *sd, ezxml_t wtl)
 			if ((rval = load_bsdf_data(sd, wdb, rowIn)) < 0)
 				return convert_errcode(rval);
 	}
-				/* separate diffuse components */
+					/* separate diffuse components */
 	extract_diffuse(&sd->rLambFront, sd->rf);
 	extract_diffuse(&sd->rLambBack, sd->rb);
-	extract_diffuse(&sd->tLamb, sd->tf);
-				/* return success */
+	if (sd->tf != NULL)
+		extract_diffuse(&sd->tLamb, sd->tf);
+	if (sd->tb != NULL)
+		extract_diffuse(&sd->tLamb, sd->tb);
+					/* return success */
 	return SDEnone;
 }
 
 /* Get Matrix BSDF value */
 static int
 SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec,
-				const FVECT inVec, const void *dist)
+				const FVECT inVec, SDComponent *sdc)
 {
-	const SDMat	*dp = (const SDMat *)dist;
+	const SDMat	*dp;
 	int		i_ndx, o_ndx;
+					/* check arguments */
+	if ((coef == NULL) | (outVec == NULL) | (inVec == NULL) | (sdc == NULL)
+			|| (dp = (SDMat *)sdc->dist) == NULL)
+		return 0;
 					/* get angle indices */
 	i_ndx = mBSDF_incndx(dp, inVec);
 	o_ndx = mBSDF_outndx(dp, outVec);
@@ -646,23 +638,28 @@ SDgetMtxBSDF(float coef[SDmaxCh], const FVECT outVec,
 	return 1;			/* XXX monochrome for now */
 }
 
-/* Query solid angle for vector */
+/* Query solid angle for vector(s) */
 static SDError
-SDqueryMtxProjSA(double *psa, const FVECT vec, int qflags, const void *dist)
+SDqueryMtxProjSA(double *psa, const FVECT v1, const RREAL *v2,
+					int qflags, SDComponent *sdc)
 {
-	const SDMat	*dp = (const SDMat *)dist;
+	const SDMat	*dp;
 	double		inc_psa, out_psa;
 					/* check arguments */
-	if ((psa == NULL) | (vec == NULL) | (dp == NULL))
+	if ((psa == NULL) | (v1 == NULL) | (sdc == NULL) ||
+			(dp = (SDMat *)sdc->dist) == NULL)
 		return SDEargument;
+	if (v2 == NULL)
+		v2 = v1;
 					/* get projected solid angles */
-	inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, vec));
-	out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, vec));
+	out_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v1));
+	inc_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v2));
+	if ((v1 != v2) & (out_psa <= 0) & (inc_psa <= 0)) {
+		inc_psa = mBSDF_outohm(dp, mBSDF_outndx(dp, v2));
+		out_psa = mBSDF_incohm(dp, mBSDF_incndx(dp, v1));
+	}
 
 	switch (qflags) {		/* record based on flag settings */
-	case SDqueryVal:
-		psa[0] = .0;
-		/* fall through */
 	case SDqueryMax:
 		if (inc_psa > psa[0])
 			psa[0] = inc_psa;
@@ -670,20 +667,24 @@ SDqueryMtxProjSA(double *psa, const FVECT vec, int qfl
 			psa[0] = out_psa;
 		break;
 	case SDqueryMin+SDqueryMax:
-		if (inc_psa > psa[0])
+		if (inc_psa > psa[1])
 			psa[1] = inc_psa;
-		if (out_psa > psa[0])
+		if (out_psa > psa[1])
 			psa[1] = out_psa;
 		/* fall through */
+	case SDqueryVal:
+		if (qflags == SDqueryVal)
+			psa[0] = M_PI;
+		/* fall through */
 	case SDqueryMin:
-		if ((inc_psa > .0) & (inc_psa < psa[0]))
+		if ((inc_psa > 0) & (inc_psa < psa[0]))
 			psa[0] = inc_psa;
-		if ((out_psa > .0) & (out_psa < psa[0]))
+		if ((out_psa > 0) & (out_psa < psa[0]))
 			psa[0] = out_psa;
 		break;
 	}
 					/* make sure it's legal */
-	return (psa[0] <= .0) ? SDEinternal : SDEnone;
+	return (psa[0] <= 0) ? SDEinternal : SDEnone;
 }
 
 /* Compute new cumulative distribution from BSDF */
@@ -721,12 +722,13 @@ make_cdist(SDMatCDst *cd, const FVECT inVec, SDMat *dp
 static const SDCDst *
 SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
 {
-	SDMat		*dp = (SDMat *)sdc->dist;
+	SDMat		*dp;
 	int		reverse;
 	SDMatCDst	myCD;
 	SDMatCDst	*cd, *cdlast;
 					/* check arguments */
-	if ((inVec == NULL) | (dp == NULL))
+	if ((inVec == NULL) | (sdc == NULL) ||
+			(dp = (SDMat *)sdc->dist) == NULL)
 		return NULL;
 	memset(&myCD, 0, sizeof(myCD));
 	myCD.indx = mBSDF_incndx(dp, inVec);
@@ -745,17 +747,15 @@ SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
 		reverse = 1;
 	}
 	cdlast = NULL;			/* check for it in cache list */
-	for (cd = (SDMatCDst *)sdc->cdList;
-				cd != NULL; cd = (SDMatCDst *)cd->next) {
+	for (cd = (SDMatCDst *)sdc->cdList; cd != NULL;
+					cdlast = cd, cd = cd->next)
 		if (cd->indx == myCD.indx && (cd->calen == myCD.calen) &
 					(cd->ob_priv == myCD.ob_priv) &
 					(cd->ob_vec == myCD.ob_vec))
 			break;
-		cdlast = cd;
-	}
 	if (cd == NULL) {		/* need to allocate new entry */
 		cd = (SDMatCDst *)malloc(sizeof(SDMatCDst) +
-					myCD.calen*sizeof(myCD.carr[0]));
+					sizeof(myCD.carr[0])*myCD.calen);
 		if (cd == NULL)
 			return NULL;
 		*cd = myCD;		/* compute cumulative distribution */
@@ -767,7 +767,7 @@ SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
 	}
 	if (cdlast != NULL) {		/* move entry to head of cache list */
 		cdlast->next = cd->next;
-		cd->next = sdc->cdList;
+		cd->next = (SDMatCDst *)sdc->cdList;
 		sdc->cdList = (SDCDst *)cd;
 	}
 	return (SDCDst *)cd;		/* ready to go */
@@ -775,19 +775,19 @@ SDgetMtxCDist(const FVECT inVec, SDComponent *sdc)
 
 /* Sample cumulative distribution */
 static SDError
-SDsampMtxCDist(FVECT outVec, double randX, const SDCDst *cdp)
+SDsampMtxCDist(FVECT ioVec, double randX, const SDCDst *cdp)
 {
 	const unsigned	maxval = ~0;
 	const SDMatCDst	*mcd = (const SDMatCDst *)cdp;
 	const unsigned	target = randX*maxval;
 	int		i, iupper, ilower;
 					/* check arguments */
-	if ((outVec == NULL) | (mcd == NULL))
+	if ((ioVec == NULL) | (mcd == NULL))
 		return SDEargument;
 					/* binary search to find index */
 	ilower = 0; iupper = mcd->calen;
 	while ((i = (iupper + ilower) >> 1) != ilower)
-		if ((long)target >= (long)mcd->carr[i])
+		if (target >= mcd->carr[i])
 			ilower = i;
 		else
 			iupper = i;
@@ -795,9 +795,9 @@ SDsampMtxCDist(FVECT outVec, double randX, const SDCDs
 	randX = (randX*maxval - mcd->carr[ilower]) /
 			(double)(mcd->carr[iupper] - mcd->carr[ilower]);
 					/* convert index to vector */
-	if ((*mcd->ob_vec)(outVec, i, randX, mcd->ob_priv))
+	if ((*mcd->ob_vec)(ioVec, i+randX, mcd->ob_priv))
 		return SDEnone;
-	strcpy(SDerrorDetail, "BSDF sampling fault");
+	strcpy(SDerrorDetail, "Matrix BSDF sampling fault");
 	return SDEinternal;
 }