src/util/glarendx.c

#ifndef lint
static const char       RCSid[] = "$Id: glarendx.c,v 2.9 2003/11/19 16:21:28 greg Exp $";
#endif
/*
 * Compute Glare Index given by program name or -t option:
 *
 *      dgi -           Daylight Glare Index
 *      brs_gi -        Building Research Station Glare Index (Petherbridge
 *                                                             & Hopkinson)
 *      ugr -           Unified Glare Rating System (Fischer)
 *      guth_dgr -      Guth discomfort glare rating
 *      guth_vcp -      Guth visual comfort probability
 *      cie_cgi -       CIE Glare Index (1983, due to Einhorn)
 *      vert_dir -      Direct vertical illuminance
 *      vert_ind -      Indirect vertical illuminance (from input)
 *      vert_ill -      Total vertical illuminance
 *
 *              12 April 1991   Greg Ward       EPFL
 *              19 April 1993   R. Compagnon    EPFL (added dgi, brs_gi, ugr)
 */

#include <string.h>

#include "standard.h"
#include "view.h"


struct glare_src {
        FVECT   dir;            /* source direction */
        double  dom;            /* solid angle */
        double  lum;            /* average luminance */
        struct glare_src        *next;
} *all_srcs = NULL;

struct glare_dir {
        double  ang;            /* angle (in radians) */
        double  indirect;       /* indirect illuminance */
        struct glare_dir        *next;
} *all_dirs = NULL;

typedef double gdfun(struct glare_dir *gd);

static gdfun dgi;
static gdfun brs_gi;
static gdfun ugr;
static gdfun guth_dgr;
static gdfun guth_vcp;
static gdfun cie_cgi;
static gdfun direct;
static gdfun indirect;
static gdfun total;

static gethfunc headline;
static void init(void);
static void read_input(void);
static void print_values(gdfun *func);
static double posindex(FVECT sd, FVECT vd, FVECT vu);


struct named_func {
        char    *name;
        gdfun   *func;
        char    *descrip;
} all_funcs[] = {
        {"dgi", dgi, "Daylight Glare Index"},
        {"brs_gi", brs_gi, "BRS Glare Index"},
        {"ugr", ugr, "Unified Glare Rating"},
        {"guth_vcp", guth_vcp, "Guth Visual Comfort Probability"},
        {"cie_cgi", cie_cgi, "CIE Glare Index (Einhorn)"},
        {"guth_dgr", guth_dgr, "Guth Disability Glare Rating"},
        {"vert_dir", direct, "Direct Vertical Illuminance"},
        {"vert_ill", total, "Total Vertical Illuminance"},
        {"vert_ind", indirect, "Indirect Vertical Illuminance"},
        {NULL}
};

#define newp(type)      (type *)malloc(sizeof(type))

char    *progname;
int     print_header = 1;

VIEW    midview = STDVIEW;

int     wrongformat = 0;


int
main(
        int     argc,
        char    *argv[]
)
{
        struct named_func       *funp;
        char    *progtail;
        int     i;
                                        /* get program name */
        progname = argv[0];
        progtail = strrchr(progname, '/');      /* final component */
        if (progtail == NULL)
                progtail = progname;
        else
                progtail++;
                                        /* get options */
        for (i = 1; i < argc && argv[i][0] == '-'; i++)
                switch (argv[i][1]) {
                case 't':
                        progtail = argv[++i];
                        break;
                case 'h':
                        print_header = 0;
                        break;
                default:
                        goto userr;
                }
        if (i < argc-1)
                goto userr;
        if (i == argc-1)                /* open file */
                if (freopen(argv[i], "r", stdin) == NULL) {
                        perror(argv[i]);
                        exit(1);
                }
                                        /* find and run calculation */
        for (funp = all_funcs; funp->name != NULL; funp++)
                if (!strcmp(funp->name, progtail)) {
                        init();
                        read_input();
                        if (print_header) {
                                printargs(i, argv, stdout);
                                putchar('\n');
                        }
                        print_values(funp->func);
                        exit(0);                /* we're done */
                }
                                        /* invalid function */
userr:
        fprintf(stderr, "Usage: %s -t type [-h] [input]\n", progname);
        fprintf(stderr, "\twhere 'type' is one of the following:\n");
        for (funp = all_funcs; funp->name != NULL; funp++)
                fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
        exit(1);
}


static int
headline(                       /* get line from header */
        char    *s,
        void    *p
)
{
        char    fmt[32];

        if (print_header)               /* copy to output */
                fputs(s, stdout);
        if (isview(s))
                sscanview(&midview, s);
        else if (isformat(s)) {
                formatval(fmt, s);
                wrongformat = strcmp(fmt, "ascii");
        }
        return(0);
}


static void
init(void)                              /* initialize calculation */
{
                                        /* read header */
        getheader(stdin, headline, NULL);
        if (wrongformat) {
                fprintf(stderr, "%s: bad input format\n", progname);
                exit(1);
        }
                                        /* set view */
        if (setview(&midview) != NULL) {
                fprintf(stderr, "%s: bad view information in input\n",
                                progname);
                exit(1);
        }
}


static void
read_input(void)                        /* read glare sources from stdin */
{
#define S_SEARCH        0
#define S_SOURCE        1
#define S_DIREC         2
        int     state = S_SEARCH;
        char    buf[128];
        register struct glare_src       *gs;
        register struct glare_dir       *gd;

        while (fgets(buf, sizeof(buf), stdin) != NULL)
                switch (state) {
                case S_SEARCH:
                        if (!strcmp(buf, "BEGIN glare source\n"))
                                state = S_SOURCE;
                        else if (!strcmp(buf, "BEGIN indirect illuminance\n"))
                                state = S_DIREC;
                        break;
                case S_SOURCE:
                        if (!strncmp(buf, "END", 3)) {
                                state = S_SEARCH;
                                break;
                        }
                        if ((gs = newp(struct glare_src)) == NULL)
                                goto memerr;
                        if (sscanf(buf, "%lf %lf %lf %lf %lf",
                                        &gs->dir[0], &gs->dir[1], &gs->dir[2],
                                        &gs->dom, &gs->lum) != 5)
                                goto readerr;
                        normalize(gs->dir);
                        gs->next = all_srcs;
                        all_srcs = gs;
                        break;
                case S_DIREC:
                        if (!strncmp(buf, "END", 3)) {
                                state = S_SEARCH;
                                break;
                        }
                        if ((gd = newp(struct glare_dir)) == NULL)
                                goto memerr;
                        if (sscanf(buf, "%lf %lf",
                                        &gd->ang, &gd->indirect) != 2)
                                goto readerr;
                        gd->ang *= PI/180.0;    /* convert to radians */
                        gd->next = all_dirs;
                        all_dirs = gd;
                        break;
                }
        return;
memerr:
        perror(progname);
        exit(1);
readerr:
        fprintf(stderr, "%s: read error on input\n", progname);
        exit(1);
#undef S_SEARCH
#undef S_SOURCE
#undef S_DIREC
}


static void
print_values(           /* print out calculations */
        gdfun *func
)
{
        register struct glare_dir       *gd;

        for (gd = all_dirs; gd != NULL; gd = gd->next)
                printf("%f\t%f\n", gd->ang*(180.0/PI), (*func)(gd));
}


static double
direct(                 /* compute direct vertical illuminance */
        struct glare_dir        *gd
)
{
        FVECT   mydir;
        double  d, dval;
        register struct glare_src       *gs;

        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
        dval = 0.0;
        for (gs = all_srcs; gs != NULL; gs = gs->next) {
                d = DOT(mydir,gs->dir);
                if (d > FTINY)
                        dval += d * gs->dom * gs->lum;
        }
        return(dval);
}


static double
indirect(                       /* return indirect vertical illuminance */
        struct glare_dir        *gd
)
{
        return(gd->indirect);
}


static double
total(                  /* return total vertical illuminance */
        struct glare_dir        *gd
)
{
        return(direct(gd)+gd->indirect);
}


/*
 * posindex -   compute glare position index from:
 *
 *      Source Direction
 *      View Direction
 *      View Up Direction
 *
 * All vectors are assumed to be normalized.
 * This function is an implementation of the method proposed by
 * Robert Levin in his 1975 JIES article.
 * This calculation presumes the view direction and up vectors perpendicular.
 * We return a value less than zero for improper positions.
 */

static double
posindex(                       /* compute position index */
        FVECT   sd,
        FVECT   vd,
        FVECT   vu
)
{
        double  sigma, tau;
        double  d;

        d = DOT(sd,vd);
        if (d <= 0.0)
                return(-1.0);
        if (d >= 1.0)
                return(1.0);
        sigma = acos(d) * (180./PI);
        d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
        if (d >= 1.0)
                tau = 0.0;
        else
                tau = acos(d) * (180./PI);
        return( exp( sigma*( (35.2 - tau*.31889 - 1.22*exp(-.22222*tau))*1e-3
                        + sigma*(21. + tau*(.26667 + tau*-.002963))*1e-5 )
                ) );
}


static double
dgi(            /* compute Daylight Glare Index */
        struct glare_dir        *gd
)
{
        register struct glare_src       *gs;
        FVECT   mydir,testdir[7],vhor;
        double  r,posn,omega,p[7],sum;
        int     i,n;

        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
        sum = 0.0; n = 0;
        for (gs = all_srcs; gs != NULL; gs = gs->next) {

                /* compute 1/p^2 weighted solid angle of the source */
                r = sqrt(1 - pow(1.-gs->dom/2./PI,2.));
                fcross(vhor,gs->dir,midview.vup);
                normalize(vhor);
                VCOPY(testdir[0],gs->dir);
                fvsum(testdir[1],gs->dir,vhor,r);
                fvsum(testdir[2],gs->dir,vhor,0.5*r);
                fvsum(testdir[5],testdir[2],midview.vup,-0.866*r);
                fvsum(testdir[2],testdir[2],midview.vup,0.866*r);
                fvsum(testdir[3],gs->dir,vhor,-r);
                fvsum(testdir[4],gs->dir,vhor,-0.5*r);
                fvsum(testdir[6],testdir[4],midview.vup,0.866*r);
                fvsum(testdir[4],testdir[4],midview.vup,-0.866*r);
                for (i = 0; i < 7; i++) {
                        normalize(testdir[i]);
                        posn = posindex(testdir[i],mydir,midview.vup);
                        if (posn <= FTINY)
                                p[i] = 0.0;
                        else
                                p[i] = 1./(posn*posn);
                }
                r = 1-gs->dom/2./PI;
                omega = gs->dom*p[0];
                omega += (r*PI*(1+1/r/r)-2*PI)*(-p[0]+(p[1]+p[2])*0.5);
                omega += (2*PI-r*PI*(1+1/r/r))*(-p[0]-0.1667*(p[1]+p[3])
                          +0.3334*(p[2]+p[4]+p[5]+p[6]));

                sum += pow(gs->lum,1.6) * pow(omega,0.8) /
                       (gd->indirect/PI + 0.07*sqrt(gs->dom)*gs->lum);
                n++;
        }
        if (n == 0)
                return(0.0);
        return( 10*log10(0.478*sum) );
}


static double
brs_gi(         /* compute BRS Glare Index */
        struct glare_dir        *gd
)
{
        register struct glare_src       *gs;
        FVECT   mydir;
        double  p;
        double  sum;

        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
        sum = 0.0;
        for (gs = all_srcs; gs != NULL; gs = gs->next) {
                p = posindex(gs->dir, mydir, midview.vup);
                if (p <= FTINY)
                        continue;
                sum += pow(gs->lum/p,1.6) * pow(gs->dom,0.8);
        }
        if (sum <= FTINY)
                return(0.0);
        sum /= gd->indirect/PI;
        return(10*log10(0.478*sum));
}


static double
guth_dgr(               /* compute Guth discomfort glare rating */
        struct glare_dir        *gd
)
{
#define q(w)    (20.4*w+1.52*pow(w,.2)-.075)
        register struct glare_src       *gs;
        FVECT   mydir;
        double  p;
        double  sum;
        double  wtot, brsum;
        int     n;

        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
        sum = wtot = brsum = 0.0; n = 0;
        for (gs = all_srcs; gs != NULL; gs = gs->next) {
                p = posindex(gs->dir, mydir, midview.vup);
                if (p <= FTINY)
                        continue;
                sum += gs->lum * q(gs->dom) / p;
                brsum += gs->lum * gs->dom;
                wtot += gs->dom;
                n++;
        }
        if (n == 0)
                return(0.0);
        return( pow(.5*sum/pow((brsum+(5.-wtot)*gd->indirect/PI)/5.,.44),
                        pow((double)n, -.0914) ) );
#undef q
}


#ifndef M_SQRT2
#define M_SQRT2 1.41421356237309504880
#endif

#define norm_integral(z)        (1.-.5*erfc((z)/M_SQRT2))


static double
guth_vcp(               /* compute Guth visual comfort probability */
        struct glare_dir        *gd
)
{
        extern double   erfc();
        double  dgr;

        dgr = guth_dgr(gd);
        if (dgr <= FTINY)
                return(100.0);
        return(100.*norm_integral(6.374-1.3227*log(dgr)));
}


static double
cie_cgi(                /* compute CIE Glare Index */
        struct glare_dir        *gd
)
{
        register struct glare_src       *gs;
        FVECT   mydir;
        double  dillum;
        double  p;
        double  sum;

        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
        sum = 0.0;
        for (gs = all_srcs; gs != NULL; gs = gs->next) {
                p = posindex(gs->dir, mydir, midview.vup);
                if (p <= FTINY)
                        continue;
                sum += gs->lum*gs->lum * gs->dom / (p*p);
        }
        if (sum <= FTINY)
                return(0.0);
        dillum = direct(gd);
        return(8.*log10(2.*sum*(1.+dillum/500.)/(dillum+gd->indirect)));
}


static double
ugr(            /* compute Unified Glare Rating */
        struct glare_dir        *gd
)
{
        register struct glare_src       *gs;
        FVECT   mydir;
        double  p;
        double  sum;

        spinvector(mydir, midview.vdir, midview.vup, gd->ang);
        sum = 0.0;
        for (gs = all_srcs; gs != NULL; gs = gs->next) {
                p = posindex(gs->dir, mydir, midview.vup);
                if (p <= FTINY)
                        continue;
                sum += gs->lum*gs->lum * gs->dom / (p*p);
        }
        if (sum <= FTINY)
                return(0.0);
        return(8.*log10(0.25*sum*PI/gd->indirect));
}
Revision:	2.10
Committed:	Fri Jan 2 12:48:36 2004 UTC (20 years, 3 months ago) by schorsch
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad3R7P2, rad3R7P1, rad4R1, rad4R0, rad3R6, rad3R6P1, rad3R8, rad3R9
Changes since 2.9:	+143 -109 lines
Log Message:	Ansification and fixed typing/prototype of getheader() and its callback.
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	schorsch	2.10	static const char RCSid[] = "$Id: glarendx.c,v 2.9 2003/11/19 16:21:28 greg Exp $";
3	greg	1.1	#endif
4			/*
5	greg	1.9	* Compute Glare Index given by program name or -t option:
6	greg	1.1	*
7	greg	2.4	* dgi - Daylight Glare Index
8			* brs_gi - Building Research Station Glare Index (Petherbridge
9			* & Hopkinson)
10			* ugr - Unified Glare Rating System (Fischer)
11	greg	1.2	* guth_dgr - Guth discomfort glare rating
12			* guth_vcp - Guth visual comfort probability
13	greg	1.9	* cie_cgi - CIE Glare Index (1983, due to Einhorn)
14			* vert_dir - Direct vertical illuminance
15			* vert_ind - Indirect vertical illuminance (from input)
16			* vert_ill - Total vertical illuminance
17	greg	1.1	*
18			* 12 April 1991 Greg Ward EPFL
19	greg	2.4	* 19 April 1993 R. Compagnon EPFL (added dgi, brs_gi, ugr)
20	greg	1.1	*/
21	schorsch	2.10
22	schorsch	2.8	#include <string.h>
23
24	greg	1.1	#include "standard.h"
25			#include "view.h"
26	schorsch	2.10
27
28			struct glare_src {
29			FVECT dir; /* source direction */
30			double dom; /* solid angle */
31			double lum; /* average luminance */
32			struct glare_src *next;
33			} *all_srcs = NULL;
34
35			struct glare_dir {
36			double ang; /* angle (in radians) */
37			double indirect; /* indirect illuminance */
38			struct glare_dir *next;
39			} *all_dirs = NULL;
40
41			typedef double gdfun(struct glare_dir *gd);
42
43			static gdfun dgi;
44			static gdfun brs_gi;
45			static gdfun ugr;
46			static gdfun guth_dgr;
47			static gdfun guth_vcp;
48			static gdfun cie_cgi;
49			static gdfun direct;
50			static gdfun indirect;
51			static gdfun total;
52
53			static gethfunc headline;
54			static void init(void);
55			static void read_input(void);
56			static void print_values(gdfun *func);
57			static double posindex(FVECT sd, FVECT vd, FVECT vu);
58
59
60	greg	1.1	struct named_func {
61			char *name;
62	schorsch	2.10	gdfun *func;
63	greg	1.9	char *descrip;
64	greg	1.1	} all_funcs[] = {
65	greg	2.4	{"dgi", dgi, "Daylight Glare Index"},
66			{"brs_gi", brs_gi, "BRS Glare Index"},
67			{"ugr", ugr, "Unified Glare Rating"},
68	greg	1.9	{"guth_vcp", guth_vcp, "Guth Visual Comfort Probability"},
69			{"cie_cgi", cie_cgi, "CIE Glare Index (Einhorn)"},
70			{"guth_dgr", guth_dgr, "Guth Disability Glare Rating"},
71			{"vert_dir", direct, "Direct Vertical Illuminance"},
72			{"vert_ill", total, "Total Vertical Illuminance"},
73			{"vert_ind", indirect, "Indirect Vertical Illuminance"},
74	greg	1.1	{NULL}
75			};
76	schorsch	2.10
77	greg	1.1	#define newp(type) (type *)malloc(sizeof(type))
78	schorsch	2.10
79	greg	1.1	char *progname;
80			int print_header = 1;
81	schorsch	2.10
82	greg	1.1	VIEW midview = STDVIEW;
83	schorsch	2.10
84	greg	1.4	int wrongformat = 0;
85	schorsch	2.10
86
87			int
88			main(
89			int argc,
90			char *argv[]
91			)
92	greg	1.1	{
93			struct named_func *funp;
94			char *progtail;
95			int i;
96			/* get program name */
97			progname = argv[0];
98	schorsch	2.8	progtail = strrchr(progname, '/'); /* final component */
99	greg	1.1	if (progtail == NULL)
100			progtail = progname;
101			else
102			progtail++;
103			/* get options */
104			for (i = 1; i < argc && argv[i][0] == '-'; i++)
105			switch (argv[i][1]) {
106			case 't':
107			progtail = argv[++i];
108			break;
109			case 'h':
110			print_header = 0;
111			break;
112			default:
113			goto userr;
114			}
115			if (i < argc-1)
116			goto userr;
117			if (i == argc-1) /* open file */
118			if (freopen(argv[i], "r", stdin) == NULL) {
119			perror(argv[i]);
120			exit(1);
121			}
122	greg	1.9	/* find and run calculation */
123	greg	1.1	for (funp = all_funcs; funp->name != NULL; funp++)
124			if (!strcmp(funp->name, progtail)) {
125	greg	1.9	init();
126			read_input();
127			if (print_header) {
128			printargs(i, argv, stdout);
129			putchar('\n');
130			}
131	greg	1.1	print_values(funp->func);
132			exit(0); /* we're done */
133			}
134			/* invalid function */
135			userr:
136	greg	1.9	fprintf(stderr, "Usage: %s -t type [-h] [input]\n", progname);
137			fprintf(stderr, "\twhere 'type' is one of the following:\n");
138			for (funp = all_funcs; funp->name != NULL; funp++)
139			fprintf(stderr, "\t%12s\t%s\n", funp->name, funp->descrip);
140	greg	1.1	exit(1);
141			}
142	schorsch	2.10
143
144			static int
145			headline( /* get line from header */
146			char *s,
147			void *p
148			)
149	greg	1.1	{
150	greg	1.4	char fmt[32];
151	schorsch	2.10
152	greg	1.1	if (print_header) /* copy to output */
153			fputs(s, stdout);
154	greg	2.2	if (isview(s))
155			sscanview(&midview, s);
156	greg	1.4	else if (isformat(s)) {
157			formatval(fmt, s);
158	greg	1.6	wrongformat = strcmp(fmt, "ascii");
159	greg	1.4	}
160	gwlarson	2.6	return(0);
161	greg	1.1	}
162	schorsch	2.10
163
164			static void
165			init(void) /* initialize calculation */
166	greg	1.9	{
167			/* read header */
168			getheader(stdin, headline, NULL);
169			if (wrongformat) {
170			fprintf(stderr, "%s: bad input format\n", progname);
171			exit(1);
172			}
173			/* set view */
174			if (setview(&midview) != NULL) {
175			fprintf(stderr, "%s: bad view information in input\n",
176			progname);
177			exit(1);
178			}
179			}
180	schorsch	2.10
181
182			static void
183			read_input(void) /* read glare sources from stdin */
184	greg	1.1	{
185			#define S_SEARCH 0
186			#define S_SOURCE 1
187			#define S_DIREC 2
188			int state = S_SEARCH;
189			char buf[128];
190			register struct glare_src *gs;
191			register struct glare_dir *gd;
192	schorsch	2.10
193	greg	1.1	while (fgets(buf, sizeof(buf), stdin) != NULL)
194			switch (state) {
195			case S_SEARCH:
196			if (!strcmp(buf, "BEGIN glare source\n"))
197			state = S_SOURCE;
198			else if (!strcmp(buf, "BEGIN indirect illuminance\n"))
199			state = S_DIREC;
200			break;
201			case S_SOURCE:
202			if (!strncmp(buf, "END", 3)) {
203			state = S_SEARCH;
204			break;
205			}
206			if ((gs = newp(struct glare_src)) == NULL)
207			goto memerr;
208			if (sscanf(buf, "%lf %lf %lf %lf %lf",
209			&gs->dir[0], &gs->dir[1], &gs->dir[2],
210			&gs->dom, &gs->lum) != 5)
211			goto readerr;
212	greg	1.2	normalize(gs->dir);
213	greg	1.1	gs->next = all_srcs;
214			all_srcs = gs;
215			break;
216			case S_DIREC:
217			if (!strncmp(buf, "END", 3)) {
218			state = S_SEARCH;
219			break;
220			}
221			if ((gd = newp(struct glare_dir)) == NULL)
222			goto memerr;
223			if (sscanf(buf, "%lf %lf",
224			&gd->ang, &gd->indirect) != 2)
225			goto readerr;
226			gd->ang = PI/180.0; / convert to radians */
227			gd->next = all_dirs;
228			all_dirs = gd;
229			break;
230			}
231			return;
232			memerr:
233			perror(progname);
234			exit(1);
235			readerr:
236			fprintf(stderr, "%s: read error on input\n", progname);
237			exit(1);
238			#undef S_SEARCH
239			#undef S_SOURCE
240			#undef S_DIREC
241			}
242	schorsch	2.10
243
244			static void
245			print_values( /* print out calculations */
246			gdfun *func
247			)
248	greg	1.1	{
249			register struct glare_dir *gd;
250	schorsch	2.10
251	greg	1.1	for (gd = all_dirs; gd != NULL; gd = gd->next)
252	schorsch	2.10	printf("%f\t%f\n", gd->ang(180.0/PI), (func)(gd));
253	greg	1.1	}
254	schorsch	2.10
255
256			static double
257			direct( /* compute direct vertical illuminance */
258			struct glare_dir *gd
259			)
260	greg	1.1	{
261			FVECT mydir;
262			double d, dval;
263			register struct glare_src *gs;
264	schorsch	2.10
265	greg	1.1	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
266			dval = 0.0;
267			for (gs = all_srcs; gs != NULL; gs = gs->next) {
268			d = DOT(mydir,gs->dir);
269			if (d > FTINY)
270			dval += d * gs->dom * gs->lum;
271			}
272			return(dval);
273	greg	1.9	}
274	schorsch	2.10
275
276			static double
277			indirect( /* return indirect vertical illuminance */
278			struct glare_dir *gd
279			)
280	greg	1.9	{
281			return(gd->indirect);
282			}
283	schorsch	2.10
284
285			static double
286			total( /* return total vertical illuminance */
287			struct glare_dir *gd
288			)
289	greg	1.9	{
290			return(direct(gd)+gd->indirect);
291	greg	1.1	}
292	schorsch	2.10
293
294	greg	1.1	/*
295			* posindex - compute glare position index from:
296			*
297			* Source Direction
298			* View Direction
299			* View Up Direction
300			*
301			* All vectors are assumed to be normalized.
302			* This function is an implementation of the method proposed by
303			* Robert Levin in his 1975 JIES article.
304	greg	1.5	* This calculation presumes the view direction and up vectors perpendicular.
305	greg	1.1	* We return a value less than zero for improper positions.
306			*/
307	schorsch	2.10
308			static double
309			posindex( /* compute position index */
310			FVECT sd,
311			FVECT vd,
312			FVECT vu
313			)
314	greg	1.1	{
315			double sigma, tau;
316			double d;
317	schorsch	2.10
318	greg	1.1	d = DOT(sd,vd);
319			if (d <= 0.0)
320			return(-1.0);
321	greg	1.5	if (d >= 1.0)
322			return(1.0);
323	greg	1.1	sigma = acos(d) * (180./PI);
324	greg	2.4	d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
325	greg	1.5	if (d >= 1.0)
326			tau = 0.0;
327			else
328			tau = acos(d) * (180./PI);
329	greg	1.1	return( exp( sigma( (35.2 - tau.31889 - 1.22exp(-.22222tau))*1e-3
330			+ sigma(21. + tau(.26667 + tau-.002963))1e-5 )
331			) );
332			}
333	schorsch	2.10
334
335			static double
336			dgi( /* compute Daylight Glare Index */
337			struct glare_dir *gd
338			)
339	greg	2.4	{
340			register struct glare_src *gs;
341			FVECT mydir,testdir[7],vhor;
342	greg	2.9	double r,posn,omega,p[7],sum;
343	greg	2.4	int i,n;
344	schorsch	2.10
345	greg	2.4	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
346			sum = 0.0; n = 0;
347			for (gs = all_srcs; gs != NULL; gs = gs->next) {
348	greg	1.1
349	greg	2.4	/* compute 1/p^2 weighted solid angle of the source */
350			r = sqrt(1 - pow(1.-gs->dom/2./PI,2.));
351			fcross(vhor,gs->dir,midview.vup);
352			normalize(vhor);
353			VCOPY(testdir[0],gs->dir);
354			fvsum(testdir[1],gs->dir,vhor,r);
355			fvsum(testdir[2],gs->dir,vhor,0.5*r);
356			fvsum(testdir[5],testdir[2],midview.vup,-0.866*r);
357			fvsum(testdir[2],testdir[2],midview.vup,0.866*r);
358			fvsum(testdir[3],gs->dir,vhor,-r);
359			fvsum(testdir[4],gs->dir,vhor,-0.5*r);
360			fvsum(testdir[6],testdir[4],midview.vup,0.866*r);
361			fvsum(testdir[4],testdir[4],midview.vup,-0.866*r);
362			for (i = 0; i < 7; i++) {
363			normalize(testdir[i]);
364	greg	2.9	posn = posindex(testdir[i],mydir,midview.vup);
365			if (posn <= FTINY)
366			p[i] = 0.0;
367			else
368			p[i] = 1./(posn*posn);
369	greg	2.4	}
370			r = 1-gs->dom/2./PI;
371			omega = gs->dom*p[0];
372			omega += (rPI(1+1/r/r)-2PI)(-p[0]+(p[1]+p[2])*0.5);
373			omega += (2PI-rPI(1+1/r/r))(-p[0]-0.1667*(p[1]+p[3])
374			+0.3334*(p[2]+p[4]+p[5]+p[6]));
375	greg	1.1
376	schorsch	2.10	sum += pow(gs->lum,1.6) * pow(omega,0.8) /
377	greg	2.4	(gd->indirect/PI + 0.07sqrt(gs->dom)gs->lum);
378			n++;
379			}
380			if (n == 0)
381			return(0.0);
382			return( 10log10(0.478sum) );
383			}
384
385	schorsch	2.10
386			static double
387			brs_gi( /* compute BRS Glare Index */
388			struct glare_dir *gd
389			)
390	greg	2.4	{
391			register struct glare_src *gs;
392			FVECT mydir;
393			double p;
394			double sum;
395	schorsch	2.10
396	greg	2.4	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
397			sum = 0.0;
398			for (gs = all_srcs; gs != NULL; gs = gs->next) {
399			p = posindex(gs->dir, mydir, midview.vup);
400			if (p <= FTINY)
401			continue;
402			sum += pow(gs->lum/p,1.6) * pow(gs->dom,0.8);
403			}
404			if (sum <= FTINY)
405			return(0.0);
406			sum /= gd->indirect/PI;
407			return(10log10(0.478sum));
408			}
409
410
411	schorsch	2.10	static double
412			guth_dgr( /* compute Guth discomfort glare rating */
413			struct glare_dir *gd
414			)
415	greg	1.1	{
416			#define q(w) (20.4w+1.52pow(w,.2)-.075)
417			register struct glare_src *gs;
418	greg	1.7	FVECT mydir;
419	greg	1.1	double p;
420			double sum;
421	greg	1.3	double wtot, brsum;
422	greg	1.1	int n;
423	schorsch	2.10
424	greg	1.7	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
425	greg	1.3	sum = wtot = brsum = 0.0; n = 0;
426	greg	1.1	for (gs = all_srcs; gs != NULL; gs = gs->next) {
427	greg	1.7	p = posindex(gs->dir, mydir, midview.vup);
428	greg	1.1	if (p <= FTINY)
429			continue;
430			sum += gs->lum * q(gs->dom) / p;
431	greg	1.3	brsum += gs->lum * gs->dom;
432			wtot += gs->dom;
433	greg	1.1	n++;
434			}
435			if (n == 0)
436			return(0.0);
437	greg	1.8	return( pow(.5sum/pow((brsum+(5.-wtot)gd->indirect/PI)/5.,.44),
438	greg	1.1	pow((double)n, -.0914) ) );
439			#undef q
440			}
441	schorsch	2.10
442
443	greg	1.1	#ifndef M_SQRT2
444			#define M_SQRT2 1.41421356237309504880
445			#endif
446	schorsch	2.10
447	greg	1.1	#define norm_integral(z) (1.-.5*erfc((z)/M_SQRT2))
448	schorsch	2.10
449
450			static double
451			guth_vcp( /* compute Guth visual comfort probability */
452			struct glare_dir *gd
453			)
454	greg	1.1	{
455	greg	2.5	extern double erfc();
456	greg	1.7	double dgr;
457	schorsch	2.10
458	greg	1.7	dgr = guth_dgr(gd);
459			if (dgr <= FTINY)
460			return(100.0);
461			return(100.norm_integral(6.374-1.3227log(dgr)));
462	greg	1.8	}
463	schorsch	2.10
464
465			static double
466			cie_cgi( /* compute CIE Glare Index */
467			struct glare_dir *gd
468			)
469	greg	1.8	{
470			register struct glare_src *gs;
471			FVECT mydir;
472			double dillum;
473			double p;
474			double sum;
475	schorsch	2.10
476	greg	1.8	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
477			sum = 0.0;
478			for (gs = all_srcs; gs != NULL; gs = gs->next) {
479			p = posindex(gs->dir, mydir, midview.vup);
480			if (p <= FTINY)
481			continue;
482			sum += gs->lumgs->lum gs->dom / (p*p);
483			}
484			if (sum <= FTINY)
485			return(0.0);
486			dillum = direct(gd);
487			return(8.log10(2.sum*(1.+dillum/500.)/(dillum+gd->indirect)));
488	greg	2.4	}
489
490
491	schorsch	2.10	static double
492			ugr( /* compute Unified Glare Rating */
493			struct glare_dir *gd
494			)
495	greg	2.4	{
496			register struct glare_src *gs;
497			FVECT mydir;
498			double p;
499			double sum;
500	schorsch	2.10
501	greg	2.4	spinvector(mydir, midview.vdir, midview.vup, gd->ang);
502			sum = 0.0;
503			for (gs = all_srcs; gs != NULL; gs = gs->next) {
504			p = posindex(gs->dir, mydir, midview.vup);
505			if (p <= FTINY)
506			continue;
507			sum += gs->lumgs->lum gs->dom / (p*p);
508			}
509			if (sum <= FTINY)
510			return(0.0);
511			return(8.log10(0.25sum*PI/gd->indirect));
512	greg	1.1	}