src/common/ccolor.c

#ifndef lint
static const char RCSid[] = "$Id: ccolor.c,v 3.1 2011/02/18 00:40:25 greg Exp $";
#endif
/*
 * Spectral color handling routines
 */

#include <stdio.h>
#include <math.h>
#include "ccolor.h"


C_COLOR         c_dfcolor = C_DEFCOLOR;

                                /* CIE 1931 Standard Observer curves */
static const C_COLOR    cie_xf = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY|C_CSEFF,
                        {14,42,143,435,1344,2839,3483,3362,2908,1954,956,
                        320,49,93,633,1655,2904,4334,5945,7621,9163,10263,
                        10622,10026,8544,6424,4479,2835,1649,874,468,227,
                        114,58,29,14,7,3,2,1,0}, 106836L, .467, .368, 362.230
                        };
static const C_COLOR    cie_yf = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY|C_CSEFF,
                        {0,1,4,12,40,116,230,380,600,910,1390,2080,3230,
                        5030,7100,8620,9540,9950,9950,9520,8700,7570,6310,
                        5030,3810,2650,1750,1070,610,320,170,82,41,21,10,
                        5,2,1,1,0,0}, 106856L, .398, .542, 493.525
                        };
static const C_COLOR    cie_zf = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY|C_CSEFF,
                        {65,201,679,2074,6456,13856,17471,17721,16692,
                        12876,8130,4652,2720,1582,782,422,203,87,39,21,17,
                        11,8,3,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
                        106770L, .147, .077, 54.363
                        };
                                /* Derived CIE 1931 Primaries (imaginary) */
static const C_COLOR    cie_xp = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY,
                        {-174,-198,-195,-197,-202,-213,-235,-272,-333,
                        -444,-688,-1232,-2393,-4497,-6876,-6758,-5256,
                        -3100,-815,1320,3200,4782,5998,6861,7408,7754,
                        7980,8120,8199,8240,8271,8292,8309,8283,8469,
                        8336,8336,8336,8336,8336,8336},
                        127424L, 1., .0,
                        };
static const C_COLOR    cie_yp = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY,
                        {-451,-431,-431,-430,-427,-417,-399,-366,-312,
                        -204,57,691,2142,4990,8810,9871,9122,7321,5145,
                        3023,1123,-473,-1704,-2572,-3127,-3474,-3704,
                        -3846,-3927,-3968,-3999,-4021,-4038,-4012,-4201,
                        -4066,-4066,-4066,-4066,-4066,-4066},
                        -23035L, .0, 1.,
                        };
static const C_COLOR    cie_zp = { 1, NULL, C_CDSPEC|C_CSSPEC|C_CSXY,
                        {4051,4054,4052,4053,4054,4056,4059,4064,4071,
                        4074,4056,3967,3677,2933,1492,313,-440,-795,
                        -904,-918,-898,-884,-869,-863,-855,-855,-851,
                        -848,-847,-846,-846,-846,-845,-846,-843,-845,
                        -845,-845,-845,-845,-845},
                        36057L, .0, .0,
                        };

                                        /* Sharp primary matrix */
static const float      toSharp[3][3] = {
        { 1.2694, -0.0988, -0.1706},
        {-0.8364,  1.8006,  0.0357},
        { 0.0297, -0.0315,  1.0018}
};
                                        /* inverse Sharp primary matrix */
static const float      fromSharp[3][3] = {
        { 0.8156,  0.0472,  0.1372},
        { 0.3791,  0.5769,  0.0440},
        {-0.0123,  0.0167,  0.9955}
};

static void
c_toSharpRGB(C_COLOR *cin, double cieY, float cout[3])
{
        double  xyz[3];

        c_ccvt(cin, C_CSXY);

        xyz[0] = cin->cx/cin->cy * cieY;
        xyz[1] = cieY;
        xyz[2] = (1. - cin->cx - cin->cy)/cin->cy * cieY;

        cout[0] = toSharp[0][0]*xyz[0] + toSharp[0][1]*xyz[1] +
                                toSharp[0][2]*xyz[2];
        cout[1] = toSharp[1][0]*xyz[0] + toSharp[1][1]*xyz[1] +
                                toSharp[1][2]*xyz[2];
        cout[2] = toSharp[2][0]*xyz[0] + toSharp[2][1]*xyz[1] +
                                toSharp[2][2]*xyz[2];
}

static double
c_fromSharpRGB(float cin[3], C_COLOR *cout)
{
        double  xyz[3], sf;

        xyz[0] = fromSharp[0][0]*cin[0] + fromSharp[0][1]*cin[1] +
                                fromSharp[0][2]*cin[2];
        xyz[1] = fromSharp[1][0]*cin[0] + fromSharp[1][1]*cin[1] +
                                fromSharp[1][2]*cin[2];
        xyz[2] = fromSharp[2][0]*cin[0] + fromSharp[2][1]*cin[1] +
                                fromSharp[2][2]*cin[2];
                                
        sf = 1./(xyz[0] + xyz[1] + xyz[2]);

        cout->cx = xyz[0] * sf;
        cout->cy = xyz[1] * sf;
        cout->flags = C_CDXY|C_CSXY;
        
        return(xyz[1]);
}

/* check if color is grey */
int
c_isgrey(C_COLOR *clr)
{
        if (!(clr->flags & (C_CSXY|C_CSSPEC)))
                return(1);              /* no settings == grey */

        c_ccvt(clr, C_CSXY);

        return(clr->cx >= .323 && clr->cx <= .343 &&
                        clr->cy >= .323 && clr->cy <= .343);
}

/* convert color representations */
void
c_ccvt(C_COLOR *clr, int fl)
{
        double  x, y, z;
        int     i;

        fl &= ~clr->flags;                      /* ignore what's done */
        if (!fl)                                /* everything's done! */
                return;
        if (!(clr->flags & (C_CSXY|C_CSSPEC))) {
                *clr = c_dfcolor;               /* nothing was set! */
                return;
        }
        if (fl & C_CSXY) {              /* cspec -> cxy */
                x = y = z = 0.;
                for (i = 0; i < C_CNSS; i++) {
                        x += cie_xf.ssamp[i] * clr->ssamp[i];
                        y += cie_yf.ssamp[i] * clr->ssamp[i];
                        z += cie_zf.ssamp[i] * clr->ssamp[i];
                }
                x /= (double)cie_xf.ssum;
                y /= (double)cie_yf.ssum;
                z /= (double)cie_zf.ssum;
                z += x + y;
                clr->cx = x / z;
                clr->cy = y / z;
                clr->flags |= C_CSXY;
        } else if (fl & C_CSSPEC) {     /* cxy -> cspec */
                x = clr->cx;
                y = clr->cy;
                z = 1. - x - y;
                clr->ssum = 0;
                for (i = 0; i < C_CNSS; i++) {
                        clr->ssamp[i] = x*cie_xp.ssamp[i] + y*cie_yp.ssamp[i]
                                        + z*cie_zp.ssamp[i] + .5;
                        if (clr->ssamp[i] < 0)          /* out of gamut! */
                                clr->ssamp[i] = 0;
                        else
                                clr->ssum += clr->ssamp[i];
                }
                clr->flags |= C_CSSPEC;
        }
        if (fl & C_CSEFF) {             /* compute efficacy */
                if (clr->flags & C_CSSPEC) {            /* from spectrum */
                        y = 0.;
                        for (i = 0; i < C_CNSS; i++)
                                y += cie_yf.ssamp[i] * clr->ssamp[i];
                        clr->eff = C_CLPWM * y / clr->ssum;
                } else /* clr->flags & C_CSXY */ {      /* from (x,y) */
                        clr->eff = clr->cx*cie_xf.eff + clr->cy*cie_yf.eff +
                                        (1. - clr->cx - clr->cy)*cie_zf.eff;
                }
                clr->flags |= C_CSEFF;
        }
}

/* mix two colors according to weights given -- cres may be c1 or c2 */
void
c_cmix(C_COLOR *cres, double w1, C_COLOR *c1, double w2, C_COLOR *c2)
{
        double  scale;
        int     i;

        if ((c1->flags|c2->flags) & C_CDSPEC) {         /* spectral mixing */
                float   cmix[C_CNSS];

                c_ccvt(c1, C_CSSPEC|C_CSEFF);
                c_ccvt(c2, C_CSSPEC|C_CSEFF);
                w1 /= c1->eff*c1->ssum;
                w2 /= c2->eff*c2->ssum;
                scale = 0.;
                for (i = 0; i < C_CNSS; i++) {
                        cmix[i] = w1*c1->ssamp[i] + w2*c2->ssamp[i];
                        if (cmix[i] > scale)
                                scale = cmix[i];
                }
                scale = C_CMAXV / scale;
                cres->ssum = 0;
                for (i = 0; i < C_CNSS; i++)
                        cres->ssum += cres->ssamp[i] = scale*cmix[i] + .5;
                cres->flags = C_CDSPEC|C_CSSPEC;
        } else {                                        /* CIE xy mixing */
                c_ccvt(c1, C_CSXY);
                c_ccvt(c2, C_CSXY);
                w1 *= c2->cy;
                w2 *= c1->cy;
                scale = w1 + w2;
                if (scale == 0.) {
                        *cres = c_dfcolor;
                        return;
                }
                scale = 1. / scale;
                cres->cx = (w1*c1->cx + w2*c2->cx) * scale;
                cres->cy = (w1*c1->cy + w2*c2->cy) * scale;
                cres->flags = C_CDXY|C_CSXY;
        }
}

/* multiply two colors -- cres may be c1 or c2 */
double
c_cmult(C_COLOR *cres, C_COLOR *c1, double y1, C_COLOR *c2, double y2)
{
        double  yres;
        int     i;

        if ((c1->flags|c2->flags) & C_CDSPEC) {
                long    cmix[C_CNSS], cmax = 0;         /* spectral multiply */

                c_ccvt(c1, C_CSSPEC|C_CSEFF);
                c_ccvt(c2, C_CSSPEC|C_CSEFF);
                for (i = 0; i < C_CNSS; i++) {
                        cmix[i] = c1->ssamp[i] * c2->ssamp[i];
                        if (cmix[i] > cmax)
                                cmax = cmix[i];
                }
                cmax /= C_CMAXV;
                if (!cmax) {
                        *cres = c_dfcolor;
                        return(0.);
                }
                cres->ssum = 0;
                for (i = 0; i < C_CNSS; i++)
                        cres->ssum += cres->ssamp[i] = cmix[i] / cmax;
                cres->flags = C_CDSPEC|C_CSSPEC;

                c_ccvt(cres, C_CSEFF);                  /* nasty, but true */
                yres = (y1 * y2 * cie_yf.ssum * C_CLPWM) /
                        (c1->eff*c1->ssum * c2->eff*c2->ssum) *
                        cres->eff*( cres->ssum*(double)cmax +
                                                C_CNSS/2.0*(cmax-1) );
        } else {
                float   rgb1[3], rgb2[3];               /* CIE xy multiply */

                c_toSharpRGB(c1, y1, rgb1);
                c_toSharpRGB(c2, y2, rgb2);
                
                rgb2[0] *= rgb1[0]; rgb2[1] *= rgb1[1]; rgb2[2] *= rgb1[2];

                yres = c_fromSharpRGB(rgb2, cres);
        }
        return(yres);
}

#define C1              3.741832e-16    /* W-m^2 */
#define C2              1.4388e-2       /* m-K */

#define bbsp(l,t)       (C1/((l)*(l)*(l)*(l)*(l)*(exp(C2/((t)*(l)))-1.)))
#define bblm(t)         (C2*0.2/(t))

/* set black body spectrum */
int
c_bbtemp(C_COLOR *clr, double tk)
{
        double  sf, wl;
        int     i;

        if (tk < 1000)
                return(0);
        wl = bblm(tk);                  /* scalefactor based on peak */
        if (wl < C_CMINWL*1e-9)
                wl = C_CMINWL*1e-9;
        else if (wl > C_CMAXWL*1e-9)
                wl = C_CMAXWL*1e-9;
        sf = C_CMAXV/bbsp(wl,tk);
        clr->ssum = 0;
        for (i = 0; i < C_CNSS; i++) {
                wl = (C_CMINWL + i*C_CWLI)*1e-9;
                clr->ssum += clr->ssamp[i] = sf*bbsp(wl,tk) + .5;
        }
        clr->flags = C_CDSPEC|C_CSSPEC;
        return(1);
}

#undef  C1
#undef  C2
#undef  bbsp
#undef  bblm
Revision:	3.2
Committed:	Thu May 17 05:47:59 2012 UTC (11 years, 11 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 3.1:	+109 -10 lines
Log Message:	Added spectral color multiplication
#	User	Rev	Content
1	greg	3.1	#ifndef lint
2	greg	3.2	static const char RCSid[] = "$Id: ccolor.c,v 3.1 2011/02/18 00:40:25 greg Exp $";
3	greg	3.1	#endif
4			/*
5			* Spectral color handling routines
6			*/
7
8			#include <stdio.h>
9			#include <math.h>
10			#include "ccolor.h"
11
12
13			C_COLOR c_dfcolor = C_DEFCOLOR;
14
15			/* CIE 1931 Standard Observer curves */
16	greg	3.2	static const C_COLOR cie_xf = { 1, NULL, C_CDSPEC\|C_CSSPEC\|C_CSXY\|C_CSEFF,
17	greg	3.1	{14,42,143,435,1344,2839,3483,3362,2908,1954,956,
18			320,49,93,633,1655,2904,4334,5945,7621,9163,10263,
19			10622,10026,8544,6424,4479,2835,1649,874,468,227,
20			114,58,29,14,7,3,2,1,0}, 106836L, .467, .368, 362.230
21			};
22	greg	3.2	static const C_COLOR cie_yf = { 1, NULL, C_CDSPEC\|C_CSSPEC\|C_CSXY\|C_CSEFF,
23	greg	3.1	{0,1,4,12,40,116,230,380,600,910,1390,2080,3230,
24			5030,7100,8620,9540,9950,9950,9520,8700,7570,6310,
25			5030,3810,2650,1750,1070,610,320,170,82,41,21,10,
26			5,2,1,1,0,0}, 106856L, .398, .542, 493.525
27			};
28	greg	3.2	static const C_COLOR cie_zf = { 1, NULL, C_CDSPEC\|C_CSSPEC\|C_CSXY\|C_CSEFF,
29	greg	3.1	{65,201,679,2074,6456,13856,17471,17721,16692,
30			12876,8130,4652,2720,1582,782,422,203,87,39,21,17,
31			11,8,3,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
32			106770L, .147, .077, 54.363
33			};
34			/* Derived CIE 1931 Primaries (imaginary) */
35	greg	3.2	static const C_COLOR cie_xp = { 1, NULL, C_CDSPEC\|C_CSSPEC\|C_CSXY,
36	greg	3.1	{-174,-198,-195,-197,-202,-213,-235,-272,-333,
37			-444,-688,-1232,-2393,-4497,-6876,-6758,-5256,
38			-3100,-815,1320,3200,4782,5998,6861,7408,7754,
39			7980,8120,8199,8240,8271,8292,8309,8283,8469,
40			8336,8336,8336,8336,8336,8336},
41			127424L, 1., .0,
42			};
43	greg	3.2	static const C_COLOR cie_yp = { 1, NULL, C_CDSPEC\|C_CSSPEC\|C_CSXY,
44	greg	3.1	{-451,-431,-431,-430,-427,-417,-399,-366,-312,
45			-204,57,691,2142,4990,8810,9871,9122,7321,5145,
46			3023,1123,-473,-1704,-2572,-3127,-3474,-3704,
47			-3846,-3927,-3968,-3999,-4021,-4038,-4012,-4201,
48			-4066,-4066,-4066,-4066,-4066,-4066},
49			-23035L, .0, 1.,
50			};
51	greg	3.2	static const C_COLOR cie_zp = { 1, NULL, C_CDSPEC\|C_CSSPEC\|C_CSXY,
52	greg	3.1	{4051,4054,4052,4053,4054,4056,4059,4064,4071,
53			4074,4056,3967,3677,2933,1492,313,-440,-795,
54			-904,-918,-898,-884,-869,-863,-855,-855,-851,
55			-848,-847,-846,-846,-846,-845,-846,-843,-845,
56			-845,-845,-845,-845,-845},
57			36057L, .0, .0,
58			};
59
60	greg	3.2	/* Sharp primary matrix */
61			static const float toSharp[3][3] = {
62			{ 1.2694, -0.0988, -0.1706},
63			{-0.8364, 1.8006, 0.0357},
64			{ 0.0297, -0.0315, 1.0018}
65			};
66			/* inverse Sharp primary matrix */
67			static const float fromSharp[3][3] = {
68			{ 0.8156, 0.0472, 0.1372},
69			{ 0.3791, 0.5769, 0.0440},
70			{-0.0123, 0.0167, 0.9955}
71			};
72
73			static void
74			c_toSharpRGB(C_COLOR *cin, double cieY, float cout[3])
75			{
76			double xyz[3];
77
78			c_ccvt(cin, C_CSXY);
79
80			xyz[0] = cin->cx/cin->cy * cieY;
81			xyz[1] = cieY;
82			xyz[2] = (1. - cin->cx - cin->cy)/cin->cy * cieY;
83
84			cout[0] = toSharp[0][0]xyz[0] + toSharp[0][1]xyz[1] +
85			toSharp[0][2]*xyz[2];
86			cout[1] = toSharp[1][0]xyz[0] + toSharp[1][1]xyz[1] +
87			toSharp[1][2]*xyz[2];
88			cout[2] = toSharp[2][0]xyz[0] + toSharp[2][1]xyz[1] +
89			toSharp[2][2]*xyz[2];
90			}
91
92			static double
93			c_fromSharpRGB(float cin[3], C_COLOR *cout)
94			{
95			double xyz[3], sf;
96
97			xyz[0] = fromSharp[0][0]cin[0] + fromSharp[0][1]cin[1] +
98			fromSharp[0][2]*cin[2];
99			xyz[1] = fromSharp[1][0]cin[0] + fromSharp[1][1]cin[1] +
100			fromSharp[1][2]*cin[2];
101			xyz[2] = fromSharp[2][0]cin[0] + fromSharp[2][1]cin[1] +
102			fromSharp[2][2]*cin[2];
103
104			sf = 1./(xyz[0] + xyz[1] + xyz[2]);
105
106			cout->cx = xyz[0] * sf;
107			cout->cy = xyz[1] * sf;
108			cout->flags = C_CDXY\|C_CSXY;
109
110			return(xyz[1]);
111			}
112
113	greg	3.1	/* check if color is grey */
114			int
115			c_isgrey(C_COLOR *clr)
116			{
117			if (!(clr->flags & (C_CSXY\|C_CSSPEC)))
118			return(1); /* no settings == grey */
119
120			c_ccvt(clr, C_CSXY);
121
122			return(clr->cx >= .323 && clr->cx <= .343 &&
123			clr->cy >= .323 && clr->cy <= .343);
124			}
125
126			/* convert color representations */
127			void
128			c_ccvt(C_COLOR *clr, int fl)
129			{
130			double x, y, z;
131			int i;
132
133			fl &= ~clr->flags; /* ignore what's done */
134			if (!fl) /* everything's done! */
135			return;
136			if (!(clr->flags & (C_CSXY\|C_CSSPEC))) {
137			clr = c_dfcolor; / nothing was set! */
138			return;
139			}
140			if (fl & C_CSXY) { /* cspec -> cxy */
141			x = y = z = 0.;
142			for (i = 0; i < C_CNSS; i++) {
143			x += cie_xf.ssamp[i] * clr->ssamp[i];
144			y += cie_yf.ssamp[i] * clr->ssamp[i];
145			z += cie_zf.ssamp[i] * clr->ssamp[i];
146			}
147			x /= (double)cie_xf.ssum;
148			y /= (double)cie_yf.ssum;
149			z /= (double)cie_zf.ssum;
150			z += x + y;
151			clr->cx = x / z;
152			clr->cy = y / z;
153			clr->flags \|= C_CSXY;
154			} else if (fl & C_CSSPEC) { /* cxy -> cspec */
155			x = clr->cx;
156			y = clr->cy;
157			z = 1. - x - y;
158			clr->ssum = 0;
159			for (i = 0; i < C_CNSS; i++) {
160			clr->ssamp[i] = xcie_xp.ssamp[i] + ycie_yp.ssamp[i]
161			+ z*cie_zp.ssamp[i] + .5;
162			if (clr->ssamp[i] < 0) /* out of gamut! */
163			clr->ssamp[i] = 0;
164			else
165			clr->ssum += clr->ssamp[i];
166			}
167			clr->flags \|= C_CSSPEC;
168			}
169			if (fl & C_CSEFF) { /* compute efficacy */
170			if (clr->flags & C_CSSPEC) { /* from spectrum */
171			y = 0.;
172			for (i = 0; i < C_CNSS; i++)
173			y += cie_yf.ssamp[i] * clr->ssamp[i];
174			clr->eff = C_CLPWM * y / clr->ssum;
175			} else /* clr->flags & C_CSXY / { / from (x,y) */
176			clr->eff = clr->cxcie_xf.eff + clr->cycie_yf.eff +
177			(1. - clr->cx - clr->cy)*cie_zf.eff;
178			}
179			clr->flags \|= C_CSEFF;
180			}
181			}
182
183	greg	3.2	/* mix two colors according to weights given -- cres may be c1 or c2 */
184	greg	3.1	void
185			c_cmix(C_COLOR cres, double w1, C_COLOR c1, double w2, C_COLOR *c2)
186			{
187			double scale;
188			int i;
189
190			if ((c1->flags\|c2->flags) & C_CDSPEC) { /* spectral mixing */
191	greg	3.2	float cmix[C_CNSS];
192
193	greg	3.1	c_ccvt(c1, C_CSSPEC\|C_CSEFF);
194			c_ccvt(c2, C_CSSPEC\|C_CSEFF);
195			w1 /= c1->eff*c1->ssum;
196			w2 /= c2->eff*c2->ssum;
197			scale = 0.;
198			for (i = 0; i < C_CNSS; i++) {
199			cmix[i] = w1c1->ssamp[i] + w2c2->ssamp[i];
200			if (cmix[i] > scale)
201			scale = cmix[i];
202			}
203			scale = C_CMAXV / scale;
204			cres->ssum = 0;
205			for (i = 0; i < C_CNSS; i++)
206			cres->ssum += cres->ssamp[i] = scale*cmix[i] + .5;
207			cres->flags = C_CDSPEC\|C_CSSPEC;
208			} else { /* CIE xy mixing */
209			c_ccvt(c1, C_CSXY);
210			c_ccvt(c2, C_CSXY);
211			w1 *= c2->cy;
212			w2 *= c1->cy;
213			scale = w1 + w2;
214			if (scale == 0.) {
215			*cres = c_dfcolor;
216			return;
217			}
218			scale = 1. / scale;
219			cres->cx = (w1c1->cx + w2c2->cx) * scale;
220			cres->cy = (w1c1->cy + w2c2->cy) * scale;
221			cres->flags = C_CDXY\|C_CSXY;
222			}
223			}
224
225	greg	3.2	/* multiply two colors -- cres may be c1 or c2 */
226			double
227			c_cmult(C_COLOR cres, C_COLOR c1, double y1, C_COLOR *c2, double y2)
228			{
229			double yres;
230			int i;
231
232			if ((c1->flags\|c2->flags) & C_CDSPEC) {
233			long cmix[C_CNSS], cmax = 0; /* spectral multiply */
234
235			c_ccvt(c1, C_CSSPEC\|C_CSEFF);
236			c_ccvt(c2, C_CSSPEC\|C_CSEFF);
237			for (i = 0; i < C_CNSS; i++) {
238			cmix[i] = c1->ssamp[i] * c2->ssamp[i];
239			if (cmix[i] > cmax)
240			cmax = cmix[i];
241			}
242			cmax /= C_CMAXV;
243			if (!cmax) {
244			*cres = c_dfcolor;
245			return(0.);
246			}
247			cres->ssum = 0;
248			for (i = 0; i < C_CNSS; i++)
249			cres->ssum += cres->ssamp[i] = cmix[i] / cmax;
250			cres->flags = C_CDSPEC\|C_CSSPEC;
251
252			c_ccvt(cres, C_CSEFF); /* nasty, but true */
253			yres = (y1 * y2 * cie_yf.ssum * C_CLPWM) /
254			(c1->effc1->ssum c2->effc2->ssum)
255			cres->eff( cres->ssum(double)cmax +
256			C_CNSS/2.0*(cmax-1) );
257			} else {
258			float rgb1[3], rgb2[3]; /* CIE xy multiply */
259
260			c_toSharpRGB(c1, y1, rgb1);
261			c_toSharpRGB(c2, y2, rgb2);
262
263			rgb2[0] = rgb1[0]; rgb2[1] = rgb1[1]; rgb2[2] *= rgb1[2];
264
265			yres = c_fromSharpRGB(rgb2, cres);
266			}
267			return(yres);
268			}
269
270	greg	3.1	#define C1 3.741832e-16 /* W-m^2 */
271			#define C2 1.4388e-2 /* m-K */
272
273			#define bbsp(l,t) (C1/((l)(l)(l)(l)(l)(exp(C2/((t)(l)))-1.)))
274	greg	3.2	#define bblm(t) (C2*0.2/(t))
275	greg	3.1
276			/* set black body spectrum */
277			int
278			c_bbtemp(C_COLOR *clr, double tk)
279			{
280			double sf, wl;
281			int i;
282
283			if (tk < 1000)
284			return(0);
285			wl = bblm(tk); /* scalefactor based on peak */
286			if (wl < C_CMINWL*1e-9)
287			wl = C_CMINWL*1e-9;
288			else if (wl > C_CMAXWL*1e-9)
289			wl = C_CMAXWL*1e-9;
290			sf = C_CMAXV/bbsp(wl,tk);
291			clr->ssum = 0;
292			for (i = 0; i < C_CNSS; i++) {
293			wl = (C_CMINWL + iC_CWLI)1e-9;
294			clr->ssum += clr->ssamp[i] = sf*bbsp(wl,tk) + .5;
295			}
296			clr->flags = C_CDSPEC\|C_CSSPEC;
297			return(1);
298			}
299
300			#undef C1
301			#undef C2
302			#undef bbsp
303			#undef bblm