src/util/glarendx.c

#ifndef lint
static const char       RCSid[] = "$Id: glarendx.c,v 2.16 2020/07/24 16:58:16 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 "paths.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))

int     print_header = 1;

VIEW    midview = STDVIEW;

int     wrongformat = 0;


int
main(
        int     argc,
        char    *argv[]
)
{
        struct named_func       *funp;
        int     i;
                                        /* get program name */
        progname = fixargv0(argv[0]);
                                        /* get options */
        for (i = 1; i < argc && argv[i][0] == '-'; i++)
                switch (argv[i][1]) {
                case 't':
                        progname = 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, progname)) {
                        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[MAXFMTLEN];

        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 */
{
        char    *err;
                                        /* read header */
        getheader(stdin, headline, NULL);
        if (wrongformat) {
                fprintf(stderr, "%s: bad input format\n", progname);
                exit(1);
        }
                                        /* set view */
        if ((err = setview(&midview)) != NULL) {
                fprintf(stderr, "%s: %s\n",
                                progname, err);
                exit(1);
        }
        if (fabs(DOT(midview.vdir, midview.vup)) > 1e-4) {
                fprintf(stderr,
                "%s: Input view direction not perpendicular to up vector\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];
        struct glare_src        *gs;
        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, FVFORMAT, &gs->dir[0], &gs->dir[1],
                                                &gs->dir[2]) != 3 ||
                                        sscanf(sskip2(buf, 3), "%lf %lf",
                                                &gs->dom, &gs->lum) != 2)
                                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
)
{
        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;
        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-FTINY)
                return(1.0);
        sigma = acos(d) * (180./PI);
        d = fabs(DOT(sd,vu)/sqrt(1.0-d*d));
        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
)
{
        struct glare_src        *gs;
        FVECT   mydir,testdir[7],vhor;
        double  r,posn,omega,p[7],sum;
        int     i;

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

                /* compute 1/p^2 weighted solid angle of the source */
                r = 1. - gs->dom*(.5/PI);
                r = sqrt(1. - r*r);
                fcross(vhor,gs->dir,midview.vup);
                if (normalize(vhor) == 0.)
                        continue;
                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++) {
                        p[i] = 0.0;
                        if (normalize(testdir[i]) == 0.)
                                continue;
                        posn = posindex(testdir[i],mydir,midview.vup);
                        if (posn > FTINY)
                                p[i] = 1./(posn*posn);
                }
                r = 1. - gs->dom*(.5/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]));
                if (omega <= 0.)
                        continue;
                sum += pow(gs->lum,1.6) * pow(omega,0.8) /
                       (gd->indirect*(1./PI) + 0.07*sqrt(gs->dom)*gs->lum);
        }
        if (sum <= FTINY)
                return(0.0);
        return( 10.*log10(0.478*sum) );
}


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