src/common/bsdf.c

#ifndef lint
static const char RCSid[] = "$Id: bsdf.c,v 2.11 2011/01/06 04:39:08 greg Exp $";
#endif
/*
 * Routines for handling BSDF data
 */

#include "standard.h"
#include "bsdf.h"
#include "paths.h"
#include "ezxml.h"
#include <ctype.h>

#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 */
                short   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] = {
        {
                "LBNL/Klems Full", 145,
                { {-5., 1},
                {5., 8},
                {15., 16},
                {25., 20},
                {35., 24},
                {45., 24},
                {55., 24},
                {65., 16},
                {75., 12},
                {90., 0} }
        }, {
                "LBNL/Klems Half", 73,
                { {-6.5, 1},
                {6.5, 8},
                {19.5, 12},
                {32.5, 16},
                {46.5, 20},
                {61.5, 12},
                {76.5, 4},
                {90., 0} }
        }, {
                "LBNL/Klems Quarter", 41,
                { {-9., 1},
                {9., 8},
                {27., 12},
                {46., 12},
                {66., 8},
                {90., 0} }
        }
};

static int      nabases = 3;    /* current number of defined bases */

#define  FEQ(a,b)       ((a)-(b) <= 1e-6 && (b)-(a) <= 1e-6)

static int
fequal(double a, double b)
{
        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
#ifdef ezxml_txt
#undef ezxml_txt
static char *
ezxml_txt(ezxml_t xml)
{
        if (xml == NULL)
                return("");
        return(xml->txt);
}
#endif


static int
ab_getvec(              /* get vector for this angle basis index */
        FVECT v,
        int ndx,
        void *p
)
{
        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
        int     li;
        double  pol, azi, d;
        
        if ((ndx < 0) | (ndx >= ab->nangles))
                return(0);
        for (li = 0; ndx >= ab->lat[li].nphis; li++)
                ndx -= ab->lat[li].nphis;
        pol = PI/180.*0.5*(ab->lat[li].tmin + ab->lat[li+1].tmin);
        azi = 2.*PI*ndx/ab->lat[li].nphis;
        v[2] = d = cos(pol);
        d = sqrt(1. - d*d);     /* sin(pol) */
        v[0] = cos(azi)*d;
        v[1] = sin(azi)*d;
        return(1);
}


static int
ab_getndx(              /* get index corresponding to the given vector */
        FVECT v,
        void *p
)
{
        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
        int     li, ndx;
        double  pol, azi, d;

        if ((v[2] < -1.0) | (v[2] > 1.0))
                return(-1);
        pol = 180.0/PI*acos(v[2]);
        azi = 180.0/PI*atan2(v[1], v[0]);
        if (azi < 0.0) azi += 360.0;
        for (li = 1; ab->lat[li].tmin <= pol; li++)
                if (!ab->lat[li].nphis)
                        return(-1);
        --li;
        ndx = (int)((1./360.)*azi*ab->lat[li].nphis + 0.5);
        if (ndx >= ab->lat[li].nphis) ndx = 0;
        while (li--)
                ndx += ab->lat[li].nphis;
        return(ndx);
}


static double
ab_getohm(              /* get solid angle for this angle basis index */
        int ndx,
        void *p
)
{
        ANGLE_BASIS  *ab = (ANGLE_BASIS *)p;
        int     li;
        double  theta, theta1;
        
        if ((ndx < 0) | (ndx >= ab->nangles))
                return(0);
        for (li = 0; ndx >= ab->lat[li].nphis; li++)
                ndx -= ab->lat[li].nphis;
        theta1 = PI/180. * ab->lat[li+1].tmin;
        if (ab->lat[li].nphis == 1) {           /* special case */
                if (ab->lat[li].tmin > FTINY)
                        error(USER, "unsupported BSDF coordinate system");
                return(2.*PI*(1. - cos(theta1)));
        }
        theta = PI/180. * ab->lat[li].tmin;
        return(2.*PI*(cos(theta) - cos(theta1))/(double)ab->lat[li].nphis);
}


static int
ab_getvecR(             /* get reverse vector for this angle basis index */
        FVECT v,
        int ndx,
        void *p
)
{
        if (!ab_getvec(v, ndx, p))
                return(0);

        v[0] = -v[0];
        v[1] = -v[1];
        v[2] = -v[2];

        return(1);
}


static int
ab_getndxR(             /* get index corresponding to the reverse vector */
        FVECT v,
        void *p
)
{
        FVECT  v2;
        
        v2[0] = -v[0];
        v2[1] = -v[1];
        v2[2] = -v[2];

        return ab_getndx(v2, p);
}


static void
load_angle_basis(       /* load custom BSDF angle basis */
        ezxml_t wab
)
{
        char    *abname = ezxml_txt(ezxml_child(wab, "AngleBasisName"));
        ezxml_t wbb;
        int     i;
        
        if (!abname || !*abname)
                return;
        for (i = nabases; i--; )
                if (!strcasecmp(abname, abase_list[i].name))
                        return;         /* assume it's the same */
        if (nabases >= MAXABASES)
                error(INTERNAL, "too many angle bases");
        strcpy(abase_list[nabases].name, abname);
        abase_list[nabases].nangles = 0;
        for (i = 0, wbb = ezxml_child(wab, "AngleBasisBlock");
                        wbb != NULL; i++, wbb = wbb->next) {
                if (i >= MAXLATS)
                        error(INTERNAL, "too many latitudes in custom basis");
                abase_list[nabases].lat[i+1].tmin = atof(ezxml_txt(
                                ezxml_child(ezxml_child(wbb,
                                        "ThetaBounds"), "UpperTheta")));
                if (!i)
                        abase_list[nabases].lat[i].tmin =
                                        -abase_list[nabases].lat[i+1].tmin;
                else if (!fequal(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
                                        "ThetaBounds"), "LowerTheta"))),
                                abase_list[nabases].lat[i].tmin))
                        error(WARNING, "theta values disagree in custom basis");
                abase_list[nabases].nangles +=
                        abase_list[nabases].lat[i].nphis =
                                atoi(ezxml_txt(ezxml_child(wbb, "nPhis")));
        }
        abase_list[nabases++].lat[i].nphis = 0;
}


static double
to_meters(              /* return factor to convert given unit to meters */
        const char *unit
)
{
        if (unit == NULL) return(1.);           /* safe assumption? */
        if (!strcasecmp(unit, "Meter")) return(1.);
        if (!strcasecmp(unit, "Foot")) return(.3048);
        if (!strcasecmp(unit, "Inch")) return(.0254);
        if (!strcasecmp(unit, "Centimeter")) return(.01);
        if (!strcasecmp(unit, "Millimeter")) return(.001);
        sprintf(errmsg, "unknown dimensional unit '%s'", unit);
        error(USER, errmsg);
}


static void
load_geometry(          /* load geometric dimensions and description (if any) */
        struct BSDF_data *dp,
        ezxml_t wdb
)
{
        ezxml_t         geom;
        double          cfact;
        const char      *fmt, *mgfstr;

        dp->dim[0] = dp->dim[1] = dp->dim[2] = 0;
        dp->mgf = NULL;
        if ((geom = ezxml_child(wdb, "Width")) != NULL)
                dp->dim[0] = atof(ezxml_txt(geom)) *
                                to_meters(ezxml_attr(geom, "unit"));
        if ((geom = ezxml_child(wdb, "Height")) != NULL)
                dp->dim[1] = atof(ezxml_txt(geom)) *
                                to_meters(ezxml_attr(geom, "unit"));
        if ((geom = ezxml_child(wdb, "Thickness")) != NULL)
                dp->dim[2] = atof(ezxml_txt(geom)) *
                                to_meters(ezxml_attr(geom, "unit"));
        if ((geom = ezxml_child(wdb, "Geometry")) == NULL ||
                        (mgfstr = ezxml_txt(geom)) == NULL)
                return;
        if ((fmt = ezxml_attr(geom, "format")) != NULL &&
                        strcasecmp(fmt, "MGF")) {
                sprintf(errmsg, "unrecognized geometry format '%s'", fmt);
                error(WARNING, errmsg);
                return;
        }
        cfact = to_meters(ezxml_attr(geom, "unit"));
        dp->mgf = (char *)malloc(strlen(mgfstr)+32);
        if (dp->mgf == NULL)
                error(SYSTEM, "out of memory in load_geometry");
        if (cfact < 0.99 || cfact > 1.01)
                sprintf(dp->mgf, "xf -s %.5f\n%s\nxf\n", cfact, mgfstr);
        else
                strcpy(dp->mgf, mgfstr);
}


static void
load_bsdf_data(         /* load BSDF distribution for this wavelength */
        struct BSDF_data *dp,
        ezxml_t wdb
)
{
        char  *cbasis = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
        char  *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis"));
        char  *sdata;
        int  i;
        
        if ((!cbasis || !*cbasis) | (!rbasis || !*rbasis)) {
                error(WARNING, "missing column/row basis for BSDF");
                return;
        }
        for (i = nabases; i--; )
                if (!strcasecmp(cbasis, abase_list[i].name)) {
                        dp->ninc = abase_list[i].nangles;
                        dp->ib_priv = (void *)&abase_list[i];
                        dp->ib_vec = ab_getvecR;
                        dp->ib_ndx = ab_getndxR;
                        dp->ib_ohm = ab_getohm;
                        break;
                }
        if (i < 0) {
                sprintf(errmsg, "undefined ColumnAngleBasis '%s'", cbasis);
                error(WARNING, errmsg);
                return;
        }
        for (i = nabases; i--; )
                if (!strcasecmp(rbasis, abase_list[i].name)) {
                        dp->nout = abase_list[i].nangles;
                        dp->ob_priv = (void *)&abase_list[i];
                        dp->ob_vec = ab_getvec;
                        dp->ob_ndx = ab_getndx;
                        dp->ob_ohm = ab_getohm;
                        break;
                }
        if (i < 0) {
                sprintf(errmsg, "undefined RowAngleBasis '%s'", cbasis);
                error(WARNING, errmsg);
                return;
        }
                                /* read BSDF data */
        sdata  = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
        if (!sdata || !*sdata) {
                error(WARNING, "missing BSDF ScatteringData");
                return;
        }
        dp->bsdf = (float *)malloc(sizeof(float)*dp->ninc*dp->nout);
        if (dp->bsdf == NULL)
                error(SYSTEM, "out of memory in load_bsdf_data");
        for (i = 0; i < dp->ninc*dp->nout; i++) {
                char  *sdnext = fskip(sdata);
                if (sdnext == NULL) {
                        error(WARNING, "bad/missing BSDF ScatteringData");
                        free(dp->bsdf); dp->bsdf = NULL;
                        return;
                }
                while (*sdnext && isspace(*sdnext))
                        sdnext++;
                if (*sdnext == ',') sdnext++;
                dp->bsdf[i] = atof(sdata);
                sdata = sdnext;
        }
        while (isspace(*sdata))
                sdata++;
        if (*sdata) {
                sprintf(errmsg, "%d extra characters after BSDF ScatteringData",
                                (int)strlen(sdata));
                error(WARNING, errmsg);
        }
}


static int
check_bsdf_data(        /* check that BSDF data is sane */
        struct BSDF_data *dp
)
{
        double          *omega_iarr, *omega_oarr;
        double          dom, contrib, hemi_total, full_total;
        int             nneg;
        FVECT           v;
        int             i, o;

        if (dp == NULL || dp->bsdf == NULL)
                return(0);
        omega_iarr = (double *)calloc(dp->ninc, sizeof(double));
        omega_oarr = (double *)calloc(dp->nout, sizeof(double));
        if ((omega_iarr == NULL) | (omega_oarr == NULL))
                error(SYSTEM, "out of memory in check_bsdf_data");
                                        /* incoming projected solid angles */
        hemi_total = .0;
        for (i = dp->ninc; i--; ) {
                dom = getBSDF_incohm(dp,i);
                if (dom <= .0) {
                        error(WARNING, "zero/negative incoming solid angle");
                        continue;
                }
                if (!getBSDF_incvec(v,dp,i) || v[2] > FTINY) {
                        error(WARNING, "illegal incoming BSDF direction");
                        free(omega_iarr); free(omega_oarr);
                        return(0);
                }
                hemi_total += omega_iarr[i] = dom * -v[2];
        }
        if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
                sprintf(errmsg, "incoming BSDF hemisphere off by %.1f%%",
                                100.*(hemi_total/PI - 1.));
                error(WARNING, errmsg);
        }
        dom = PI / hemi_total;          /* fix normalization */
        for (i = dp->ninc; i--; )
                omega_iarr[i] *= dom;
                                        /* outgoing projected solid angles */
        hemi_total = .0;
        for (o = dp->nout; o--; ) {
                dom = getBSDF_outohm(dp,o);
                if (dom <= .0) {
                        error(WARNING, "zero/negative outgoing solid angle");
                        continue;
                }
                if (!getBSDF_outvec(v,dp,o) || v[2] < -FTINY) {
                        error(WARNING, "illegal outgoing BSDF direction");
                        free(omega_iarr); free(omega_oarr);
                        return(0);
                }
                hemi_total += omega_oarr[o] = dom * v[2];
        }
        if ((hemi_total > 1.02*PI) | (hemi_total < 0.98*PI)) {
                sprintf(errmsg, "outgoing BSDF hemisphere off by %.1f%%",
                                100.*(hemi_total/PI - 1.));
                error(WARNING, errmsg);
        }
        dom = PI / hemi_total;          /* fix normalization */
        for (o = dp->nout; o--; )
                omega_oarr[o] *= dom;
        nneg = 0;                       /* check outgoing totals */
        for (i = 0; i < dp->ninc; i++) {
                hemi_total = .0;
                for (o = dp->nout; o--; ) {
                        double  f = BSDF_value(dp,i,o);
                        if (f >= .0)
                                hemi_total += f*omega_oarr[o];
                        else {
                                nneg += (f < -FTINY);
                                BSDF_value(dp,i,o) = .0f;
                        }
                }
                if (hemi_total > 1.01) {
                        sprintf(errmsg,
                        "incoming BSDF direction %d passes %.1f%% of light",
                                        i, 100.*hemi_total);
                        error(WARNING, errmsg);
                }
        }
        if (nneg) {
                sprintf(errmsg, "%d negative BSDF values (ignored)", nneg);
                error(WARNING, errmsg);
        }
        full_total = .0;                /* reverse roles and check again */
        for (o = 0; o < dp->nout; o++) {
                hemi_total = .0;
                for (i = dp->ninc; i--; )
                        hemi_total += BSDF_value(dp,i,o) * omega_iarr[i];

                if (hemi_total > 1.01) {
                        sprintf(errmsg,
                        "outgoing BSDF direction %d collects %.1f%% of light",
                                        o, 100.*hemi_total);
                        error(WARNING, errmsg);
                }
                full_total += hemi_total*omega_oarr[o];
        }
        full_total /= PI;
        if (full_total > 1.00001) {
                sprintf(errmsg, "BSDF transfers %.4f%% of light",
                                100.*full_total);
                error(WARNING, errmsg);
        }
        free(omega_iarr); free(omega_oarr);
        return(1);
}


struct BSDF_data *
load_BSDF(              /* load BSDF data from file */
        char *fname
)
{
        char                    *path;
        ezxml_t                 fl, wtl, wld, wdb;
        struct BSDF_data        *dp;
        
        path = getpath(fname, getrlibpath(), R_OK);
        if (path == NULL) {
                sprintf(errmsg, "cannot find BSDF file \"%s\"", fname);
                error(WARNING, errmsg);
                return(NULL);
        }
        fl = ezxml_parse_file(path);
        if (fl == NULL) {
                sprintf(errmsg, "cannot open BSDF \"%s\"", path);
                error(WARNING, errmsg);
                return(NULL);
        }
        if (ezxml_error(fl)[0]) {
                sprintf(errmsg, "BSDF \"%s\" %s", path, ezxml_error(fl));
                error(WARNING, errmsg);
                ezxml_free(fl);
                return(NULL);
        }
        if (strcmp(ezxml_name(fl), "WindowElement")) {
                sprintf(errmsg,
                        "BSDF \"%s\": top level node not 'WindowElement'",
                                path);
                error(WARNING, errmsg);
                ezxml_free(fl);
                return(NULL);
        }
        wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
        if (strcasecmp(ezxml_txt(ezxml_child(ezxml_child(wtl,
                        "DataDefinition"), "IncidentDataStructure")),
                        "Columns")) {
                sprintf(errmsg,
                        "BSDF \"%s\": unsupported IncidentDataStructure",
                                path);
                error(WARNING, errmsg);
                ezxml_free(fl);
                return(NULL);
        }               
        load_angle_basis(ezxml_child(ezxml_child(wtl,
                                "DataDefinition"), "AngleBasis"));
        dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
        load_geometry(dp, ezxml_child(wtl, "Material"));
        for (wld = ezxml_child(wtl, "WavelengthData");
                                wld != NULL; wld = wld->next) {
                if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible"))
                        continue;
                wdb = ezxml_child(wld, "WavelengthDataBlock");
                if (wdb == NULL) continue;
                if (strcasecmp(ezxml_txt(ezxml_child(wdb,"WavelengthDataDirection")),
                                        "Transmission Front"))
                        continue;
                load_bsdf_data(dp, wdb);        /* load front BTDF */
                break;                          /* ignore the rest */
        }
        ezxml_free(fl);                         /* done with XML file */
        if (!check_bsdf_data(dp)) {
                sprintf(errmsg, "bad/missing BTDF data in \"%s\"", path);
                error(WARNING, errmsg);
                free_BSDF(dp);
                dp = NULL;
        }
        return(dp);
}


void
free_BSDF(              /* free BSDF data structure */
        struct BSDF_data *b
)
{
        if (b == NULL)
                return;
        if (b->mgf != NULL)
                free(b->mgf);
        if (b->bsdf != NULL)
                free(b->bsdf);
        free(b);
}


int
r_BSDF_incvec(          /* compute random input vector at given location */
        FVECT v,
        struct BSDF_data *b,
        int i,
        double rv,
        MAT4 xm
)
{
        FVECT   pert;
        double  rad;
        int     j;
        
        if (!getBSDF_incvec(v, b, i))
                return(0);
        rad = sqrt(getBSDF_incohm(b, i) / PI);
        multisamp(pert, 3, rv);
        for (j = 0; j < 3; j++)
                v[j] += rad*(2.*pert[j] - 1.);
        if (xm != NULL)
                multv3(v, v, xm);
        return(normalize(v) != 0.0);
}


int
r_BSDF_outvec(          /* compute random output vector at given location */
        FVECT v,
        struct BSDF_data *b,
        int o,
        double rv,
        MAT4 xm
)
{
        FVECT   pert;
        double  rad;
        int     j;
        
        if (!getBSDF_outvec(v, b, o))
                return(0);
        rad = sqrt(getBSDF_outohm(b, o) / PI);
        multisamp(pert, 3, rv);
        for (j = 0; j < 3; j++)
                v[j] += rad*(2.*pert[j] - 1.);
        if (xm != NULL)
                multv3(v, v, xm);
        return(normalize(v) != 0.0);
}


static int
addrot(                 /* compute rotation (x,y,z) => (xp,yp,zp) */
        char *xfarg[],
        FVECT xp,
        FVECT yp,
        FVECT zp
)
{
        static char     bufs[3][16];
        int     bn = 0;
        char    **xfp = xfarg;
        double  theta;

        if (yp[2]*yp[2] + zp[2]*zp[2] < 2.*FTINY*FTINY) {
                /* Special case for X' along Z-axis */
                theta = -atan2(yp[0], yp[1]);
                *xfp++ = "-ry";
                *xfp++ = xp[2] < 0.0 ? "90" : "-90";
                *xfp++ = "-rz";
                sprintf(bufs[bn], "%f", theta*(180./PI));
                *xfp++ = bufs[bn++];
                return(xfp - xfarg);
        }
        theta = atan2(yp[2], zp[2]);
        if (!FEQ(theta,0.0)) {
                *xfp++ = "-rx";
                sprintf(bufs[bn], "%f", theta*(180./PI));
                *xfp++ = bufs[bn++];
        }
        theta = asin(-xp[2]);
        if (!FEQ(theta,0.0)) {
                *xfp++ = "-ry";
                sprintf(bufs[bn], " %f", theta*(180./PI));
                *xfp++ = bufs[bn++];
        }
        theta = atan2(xp[1], xp[0]);
        if (!FEQ(theta,0.0)) {
                *xfp++ = "-rz";
                sprintf(bufs[bn], "%f", theta*(180./PI));
                *xfp++ = bufs[bn++];
        }
        *xfp = NULL;
        return(xfp - xfarg);
}


int
getBSDF_xfm(            /* compute BSDF orient. -> world orient. transform */
        MAT4 xm,
        FVECT nrm,
        UpDir ud,
        char *xfbuf
)
{
        char    *xfargs[7];
        XF      myxf;
        FVECT   updir, xdest, ydest;
        int     i;

        updir[0] = updir[1] = updir[2] = 0.;
        switch (ud) {
        case UDzneg:
                updir[2] = -1.;
                break;
        case UDyneg:
                updir[1] = -1.;
                break;
        case UDxneg:
                updir[0] = -1.;
                break;
        case UDxpos:
                updir[0] = 1.;
                break;
        case UDypos:
                updir[1] = 1.;
                break;
        case UDzpos:
                updir[2] = 1.;
                break;
        case UDunknown:
                return(0);
        }
        fcross(xdest, updir, nrm);
        if (normalize(xdest) == 0.0)
                return(0);
        fcross(ydest, nrm, xdest);
        xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
        copymat4(xm, myxf.xfm);
        if (xfbuf == NULL)
                return(1);
                                /* return xf arguments as well */
        for (i = 0; xfargs[i] != NULL; i++) {
                *xfbuf++ = ' ';
                strcpy(xfbuf, xfargs[i]);
                while (*xfbuf) ++xfbuf;
        }
        return(1);
}
Revision:	2.12
Committed:	Thu Jan 6 04:40:22 2011 UTC (14 years, 9 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.11:	+6 -6 lines
Log Message:	Typo!
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.12	static const char RCSid[] = "$Id: bsdf.c,v 2.11 2011/01/06 04:39:08 greg Exp $";
3	greg	2.1	#endif
4			/*
5			* Routines for handling BSDF data
6			*/
7
8			#include "standard.h"
9			#include "bsdf.h"
10			#include "paths.h"
11			#include "ezxml.h"
12			#include <ctype.h>
13
14			#define MAXLATS 46 /* maximum number of latitudes */
15
16			/* BSDF angle specification */
17			typedef struct {
18			char name[64]; /* basis name */
19			int nangles; /* total number of directions */
20			struct {
21			float tmin; /* starting theta */
22			short nphis; /* number of phis (0 term) */
23			} lat[MAXLATS+1]; /* latitudes */
24			} ANGLE_BASIS;
25
26	greg	2.4	#define MAXABASES 7 /* limit on defined bases */
27	greg	2.1
28			static ANGLE_BASIS abase_list[MAXABASES] = {
29			{
30			"LBNL/Klems Full", 145,
31			{ {-5., 1},
32			{5., 8},
33			{15., 16},
34			{25., 20},
35			{35., 24},
36			{45., 24},
37			{55., 24},
38			{65., 16},
39			{75., 12},
40			{90., 0} }
41			}, {
42			"LBNL/Klems Half", 73,
43			{ {-6.5, 1},
44			{6.5, 8},
45			{19.5, 12},
46			{32.5, 16},
47			{46.5, 20},
48			{61.5, 12},
49			{76.5, 4},
50			{90., 0} }
51			}, {
52			"LBNL/Klems Quarter", 41,
53			{ {-9., 1},
54			{9., 8},
55			{27., 12},
56			{46., 12},
57			{66., 8},
58			{90., 0} }
59			}
60			};
61
62			static int nabases = 3; /* current number of defined bases */
63
64	greg	2.9	#define FEQ(a,b) ((a)-(b) <= 1e-6 && (b)-(a) <= 1e-6)
65
66			static int
67			fequal(double a, double b)
68			{
69			if (b != .0)
70			a = a/b - 1.;
71			return((a <= 1e-6) & (a >= -1e-6));
72			}
73	greg	2.3
74			// returns the name of the given tag
75			#ifdef ezxml_name
76			#undef ezxml_name
77			static char *
78			ezxml_name(ezxml_t xml)
79			{
80			if (xml == NULL)
81			return(NULL);
82			return(xml->name);
83			}
84			#endif
85
86			// returns the given tag's character content or empty string if none
87			#ifdef ezxml_txt
88			#undef ezxml_txt
89			static char *
90			ezxml_txt(ezxml_t xml)
91			{
92			if (xml == NULL)
93			return("");
94			return(xml->txt);
95			}
96			#endif
97
98	greg	2.1
99			static int
100			ab_getvec( /* get vector for this angle basis index */
101			FVECT v,
102			int ndx,
103			void *p
104			)
105			{
106			ANGLE_BASIS ab = (ANGLE_BASIS )p;
107			int li;
108	greg	2.2	double pol, azi, d;
109	greg	2.1
110			if ((ndx < 0) \| (ndx >= ab->nangles))
111			return(0);
112			for (li = 0; ndx >= ab->lat[li].nphis; li++)
113			ndx -= ab->lat[li].nphis;
114	greg	2.2	pol = PI/180.0.5(ab->lat[li].tmin + ab->lat[li+1].tmin);
115	greg	2.1	azi = 2.PIndx/ab->lat[li].nphis;
116	greg	2.2	v[2] = d = cos(pol);
117			d = sqrt(1. - dd); / sin(pol) */
118	greg	2.1	v[0] = cos(azi)*d;
119			v[1] = sin(azi)*d;
120			return(1);
121			}
122
123
124			static int
125			ab_getndx( /* get index corresponding to the given vector */
126			FVECT v,
127			void *p
128			)
129			{
130			ANGLE_BASIS ab = (ANGLE_BASIS )p;
131			int li, ndx;
132	greg	2.2	double pol, azi, d;
133	greg	2.1
134			if ((v[2] < -1.0) \| (v[2] > 1.0))
135			return(-1);
136	greg	2.2	pol = 180.0/PI*acos(v[2]);
137	greg	2.1	azi = 180.0/PI*atan2(v[1], v[0]);
138			if (azi < 0.0) azi += 360.0;
139	greg	2.2	for (li = 1; ab->lat[li].tmin <= pol; li++)
140	greg	2.1	if (!ab->lat[li].nphis)
141			return(-1);
142			--li;
143			ndx = (int)((1./360.)aziab->lat[li].nphis + 0.5);
144			if (ndx >= ab->lat[li].nphis) ndx = 0;
145			while (li--)
146			ndx += ab->lat[li].nphis;
147			return(ndx);
148			}
149
150
151			static double
152			ab_getohm( /* get solid angle for this angle basis index */
153			int ndx,
154			void *p
155			)
156			{
157			ANGLE_BASIS ab = (ANGLE_BASIS )p;
158			int li;
159			double theta, theta1;
160
161			if ((ndx < 0) \| (ndx >= ab->nangles))
162			return(0);
163			for (li = 0; ndx >= ab->lat[li].nphis; li++)
164			ndx -= ab->lat[li].nphis;
165			theta1 = PI/180. * ab->lat[li+1].tmin;
166			if (ab->lat[li].nphis == 1) { /* special case */
167			if (ab->lat[li].tmin > FTINY)
168			error(USER, "unsupported BSDF coordinate system");
169			return(2.PI(1. - cos(theta1)));
170			}
171			theta = PI/180. * ab->lat[li].tmin;
172			return(2.PI(cos(theta) - cos(theta1))/(double)ab->lat[li].nphis);
173			}
174
175
176			static int
177			ab_getvecR( /* get reverse vector for this angle basis index */
178			FVECT v,
179			int ndx,
180			void *p
181			)
182			{
183			if (!ab_getvec(v, ndx, p))
184			return(0);
185
186			v[0] = -v[0];
187			v[1] = -v[1];
188			v[2] = -v[2];
189
190			return(1);
191			}
192
193
194			static int
195			ab_getndxR( /* get index corresponding to the reverse vector */
196			FVECT v,
197			void *p
198			)
199			{
200			FVECT v2;
201
202			v2[0] = -v[0];
203			v2[1] = -v[1];
204			v2[2] = -v[2];
205
206			return ab_getndx(v2, p);
207			}
208
209
210			static void
211	greg	2.4	load_angle_basis( /* load custom BSDF angle basis */
212	greg	2.3	ezxml_t wab
213			)
214			{
215			char *abname = ezxml_txt(ezxml_child(wab, "AngleBasisName"));
216			ezxml_t wbb;
217			int i;
218
219	greg	2.4	if (!abname \|\| !*abname)
220	greg	2.3	return;
221			for (i = nabases; i--; )
222	greg	2.12	if (!strcasecmp(abname, abase_list[i].name))
223	greg	2.4	return; /* assume it's the same */
224	greg	2.3	if (nabases >= MAXABASES)
225			error(INTERNAL, "too many angle bases");
226			strcpy(abase_list[nabases].name, abname);
227			abase_list[nabases].nangles = 0;
228			for (i = 0, wbb = ezxml_child(wab, "AngleBasisBlock");
229			wbb != NULL; i++, wbb = wbb->next) {
230			if (i >= MAXLATS)
231			error(INTERNAL, "too many latitudes in custom basis");
232			abase_list[nabases].lat[i+1].tmin = atof(ezxml_txt(
233			ezxml_child(ezxml_child(wbb,
234			"ThetaBounds"), "UpperTheta")));
235			if (!i)
236			abase_list[nabases].lat[i].tmin =
237			-abase_list[nabases].lat[i+1].tmin;
238	greg	2.9	else if (!fequal(atof(ezxml_txt(ezxml_child(ezxml_child(wbb,
239	greg	2.3	"ThetaBounds"), "LowerTheta"))),
240			abase_list[nabases].lat[i].tmin))
241			error(WARNING, "theta values disagree in custom basis");
242			abase_list[nabases].nangles +=
243			abase_list[nabases].lat[i].nphis =
244			atoi(ezxml_txt(ezxml_child(wbb, "nPhis")));
245			}
246			abase_list[nabases++].lat[i].nphis = 0;
247			}
248
249
250	greg	2.6	static double
251			to_meters( /* return factor to convert given unit to meters */
252			const char *unit
253			)
254			{
255			if (unit == NULL) return(1.); /* safe assumption? */
256			if (!strcasecmp(unit, "Meter")) return(1.);
257			if (!strcasecmp(unit, "Foot")) return(.3048);
258			if (!strcasecmp(unit, "Inch")) return(.0254);
259			if (!strcasecmp(unit, "Centimeter")) return(.01);
260	greg	2.8	if (!strcasecmp(unit, "Millimeter")) return(.001);
261	greg	2.6	sprintf(errmsg, "unknown dimensional unit '%s'", unit);
262			error(USER, errmsg);
263			}
264
265
266			static void
267			load_geometry( /* load geometric dimensions and description (if any) */
268			struct BSDF_data *dp,
269			ezxml_t wdb
270			)
271			{
272			ezxml_t geom;
273			double cfact;
274			const char fmt, mgfstr;
275
276			dp->dim[0] = dp->dim[1] = dp->dim[2] = 0;
277			dp->mgf = NULL;
278			if ((geom = ezxml_child(wdb, "Width")) != NULL)
279			dp->dim[0] = atof(ezxml_txt(geom)) *
280			to_meters(ezxml_attr(geom, "unit"));
281			if ((geom = ezxml_child(wdb, "Height")) != NULL)
282			dp->dim[1] = atof(ezxml_txt(geom)) *
283			to_meters(ezxml_attr(geom, "unit"));
284			if ((geom = ezxml_child(wdb, "Thickness")) != NULL)
285			dp->dim[2] = atof(ezxml_txt(geom)) *
286			to_meters(ezxml_attr(geom, "unit"));
287			if ((geom = ezxml_child(wdb, "Geometry")) == NULL \|\|
288			(mgfstr = ezxml_txt(geom)) == NULL)
289			return;
290			if ((fmt = ezxml_attr(geom, "format")) != NULL &&
291			strcasecmp(fmt, "MGF")) {
292			sprintf(errmsg, "unrecognized geometry format '%s'", fmt);
293			error(WARNING, errmsg);
294			return;
295			}
296			cfact = to_meters(ezxml_attr(geom, "unit"));
297			dp->mgf = (char *)malloc(strlen(mgfstr)+32);
298			if (dp->mgf == NULL)
299			error(SYSTEM, "out of memory in load_geometry");
300			if (cfact < 0.99 \|\| cfact > 1.01)
301			sprintf(dp->mgf, "xf -s %.5f\n%s\nxf\n", cfact, mgfstr);
302			else
303			strcpy(dp->mgf, mgfstr);
304			}
305
306
307	greg	2.3	static void
308	greg	2.1	load_bsdf_data( /* load BSDF distribution for this wavelength */
309			struct BSDF_data *dp,
310			ezxml_t wdb
311			)
312			{
313			char *cbasis = ezxml_txt(ezxml_child(wdb,"ColumnAngleBasis"));
314			char *rbasis = ezxml_txt(ezxml_child(wdb,"RowAngleBasis"));
315			char *sdata;
316			int i;
317
318	greg	2.4	if ((!cbasis \|\| !cbasis) \| (!rbasis \|\| !rbasis)) {
319	greg	2.1	error(WARNING, "missing column/row basis for BSDF");
320			return;
321			}
322			for (i = nabases; i--; )
323	greg	2.12	if (!strcasecmp(cbasis, abase_list[i].name)) {
324	greg	2.1	dp->ninc = abase_list[i].nangles;
325			dp->ib_priv = (void *)&abase_list[i];
326			dp->ib_vec = ab_getvecR;
327			dp->ib_ndx = ab_getndxR;
328			dp->ib_ohm = ab_getohm;
329			break;
330			}
331			if (i < 0) {
332	greg	2.4	sprintf(errmsg, "undefined ColumnAngleBasis '%s'", cbasis);
333	greg	2.1	error(WARNING, errmsg);
334			return;
335			}
336			for (i = nabases; i--; )
337	greg	2.12	if (!strcasecmp(rbasis, abase_list[i].name)) {
338	greg	2.1	dp->nout = abase_list[i].nangles;
339			dp->ob_priv = (void *)&abase_list[i];
340			dp->ob_vec = ab_getvec;
341			dp->ob_ndx = ab_getndx;
342			dp->ob_ohm = ab_getohm;
343			break;
344			}
345			if (i < 0) {
346	greg	2.4	sprintf(errmsg, "undefined RowAngleBasis '%s'", cbasis);
347	greg	2.1	error(WARNING, errmsg);
348			return;
349			}
350			/* read BSDF data */
351			sdata = ezxml_txt(ezxml_child(wdb,"ScatteringData"));
352	greg	2.4	if (!sdata \|\| !*sdata) {
353	greg	2.1	error(WARNING, "missing BSDF ScatteringData");
354			return;
355			}
356			dp->bsdf = (float )malloc(sizeof(float)dp->ninc*dp->nout);
357			if (dp->bsdf == NULL)
358			error(SYSTEM, "out of memory in load_bsdf_data");
359			for (i = 0; i < dp->ninc*dp->nout; i++) {
360			char *sdnext = fskip(sdata);
361			if (sdnext == NULL) {
362			error(WARNING, "bad/missing BSDF ScatteringData");
363			free(dp->bsdf); dp->bsdf = NULL;
364			return;
365			}
366			while (sdnext && isspace(sdnext))
367			sdnext++;
368			if (*sdnext == ',') sdnext++;
369			dp->bsdf[i] = atof(sdata);
370			sdata = sdnext;
371			}
372			while (isspace(*sdata))
373			sdata++;
374			if (*sdata) {
375			sprintf(errmsg, "%d extra characters after BSDF ScatteringData",
376	greg	2.5	(int)strlen(sdata));
377	greg	2.1	error(WARNING, errmsg);
378			}
379			}
380
381
382			static int
383			check_bsdf_data( /* check that BSDF data is sane */
384			struct BSDF_data *dp
385			)
386			{
387	greg	2.2	double omega_iarr, omega_oarr;
388	greg	2.7	double dom, contrib, hemi_total, full_total;
389	greg	2.1	int nneg;
390	greg	2.2	FVECT v;
391	greg	2.1	int i, o;
392
393			if (dp == NULL \|\| dp->bsdf == NULL)
394			return(0);
395	greg	2.2	omega_iarr = (double *)calloc(dp->ninc, sizeof(double));
396			omega_oarr = (double *)calloc(dp->nout, sizeof(double));
397			if ((omega_iarr == NULL) \| (omega_oarr == NULL))
398	greg	2.1	error(SYSTEM, "out of memory in check_bsdf_data");
399	greg	2.2	/* incoming projected solid angles */
400			hemi_total = .0;
401			for (i = dp->ninc; i--; ) {
402			dom = getBSDF_incohm(dp,i);
403			if (dom <= .0) {
404			error(WARNING, "zero/negative incoming solid angle");
405			continue;
406			}
407			if (!getBSDF_incvec(v,dp,i) \|\| v[2] > FTINY) {
408			error(WARNING, "illegal incoming BSDF direction");
409			free(omega_iarr); free(omega_oarr);
410			return(0);
411			}
412			hemi_total += omega_iarr[i] = dom * -v[2];
413			}
414			if ((hemi_total > 1.02PI) \| (hemi_total < 0.98PI)) {
415			sprintf(errmsg, "incoming BSDF hemisphere off by %.1f%%",
416			100.*(hemi_total/PI - 1.));
417			error(WARNING, errmsg);
418			}
419			dom = PI / hemi_total; /* fix normalization */
420			for (i = dp->ninc; i--; )
421			omega_iarr[i] *= dom;
422			/* outgoing projected solid angles */
423	greg	2.1	hemi_total = .0;
424			for (o = dp->nout; o--; ) {
425			dom = getBSDF_outohm(dp,o);
426			if (dom <= .0) {
427	greg	2.2	error(WARNING, "zero/negative outgoing solid angle");
428	greg	2.1	continue;
429			}
430			if (!getBSDF_outvec(v,dp,o) \|\| v[2] < -FTINY) {
431			error(WARNING, "illegal outgoing BSDF direction");
432	greg	2.2	free(omega_iarr); free(omega_oarr);
433	greg	2.1	return(0);
434			}
435	greg	2.2	hemi_total += omega_oarr[o] = dom * v[2];
436	greg	2.1	}
437			if ((hemi_total > 1.02PI) \| (hemi_total < 0.98PI)) {
438			sprintf(errmsg, "outgoing BSDF hemisphere off by %.1f%%",
439			100.*(hemi_total/PI - 1.));
440			error(WARNING, errmsg);
441			}
442	greg	2.2	dom = PI / hemi_total; /* fix normalization */
443	greg	2.1	for (o = dp->nout; o--; )
444	greg	2.2	omega_oarr[o] *= dom;
445			nneg = 0; /* check outgoing totals */
446			for (i = 0; i < dp->ninc; i++) {
447	greg	2.1	hemi_total = .0;
448			for (o = dp->nout; o--; ) {
449			double f = BSDF_value(dp,i,o);
450	greg	2.2	if (f >= .0)
451			hemi_total += f*omega_oarr[o];
452			else {
453			nneg += (f < -FTINY);
454			BSDF_value(dp,i,o) = .0f;
455			}
456	greg	2.1	}
457	greg	2.8	if (hemi_total > 1.01) {
458	greg	2.2	sprintf(errmsg,
459			"incoming BSDF direction %d passes %.1f%% of light",
460			i, 100.*hemi_total);
461	greg	2.1	error(WARNING, errmsg);
462			}
463			}
464	greg	2.2	if (nneg) {
465			sprintf(errmsg, "%d negative BSDF values (ignored)", nneg);
466	greg	2.1	error(WARNING, errmsg);
467			}
468	greg	2.7	full_total = .0; /* reverse roles and check again */
469	greg	2.2	for (o = 0; o < dp->nout; o++) {
470			hemi_total = .0;
471			for (i = dp->ninc; i--; )
472			hemi_total += BSDF_value(dp,i,o) * omega_iarr[i];
473
474	greg	2.8	if (hemi_total > 1.01) {
475	greg	2.2	sprintf(errmsg,
476			"outgoing BSDF direction %d collects %.1f%% of light",
477			o, 100.*hemi_total);
478			error(WARNING, errmsg);
479			}
480	greg	2.7	full_total += hemi_total*omega_oarr[o];
481			}
482			full_total /= PI;
483	greg	2.8	if (full_total > 1.00001) {
484			sprintf(errmsg, "BSDF transfers %.4f%% of light",
485	greg	2.7	100.*full_total);
486			error(WARNING, errmsg);
487	greg	2.2	}
488			free(omega_iarr); free(omega_oarr);
489	greg	2.1	return(1);
490			}
491
492	greg	2.6
493	greg	2.1	struct BSDF_data *
494			load_BSDF( /* load BSDF data from file */
495			char *fname
496			)
497			{
498			char *path;
499			ezxml_t fl, wtl, wld, wdb;
500			struct BSDF_data *dp;
501
502			path = getpath(fname, getrlibpath(), R_OK);
503			if (path == NULL) {
504			sprintf(errmsg, "cannot find BSDF file \"%s\"", fname);
505			error(WARNING, errmsg);
506			return(NULL);
507			}
508			fl = ezxml_parse_file(path);
509			if (fl == NULL) {
510			sprintf(errmsg, "cannot open BSDF \"%s\"", path);
511			error(WARNING, errmsg);
512			return(NULL);
513			}
514			if (ezxml_error(fl)[0]) {
515			sprintf(errmsg, "BSDF \"%s\" %s", path, ezxml_error(fl));
516			error(WARNING, errmsg);
517			ezxml_free(fl);
518			return(NULL);
519			}
520			if (strcmp(ezxml_name(fl), "WindowElement")) {
521			sprintf(errmsg,
522			"BSDF \"%s\": top level node not 'WindowElement'",
523			path);
524			error(WARNING, errmsg);
525			ezxml_free(fl);
526			return(NULL);
527			}
528			wtl = ezxml_child(ezxml_child(fl, "Optical"), "Layer");
529	greg	2.10	if (strcasecmp(ezxml_txt(ezxml_child(ezxml_child(wtl,
530			"DataDefinition"), "IncidentDataStructure")),
531			"Columns")) {
532			sprintf(errmsg,
533			"BSDF \"%s\": unsupported IncidentDataStructure",
534			path);
535			error(WARNING, errmsg);
536			ezxml_free(fl);
537			return(NULL);
538			}
539	greg	2.3	load_angle_basis(ezxml_child(ezxml_child(wtl,
540			"DataDefinition"), "AngleBasis"));
541	greg	2.1	dp = (struct BSDF_data *)calloc(1, sizeof(struct BSDF_data));
542	greg	2.6	load_geometry(dp, ezxml_child(wtl, "Material"));
543	greg	2.1	for (wld = ezxml_child(wtl, "WavelengthData");
544			wld != NULL; wld = wld->next) {
545	greg	2.12	if (strcasecmp(ezxml_txt(ezxml_child(wld,"Wavelength")), "Visible"))
546	greg	2.1	continue;
547			wdb = ezxml_child(wld, "WavelengthDataBlock");
548			if (wdb == NULL) continue;
549	greg	2.12	if (strcasecmp(ezxml_txt(ezxml_child(wdb,"WavelengthDataDirection")),
550	greg	2.1	"Transmission Front"))
551			continue;
552			load_bsdf_data(dp, wdb); /* load front BTDF */
553			break; /* ignore the rest */
554			}
555			ezxml_free(fl); /* done with XML file */
556			if (!check_bsdf_data(dp)) {
557			sprintf(errmsg, "bad/missing BTDF data in \"%s\"", path);
558			error(WARNING, errmsg);
559			free_BSDF(dp);
560			dp = NULL;
561			}
562			return(dp);
563			}
564
565
566			void
567			free_BSDF( /* free BSDF data structure */
568			struct BSDF_data *b
569			)
570			{
571			if (b == NULL)
572			return;
573	greg	2.6	if (b->mgf != NULL)
574			free(b->mgf);
575	greg	2.1	if (b->bsdf != NULL)
576			free(b->bsdf);
577			free(b);
578			}
579
580
581			int
582			r_BSDF_incvec( /* compute random input vector at given location */
583			FVECT v,
584			struct BSDF_data *b,
585			int i,
586			double rv,
587			MAT4 xm
588			)
589			{
590			FVECT pert;
591			double rad;
592			int j;
593
594			if (!getBSDF_incvec(v, b, i))
595			return(0);
596			rad = sqrt(getBSDF_incohm(b, i) / PI);
597			multisamp(pert, 3, rv);
598			for (j = 0; j < 3; j++)
599			v[j] += rad(2.pert[j] - 1.);
600			if (xm != NULL)
601			multv3(v, v, xm);
602			return(normalize(v) != 0.0);
603			}
604
605
606			int
607			r_BSDF_outvec( /* compute random output vector at given location */
608			FVECT v,
609			struct BSDF_data *b,
610			int o,
611			double rv,
612			MAT4 xm
613			)
614			{
615			FVECT pert;
616			double rad;
617			int j;
618
619			if (!getBSDF_outvec(v, b, o))
620			return(0);
621			rad = sqrt(getBSDF_outohm(b, o) / PI);
622			multisamp(pert, 3, rv);
623			for (j = 0; j < 3; j++)
624			v[j] += rad(2.pert[j] - 1.);
625			if (xm != NULL)
626			multv3(v, v, xm);
627			return(normalize(v) != 0.0);
628			}
629
630
631			static int
632			addrot( /* compute rotation (x,y,z) => (xp,yp,zp) */
633			char *xfarg[],
634			FVECT xp,
635			FVECT yp,
636			FVECT zp
637			)
638			{
639			static char bufs[3][16];
640			int bn = 0;
641			char **xfp = xfarg;
642			double theta;
643
644			if (yp[2]yp[2] + zp[2]zp[2] < 2.FTINYFTINY) {
645			/* Special case for X' along Z-axis */
646			theta = -atan2(yp[0], yp[1]);
647			*xfp++ = "-ry";
648			*xfp++ = xp[2] < 0.0 ? "90" : "-90";
649			*xfp++ = "-rz";
650			sprintf(bufs[bn], "%f", theta*(180./PI));
651			*xfp++ = bufs[bn++];
652			return(xfp - xfarg);
653			}
654			theta = atan2(yp[2], zp[2]);
655			if (!FEQ(theta,0.0)) {
656			*xfp++ = "-rx";
657			sprintf(bufs[bn], "%f", theta*(180./PI));
658			*xfp++ = bufs[bn++];
659			}
660			theta = asin(-xp[2]);
661			if (!FEQ(theta,0.0)) {
662			*xfp++ = "-ry";
663			sprintf(bufs[bn], " %f", theta*(180./PI));
664			*xfp++ = bufs[bn++];
665			}
666			theta = atan2(xp[1], xp[0]);
667			if (!FEQ(theta,0.0)) {
668			*xfp++ = "-rz";
669			sprintf(bufs[bn], "%f", theta*(180./PI));
670			*xfp++ = bufs[bn++];
671			}
672			*xfp = NULL;
673			return(xfp - xfarg);
674			}
675
676
677			int
678			getBSDF_xfm( /* compute BSDF orient. -> world orient. transform */
679			MAT4 xm,
680			FVECT nrm,
681	greg	2.6	UpDir ud,
682			char *xfbuf
683	greg	2.1	)
684			{
685			char *xfargs[7];
686			XF myxf;
687			FVECT updir, xdest, ydest;
688	greg	2.6	int i;
689	greg	2.1
690			updir[0] = updir[1] = updir[2] = 0.;
691			switch (ud) {
692			case UDzneg:
693			updir[2] = -1.;
694			break;
695			case UDyneg:
696			updir[1] = -1.;
697			break;
698			case UDxneg:
699			updir[0] = -1.;
700			break;
701			case UDxpos:
702			updir[0] = 1.;
703			break;
704			case UDypos:
705			updir[1] = 1.;
706			break;
707			case UDzpos:
708			updir[2] = 1.;
709			break;
710			case UDunknown:
711			return(0);
712			}
713			fcross(xdest, updir, nrm);
714			if (normalize(xdest) == 0.0)
715			return(0);
716			fcross(ydest, nrm, xdest);
717			xf(&myxf, addrot(xfargs, xdest, ydest, nrm), xfargs);
718			copymat4(xm, myxf.xfm);
719	greg	2.6	if (xfbuf == NULL)
720			return(1);
721			/* return xf arguments as well */
722			for (i = 0; xfargs[i] != NULL; i++) {
723			*xfbuf++ = ' ';
724			strcpy(xfbuf, xfargs[i]);
725			while (*xfbuf) ++xfbuf;
726			}
727	greg	2.1	return(1);
728			}