src/util/cmglare.c

/*
 *  cmglare.c - routines for calculating glare autonomy.
 *
 *  N. Jones
 */

/*
 * Copyright (c) 2017-2019 Nathaniel Jones
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in all
 * copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 * SOFTWARE.
 */

#include "rtmath.h"
#include "cmglare.h"

#define LUMINOUS_EFFICACY       179     /* lumens per Watt */
#define LUMINANCE_THRESHOLD     100 /* minimum threshold that will be interpreted as luminance rather than ratio between Ev and glare source */

#define ANGLE(u,v)      acos(DOT((u),(v)))      // TODO need to normalize?

typedef struct reinhart_sky {
        //int mf;                                       /* Linear divisions per Tregenza patch. */
        int rings;                                      /* Number of rings of sky patches. */
        //int patches;                          /* Number of sky patches. */
        double ringElevationAngle;      /* Angle in radians between rings. */
        int *patchesPerRow;                     /* Sky patches per row. */
        int *firstPatchIndex;           /* Index of first sky patch in each row. */
        double *solidAngle;                     /* Solid angle of each patch in each row. */
} ReinhartSky;

static const int tnaz[] = { 30, 30, 24, 24, 18, 12, 6 };        /* Number of patches per row */

extern char* progname;

static ReinhartSky* make_sky(const CMATRIX *smx)
{
        int mf, count, i, j;
        double remaining = PI / 2;

        /* Check sky partitionss */
        switch (smx->nrows) { // nrows is Reinhart sky subdivisoins plus one for miss
        case 2: mf = 0; break;
        case 146: mf = 1; break;
        case 578: mf = 2; break;
        case 1298: mf = 3; break;
        case 2306: mf = 4; break;
        case 3602: mf = 5; break;
        case 5186: mf = 6; break;
        case 7058: mf = 7; break;
        case 9218: mf = 8; break;
        case 11666: mf = 9; break;
        case 14402: mf = 10; break;
        case 17426: mf = 11; break;
        case 20738: mf = 12; break;
        default:
                fprintf(stderr,
                        "%s: unknown number of sky patches %d\n",
                        progname, smx->nrows);
                return NULL;
        }

        /* Allocate sky */
        ReinhartSky *sky = (ReinhartSky*)malloc(sizeof(ReinhartSky));
        if (!sky) goto memerr;

        /* Calculate patches per row */
        sky->rings = 7 * mf + 1;
        sky->patchesPerRow = (int*)malloc(sky->rings * sizeof(int));
        sky->firstPatchIndex = (int*)malloc(sky->rings * sizeof(int));
        if (!sky->patchesPerRow || !sky->firstPatchIndex) goto memerr;
        count = 1; // The below horizon patch
        for (i = 0; i < 7; i++)
                for (j = 0; j < mf; j++) {
                        sky->firstPatchIndex[i * mf + j] = count;
                        count += sky->patchesPerRow[i * mf + j] = tnaz[i] * mf;
                }
        sky->firstPatchIndex[7 * mf] = count;
        sky->patchesPerRow[7 * mf] = 1;
        //sky->patches = count + 1;

        /* Calculate solid angle of patches in each row */
        sky->solidAngle = (double*)malloc(sky->rings * sizeof(double));
        if (!sky->solidAngle) goto memerr;
        sky->ringElevationAngle = PI / (2 * sky->rings - 1);
        sky->solidAngle[0] = 2 * PI;
        for (i = 1; i < sky->rings; i++) {
                remaining -= sky->ringElevationAngle;
                sky->solidAngle[i] = 2 * PI * (1 - cos(remaining)); // solid angle of cap
                sky->solidAngle[i - 1] -= sky->solidAngle[i];
                sky->solidAngle[i - 1] /= sky->patchesPerRow[i - 1];
        }

        return sky;

memerr:
        fprintf(stderr,
                "%s: out of memory\n",
                progname);
        return NULL;
}

static void free_sky(ReinhartSky *sky)
{
        if (sky) {
                free(sky->patchesPerRow);
                free(sky->firstPatchIndex);
                free(sky->solidAngle);
                free(sky);
        }
}

static void get_patch_direction(const ReinhartSky *sky, const int patch, FVECT target)
{
        //if (patch >= sky->patches || patch < 0) throw new RuntimeException("Illegal patch " + patch);
        if (!patch) {
                target[0] = target[1] = 0;
                target[2] = -1;
                return; // Ignore below horizon?
        }
        int row = sky->rings - 1;
        while (patch < sky->firstPatchIndex[row]) row--;
        const double alt = PI / 2 - sky->ringElevationAngle * (row - (sky->rings - 1));
        const double azi = 2 * PI * (patch - sky->firstPatchIndex[row]) / sky->patchesPerRow[row];
        const double cos_alt = cos(alt);
        target[0] = cos_alt * -sin(azi);
        target[1] = cos_alt * -cos(azi);
        target[2] = sin(alt);
}

static double get_patch_solid_angle(const ReinhartSky *sky, const int patch, const double cos_theta)
{
        //if (patch >= sky->patches || patch < 0) throw new RuntimeException("Illegal patch " + patch);
        if (!patch) return 2 * (PI - 2 * acos(cos_theta)); // Solid angle overlap between visible hemisphere and ground hemisphere
        int row = sky->rings - 1;
        while (patch < sky->firstPatchIndex[row]) row--;
        return sky->solidAngle[row];
}

static double get_guth(const FVECT dir, const FVECT forward, const FVECT up)
{
        double posindex;
        FVECT hv, temp;
        VCROSS(hv, forward, up);
        normalize(hv);
        VCROSS(temp, forward, hv);
        double phi = ANGLE(dir, temp) - PI / 2;

        /* Guth model, equation from IES lighting handbook */
        if (phi >= 0) {
                double sigma = ANGLE(dir, forward);
                VSUM(hv, forward, dir, 1 / DOT(dir, forward));
                normalize(hv);
                double tau = ANGLE(hv, up);
                tau *= 180.0 / PI;
                sigma *= 180.0 / PI;

                if (sigma <= 0)
                        sigma = -sigma;

                posindex = exp((35.2 - 0.31889 * tau - 1.22 * exp(-2.0 * tau / 9.0)) / 1000.0 * sigma + (21.0 + 0.26667 * tau - 0.002963 * tau * tau) / 100000.0 * sigma * sigma);
        }
        /* below line of sight, using Iwata model */
        else {
                double teta = PI / 2 - ANGLE(dir, hv);

                if (teta == 0.0)
                        teta = FTINY;

                double fact = 0.8;
                double d = 1 / tan(phi);
                double s = tan(teta) / tan(phi);
                double r = sqrt((1 + s * s) / (d * d));
                if (r > 0.6)
                        fact = 1.2;
                if (r > 3)
                        r = 3.0;

                posindex = 1.0 + fact * r;
        }
        if (posindex > 16)
                posindex = 16.0;

        return posindex;
}

float* cm_glare(const CMATRIX *dcmx, const CMATRIX *evmx, const CMATRIX *smx, const int *occupied, const double dgp_limit, const double dgp_threshold, const FVECT *views, const FVECT dir, const FVECT up)
{
        int p, t, c;
        int hourly_output = dgp_limit < 0;
        float *dgp_list;
        ReinhartSky *sky;
        FVECT vdir;

        /* Check consistancy */
        if ((dcmx->nrows != evmx->nrows) | (dcmx->ncols != smx->nrows) | (evmx->ncols != smx->ncols)) {
                fprintf(stderr,
                        "%s: inconsistant matrix dimensions: dc(%d, %d) ev(%d, %d) s(%d, %d)\n",
                        progname, dcmx->nrows, dcmx->ncols, evmx->nrows, evmx->ncols, smx->nrows, smx->ncols);
                return NULL;
        }

        /* Create output buffer */
        dgp_list = (float*)malloc(evmx->nrows * (hourly_output ? evmx->ncols : 1) * sizeof(float));
        if (!dgp_list) {
                fprintf(stderr,
                        "%s: out of memory in cm_glare()\n",
                        progname);
                return NULL;
        }

        /* Create sky */
        sky = make_sky(smx);
        if (!sky) return NULL;

        /* Calculate glare limit */
        double ev_max = -1;
        if (!hourly_output) {
                ev_max = (dgp_limit - 0.16) / 5.87e-5;
                if (ev_max < 0) ev_max = 0;
        }

        /* For each position and direction */
        if (!views) VCOPY(vdir, dir);
        for (p = 0; p < evmx->nrows; p++) {
                /* For each time step */
                int occupied_hours = 0;
                int glare_hours = 0;
                if (views) VCOPY(vdir, views[p]);
                for (t = 0; t < evmx->ncols; t++) {
                        if (!occupied[t]) {
                                /* Not occupied */
                                if (hourly_output) dgp_list[p * evmx->ncols + t] = 0.0f;
                        }
                        else {
                                /* Occupied */
                                double illum = LUMINOUS_EFFICACY * bright(cm_lval(evmx, p, t));
                                occupied_hours++;

                                if (illum <= FTINY) {
                                        /* No light, so no glare */
                                        if (hourly_output) dgp_list[p * evmx->ncols + t] = 0.0f;
                                }
                                else if ((illum >= ev_max) & (!hourly_output)) {
                                        /* Guarangeed glare */
                                        glare_hours++;
                                }
                                else {
                                        /* Calculate enhanced simplified daylight glare probability */
                                        double sum = 0.0;
                                        FVECT patch_normal;
                                        for (c = 0; c < smx->nrows; c++) {
                                                const double dc = bright(cm_lval(dcmx, p, c));
                                                if (dc > 0) {
                                                        get_patch_direction(sky, c, patch_normal);
                                                        if (!c) {
                                                                /* Direction toward the center of the visible ground */
                                                                VADD(patch_normal, patch_normal, vdir);
                                                                if (normalize(patch_normal) == 0) continue;
                                                        }
                                                        const double cos_theta = DOT(vdir, patch_normal);
                                                        if (cos_theta <= FTINY) continue;
                                                        const double s = LUMINOUS_EFFICACY * bright(cm_lval(smx, c, t));
                                                        const double omega = get_patch_solid_angle(sky, c, cos_theta);
                                                        const double patch_luminance = (dc * s) / (omega * cos_theta);

                                                        double min_patch_luminance = dgp_threshold;
                                                        if (dgp_threshold < LUMINANCE_THRESHOLD)
                                                                min_patch_luminance *= illum / PI; // TODO should use average luminance, not illuminance
                                                        if (patch_luminance < min_patch_luminance) continue;
                                                        const double P = get_guth(patch_normal, vdir, up);
                                                        sum += (patch_luminance * patch_luminance * omega) / (P * P);
                                                }
                                        }

                                        //double dgp = 5.87e-5 * illum + 0.092 * log10(1 + dgp / pow(illum, 1.87)) + 0.159;
                                        double eDGPs = 5.87e-5 * illum + 0.0918 * log10(1 + sum / pow(illum, 1.87)) + 0.16;
                                        if (illum < 1000) /* low light correction */
                                                eDGPs *= exp(0.024 * illum - 4) / (1 + exp(0.024 * illum - 4));
                                        //eDGPs /= 1.1 - 0.5 * age / 100.0; /* age correction */
                                        if (eDGPs > 1.0) eDGPs = 1.0;

                                        if (hourly_output)
                                                dgp_list[p * evmx->ncols + t] = (float)eDGPs;
                                        else if (eDGPs >= dgp_limit)
                                                glare_hours++;
                                }
                        }
                }
                if (!hourly_output) {
                        /* Save glare autonomy */
                        dgp_list[p] = (float)(occupied_hours - glare_hours) / occupied_hours;
                }
        }

        free_sky(sky);

        return dgp_list;
}

static int getvec(FVECT vec, const int dtype, FILE *fp)         /* get a vector from fp */
{
        static float  vf[3];
        static double  vd[3];
        char  buf[32];
        int  i;

        switch (dtype) {
        case DTascii:
                for (i = 0; i < 3; i++) {
                        if (fgetword(buf, sizeof(buf), fp) == NULL ||
                                !isflt(buf))
                                return(-1);
                        vec[i] = atof(buf);
                }
                break;
        case DTfloat:
                if (getbinary(vf, sizeof(float), 3, fp) != 3)
                        return(-1);
                VCOPY(vec, vf);
                break;
        case DTdouble:
                if (getbinary(vd, sizeof(double), 3, fp) != 3)
                        return(-1);
                VCOPY(vec, vd);
                break;
        default:
                fprintf(stderr,
                        "%s: botched input format\n",
                        progname);
                return(-1);
        }
        return(0);
}

int cm_load_schedule(const int count, int* schedule, FILE *fp)
{
        char buf[512];
        char *cp;
        char *comma;
        double val;
        int i = 0;

        while (fgetline(buf, sizeof(buf), fp) != NULL) {
                if (buf[0] == '#') continue; // Comment line
                comma = NULL;
                for (cp = buf; *cp; cp++) {
                        /* If there are multiple commas, assume the value is after the last comma */
                        if (*cp == ',') {
                                comma = cp; /* Record position of last comma */
                        }
                }
                if (comma)
                        val = atof(comma + 1);
                else
                        val = atof(buf);

                if (i < count) {
                        /* Add the value to the schedule */
                        schedule[i++] = (val > 0);
                }
                else {
                        fprintf(stderr,
                                "%s: too many schedule entries\n",
                                progname);
                        return(1);
                }
        }
        fclose(fp);

        if (i < count) {
                fprintf(stderr,
                        "%s: too few schedule entries\n",
                        progname);
                return(-1);
        }

        return 0;
}

FVECT* cm_load_views(const int nrows, const int dtype, FILE *fp)
{
        int i;
        double d;
        FVECT orig;

        FVECT *views = (FVECT*)malloc(nrows * sizeof(FVECT));
        if (!views) {
                fprintf(stderr,
                        "%s: out of memory in cm_load_views()\n",
                        progname);
                return NULL;
        }

        for (i = 0; i < nrows; i++) {
                if (getvec(orig, dtype, fp) | getvec(views[i], dtype, fp)) {
                        fprintf(stderr,
                                "%s: unexpected end of input, missing %d entries\n",
                                progname, i);
                        return NULL;
                }

                d = normalize(views[i]);
                if (d == 0.0) {                         /* zero ==> flush */
                        fprintf(stderr,
                                "%s: zero length direction detected\n",
                                progname);
                        return NULL;
                }
        }

        return views;
}

int cm_write_glare(const float *mp, const int nrows, const int ncols, const int dtype, FILE *fp)
{
        static const char       tabEOL[2] = { '\t', '\n' };
        int                     r, c;
        double  dc[1];

        switch (dtype) {
        case DTascii:
                for (r = 0; r < nrows; r++)
                        for (c = 0; c < ncols; c++, mp++)
                                fprintf(fp, "%.6e%c",
                                mp[0],
                                tabEOL[c >= ncols - 1]);
                break;
        case DTfloat:
                r = ncols*nrows;
                while (r > 0) {
                        c = putbinary(mp, sizeof(float), r, fp);
                        if (c <= 0)
                                return(0);
                        mp += c;
                        r -= c;
                }
                break;
        case DTdouble:
                r = ncols*nrows;
                while (r--) {
                        dc[0] = mp[0];
                        if (putbinary(dc, sizeof(double), 1, fp) != 1)
                                return(0);
                        mp++;
                }
                break;
        default:
                fputs("Unsupported data type in cm_write_glare()!\n", stderr);
                return(0);
        }
        return(fflush(fp) == 0);
}
Revision:	2.1
Committed:	Mon Sep 9 17:19:51 2019 UTC (4 years, 7 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad5R4, rad5R3, HEAD
Log Message:	Added Nathaniel Jones' dcglare utility
#	User	Rev	Content
1	greg	2.1	/*
2			* cmglare.c - routines for calculating glare autonomy.
3			*
4			* N. Jones
5			*/
6
7			/*
8			* Copyright (c) 2017-2019 Nathaniel Jones
9			*
10			* Permission is hereby granted, free of charge, to any person obtaining a copy
11			* of this software and associated documentation files (the "Software"), to deal
12			* in the Software without restriction, including without limitation the rights
13			* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
14			* copies of the Software, and to permit persons to whom the Software is
15			* furnished to do so, subject to the following conditions:
16			*
17			* The above copyright notice and this permission notice shall be included in all
18			* copies or substantial portions of the Software.
19			*
20			* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
21			* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22			* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
23			* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24			* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
25			* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
26			* SOFTWARE.
27			*/
28
29			#include "rtmath.h"
30			#include "cmglare.h"
31
32			#define LUMINOUS_EFFICACY 179 /* lumens per Watt */
33			#define LUMINANCE_THRESHOLD 100 /* minimum threshold that will be interpreted as luminance rather than ratio between Ev and glare source */
34
35			#define ANGLE(u,v) acos(DOT((u),(v))) // TODO need to normalize?
36
37			typedef struct reinhart_sky {
38			//int mf; /* Linear divisions per Tregenza patch. */
39			int rings; /* Number of rings of sky patches. */
40			//int patches; /* Number of sky patches. */
41			double ringElevationAngle; /* Angle in radians between rings. */
42			int patchesPerRow; / Sky patches per row. */
43			int firstPatchIndex; / Index of first sky patch in each row. */
44			double solidAngle; / Solid angle of each patch in each row. */
45			} ReinhartSky;
46
47			static const int tnaz[] = { 30, 30, 24, 24, 18, 12, 6 }; /* Number of patches per row */
48
49			extern char* progname;
50
51			static ReinhartSky* make_sky(const CMATRIX *smx)
52			{
53			int mf, count, i, j;
54			double remaining = PI / 2;
55
56			/* Check sky partitionss */
57			switch (smx->nrows) { // nrows is Reinhart sky subdivisoins plus one for miss
58			case 2: mf = 0; break;
59			case 146: mf = 1; break;
60			case 578: mf = 2; break;
61			case 1298: mf = 3; break;
62			case 2306: mf = 4; break;
63			case 3602: mf = 5; break;
64			case 5186: mf = 6; break;
65			case 7058: mf = 7; break;
66			case 9218: mf = 8; break;
67			case 11666: mf = 9; break;
68			case 14402: mf = 10; break;
69			case 17426: mf = 11; break;
70			case 20738: mf = 12; break;
71			default:
72			fprintf(stderr,
73			"%s: unknown number of sky patches %d\n",
74			progname, smx->nrows);
75			return NULL;
76			}
77
78			/* Allocate sky */
79			ReinhartSky sky = (ReinhartSky)malloc(sizeof(ReinhartSky));
80			if (!sky) goto memerr;
81
82			/* Calculate patches per row */
83			sky->rings = 7 * mf + 1;
84			sky->patchesPerRow = (int)malloc(sky->rings sizeof(int));
85			sky->firstPatchIndex = (int)malloc(sky->rings sizeof(int));
86			if (!sky->patchesPerRow \|\| !sky->firstPatchIndex) goto memerr;
87			count = 1; // The below horizon patch
88			for (i = 0; i < 7; i++)
89			for (j = 0; j < mf; j++) {
90			sky->firstPatchIndex[i * mf + j] = count;
91			count += sky->patchesPerRow[i * mf + j] = tnaz[i] * mf;
92			}
93			sky->firstPatchIndex[7 * mf] = count;
94			sky->patchesPerRow[7 * mf] = 1;
95			//sky->patches = count + 1;
96
97			/* Calculate solid angle of patches in each row */
98			sky->solidAngle = (double)malloc(sky->rings sizeof(double));
99			if (!sky->solidAngle) goto memerr;
100			sky->ringElevationAngle = PI / (2 * sky->rings - 1);
101			sky->solidAngle[0] = 2 * PI;
102			for (i = 1; i < sky->rings; i++) {
103			remaining -= sky->ringElevationAngle;
104			sky->solidAngle[i] = 2 * PI * (1 - cos(remaining)); // solid angle of cap
105			sky->solidAngle[i - 1] -= sky->solidAngle[i];
106			sky->solidAngle[i - 1] /= sky->patchesPerRow[i - 1];
107			}
108
109			return sky;
110
111			memerr:
112			fprintf(stderr,
113			"%s: out of memory\n",
114			progname);
115			return NULL;
116			}
117
118			static void free_sky(ReinhartSky *sky)
119			{
120			if (sky) {
121			free(sky->patchesPerRow);
122			free(sky->firstPatchIndex);
123			free(sky->solidAngle);
124			free(sky);
125			}
126			}
127
128			static void get_patch_direction(const ReinhartSky *sky, const int patch, FVECT target)
129			{
130			//if (patch >= sky->patches \|\| patch < 0) throw new RuntimeException("Illegal patch " + patch);
131			if (!patch) {
132			target[0] = target[1] = 0;
133			target[2] = -1;
134			return; // Ignore below horizon?
135			}
136			int row = sky->rings - 1;
137			while (patch < sky->firstPatchIndex[row]) row--;
138			const double alt = PI / 2 - sky->ringElevationAngle * (row - (sky->rings - 1));
139			const double azi = 2 * PI * (patch - sky->firstPatchIndex[row]) / sky->patchesPerRow[row];
140			const double cos_alt = cos(alt);
141			target[0] = cos_alt * -sin(azi);
142			target[1] = cos_alt * -cos(azi);
143			target[2] = sin(alt);
144			}
145
146			static double get_patch_solid_angle(const ReinhartSky *sky, const int patch, const double cos_theta)
147			{
148			//if (patch >= sky->patches \|\| patch < 0) throw new RuntimeException("Illegal patch " + patch);
149			if (!patch) return 2 * (PI - 2 * acos(cos_theta)); // Solid angle overlap between visible hemisphere and ground hemisphere
150			int row = sky->rings - 1;
151			while (patch < sky->firstPatchIndex[row]) row--;
152			return sky->solidAngle[row];
153			}
154
155			static double get_guth(const FVECT dir, const FVECT forward, const FVECT up)
156			{
157			double posindex;
158			FVECT hv, temp;
159			VCROSS(hv, forward, up);
160			normalize(hv);
161			VCROSS(temp, forward, hv);
162			double phi = ANGLE(dir, temp) - PI / 2;
163
164			/* Guth model, equation from IES lighting handbook */
165			if (phi >= 0) {
166			double sigma = ANGLE(dir, forward);
167			VSUM(hv, forward, dir, 1 / DOT(dir, forward));
168			normalize(hv);
169			double tau = ANGLE(hv, up);
170			tau *= 180.0 / PI;
171			sigma *= 180.0 / PI;
172
173			if (sigma <= 0)
174			sigma = -sigma;
175
176			posindex = exp((35.2 - 0.31889 * tau - 1.22 * exp(-2.0 * tau / 9.0)) / 1000.0 * sigma + (21.0 + 0.26667 * tau - 0.002963 * tau * tau) / 100000.0 * sigma * sigma);
177			}
178			/* below line of sight, using Iwata model */
179			else {
180			double teta = PI / 2 - ANGLE(dir, hv);
181
182			if (teta == 0.0)
183			teta = FTINY;
184
185			double fact = 0.8;
186			double d = 1 / tan(phi);
187			double s = tan(teta) / tan(phi);
188			double r = sqrt((1 + s * s) / (d * d));
189			if (r > 0.6)
190			fact = 1.2;
191			if (r > 3)
192			r = 3.0;
193
194			posindex = 1.0 + fact * r;
195			}
196			if (posindex > 16)
197			posindex = 16.0;
198
199			return posindex;
200			}
201
202			float* cm_glare(const CMATRIX dcmx, const CMATRIX evmx, const CMATRIX smx, const int occupied, const double dgp_limit, const double dgp_threshold, const FVECT *views, const FVECT dir, const FVECT up)
203			{
204			int p, t, c;
205			int hourly_output = dgp_limit < 0;
206			float *dgp_list;
207			ReinhartSky *sky;
208			FVECT vdir;
209
210			/* Check consistancy */
211			if ((dcmx->nrows != evmx->nrows) \| (dcmx->ncols != smx->nrows) \| (evmx->ncols != smx->ncols)) {
212			fprintf(stderr,
213			"%s: inconsistant matrix dimensions: dc(%d, %d) ev(%d, %d) s(%d, %d)\n",
214			progname, dcmx->nrows, dcmx->ncols, evmx->nrows, evmx->ncols, smx->nrows, smx->ncols);
215			return NULL;
216			}
217
218			/* Create output buffer */
219			dgp_list = (float)malloc(evmx->nrows (hourly_output ? evmx->ncols : 1) * sizeof(float));
220			if (!dgp_list) {
221			fprintf(stderr,
222			"%s: out of memory in cm_glare()\n",
223			progname);
224			return NULL;
225			}
226
227			/* Create sky */
228			sky = make_sky(smx);
229			if (!sky) return NULL;
230
231			/* Calculate glare limit */
232			double ev_max = -1;
233			if (!hourly_output) {
234			ev_max = (dgp_limit - 0.16) / 5.87e-5;
235			if (ev_max < 0) ev_max = 0;
236			}
237
238			/* For each position and direction */
239			if (!views) VCOPY(vdir, dir);
240			for (p = 0; p < evmx->nrows; p++) {
241			/* For each time step */
242			int occupied_hours = 0;
243			int glare_hours = 0;
244			if (views) VCOPY(vdir, views[p]);
245			for (t = 0; t < evmx->ncols; t++) {
246			if (!occupied[t]) {
247			/* Not occupied */
248			if (hourly_output) dgp_list[p * evmx->ncols + t] = 0.0f;
249			}
250			else {
251			/* Occupied */
252			double illum = LUMINOUS_EFFICACY * bright(cm_lval(evmx, p, t));
253			occupied_hours++;
254
255			if (illum <= FTINY) {
256			/* No light, so no glare */
257			if (hourly_output) dgp_list[p * evmx->ncols + t] = 0.0f;
258			}
259			else if ((illum >= ev_max) & (!hourly_output)) {
260			/* Guarangeed glare */
261			glare_hours++;
262			}
263			else {
264			/* Calculate enhanced simplified daylight glare probability */
265			double sum = 0.0;
266			FVECT patch_normal;
267			for (c = 0; c < smx->nrows; c++) {
268			const double dc = bright(cm_lval(dcmx, p, c));
269			if (dc > 0) {
270			get_patch_direction(sky, c, patch_normal);
271			if (!c) {
272			/* Direction toward the center of the visible ground */
273			VADD(patch_normal, patch_normal, vdir);
274			if (normalize(patch_normal) == 0) continue;
275			}
276			const double cos_theta = DOT(vdir, patch_normal);
277			if (cos_theta <= FTINY) continue;
278			const double s = LUMINOUS_EFFICACY * bright(cm_lval(smx, c, t));
279			const double omega = get_patch_solid_angle(sky, c, cos_theta);
280			const double patch_luminance = (dc * s) / (omega * cos_theta);
281
282			double min_patch_luminance = dgp_threshold;
283			if (dgp_threshold < LUMINANCE_THRESHOLD)
284			min_patch_luminance *= illum / PI; // TODO should use average luminance, not illuminance
285			if (patch_luminance < min_patch_luminance) continue;
286			const double P = get_guth(patch_normal, vdir, up);
287			sum += (patch_luminance * patch_luminance * omega) / (P * P);
288			}
289			}
290
291			//double dgp = 5.87e-5 * illum + 0.092 * log10(1 + dgp / pow(illum, 1.87)) + 0.159;
292			double eDGPs = 5.87e-5 * illum + 0.0918 * log10(1 + sum / pow(illum, 1.87)) + 0.16;
293			if (illum < 1000) /* low light correction */
294			eDGPs = exp(0.024 illum - 4) / (1 + exp(0.024 * illum - 4));
295			//eDGPs /= 1.1 - 0.5 * age / 100.0; /* age correction */
296			if (eDGPs > 1.0) eDGPs = 1.0;
297
298			if (hourly_output)
299			dgp_list[p * evmx->ncols + t] = (float)eDGPs;
300			else if (eDGPs >= dgp_limit)
301			glare_hours++;
302			}
303			}
304			}
305			if (!hourly_output) {
306			/* Save glare autonomy */
307			dgp_list[p] = (float)(occupied_hours - glare_hours) / occupied_hours;
308			}
309			}
310
311			free_sky(sky);
312
313			return dgp_list;
314			}
315
316			static int getvec(FVECT vec, const int dtype, FILE fp) / get a vector from fp */
317			{
318			static float vf[3];
319			static double vd[3];
320			char buf[32];
321			int i;
322
323			switch (dtype) {
324			case DTascii:
325			for (i = 0; i < 3; i++) {
326			if (fgetword(buf, sizeof(buf), fp) == NULL \|\|
327			!isflt(buf))
328			return(-1);
329			vec[i] = atof(buf);
330			}
331			break;
332			case DTfloat:
333			if (getbinary(vf, sizeof(float), 3, fp) != 3)
334			return(-1);
335			VCOPY(vec, vf);
336			break;
337			case DTdouble:
338			if (getbinary(vd, sizeof(double), 3, fp) != 3)
339			return(-1);
340			VCOPY(vec, vd);
341			break;
342			default:
343			fprintf(stderr,
344			"%s: botched input format\n",
345			progname);
346			return(-1);
347			}
348			return(0);
349			}
350
351			int cm_load_schedule(const int count, int* schedule, FILE *fp)
352			{
353			char buf[512];
354			char *cp;
355			char *comma;
356			double val;
357			int i = 0;
358
359			while (fgetline(buf, sizeof(buf), fp) != NULL) {
360			if (buf[0] == '#') continue; // Comment line
361			comma = NULL;
362			for (cp = buf; *cp; cp++) {
363			/* If there are multiple commas, assume the value is after the last comma */
364			if (*cp == ',') {
365			comma = cp; /* Record position of last comma */
366			}
367			}
368			if (comma)
369			val = atof(comma + 1);
370			else
371			val = atof(buf);
372
373			if (i < count) {
374			/* Add the value to the schedule */
375			schedule[i++] = (val > 0);
376			}
377			else {
378			fprintf(stderr,
379			"%s: too many schedule entries\n",
380			progname);
381			return(1);
382			}
383			}
384			fclose(fp);
385
386			if (i < count) {
387			fprintf(stderr,
388			"%s: too few schedule entries\n",
389			progname);
390			return(-1);
391			}
392
393			return 0;
394			}
395
396			FVECT* cm_load_views(const int nrows, const int dtype, FILE *fp)
397			{
398			int i;
399			double d;
400			FVECT orig;
401
402			FVECT views = (FVECT)malloc(nrows * sizeof(FVECT));
403			if (!views) {
404			fprintf(stderr,
405			"%s: out of memory in cm_load_views()\n",
406			progname);
407			return NULL;
408			}
409
410			for (i = 0; i < nrows; i++) {
411			if (getvec(orig, dtype, fp) \| getvec(views[i], dtype, fp)) {
412			fprintf(stderr,
413			"%s: unexpected end of input, missing %d entries\n",
414			progname, i);
415			return NULL;
416			}
417
418			d = normalize(views[i]);
419			if (d == 0.0) { /* zero ==> flush */
420			fprintf(stderr,
421			"%s: zero length direction detected\n",
422			progname);
423			return NULL;
424			}
425			}
426
427			return views;
428			}
429
430			int cm_write_glare(const float mp, const int nrows, const int ncols, const int dtype, FILE fp)
431			{
432			static const char tabEOL[2] = { '\t', '\n' };
433			int r, c;
434			double dc[1];
435
436			switch (dtype) {
437			case DTascii:
438			for (r = 0; r < nrows; r++)
439			for (c = 0; c < ncols; c++, mp++)
440			fprintf(fp, "%.6e%c",
441			mp[0],
442			tabEOL[c >= ncols - 1]);
443			break;
444			case DTfloat:
445			r = ncols*nrows;
446			while (r > 0) {
447			c = putbinary(mp, sizeof(float), r, fp);
448			if (c <= 0)
449			return(0);
450			mp += c;
451			r -= c;
452			}
453			break;
454			case DTdouble:
455			r = ncols*nrows;
456			while (r--) {
457			dc[0] = mp[0];
458			if (putbinary(dc, sizeof(double), 1, fp) != 1)
459			return(0);
460			mp++;
461			}
462			break;
463			default:
464			fputs("Unsupported data type in cm_write_glare()!\n", stderr);
465			return(0);
466			}
467			return(fflush(fp) == 0);
468			}