CIE.inc

/*
  Persistence of Vision Ray Tracer Include File
  
  COMMISSION INTERNATIONALE DE L'ECLAIRAGE

  CIE XYZ Color Model for PoV-Ray, Version 1.8

  This started as a simple replacement for the spectrum sample calculation
  for Jaimes' lightsys but meanwhile it's a kind of multi purpose CIE color
  transformation system.

  -Ive 2003  email: ive@lilysoft.org
  
  -------------------------------------------------------------------------

  update (version 1.71)
    - corrected some pass-by-reference issue
    - email address updated
    
    -Ive 2012 

  -------------------------------------------------------------------------

  update (version 1.8)
    - added the XYZ2K() macro to calculate the correlated color temperature
    from a given XYZ coordinate.
    - corrected some errors in the comments.
    
    -Ive 2016

  -------------------------------------------------------------------------

  Main features:

  CIE observer
    CIE 1931 2°observer, CIE 1964 10°observer, CIE 1978 Judd/Vos

  CIE reference white
    Illuminant A, B, C, D50, D55, D65, D75, E and all F-types are supported.  

  CIE color system and default illuminant
    with 13 predifined color systems together with the illuminant and the
    possibility to select also 'custom-made' systems.

  CIE color system with custum illuminant
    use any non-default illuminant for the predefined color system
    already 20 CIE standard illuminants are defined and in addition you 
    can convert temperatures to whitepoints.

  CIE color order systems and conversion utilities between 
    RGB, XYZ, xyY, Lab, Lch

  Wavelength, blackbody, daylight(D illuminant), reflective-, emissive-
  and line-spectrum color conversion tools

  Gamut mapping with different functions (for xyz->RGB conversion)

  Chromatic adaption with different functions 
    (used for Lab,Lch <-> xyz and reflectance spectrum -> xyz and
     'reference rgb' conversion)
     
  Grayscale feature 

*/   

// Standard include file beginning stuff 
#ifndef (CIE_Inc_Temp)
#declare CIE_Inc_Temp = version;


#ifdef(View_POV_Include_Stack)
  #debug "Including CIE.inc\n"
#end


// additional CIE.inc settings
#ifndef (CIE_Debug)
  #declare CIE_Debug = false;
#end  

#ifndef (CIE_MultiObserver)
  #declare CIE_MultiObserver = false;
#end

#ifndef (CIE_IntegralStep)
  #declare CIE_IntegralStep = 5; // step width for integral over spectra
#end


//=====================================================================
// first all the lookup tables and pre-definitions
//
/*  CIE observer color matching function tables
 *  xyz spectral coordinates in steps of 5 nm.
 *  see also Fig. 13.24 in [Foley96 p.581]     
 *  taken from http://cvision.ucsd.edu/cie.htm
 *  For consistence and easy access all tables are reduced to the 
 *  range of 380 to 825 nanometers.
 *  The first column is the original 1931 2° observer CMF, the second
 *  the 1964 10° and the third is the 2° observer modified by Judd/Vos
 */
#if (CIE_MultiObserver)
#declare CMF_table = array[90][3] {
{<0.001368000000,0.000039000000,0.006450001000>,<0.000159952000,0.000017364000,0.000704776000>,<0.002689900000,0.000200000000,0.012260000000>},//380
{<0.002236000000,0.000064000000,0.010549990000>,<0.000662440000,0.000071560000,0.002927800000>,<0.005310500000,0.000395560000,0.024222000000>},//385
{<0.004243000000,0.000120000000,0.020050010000>,<0.002361600000,0.000253400000,0.010482200000>,<0.010781000000,0.000800000000,0.049250000000>},//390
{<0.007650000000,0.000217000000,0.036210000000>,<0.007242300000,0.000768500000,0.032344000000>,<0.020792000000,0.001545700000,0.095135000000>},//395
{<0.014310000000,0.000396000000,0.067850010000>,<0.019109700000,0.002004400000,0.086010900000>,<0.037981000000,0.002800000000,0.174090000000>},//400
{<0.023190000000,0.000640000000,0.110200000000>,<0.043400000000,0.004509000000,0.197120000000>,<0.063157000000,0.004656200000,0.290130000000>},//405
{<0.043510000000,0.001210000000,0.207400000000>,<0.084736000000,0.008756000000,0.389366000000>,<0.099941000000,0.007400000000,0.460530000000>},//410
{<0.077630000000,0.002180000000,0.371300000000>,<0.140638000000,0.014456000000,0.656760000000>,<0.158240000000,0.011779000000,0.731660000000>},//415
{<0.134380000000,0.004000000000,0.645600000000>,<0.204492000000,0.021391000000,0.972542000000>,<0.229480000000,0.017500000000,1.065800000000>},//420
{<0.214770000000,0.007300000000,1.039050100000>,<0.264737000000,0.029497000000,1.282500000000>,<0.281080000000,0.022678000000,1.314600000000>},//425
{<0.283900000000,0.011600000000,1.385600000000>,<0.314679000000,0.038676000000,1.553480000000>,<0.310950000000,0.027300000000,1.467200000000>},//430
{<0.328500000000,0.016840000000,1.622960000000>,<0.357719000000,0.049602000000,1.798500000000>,<0.330720000000,0.032584000000,1.579600000000>},//435
{<0.348280000000,0.023000000000,1.747060000000>,<0.383734000000,0.062077000000,1.967280000000>,<0.333360000000,0.037900000000,1.616600000000>},//440
{<0.348060000000,0.029800000000,1.782600000000>,<0.386726000000,0.074704000000,2.027300000000>,<0.316720000000,0.042391000000,1.568200000000>},//445
{<0.336200000000,0.038000000000,1.772110000000>,<0.370702000000,0.089456000000,1.994800000000>,<0.288820000000,0.046800000000,1.471700000000>},//450
{<0.318700000000,0.048000000000,1.744100000000>,<0.342957000000,0.106256000000,1.900700000000>,<0.259690000000,0.052122000000,1.374000000000>},//455
{<0.290800000000,0.060000000000,1.669200000000>,<0.302273000000,0.128201000000,1.745370000000>,<0.232760000000,0.060000000000,1.291700000000>},//460
{<0.251100000000,0.073900000000,1.528100000000>,<0.254085000000,0.152761000000,1.554900000000>,<0.209990000000,0.072942000000,1.235600000000>},//465
{<0.195360000000,0.090980000000,1.287640000000>,<0.195618000000,0.185190000000,1.317560000000>,<0.174760000000,0.090980000000,1.113800000000>},//470
{<0.142100000000,0.112600000000,1.041900000000>,<0.132349000000,0.219940000000,1.030200000000>,<0.132870000000,0.112840000000,0.942200000000>},//475
{<0.095640000000,0.139020000000,0.812950100000>,<0.080507000000,0.253589000000,0.772125000000>,<0.091944000000,0.139020000000,0.755960000000>},//480
{<0.057950010000,0.169300000000,0.616200000000>,<0.041072000000,0.297665000000,0.570600000000>,<0.056985000000,0.169870000000,0.586400000000>},//485
{<0.032010000000,0.208020000000,0.465180000000>,<0.016172000000,0.339133000000,0.415254000000>,<0.031731000000,0.208020000000,0.446690000000>},//490
{<0.014700000000,0.258600000000,0.353300000000>,<0.005132000000,0.395379000000,0.302356000000>,<0.014613000000,0.258080000000,0.341160000000>},//495
{<0.004900000000,0.323000000000,0.272000000000>,<0.003816000000,0.460777000000,0.218502000000>,<0.004849100000,0.323000000000,0.264370000000>},//500
{<0.002400000000,0.407300000000,0.212300000000>,<0.015444000000,0.531360000000,0.159249000000>,<0.002321500000,0.405400000000,0.205940000000>},//505
{<0.009300000000,0.503000000000,0.158200000000>,<0.037465000000,0.606741000000,0.112044000000>,<0.009289900000,0.503000000000,0.154450000000>},//510
{<0.029100000000,0.608200000000,0.111700000000>,<0.071358000000,0.685660000000,0.082248000000>,<0.029278000000,0.608110000000,0.109180000000>},//515
{<0.063270000000,0.710000000000,0.078249990000>,<0.117749000000,0.761757000000,0.060709000000>,<0.063791000000,0.710000000000,0.076585000000>},//520
{<0.109600000000,0.793200000000,0.057250010000>,<0.172953000000,0.823330000000,0.043050000000>,<0.110810000000,0.795100000000,0.056227000000>},//525
{<0.165500000000,0.862000000000,0.042160000000>,<0.236491000000,0.875211000000,0.030451000000>,<0.166920000000,0.862000000000,0.041366000000>},//530
{<0.225749900000,0.914850100000,0.029840000000>,<0.304213000000,0.923810000000,0.020584000000>,<0.227680000000,0.915050000000,0.029353000000>},//535
{<0.290400000000,0.954000000000,0.020300000000>,<0.376772000000,0.961988000000,0.013676000000>,<0.292690000000,0.954000000000,0.020042000000>},//540
{<0.359700000000,0.980300000000,0.013400000000>,<0.451584000000,0.982200000000,0.007918000000>,<0.362250000000,0.980040000000,0.013312000000>},//545
{<0.433449900000,0.994950100000,0.008749999000>,<0.529826000000,0.991761000000,0.003988000000>,<0.436350000000,0.994950000000,0.008782300000>},//550
{<0.512050100000,1.000000000000,0.005749999000>,<0.616053000000,0.999110000000,0.001091000000>,<0.515130000000,1.000100000000,0.005857300000>},//555
{<0.594500000000,0.995000000000,0.003900000000>,<0.705224000000,0.997340000000,0.000000000000>,<0.597480000000,0.995000000000,0.004049300000>},//560
{<0.678400000000,0.978600000000,0.002749999000>,<0.793832000000,0.982380000000,0.000000000000>,<0.681210000000,0.978750000000,0.002921700000>},//565
{<0.762100000000,0.952000000000,0.002100000000>,<0.878655000000,0.955552000000,0.000000000000>,<0.764250000000,0.952000000000,0.002277100000>},//570
{<0.842500000000,0.915400000000,0.001800000000>,<0.951162000000,0.915175000000,0.000000000000>,<0.843940000000,0.915580000000,0.001970600000>},//575
{<0.916300000000,0.870000000000,0.001650001000>,<1.014160000000,0.868934000000,0.000000000000>,<0.916350000000,0.870000000000,0.001806600000>},//580
{<0.978600000000,0.816300000000,0.001400000000>,<1.074300000000,0.825623000000,0.000000000000>,<0.977030000000,0.816230000000,0.001544900000>},//585
{<1.026300000000,0.757000000000,0.001100000000>,<1.118520000000,0.777405000000,0.000000000000>,<1.023000000000,0.757000000000,0.001234800000>},//590
{<1.056700000000,0.694900000000,0.001000000000>,<1.134300000000,0.720353000000,0.000000000000>,<1.051300000000,0.694830000000,0.001117700000>},//595
{<1.062200000000,0.631000000000,0.000800000000>,<1.123990000000,0.658341000000,0.000000000000>,<1.055000000000,0.631000000000,0.000905640000>},//600
{<1.045600000000,0.566800000000,0.000600000000>,<1.089100000000,0.593878000000,0.000000000000>,<1.036200000000,0.566540000000,0.000694670000>},//605
{<1.002600000000,0.503000000000,0.000340000000>,<1.030480000000,0.527963000000,0.000000000000>,<0.992390000000,0.503000000000,0.000428850000>},//610
{<0.938400000000,0.441200000000,0.000240000000>,<0.950740000000,0.461834000000,0.000000000000>,<0.928610000000,0.441720000000,0.000318170000>},//615
{<0.854449900000,0.381000000000,0.000190000000>,<0.856297000000,0.398057000000,0.000000000000>,<0.843460000000,0.381000000000,0.000255980000>},//620
{<0.751400000000,0.321000000000,0.000100000000>,<0.754930000000,0.339554000000,0.000000000000>,<0.739830000000,0.320520000000,0.000156790000>},//625
{<0.642400000000,0.265000000000,0.000049999990>,<0.647467000000,0.283493000000,0.000000000000>,<0.632890000000,0.265000000000,0.000097694000>},//630
{<0.541900000000,0.217000000000,0.000030000000>,<0.535110000000,0.228254000000,0.000000000000>,<0.533510000000,0.217020000000,0.000068944000>},//635
{<0.447900000000,0.175000000000,0.000020000000>,<0.431567000000,0.179828000000,0.000000000000>,<0.440620000000,0.175000000000,0.000051165000>},//640
{<0.360800000000,0.138200000000,0.000010000000>,<0.343690000000,0.140211000000,0.000000000000>,<0.354530000000,0.138120000000,0.000036016000>},//645
{<0.283500000000,0.107000000000,0.000000000000>,<0.268329000000,0.107633000000,0.000000000000>,<0.278620000000,0.107000000000,0.000024238000>},//650
{<0.218700000000,0.081600000000,0.000000000000>,<0.204300000000,0.081187000000,0.000000000000>,<0.214850000000,0.081652000000,0.000016915000>},//655
{<0.164900000000,0.061000000000,0.000000000000>,<0.152568000000,0.060281000000,0.000000000000>,<0.161610000000,0.061000000000,0.000011906000>},//660
{<0.121200000000,0.044580000000,0.000000000000>,<0.112210000000,0.044096000000,0.000000000000>,<0.118200000000,0.044327000000,0.000008148900>},//665
{<0.087400000000,0.032000000000,0.000000000000>,<0.081260600000,0.031800400000,0.000000000000>,<0.085753000000,0.032000000000,0.000005600600>},//670
{<0.063600000000,0.023200000000,0.000000000000>,<0.057930000000,0.022601700000,0.000000000000>,<0.063077000000,0.023454000000,0.000003954400>},//675
{<0.046770000000,0.017000000000,0.000000000000>,<0.040850800000,0.015905100000,0.000000000000>,<0.045834000000,0.017000000000,0.000002791200>},//680
{<0.032900000000,0.011920000000,0.000000000000>,<0.028623000000,0.011130300000,0.000000000000>,<0.032057000000,0.011872000000,0.000001917600>},//685
{<0.022700000000,0.008210000000,0.000000000000>,<0.019941300000,0.007748800000,0.000000000000>,<0.022187000000,0.008210000000,0.000001313500>},//690
{<0.015840000000,0.005723000000,0.000000000000>,<0.013842000000,0.005375100000,0.000000000000>,<0.015612000000,0.005772300000,0.000000915190>},//695
{<0.011359160000,0.004102000000,0.000000000000>,<0.009576880000,0.003717740000,0.000000000000>,<0.011098000000,0.004102000000,0.000000647670>},//700
{<0.008110916000,0.002929000000,0.000000000000>,<0.006605200000,0.002564560000,0.000000000000>,<0.007923300000,0.002929100000,0.000000463520>},//705
{<0.005790346000,0.002091000000,0.000000000000>,<0.004552630000,0.001768470000,0.000000000000>,<0.005653100000,0.002091000000,0.000000333040>},//710
{<0.004106457000,0.001484000000,0.000000000000>,<0.003144700000,0.001222390000,0.000000000000>,<0.004003900000,0.001482200000,0.000000238230>},//715
{<0.002899327000,0.001047000000,0.000000000000>,<0.002174960000,0.000846190000,0.000000000000>,<0.002825300000,0.001047000000,0.000000170260>},//720
{<0.002049190000,0.000740000000,0.000000000000>,<0.001505700000,0.000586440000,0.000000000000>,<0.001994700000,0.000740150000,0.000000122070>},//725
{<0.001439971000,0.000520000000,0.000000000000>,<0.001044760000,0.000407410000,0.000000000000>,<0.001399400000,0.000520000000,0.000000087107>},//730
{<0.000999949300,0.000361100000,0.000000000000>,<0.000727450000,0.000284041000,0.000000000000>,<0.000969800000,0.000360930000,0.000000061455>},//735
{<0.000690078600,0.000249200000,0.000000000000>,<0.000508258000,0.000198730000,0.000000000000>,<0.000668470000,0.000249200000,0.000000043162>},//740
{<0.000476021300,0.000171900000,0.000000000000>,<0.000356380000,0.000139550000,0.000000000000>,<0.000461410000,0.000172310000,0.000000030379>},//745
{<0.000332301100,0.000120000000,0.000000000000>,<0.000250969000,0.000098428000,0.000000000000>,<0.000320730000,0.000120000000,0.000000021554>},//750
{<0.000234826100,0.000084800000,0.000000000000>,<0.000177730000,0.000069819000,0.000000000000>,<0.000225730000,0.000084620000,0.000000015493>},//755
{<0.000166150500,0.000060000000,0.000000000000>,<0.000126390000,0.000049737000,0.000000000000>,<0.000159730000,0.000060000000,0.000000011204>},//760
{<0.000117413000,0.000042400000,0.000000000000>,<0.000090151000,0.000035540500,0.000000000000>,<0.000112750000,0.000042446000,0.000000008087>},//765
{<0.000083075270,0.000030000000,0.000000000000>,<0.000064525800,0.000025486000,0.000000000000>,<0.000079513000,0.000030000000,0.000000005834>},//770
{<0.000058706520,0.000021200000,0.000000000000>,<0.000046339000,0.000018338400,0.000000000000>,<0.000056087000,0.000021210000,0.000000004211>},//775
{<0.000041509940,0.000014990000,0.000000000000>,<0.000033411700,0.000013249000,0.000000000000>,<0.000039541000,0.000014989000,0.000000003038>},//780
{<0.000029353260,0.000010600000,0.000000000000>,<0.000024209000,0.000009619600,0.000000000000>,<0.000027852000,0.000010584000,0.000000002191>},//785
{<0.000020673830,0.000007465700,0.000000000000>,<0.000017611500,0.000007012800,0.000000000000>,<0.000019597000,0.000007465600,0.000000001578>},//790
{<0.000014559770,0.000005257800,0.000000000000>,<0.000012855000,0.000005129800,0.000000000000>,<0.000013770000,0.000005259200,0.000000001135>},//795
{<0.000010253980,0.000003702900,0.000000000000>,<0.000009413630,0.000003764730,0.000000000000>,<0.000009670000,0.000003702800,0.000000000816>},//800
{<0.000007221456,0.000002607800,0.000000000000>,<0.000006913000,0.000002770810,0.000000000000>,<0.000006791800,0.000002607600,0.000000000586>},//805
{<0.000005085868,0.000001836600,0.000000000000>,<0.000005093470,0.000002046130,0.000000000000>,<0.000004770600,0.000001836500,0.000000000421>},//810
{<0.000003581652,0.000001293400,0.000000000000>,<0.000003767100,0.000001516770,0.000000000000>,<0.000003355000,0.000001295000,0.000000000303>},//815
{<0.000002522525,0.000000910930,0.000000000000>,<0.000002795310,0.000001128090,0.000000000000>,<0.000002353400,0.000000910920,0.000000000218>},//820
{<0.000001776509,0.000000641530,0.000000000000>,<0.000002082000,0.000000842160,0.000000000000>,<0.000001637700,0.000000635640,0.000000000155>} //825
}
#end // if MultiObserver

/*  CIE observer color matching function 
 *  initialized for the 1931 2° standard observer.
 *  This spline is automatical reinitialized if
 *  any other CMF is selected. 
 *  See the macro CIE_Observer for more details. 
 */
#declare CMF_xyz = spline { cubic_spline  
  375, <0,0,0>
  380, <0.001368000000, 0.000039000000, 0.006450001000>
  385, <0.002236000000, 0.000064000000, 0.010549990000>
  390, <0.004243000000, 0.000120000000, 0.020050010000>
  395, <0.007650000000, 0.000217000000, 0.036210000000>
  400, <0.014310000000, 0.000396000000, 0.067850010000>
  405, <0.023190000000, 0.000640000000, 0.110200000000>
  410, <0.043510000000, 0.001210000000, 0.207400000000>
  415, <0.077630000000, 0.002180000000, 0.371300000000>
  420, <0.134380000000, 0.004000000000, 0.645600000000>
  425, <0.214770000000, 0.007300000000, 1.039050100000>
  430, <0.283900000000, 0.011600000000, 1.385600000000>
  435, <0.328500000000, 0.016840000000, 1.622960000000>
  440, <0.348280000000, 0.023000000000, 1.747060000000>
  445, <0.348060000000, 0.029800000000, 1.782600000000>
  450, <0.336200000000, 0.038000000000, 1.772110000000>
  455, <0.318700000000, 0.048000000000, 1.744100000000>
  460, <0.290800000000, 0.060000000000, 1.669200000000>
  465, <0.251100000000, 0.073900000000, 1.528100000000>
  470, <0.195360000000, 0.090980000000, 1.287640000000>
  475, <0.142100000000, 0.112600000000, 1.041900000000>
  480, <0.095640000000, 0.139020000000, 0.812950100000>
  485, <0.057950010000, 0.169300000000, 0.616200000000>
  490, <0.032010000000, 0.208020000000, 0.465180000000>
  495, <0.014700000000, 0.258600000000, 0.353300000000>
  500, <0.004900000000, 0.323000000000, 0.272000000000>
  505, <0.002400000000, 0.407300000000, 0.212300000000>
  510, <0.009300000000, 0.503000000000, 0.158200000000>
  515, <0.029100000000, 0.608200000000, 0.111700000000>
  520, <0.063270000000, 0.710000000000, 0.078249990000>
  525, <0.109600000000, 0.793200000000, 0.057250010000>
  530, <0.165500000000, 0.862000000000, 0.042160000000>
  535, <0.225749900000, 0.914850100000, 0.029840000000>
  540, <0.290400000000, 0.954000000000, 0.020300000000>
  545, <0.359700000000, 0.980300000000, 0.013400000000>
  550, <0.433449900000, 0.994950100000, 0.008749999000>
  555, <0.512050100000, 1.000000000000, 0.005749999000>
  560, <0.594500000000, 0.995000000000, 0.003900000000>
  565, <0.678400000000, 0.978600000000, 0.002749999000>
  570, <0.762100000000, 0.952000000000, 0.002100000000>
  575, <0.842500000000, 0.915400000000, 0.001800000000>
  580, <0.916300000000, 0.870000000000, 0.001650001000>
  585, <0.978600000000, 0.816300000000, 0.001400000000>
  590, <1.026300000000, 0.757000000000, 0.001100000000>
  595, <1.056700000000, 0.694900000000, 0.001000000000>
  600, <1.062200000000, 0.631000000000, 0.000800000000>
  605, <1.045600000000, 0.566800000000, 0.000600000000>
  610, <1.002600000000, 0.503000000000, 0.000340000000>
  615, <0.938400000000, 0.441200000000, 0.000240000000>
  620, <0.854449900000, 0.381000000000, 0.000190000000>
  625, <0.751400000000, 0.321000000000, 0.000100000000>
  630, <0.642400000000, 0.265000000000, 0.000049999990>
  635, <0.541900000000, 0.217000000000, 0.000030000000>
  640, <0.447900000000, 0.175000000000, 0.000020000000>
  645, <0.360800000000, 0.138200000000, 0.000010000000>
  650, <0.283500000000, 0.107000000000, 0.000000000000>
  655, <0.218700000000, 0.081600000000, 0.000000000000>
  660, <0.164900000000, 0.061000000000, 0.000000000000>
  665, <0.121200000000, 0.044580000000, 0.000000000000>
  670, <0.087400000000, 0.032000000000, 0.000000000000>
  675, <0.063600000000, 0.023200000000, 0.000000000000>
  680, <0.046770000000, 0.017000000000, 0.000000000000>
  685, <0.032900000000, 0.011920000000, 0.000000000000>
  690, <0.022700000000, 0.008210000000, 0.000000000000>
  695, <0.015840000000, 0.005723000000, 0.000000000000>
  700, <0.011359160000, 0.004102000000, 0.000000000000>
  705, <0.008110916000, 0.002929000000, 0.000000000000>
  710, <0.005790346000, 0.002091000000, 0.000000000000>
  715, <0.004106457000, 0.001484000000, 0.000000000000>
  720, <0.002899327000, 0.001047000000, 0.000000000000>
  725, <0.002049190000, 0.000740000000, 0.000000000000>
  730, <0.001439971000, 0.000520000000, 0.000000000000>
  735, <0.000999949300, 0.000361100000, 0.000000000000>
  740, <0.000690078600, 0.000249200000, 0.000000000000>
  745, <0.000476021300, 0.000171900000, 0.000000000000>
  750, <0.000332301100, 0.000120000000, 0.000000000000>
  755, <0.000234826100, 0.000084800000, 0.000000000000>
  760, <0.000166150500, 0.000060000000, 0.000000000000>
  765, <0.000117413000, 0.000042400000, 0.000000000000>
  770, <0.000083075270, 0.000030000000, 0.000000000000>
  775, <0.000058706520, 0.000021200000, 0.000000000000>
  780, <0.000041509940, 0.000014990000, 0.000000000000>
  785, <0.000029353260, 0.000010600000, 0.000000000000>
  790, <0.000020673830, 0.000007465700, 0.000000000000>
  795, <0.000014559770, 0.000005257800, 0.000000000000>
  800, <0.000010253980, 0.000003702900, 0.000000000000>
  805, <0.000007221456, 0.000002607800, 0.000000000000>
  810, <0.000005085868, 0.000001836600, 0.000000000000>
  815, <0.000003581652, 0.000001293400, 0.000000000000>
  820, <0.000002522525, 0.000000910930, 0.000000000000>
  825, <0.000001776509, 0.000000641530, 0.000000000000>
  830, <0,0,0>
}
  
/* Lookup table containing S0, S1, S2 to calculate the spectra for the 
 * CIE D standard illuminants (D50, D55, D65, D75) and also custum D 
 * illuminants for color temperatures between 4000K and 25000K.
 * See the macro Daylight for the formula and more details.
 */
#declare DS012 = array[90][3] {
  { 63.40, 38.50,  3.00},// 380
  { 64.60, 36.75,  2.10},// 385
  { 65.80, 35.00,  1.20},// 390
  { 80.30, 39.20,  0.05},// 395
  { 94.80, 43.40, -1.10},// 400
  { 99.80, 44.85, -0.80},// 405
  {104.80, 46.30, -0.50},// 410
  {105.35, 45.10, -0.60},// 415
  {105.90, 43.90, -0.70},// 420
  {101.35, 40.50, -0.95},// 425
  { 96.80, 37.10, -1.20},// 430
  {105.35, 36.90, -1.90},// 435
  {113.90, 36.70, -2.60},// 440
  {119.75, 36.30, -2.75},// 445
  {125.60, 35.90, -2.90},// 450
  {125.55, 34.25, -2.85},// 455 
  {125.50, 32.60, -2.80},// 460
  {123.40, 30.25, -2.70},// 465
  {121.30, 27.90, -2.60},// 470
  {121.30, 26.10, -2.60},// 475
  {121.30, 24.30, -2.60},// 480
  {117.40, 22.20, -2.20},// 485
  {113.50, 20.10, -1.80},// 490
  {113.30, 18.15, -1.65},// 495
  {113.10, 16.20, -1.50},// 500
  {111.95, 14.70, -1.40},// 505
  {110.80, 13.20, -1.30},// 510
  {108.65, 10.90, -1.25},// 515
  {106.50,  8.60, -1.20},// 520
  {107.65,  7.35, -1.10},// 525
  {108.80,  6.10, -1.00},// 530
  {107.05,  5.15, -0.75},// 535
  {105.30,  4.20, -0.50},// 540
  {104.85,  3.05, -0.40},// 545
  {104.40,  1.90, -0.30},// 550
  {102.20,  0.95, -0.15},// 555
  {100.00,  0.00,  0.00},// 560
  { 98.00, -0.80,  0.10},// 565
  { 96.00, -1.60,  0.20},// 570
  { 95.55, -2.55,  0.35},// 575
  { 95.10, -3.50,  0.50},// 580
  { 92.10, -3.50,  1.30},// 585
  { 89.10, -3.50,  2.10},// 590
  { 89.80, -4.65,  2.65},// 595
  { 90.50, -5.80,  3.20},// 600
  { 90.40, -6.50,  3.65},// 605
  { 90.30, -7.20,  4.10},// 610
  { 89.35, -7.90,  4.40},// 615
  { 88.40, -8.60,  4.70},// 620
  { 86.20, -9.05,  4.90},// 625
  { 84.00, -9.50,  5.10},// 630
  { 84.55,-10.20,  5.90},// 635
  { 85.10,-10.90,  6.70},// 640
  { 83.50,-10.80,  7.00},// 645
  { 81.90,-10.70,  7.30},// 650
  { 82.25,-11.35,  7.95},// 655
  { 82.60,-12.00,  8.60},// 660
  { 83.75,-13.00,  9.20},// 665
  { 84.90,-14.00,  9.80},// 670
  { 83.10,-13.80, 10.00},// 675
  { 81.30,-13.60, 10.20},// 680
  { 76.60,-12.80,  9.25},// 685
  { 71.90,-12.00,  8.30},// 690
  { 73.10,-12.65,  8.95},// 695
  { 74.30,-13.30,  9.60},// 700
  { 75.35,-13.10,  9.05},// 705
  { 76.40,-12.90,  8.50},// 710
  { 69.85,-11.75,  7.75},// 715
  { 63.30,-10.60,  7.00},// 720
  { 67.50,-11.10,  7.30},// 725
  { 71.70,-11.60,  7.60},// 730
  { 74.35,-11.90,  7.80},// 735
  { 77.00,-12.20,  8.00},// 740
  { 71.10,-11.20,  7.35},// 745
  { 65.20,-10.20,  6.70},// 750
  { 56.45, -9.00,  5.95},// 755
  { 47.70, -7.80,  5.20},// 760
  { 58.15, -9.50,  6.30},// 765
  { 68.60,-11.20,  7.40},// 770
  { 66.80,-10.80,  7.10},// 775
  { 65.00,-10.40,  6.80},// 780
  { 65.50,-10.50,  6.90},// 785
  { 66.00,-10.60,  7.00},// 790
  { 63.50,-10.15,  6.70},// 795
  { 61.00, -9.70,  6.40},// 800
  { 57.15, -9.00,  5.95},// 805
  { 53.30, -8.30,  5.50},// 810
  { 56.10, -8.80,  5.80},// 815
  { 58.90, -9.30,  6.10},// 820
  { 60.40, -9.55,  6.30},// 825
}

/*  RefWhiteSP is initialized to D50 and recalculated if
 *  any other D illuminant or the A illuminant is requested.
 *  In every other case it is set equal to the the CIE standard
 *  illuminants as stored in the file "espd_cie_standard.inc".
 */ 
#declare RefWhiteSP = spline{ linear_spline
  380, 0.244992
  385, 0.271905
  390, 0.298818
  395, 0.396023
  400, 0.493228
  405, 0.529255
  410, 0.565283
  415, 0.582884
  420, 0.600485
  425, 0.589394
  430, 0.578302
  435, 0.663340
  440, 0.748379
  445, 0.810490
  450, 0.872601
  455, 0.889420
  460, 0.906239
  465, 0.910011
  470, 0.913782
  475, 0.932478
  480, 0.951174
  485, 0.935437
  490, 0.919700
  495, 0.938498
  500, 0.957296
  505, 0.961738
  510, 0.966181
  515, 0.968752
  520, 0.971323
  525, 0.996169
  530, 1.021015
  535, 1.014289
  540, 1.007563
  545, 1.015370
  550, 1.023177
  555, 1.011589
  560, 1.000000
  565, 0.988672
  570, 0.977344
  575, 0.983255
  580, 0.989167
  585, 0.962068
  590, 0.934969
  595, 0.955909
  600, 0.976848
  605, 0.984751
  610, 0.992653
  615, 0.991512
  620, 0.990371
  625, 0.973770
  630, 0.957169
  635, 0.972841
  640, 0.988513
  645, 0.972562
  650, 0.956612
  655, 0.969220
  660, 0.981829
  665, 1.005892
  670, 1.029954
  675, 1.010602
  680, 0.991250
  685, 0.932496
  690, 0.873741
  695, 0.894850
  700, 0.915958
  705, 0.922386
  710, 0.928814
  715, 0.848650
  720, 0.768485
  725, 0.816766
  730, 0.865048
  735, 0.895389
  740, 0.925730
  745, 0.853986
  750, 0.782242
  755, 0.679558
  760, 0.576875
}


/*  xy reference values for perfect diffuse reflection used as 
 *  predefined whitepoints for the 2° and the 10° observer.
 *  Whitepoint format:
 *    (x-white, y-white)   for the 1931  2° observers
 *    (x-white, y-white)   for the 1964 10° observer
 *    (x-white, y-white)   for the Judd/Vos 2° observer
 *  You can also use kelvin temperatures to create custom-whitepoints
 *  by using the Blackbody2Whitepoint or Daylight2Whitepoint macros.
 * 
 *  Due to a revision of the Boltzmann constant we have the correlated
 *  Kelvin temperature for the daylight illuminants: (1.4388/1.4380)*D(K)
 *  D50 - 5002.78 K
 *  D55 - 5503.06 K
 *  D65 - 6503.62 K
 *  D75 - 7504.17 K
 */
#if (CIE_MultiObserver) 
#declare Illuminant_A = array[3][2] {   // blackbody 2856 K
  {0.44753, 0.40744}, // 1931  2°observer
  {0.45113, 0.40593}, // 1964 10°observer
  {0.44927, 0.41273}  // 1978  2°observer
}
#declare Illuminant_B = array[3][2] {   // ca. blackbody 4871 K
  {0.34849, 0.35173}, // 1931  2°observer
  {0.34990, 0.35279}, // 1964 10°observer
  {0.35189, 0.35853}  // 1978  2°observer
}
#declare Illuminant_C = array[3][2] {   // ca. blackbody 6774 K
  {0.31006, 0.31616}, // 1931  2°observer
  {0.31038, 0.31905}, // 1964 10°observer
  {0.31373, 0.32291}  // 1978  2°observer
}
#declare Illuminant_D50 = array[3][2] { // daylight 5002.78 K
  {0.34566, 0.35850}, // 1931  2°observer
  {0.34771, 0.35951}, // 1964 10°observer
  {0.34919, 0.36560}  // 1978  2°observer
}
#declare Illuminant_D55 = array[3][2] { // daylight 5503.06 K
  {0.33242, 0.34744}, // 1931  2°observer
  {0.33411, 0.34876}, // 1964 10°observer
  {0.33587, 0.35417}  // 1978  2°observer
}
#declare Illuminant_D65 = array[3][2] { // daylight 6503.62 K
  {0.31271, 0.32902}, // 1931  2°observer
  {0.31379, 0.33096}, // 1964 10°observer
  {0.31598, 0.33502}  // 1978  2°observer
}
#declare Illuminant_D75 = array[3][2] { // daylight 7504.17 K
  {0.29903, 0.31487}, // 1931  2°observer
  {0.29967, 0.31739}, // 1964 10°observer
  {0.30216, 0.32027}  // 1978  2°observer
}
#declare Illuminant_E = array[3][2] {  // equal energy 5469 K
  {1/3, 1/3},       // 1931  2°observer
  {1/3, 1/3}        // 1964 10°observer
  {1/3, 1/3}        // 1978  2°observer
} 
#declare Illuminant_F1 = array[3][2] { // fluorescent 6430 K
  {0.31306, 0.33711}, //  2°observer
  {0.31810, 0.33549}, // 10°observer
  {0.31690, 0.34431}  //  2°observer
}
#declare Illuminant_F2 = array[3][2] { // fluorescent 4230 K
  {0.37207, 0.37512}, //  2°observer
  {0.37927, 0.36723}, // 10°observer
  {0.37463, 0.38025}  //  2°observer
}
#declare Illuminant_F3 = array[3][2] { // fluorescent 3450 K
  {0.40909, 0.39412}, //  2°observer
  {0.41764, 0.38312}, // 10°observer
  {0.41055, 0.39782}  //  2°observer
}
#declare Illuminant_F4 = array[3][2] { // fluorescent 2940 K
  {0.43987, 0.40314}, //  2°observer
  {0.44895, 0.39064}, // 10°observer
  {0.44035, 0.40579}  //  2°observer
}
#declare Illuminant_F5 = array[3][2] { // fluorescent 6350 K
  {0.31376, 0.34516}, //  2°observer
  {0.31974, 0.34236}, // 10°observer
  {0.31749, 0.35215}  //  2°observer
}
#declare Illuminant_F6 = array[3][2] { // fluorescent 4150 K
  {0.37788, 0.38819}, //  2°observer
  {0.38662, 0.37837}, // 10°observer
  {0.38016, 0.39285}  //  2°observer
}
#declare Illuminant_F7 = array[3][2] { // broad band fluorescent 6500 K
  {0.31285, 0.32917}, //  2°observer
  {0.31564, 0.32951}, // 10°observer
  {0.31646, 0.33629}  //  2°observer
}
#declare Illuminant_F8 = array[3][2] { // broad band fluorescent 5000 K
  {0.34581, 0.35862}, //  2°observer
  {0.34896, 0.35931}, // 10°observer
  {0.34975, 0.36697}  //  2°observer
}
#declare Illuminant_F9 = array[3][2] { // broad band fluorescent 4150 K
  {0.37410, 0.37267}, //  2°observer
  {0.37826, 0.37037}, // 10°observer
  {0.37738, 0.37978}  //  2°observer
}
#declare Illuminant_F10 = array[3][2] { // narrow band fluorescent 5000 K
  {0.34579, 0.35876}, //  2°observer
  {0.35061, 0.35430}, // 10°observer
  {0.35014, 0.36779}  //  2°observer
}
#declare Illuminant_F11 = array[3][2] { // narrow band fluorescent 4000 K
  {0.38054, 0.37692}, //  2°observer
  {0.38543, 0.37109}, // 10°observer
  {0.38363, 0.38392}  //  2°observer
}
#declare Illuminant_F12 = array[3][2] { // narrow band fluorescent 3000 K
  {0.43702, 0.40421}, //  2°observer
  {0.44265, 0.39706}, // 10°observer
  {0.43752, 0.40780}  //  2°observer
}

#else

#declare Illuminant_A   = array[2] {0.44753,0.40744} // blackbody 2856 K
#declare Illuminant_B   = array[2] {0.34849,0.35173} // ca. blackbody 4871 K
#declare Illuminant_C   = array[2] {0.31006,0.31616} // ca. blackbody 6774 K
#declare Illuminant_D50 = array[2] {0.34566,0.35850} // daylight 5002.78 K
#declare Illuminant_D55 = array[2] {0.33242,0.34744} // daylight 5503.06 K
#declare Illuminant_D65 = array[2] {0.31271,0.32902} // daylight 6503.62 K
#declare Illuminant_D75 = array[2] {0.29903,0.31487} // daylight 7504.17 K
#declare Illuminant_E   = array[2] {1/3, 1/3}        // equal energy 5469 K
#declare Illuminant_F1  = array[2] {0.31306,0.33711} // fluorescent 6430 K
#declare Illuminant_F2  = array[2] {0.37207,0.37512} // fluorescent 4230 K
#declare Illuminant_F3  = array[2] {0.40909,0.39412} // fluorescent 3450 K
#declare Illuminant_F4  = array[2] {0.43987,0.40314} // fluorescent 2940 K
#declare Illuminant_F5  = array[2] {0.31376,0.34516} // fluorescent 6350 K
#declare Illuminant_F6  = array[2] {0.37788,0.38819} // fluorescent 4150 K
#declare Illuminant_F7  = array[2] {0.31285,0.32917} // broad band fluorescent 6500 K
#declare Illuminant_F8  = array[2] {0.34581,0.35862} // broad band fluorescent 5000 K
#declare Illuminant_F9  = array[2] {0.37410,0.37267} // broad band fluorescent 4150 K
#declare Illuminant_F10 = array[2] {0.34579,0.35876} // narrow band fluorescent 5000 K
#declare Illuminant_F11 = array[2] {0.38054,0.37692} // narrow band fluorescent 4000 K
#declare Illuminant_F12 = array[2] {0.43702,0.40421} // narrow band fluorescent 3000 K
  
#end


/*  Color System xy coordinates and whitepoints
 *  sRGB:   CCIR Recommendation 709   (used by HP, Microsoft)
 *  CIE:    CIE 'ideal' gamut where the corners are on the contour of the horseshoe
 *  ITU:    CCIR Report 476-1         (YUV color order, Pal Secam, european tv)
 *  NTSC:   CCIR Recommendation 601-1 (YIQ color order, american tv)
 *  SMPTE:  CIE D65
 *  Adobe:  Adobe (Photoshop ICM)
 *  Match:  Adobe ColorMatch (Photoshop ICM)
 *  Apple:  Apple Monitor (Photoshop ICM)
 *  Beta:   RGB working space created by Bruce Lindbloom 
              http://www.brucelindbloom.com/BetaRGB.html
 *  Dell phosphors: Dell data sheet
 *  Short-persistence phosphors: [Foley96, p.583]
 *  Long-persistence phosphors:  [Foley96, p.583]
 *  HOT: Hydrogen, Oxygen, Thermal color space uses primaries found in deep space
 *       see http://www.nightscapes.net/techniques/HOTcolors/HOTColorSpace.html
 *
 *  Color System format:
 *   (x-red,   y-red),
 *   (x-green, y-green),
 *   (x-blue,  y-blue),
 *   (x-white, y-white), whitepoint for 1931 2°observer
 *   (x-white, y-white), whitepoint for 1964 10°observer
 *   (x-white, y-white)  whitepoint for 1978 2°observer (Judd/Vos)
 */
#if (CIE_MultiObserver)
 
#declare sRGB_ColSys = array[6][2] { // <- initialized as default
  {0.640,  0.330},  // red
  {0.300,  0.600},  // green
  {0.150,  0.060},  // blue  
  {0.31271,0.32902},// Illuminant D65  1931 observer
  {0.31379,0.33096},// Illuminant D65  1964
  {0.31598,0.33502} // Illuminant D65  1978
}
#declare CIE_ColSys = array[6][2] {
  {0.7347, 0.2653}, // red
  {0.2738, 0.7174}, // green
  {0.1666, 0.0088}, // blue
  {1/3, 1/3},       // Illuminant E  1931 observer
  {1/3, 1/3},       // Illuminant E  1964
  {1/3, 1/3}        // Illuminant E  1978
}  
#declare ITU_ColSys = array[6][2] {
  {0.640,  0.330},  // red
  {0.290,  0.600},  // green
  {0.150,  0.060},  // blue
  {0.31271,0.32902},// Illuminant D65  1931 observer
  {0.31379,0.33096},// Illuminant D65  1964
  {0.31598,0.33502} // Illuminant D65  1978
}
#declare NTSC_ColSys = array[6][2] {
  {0.670,  0.330},  // red
  {0.210,  0.710},  // green
  {0.140,  0.080},  // blue
  {0.31006,0.31616},// Illuminant C  1931 observer
  {0.31038,0.31905},// Illuminant C  1964
  {0.31373,0.32291} // Illuminant C  1978
}
#declare SMPTE_ColSys = array[6][2] {
  {0.630,  0.340},  // red
  {0.310,  0.595},  // green
  {0.155,  0.070},  // blue
  {0.31271,0.32902},// Illuminant D65  1931 observer
  {0.31379,0.33096},// Illuminant D65  1964
  {0.31598,0.33502} // Illuminant D65  1978
}
#declare Adobe_ColSys = array[6][2] {
  {0.640,  0.330},  // red
  {0.210,  0.710},  // green
  {0.150,  0.060},  // blue
  {0.31271,0.32902},// Illuminant D65  1931 observer
  {0.31379,0.33096},// Illuminant D65  1964
  {0.31598,0.33502} // Illuminant D65  1978
}
#declare Match_ColSys = array[6][2] {
  {0.630,  0.340},  // red
  {0.290,  0.610},  // green
  {0.150,  0.075},  // blue
  {0.34566,0.35850},// Illuminant D50  1931 observer
  {0.34771,0.35951},// Illuminant D50  1964
  {0.34919,0.36560} // Illuminant D50  1978
}
#declare Apple_ColSys = array[6][2] {
  {0.630, 0.340},   // red
  {0.280, 0.600},   // green
  {0.160, 0.070},   // blue
  {0.31271,0.32902},// Illuminant D65  1931 observer
  {0.31379,0.33096},// Illuminant D65  1964
  {0.31598,0.33502} // Illuminant D65  1978
}
#declare Beta_ColSys = array[6][2] {
  {0.6888, 0.3112}, // red
  {0.1986, 0.7551}, // green
  {0.1265, 0.0352}, // blue
  {0.34566,0.35850},// Illuminant D50  1931 observer
  {0.34771,0.35951},// Illuminant D50  1964
  {0.34919,0.36560} // Illuminant D50  1978
}
#declare Dell_ColSys = array[6][2] {
  {0.625,  0.340},  // red
  {0.275,  0.605},  // green
  {0.150,  0.065},  // blue
  {0.31271,0.32902},// Illuminant D65  1931 observer
  {0.31379,0.33096},// Illuminant D65  1964
  {0.31598,0.33502} // Illuminant D65  1978
}
#declare ShortPersistence_ColSys = array[6][2] {
  {0.610,  0.350},  // red
  {0.290,  0.590},  // green
  {0.150,  0.063},  // blue
  {0.31006,0.31616},// Illuminant C  1931 observer
  {0.31038,0.31905},// Illuminant C  1964
  {0.31373,0.32291} // Illuminant C  1978
}
#declare LongPersistence_ColSys = array[6][2] {
  {0.620,  0.330},  // red
  {0.210,  0.685},  // green
  {0.150,  0.063},  // blue
  {0.31006,0.31616},// Illuminant C  1931 observer
  {0.31038,0.31905},// Illuminant C  1964
  {0.31373,0.32291} // Illuminant C  1978
}
#declare HOT_ColSys = array[6][2] {
  {0.7287, 0.2713}, // red
  {0.0082, 0.5384}, // green
  {0.1491, 0.0544}, // blue
  {0.23991,0.23414},// infinite temperature  1931 observer
  {0.23855,0.23873},// infinite temperature  1964
  {0.24141,0.23253} // infinite temperature  1978
}             

#else //no MultiObserver

/*  Simplified version without the possibility to change the observer
 *  Color System format:
 *   (x-red,   y-red),
 *   (x-green, y-green),
 *   (x-blue,  y-blue),
 *   (x-white, y-white)  whitepoint for 1931 2°observer
 */
#declare sRGB_ColSys = array[4][2] { // <- initialized as default
  {0.640,  0.330},  // red
  {0.300,  0.600},  // green
  {0.150,  0.060},  // blue
  {0.31271,0.32902} // Illuminant D65
}
#declare CIE_ColSys = array[4][2] {
  {0.7347, 0.2653}, // red
  {0.2738, 0.7174}, // green
  {0.1666, 0.0088}, // blue
  {1/3, 1/3}        // Illuminant E
}
#declare ITU_ColSys = array[4][2] {
  {0.640,  0.330},  // red
  {0.290,  0.600},  // green
  {0.150,  0.060},  // blue
  {0.31271,0.32902} // Illuminant D65
}
#declare NTSC_ColSys = array[4][2] {
  {0.670,  0.330},  // red
  {0.210,  0.710},  // green
  {0.140,  0.080},  // blue
  {0.31006,0.31616} // Illuminant C
}
#declare SMPTE_ColSys = array[4][2] {
  {0.630,  0.340},  // red
  {0.310,  0.595},  // green
  {0.155,  0.070},  // blue
  {0.31271,0.32902} // Illuminant D65
}
#declare Adobe_ColSys = array[4][2] {
  {0.640,  0.330},  // red
  {0.210,  0.710},  // green
  {0.150,  0.060},  // blue
  {0.31271,0.32902} // Illuminant D65
}
#declare Match_ColSys = array[4][2] {
  {0.630,  0.340},  // red
  {0.290,  0.610},  // green
  {0.150,  0.075},  // blue
  {0.34566,0.35850} // Illuminant D50
}
#declare Apple_ColSys = array[4][2] {
  {0.630,  0.340},  // red
  {0.280,  0.600},  // green
  {0.160,  0.070},  // blue
  {0.31271,0.32902} // Illuminant D65
}
#declare Beta_ColSys = array[4][2] {
  {0.6888, 0.3112}, // red
  {0.1986, 0.7551}, // green
  {0.1265, 0.0352}, // blue
  {0.34566,0.35850} // Illuminant D50
}
#declare Dell_ColSys = array[4][2] {
  {0.625,  0.340},  // red
  {0.275,  0.605},  // green
  {0.150,  0.065},  // blue
  {0.31271,0.32902} // Illuminant D65
}
#declare ShortPersistence_ColSys = array[4][2] {
  {0.610,  0.350},  // red
  {0.290,  0.590},  // green
  {0.150,  0.063},  // blue
  {0.31006,0.31616} // Illuminant C
}
#declare LongPersistence_ColSys = array[4][2] {
  {0.620,  0.330},  // red
  {0.210,  0.685},  // green
  {0.150,  0.063},  // blue
  {0.31096,0.31616} // Illuminant C
}
#declare HOT_ColSys = array[4][2] {
  {0.7287, 0.2713}, // red
  {0.0082, 0.5384}, // green
  {0.1491, 0.0544}, // blue
  {0.23991,0.23414} // infinite temperature
}             
#end //if MultiObserver


/*  Chroma adaption matrices
 *    Scaling
 *    Von Kries
 *    Bradford
 */
#declare Scaling_MM = array[3][3] { // XYZ Scaling Chroma Match
  {1.0, 0.0, 0.0},
  {0.0, 1.0, 0.0},
  {0.0, 0.0, 1.0}
}
#declare VonKries_MM = array[3][3] { // Von Kries Chroma Match
  { 0.40024,  0.70760, -0.08081},
  {-0.22630,  1.16532,  0.04570},
  { 0.00000,  0.00000,  0.91822}
}
#declare Bradford_MM = array[3][3] { // Bradford Chroma Match
  { 0.8951,   0.2664, -0.1614},
  {-0.7502,   1.7135,  0.0367},     
  { 0.0389,  -0.0685,  1.0296}
} 


//=====================================================================
// CIE color system default settings
//=====================================================================

/*  definition for the color matching functions CMF
 */
#ifdef (CIE_MultiObserver)  
 #declare CIE_1931 = 0;
 #declare CIE_1964 = 1;
 #declare CIE_1978 = 2;
#end

/* definition for chromatic adaption functions
 */
#declare Scaling_ChromaMatch  = 1;
#declare VonKries_ChromaMatch = 2;
#declare Bradford_ChromaMatch = 3;

/*  hard coded and precalculated values
 *  these values are *local* so do not change them, if necessary they are
 *  automatical recalculated by using the appropriate macros
 */ 
#declare CIE_Spectrum = 0;                           // 1931 2° Observer
#declare xy_ColSys    = sRGB_ColSys;                 // Color system
#declare ReferenceWhiteXY  = Illuminant_D50;         // reference white
#declare ReferenceWhiteXYZ = <0.96418, 1.0, 0.82522>;// XYZ for D50 CMF 1931 2° observer
#declare ChromaMatchFunction = Bradford_ChromaMatch; // chroma adaption function
#declare GamutMapFunction = 1;                       // clipping

/*  Color system dependent matrix for rgb->xyz and xyz->rgb
 *  conversion. Hard coded and pre-calculated for the sRGB
 *  color system, illuminant D65 and CIE 1931 2°standard
 *  observer (the default setting).
 *  CCIR recommendation 709
 */
#declare csMatXYZ = array[3][3] {     
   { 0.48499, 0.34891, 0.13029},
   { 0.25007, 0.69781, 0.05211},
   { 0.02273, 0.11630, 0.68618}
}
#declare csMatRGB = array[3][3] {  
   { 3.24102,-1.53740,-0.49862},
   {-0.96922, 1.87593, 0.04155},
   { 0.05564,-0.20401, 1.05714}
}

/*  Adaption function dependent matrix for chroma matching
 *  Hard coded and pre-calculated for the Bradford match,
 *  reference white D50, 2° observer and sRGB color system
 *  with D65
 */
#declare ChromaMatS = array[3][3] {   
  { 0.955565, -0.023045,  0.063180},
  {-0.028296,  1.009942,  0.021013},
  { 0.012301, -0.020488,  1.329995}
}
#declare ChromaMatD = array[3][3] {
  { 1.047825,  0.022892, -0.050137},
  { 0.029550,  0.990484, -0.017053},
  {-0.009236,  0.015047,  0.752084}
}


//=====================================================================
// internal helper macros
//=====================================================================

// convert a color component from the active color system to xyz
// where: z = 1 - x - y
// as the whitepoint depends on the CIE observer, it is also
// checked for the requested illuminant
#macro ColSys2xyz(Idx)
  #if ((Idx>3)|(Idx<0)) #error "CIE.inc: Invalid ColSys index" #end
 #if (CIE_MultiObserver & (Idx=3))
  #local Idx=Idx+CIE_Spectrum;
 #end  
  #local X = xy_ColSys[Idx][0];
  #local Y = xy_ColSys[Idx][1];
  (<X, Y, 1.0-X-Y>)
#end


// chroma match matrix calculation:  M * [Xd/Xs Yd/Ys Zd/Zs] * Mi
#macro ApplyMatchMat(M, Vd, Vs)  
  #local V = ApplyMat(M,Vd)/ApplyMat(M,Vs);
  #local D = array[3][3] {
    {V.x, 0, 0},
    {0, V.y, 0},
    {0, 0, V.z}}
  MultMat33(InvertMat33(M), MultMat33(D,M))
#end

#macro CreateChromaMat()
  #local S = ReferenceWhiteXYZ;
  #local D = xyz2XYZ(ColSys2xyz(3));
  
  #switch (ChromaMatchFunction)
    #case (Scaling_ChromaMatch)
      #declare ChromaMatS = ApplyMatchMat(Scaling_MM, D,S);
      #declare ChromaMatD = InvertMat33(ChromaMatS);
    #break
    #case (VonKries_ChromaMatch)
      #declare ChromaMatS = ApplyMatchMat(VonKries_MM, D,S);
      #declare ChromaMatD = InvertMat33(ChromaMatS);
    #break
    #case (Bradford_ChromaMatch)
      #declare ChromaMatS = ApplyMatchMat(Bradford_MM, D,S);
      #declare ChromaMatD = InvertMat33(ChromaMatS);
    #break
    #else 
      #declare ChromaMatS = Scaling_MM; 
      #declare ChromaMatD = Scaling_MM; 
  #end
  
 #if (CIE_Debug)
  #debug "\n"
  #debug "Source White\n"
  #debug concat(str(S.x,2,5),", ",str(S.y,2,5),", ",str(S.z,2,5),"\n")
  #debug "Destination White\n"
  #debug concat(str(D.x,2,5),", ",str(D.y,2,5),", ",str(D.z,2,5),"\n")
  #debug "\n"
  #debug "Chroma Source Matrix\n"
  #debug concat(str(ChromaMatS[0][0],2,6),", ",str(ChromaMatS[0][1],2,6),", ",str(ChromaMatS[0][2],2,6),"\n")
  #debug concat(str(ChromaMatS[1][0],2,6),", ",str(ChromaMatS[1][1],2,6),", ",str(ChromaMatS[1][2],2,6),"\n")
  #debug concat(str(ChromaMatS[2][0],2,6),", ",str(ChromaMatS[2][1],2,6),", ",str(ChromaMatS[2][2],2,6),"\n")
  #debug "\n"
  #debug "Chroma Destination Matrix\n"
  #debug concat(str(ChromaMatD[0][0],2,6),", ",str(ChromaMatD[0][1],2,6),", ",str(ChromaMatD[0][2],2,6),"\n")
  #debug concat(str(ChromaMatD[1][0],2,6),", ",str(ChromaMatD[1][1],2,6),", ",str(ChromaMatD[1][2],2,6),"\n")
  #debug concat(str(ChromaMatD[2][0],2,6),", ",str(ChromaMatD[2][1],2,6),", ",str(ChromaMatD[2][2],2,6),"\n")
 #end 
#end

#macro ChromaMatchSource(XYZ)
  < ChromaMatS[0][0]*XYZ.x + ChromaMatS[0][1]*XYZ.y + ChromaMatS[0][2]*XYZ.z,
    ChromaMatS[1][0]*XYZ.x + ChromaMatS[1][1]*XYZ.y + ChromaMatS[1][2]*XYZ.z,
    ChromaMatS[2][0]*XYZ.x + ChromaMatS[2][1]*XYZ.y + ChromaMatS[2][2]*XYZ.z >
#end

#macro ChromaMatchDest(XYZ)
  < ChromaMatD[0][0]*XYZ.x + ChromaMatD[0][1]*XYZ.y + ChromaMatD[0][2]*XYZ.z,
    ChromaMatD[1][0]*XYZ.x + ChromaMatD[1][1]*XYZ.y + ChromaMatD[1][2]*XYZ.z,
    ChromaMatD[2][0]*XYZ.x + ChromaMatD[2][1]*XYZ.y + ChromaMatD[2][2]*XYZ.z >
#end

//=====================================================================
// internal matrix helper macros
//=====================================================================
#macro ApplyMat(Mat,XYZ)  // apply matrix transformation to 3d vector 
  < Mat[0][0]*XYZ.x + Mat[0][1]*XYZ.y + Mat[0][2]*XYZ.z,
    Mat[1][0]*XYZ.x + Mat[1][1]*XYZ.y + Mat[1][2]*XYZ.z,
    Mat[2][0]*XYZ.x + Mat[2][1]*XYZ.y + Mat[2][2]*XYZ.z >
#end

#macro MultMat33(M1,M2)  // multiply 2 3x3 matrices
  #local I=0; 
  #local Mat = array[3][3];
  #while (I<3) 
    #local J=0;  
    #while (J<3)   
      #declare Mat[I][J] = 0;
      #local K = 0;
      #while (K<3)
        #declare Mat[I][J] = Mat[I][J] + M1[I][K] * M2[K][J];
        #local K=K+1;
      #end  
      #local J=J+1;
    #end
    #local I=I+1;
  #end  
  Mat    
#end

#macro InvertMat33(M)  // invert a 3x3 matrix
  #local vX = <M[0][0], M[0][1], M[0][2]>;
  #local vY = <M[1][0], M[1][1], M[1][2]>;
  #local vZ = <M[2][0], M[2][1], M[2][2]>;
  #local Det = vdot(vX, vcross(vY, vZ));
  #local v0 =  vcross(vY, vZ)/Det;
  #local v1 = -vcross(vX, vZ)/Det;
  #local v2 =  vcross(vX, vY)/Det;
  array[3][3] {
    { v0.x, v1.x, v2.x }
    { v0.y, v1.y, v2.y }
    { v0.z, v1.z, v2.z }
  }
#end 


//=====================================================================
// Select the CIE observer.
// possible values are
//   CIE_1931  (2° observer)
//   CIE_1964  (10° observer)
//   CIE_1978  (2° observer modified by Judd/Vos)
//=====================================================================
#if (CIE_MultiObserver)
#macro CIE_Observer(S)
  #if (S!=CIE_Spectrum)         // need for recalculation
    #if ((S<0)|(S>2)) #error "CIE.inc: invalid observer selected" #end
    #declare CIE_Spectrum = S;
    #declare CMF_xyz = spline { cubic_spline
      375, <0,0,0>
      #local I=0;
      #while (I<90) 380+I*5, CMF_table[I][S] #local I=I+1; #end
      830, <0,0,0>
    }
    CIE_ReferenceWhite(ReferenceWhiteXY) // also ReferenceWhite has changed
  #end                                   // this forces the recalculation of
#end                                     // the chroma matrices
#end // if MultiObserver


//=====================================================================
// Change the reference whitepoint for the CIE observer.
// Usually there is no need to do this because all common applications 
// color management systems and spectrometer do work with D50 which is 
// also here the default setting.
// AND DO NOT CONFUSE THIS WITH THE COLOR SYSTEM WHITEPOINT!
//===================================================================== 

#macro CIE_ReferenceWhite(RW)
 
 #if (CIE_MultiObserver)
  #macro IlluID(I)
    I[0][0]+I[0][1]
  #end
 #else
  #macro IlluID(I)
    I[0]+I[1]
  #end
 #end

  #macro CreateRefWhite_D(K)
    #local M1=0; #local M2=0;
    DaylightM1M2((1.4388/1.4380)*K,M1,M2)
    #declare RefWhiteSP = spline{ linear_spline
      #local I=0;
      #while (I<77)
        380+I*5, (DS012[I][0]+M1*DS012[I][1]+M2*DS012[I][2])/100
        #local I=I+1;
      #end
    }
  #end
  
  #macro CreateRefWhite_A()
    #local D=PlanckBlackBody(560*1e-9, 2856);
    #declare RefWhiteSP = spline{ linear_spline
      #local I=0;
      #while (I<77)
        #local WL = 380+I*5;
        WL, PlanckBlackBody(WL*1e-9, 2856)/D
        #local I=I+1;
      #end
    }
  #end

  #switch (IlluID(RW))
    #case (IlluID(Illuminant_A))
      CreateRefWhite_A()
    #break
    #case (IlluID(Illuminant_B))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_B;
    #break
    #case (IlluID(Illuminant_C))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_C;
    #break
    #case (IlluID(Illuminant_D50))
      CreateRefWhite_D(5000)
    #break
    #case (IlluID(Illuminant_D55))
      CreateRefWhite_D(5500)
    #break
    #case (IlluID(Illuminant_D65))
      CreateRefWhite_D(6500)
    #break
    #case (IlluID(Illuminant_D75))
      CreateRefWhite_D(7500)
    #break
    #case (IlluID(Illuminant_E))
      #declare RefWhiteSP = spline{linear_spline 360, 1.0 830, 1.0}
    #break
    #case (IlluID(Illuminant_F1))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_F1;
    #break
    #case (IlluID(Illuminant_F2))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_F2;
    #break
    #case (IlluID(Illuminant_F3))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_F3;
    #break
    #case (IlluID(Illuminant_F4))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_F4;
    #break
    #case (IlluID(Illuminant_F5))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_F5;
    #break
    #case (IlluID(Illuminant_F6))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_F6;
    #break
    #case (IlluID(Illuminant_F7))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_F7;
    #break
    #case (IlluID(Illuminant_F8))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_F8;
    #break
    #case (IlluID(Illuminant_F9))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_F9;
    #break
    #case (IlluID(Illuminant_F10))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_F10;
    #break
    #case (IlluID(Illuminant_F11))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_F11;
    #break
    #case (IlluID(Illuminant_F12))
      #include "espd_cie_standard.inc"
      #declare RefWhiteSP = ES_Illuminant_F12;
    #break
    #else
      #error "CIE.inc: Invalid CIE reference white selected"
  #end
 #if (CIE_MultiObserver)
  #local X = RW[CIE_Spectrum][0];
  #local Y = RW[CIE_Spectrum][1];
 #else
  #local X = RW[0];
  #local Y = RW[1];
 #end  
  #declare ReferenceWhiteXY = RW;
  #declare ReferenceWhiteXYZ = <X/Y, 1.0, (1-X-Y)/Y>;
  CreateChromaMat()          // update the chromatic adaption matrices
  CIE_ColorSystem(xy_ColSys) // we need also a new rgb->xyz transformation matrix
#end


//=====================================================================
// select the chroma matching function
// possible values are
//   Scaling_ChromaMatch
//   VonKries_ChromaMatch
//   Bradford_ChromaMatch
//   or false/off/0 to switch it completely off 
//=====================================================================
#macro CIE_ChromaticAdaption(Match)
  #if (Match!=ChromaMatchFunction)
    #if (ChromaMatchFunction>3) 
      #error "CIE.inc: Invalid chromatic adaption selected" 
    #end
    #declare ChromaMatchFunction = Match;
    CreateChromaMat()
  #end
#end


//=====================================================================
// select the gamut mapping function for xyz->rgb conversion
//   0/off/false: no mapping
//   1: clipping
//   2: triangle intersection in rgb space
//=====================================================================
#macro CIE_GamutMapping(GM)
  #declare GamutMapFunction = GM;
#end


//=====================================================================
// Select a color/device system and create the required transformation
// matrix
//=====================================================================
#macro CIE_ColorSystem(CS)
// first check if a recalculation of the chroma matrices is required
  #local Idx=3+CIE_Spectrum;
  #if ((xy_ColSys[Idx][0]!=CS[Idx][0]) | (xy_ColSys[Idx][1]!=CS[Idx][1]))
    #declare xy_ColSys = CS; // set new color system
    CreateChromaMat()        // calculate chroma matrices
  #else
    #declare xy_ColSys = CS; // just set the active new color system
  #end
// copy into local vars
  #local XR = CS[0][0]; #local YR = CS[0][1]; #local ZR = 1-(XR+YR);
  #local XG = CS[1][0]; #local YG = CS[1][1]; #local ZG = 1-(XG+YG);
  #local XB = CS[2][0]; #local YB = CS[2][1]; #local ZB = 1-(XB+YB);
  #local XW = CS[Idx][0]; #local YW = CS[Idx][1];
// UV-W coordinates
  #local D = (XR-XB)*(YG-YB)-(YR-YB)*(XG-XB);
  #local U = (XW-XB)*(YG-YB)-(YW-YB)*(XG-XB);
  #local V = (XR-XB)*(YW-YB)-(YR-YB)*(XW-XB);
  #local U = U/D;   #local V = V/D;  #local W = 1-(U+V);
  #local U = U/YW;  #local V = V/YW; #local W = W/YW;
// rgb->xyz matrix
  #declare csMatXYZ[0][0]=U*XR; #declare csMatXYZ[0][1]=V*XG; #declare csMatXYZ[0][2]=W*XB; 
  #declare csMatXYZ[1][0]=U*YR; #declare csMatXYZ[1][1]=V*YG; #declare csMatXYZ[1][2]=W*YB; 
  #declare csMatXYZ[2][0]=U*ZR; #declare csMatXYZ[2][1]=V*ZG; #declare csMatXYZ[2][2]=W*ZB;
// the xyz-rgb matrix is the invert of the rgb-xyz matrix
  #declare csMatRGB = InvertMat33(csMatXYZ)

// the special matrix for rgb->rgb conversion
 #if (CIE_MultiObserver)
  #local XW = ReferenceWhiteXY[CIE_Spectrum][0];
  #local YW = ReferenceWhiteXY[CIE_Spectrum][1];
 #else
  #local XW = ReferenceWhiteXY[0];
  #local YW = ReferenceWhiteXY[1];
 #end
  #if ((XW!=CS[Idx][0]) & (YW!=CS[Idx][1]))
    #local U = (XW-XB)*(YG-YB)-(YW-YB)*(XG-XB);
    #local V = (XR-XB)*(YW-YB)-(YR-YB)*(XW-XB);
    #local U = U/D;   #local V = V/D;  #local W = 1-(U+V);
    #local U = U/YW;  #local V = V/YW; #local W = W/YW;
// new rgb->xyz matrix
    #declare csMatXYZ[0][0]=U*XR; #declare csMatXYZ[0][1]=V*XG; #declare csMatXYZ[0][2]=W*XB; 
    #declare csMatXYZ[1][0]=U*YR; #declare csMatXYZ[1][1]=V*YG; #declare csMatXYZ[1][2]=W*YB; 
    #declare csMatXYZ[2][0]=U*ZR; #declare csMatXYZ[2][1]=V*ZG; #declare csMatXYZ[2][2]=W*ZB;
  #end 

 #if (CIE_Debug)
  #debug "\n"
  #if (CIE_MultiObserver)
   #switch (CIE_Spectrum)
    #case (CIE_1931) #debug "CIE 1931 2 degree observer\n"  #break
    #case (CIE_1964) #debug "CIE 1964 10 degree observer\n" #break
    #case (CIE_1978) #debug "CIE 1931 2 degree observer (mod. Judd/Vos 1978)\n" #break
   #end
  #else
   #debug "CIE 1931 2 degree observer\n"
  #end
  #debug "\n"
  #debug "color system\n"
  #debug "       x       y\n"
  #debug concat("red   ",str(xy_ColSys[0][0],1,4),"  ",str(xy_ColSys[0][1],1,4),"\n")
  #debug concat("green ",str(xy_ColSys[1][0],1,4),"  ",str(xy_ColSys[1][1],1,4),"\n")
  #debug concat("blue  ",str(xy_ColSys[2][0],1,4),"  ",str(xy_ColSys[2][1],1,4),"\n")
  #debug concat("white ",str(xy_ColSys[3+CIE_Spectrum][0],1,5)," ",str(xy_ColSys[3+CIE_Spectrum][1],1,5))
  #debug "\n\n"
  #debug "rgb-xyz conversion matrix\n"
  #debug concat(str(csMatXYZ[0][0],2,5),", ",str(csMatXYZ[0][1],2,5),", ",str(csMatXYZ[0][2],2,5),"\n")
  #debug concat(str(csMatXYZ[1][0],2,5),", ",str(csMatXYZ[1][1],2,5),", ",str(csMatXYZ[1][2],2,5),"\n")
  #debug concat(str(csMatXYZ[2][0],2,5),", ",str(csMatXYZ[2][1],2,5),", ",str(csMatXYZ[2][2],2,5),"\n")
  #debug "\n"
  #debug "xyz-rgb conversion matrix\n"
  #debug concat(str(csMatRGB[0][0],2,5),", ",str(csMatRGB[0][1],2,5),", ",str(csMatRGB[0][2],2,5),"\n")
  #debug concat(str(csMatRGB[1][0],2,5),", ",str(csMatRGB[1][1],2,5),", ",str(csMatRGB[1][2],2,5),"\n")
  #debug concat(str(csMatRGB[2][0],2,5),", ",str(csMatRGB[2][1],2,5),", ",str(csMatRGB[2][2],2,5),"\n")
 #end // if debug
#end


//=====================================================================
// Select a color/device system but change also the whitepoint to any
// non default value.
// The color system with the new illuminant is also activated!
//=====================================================================
#if (CIE_MultiObserver)

#macro CIE_ColorSystemWhitepoint(CS,WP)
  #declare NewCS = CS;    // do not overwrite the original one
  #declare NewCS[3][0] = WP[0][0];
  #declare NewCS[3][1] = WP[0][1];
  #declare NewCS[4][0] = WP[1][0];
  #declare NewCS[4][1] = WP[1][1];
  #declare NewCS[5][0] = WP[2][0];
  #declare NewCS[5][1] = WP[2][1];
  CIE_ColorSystem(NewCS)    // activate
#end

#else

#macro CIE_ColorSystemWhitepoint(CS,WP)
  #declare NewCS = CS;    // do not overwrite the original one
  #declare NewCS[3][0] = WP[0];
  #declare NewCS[3][1] = WP[1];
  CIE_ColorSystem(NewCS)    // activate
#end

#end


//=========================================================                         
// xyz <-> XYZ conversion
//=========================================================

#macro xyz2XYZ(XYZ)
  <XYZ.x/XYZ.y, 1.0, XYZ.z/XYZ.y>
#end

#macro XYZ2xyz(XYZ)
  #local Y = 1/(XYZ.x+XYZ.z+1);
  <XYZ.x*Y, Y, XYZ.z*Y>
#end 


//=========================================================                         
// xyz <-> xyY conversion
//=========================================================

#macro xyz2xyY(XYZ)
  #local Sum = XYZ.x + XYZ.y + XYZ.z;
// note to myself:   
// if Sum=0 then x and y should correspondent
// to the color system whitepoint!
// Is this possible for (synthetical) colors?
// Maybe I should check this here and set x/y in the correct
// way to avoid the division by 0 ??? 
  <XYZ.x/Sum, XYZ.y/Sum, XYZ.y>
#end 

#macro xyY2xyz(xyY)  // xyY.z = Y  !!!
  #local Yy = xyY.z/xyY.y;
  <xyY.x*Yy,  xyY.z, (1.0-xyY.x-xyY.y)*Yy>
#end  

                         
//=========================================================
// xyz <-> rgb conversion  
// Note: 
// xyz->rgb transformation uses reference white
// rgb->xyz transformation uses color system whitepoint 
//=========================================================
// convert RGB to xyz by using the pre calculated transformation
// matrix csMatXYZ. 
#macro RGB2xyz(RGB)
  < csMatXYZ[0][0]*RGB.red + csMatXYZ[0][1]*RGB.green + csMatXYZ[0][2]*RGB.blue,
    csMatXYZ[1][0]*RGB.red + csMatXYZ[1][1]*RGB.green + csMatXYZ[1][2]*RGB.blue,
    csMatXYZ[2][0]*RGB.red + csMatXYZ[2][1]*RGB.green + csMatXYZ[2][2]*RGB.blue >
#end

// convert xyz to RGB by using the pre calculated transformation
// matrix csMatRGB.
// note that the resulting RGB components can be negative if they are
// not within the gamut as defined by the current color system
#macro xyz2RGB(XYZ)
  < csMatRGB[0][0]*XYZ.x + csMatRGB[0][1]*XYZ.y + csMatRGB[0][2]*XYZ.z,
    csMatRGB[1][0]*XYZ.x + csMatRGB[1][1]*XYZ.y + csMatRGB[1][2]*XYZ.z,
    csMatRGB[2][0]*XYZ.x + csMatRGB[2][1]*XYZ.y + csMatRGB[2][2]*XYZ.z >
#end


// gray mapping feature
// input from: 
//  0.0 - full rgb colors to 1.0 - grayscale 
#macro CIE_Saturation(S)
  #undef xyz2RGB
  #switch (S)
   #case (1)
     #macro xyz2RGB(XYZ)
      < csMatRGB[0][0]*XYZ.x + csMatRGB[0][1]*XYZ.y + csMatRGB[0][2]*XYZ.z,
        csMatRGB[1][0]*XYZ.x + csMatRGB[1][1]*XYZ.y + csMatRGB[1][2]*XYZ.z,
        csMatRGB[2][0]*XYZ.x + csMatRGB[2][1]*XYZ.y + csMatRGB[2][2]*XYZ.z >
     #end #break
   #case (0)
     #macro xyz2RGB(XYZ)
      < XYZ.y, XYZ.y, XYZ.y >
     #end #break
   #else 
     #if ((S>1)|(S<0)) #error "CIE.inc: saturation range is 0.0 .. 1.0" #end
     #declare CIE_SaturationValue = S;
     #macro xyz2RGB(XYZ)
      (< csMatRGB[0][0]*XYZ.x + csMatRGB[0][1]*XYZ.y + csMatRGB[0][2]*XYZ.z,
         csMatRGB[1][0]*XYZ.x + csMatRGB[1][1]*XYZ.y + csMatRGB[1][2]*XYZ.z,
         csMatRGB[2][0]*XYZ.x + csMatRGB[2][1]*XYZ.y + csMatRGB[2][2]*XYZ.z >*(CIE_SaturationValue)
         + < XYZ.y, XYZ.y, XYZ.y >*(1-CIE_SaturationValue) )
     #end    
  #end
#end

//=========================================================
// rgb->xyz->rgb conversion
// even if this sounds strange this macro can be used to
// transform a rgb color in the same way as reflective 
// spectra. The input value can be seen as reference rgb 
// and is illuminated by reference white and the output 
// rgb is modified in case the color system whitepoint
// is different.
//=========================================================
#macro ReferenceRGB(RGB)
  MapGamut(xyz2RGB(ChromaMatchSource(RGB2xyz(RGB))))
#end

#macro ReferenceRGB2RGB(RGB)
  xyz2RGB(ChromaMatchSource(RGB2xyz(RGB)))
#end


//=========================================================
// xyY <-> rgb conversion
//=========================================================
#macro RGB2xyY(RGB)
  xyz2xyY(RGB2xyz(RGB))
#end

#macro xyY2RGB(xyY)  
  xyz2RGB(xyY2xyz(xyY))
#end


//=========================================================
// Lab/Lch<->xyz<->rgb conversion 
// The luminance mapping is composed of two functions but 
// shows a discontinuity at the junction point.
// So there is a new (April 2003) draft from the CIE to use
// e = 216/24389 instead of 0.008856
// k = 24389/27  instead of 903.3
//=========================================================
#macro xyz2Lab(XYZ)
  #local E = 216/24389;
  #local J = 16/116;
  #local C = ChromaMatchDest(XYZ)/ReferenceWhiteXYZ;
  #if (C.y > E)
    #local F = pow(C.y,1/3);
    #local L = 116*F-16;
  #else
    #local F = 7.787*C.y+J;
    #local L = (24389/27)*C.y;
  #end
  < L,
    500*((C.x>E ? pow(C.x,1/3) : 7.787*C.x+J) - F),
    200*(F - (C.z>E ? pow(C.z,1/3) : 7.787*C.z+J)) >
#end

#macro Lab2xyz(Lab)
  #local E = 216/24389;
  #local J = 16/116;
  #local Y = (Lab.x +16) / 116;
  #local X = Lab.y/500 + Y;
  #local Z = Y - Lab.z/200;
  #local powX = pow(X,3);
  #local powY = pow(Y,3);
  #local powZ = pow(Z,3);
  #local X = (powX>E ? powX : (X-J) / 7.787);
  #local Y = (powY>E ? powY : (Y-J) / 7.787);
  #local Z = (powZ>E ? powZ : (Z-J) / 7.787);
  (ChromaMatchSource(<X,Y,Z>*ReferenceWhiteXYZ))
#end

#macro Lab2Lch(Lab)
  <Lab.x, vlength(Lab*<0,1,1>), (#if (abs(Lab.z)<1e-9) 0 #else atan2(Lab.z, Lab.y) #end) >
#end

#macro Lch2Lab(Lch)
  <Lch.x, cos(Lch.z)*Lch.y, sin(Lch.z)*Lch.y>
#end

#macro xyz2Lch(XYZ)
  Lab2Lch(xyz2Lab(XYZ))
#end

#macro Lch2xyz(Lch)
  Lab2xyz(Lch2Lab(Lch))
#end

#macro RGB2Lab(RGB)
  xyz2Lab(RGB2xyz(RGB))
#end

#macro Lab2RGB(Lab)
  xyz2RGB(Lab2xyz(Lab))
#end

#macro RGB2Lch(RGB)
  Lab2Lch(RGB2Lab(RGB))
#end

#macro Lch2RGB(Lch)
  Lab2RGB(Lch2Lab(Lch))
#end


//=====================================================================
// Gamut mapping is a kind of science for itself but here it's
// done in a quite simple way.
// In some cases method 2 gives better results but in some others it is
// even less accurate than just clipping the negative rgb color values.
//=====================================================================
#macro MapGamut(RGB_Input)
  #local RGB = RGB_Input;
  #if (min(RGB.red, RGB.green, RGB.blue) < 0.0)
   // Outside of the gamut, a little more work to do!
    #switch (GamutMapFunction)
    // 0 and any other: do not map
     #case (1)
      // clip the negative values
       #if (RGB.red<0)   #local RGB = <0, RGB.green, RGB.blue>; #end
       #if (RGB.green<0) #local RGB = <RGB.red, 0, RGB.blue>;   #end
       #if (RGB.blue<0)  #local RGB = <RGB.red, RGB.green, 0>;  #end
     #break  
     #case (2)
      // We need the RGB value for the whitepoint (and this point has to
      // be inside) and calculate the intersection with the edge of the 
      // maxwell triangle (the gamut as set by the active color  system) 
      // by a (fictive) line defined by our RGB value and the location 
      // of the whitepoint.  
       #local W = xyz2RGB(ColSys2xyz(3));  // RGB for the colsys whitepoint
       #if  ((RGB.red < RGB.green) & (RGB.red < RGB.blue))
         #local P = W.red / (W.red - RGB.red);
       #else
         #if  ((RGB.green < RGB.red) & (RGB.green < RGB.blue))
           #local P = W.green / (W.green - RGB.green);
         #else 
           #local P = W.blue / (W.blue - RGB.blue);
         #end  
       #end
      // finally calculate the desaturated gamut-constrained RGB weights 
       #local RGB = <W.red+P*(RGB.red-W.red), W.green+P*(RGB.green-W.green), W.blue+P*(RGB.blue-W.blue)>;
     #break
     #case (3)     
      // another try on gamut mapping by just desaturating the RGB color
       #local W = -min(0,RGB.red,RGB.green,RGB.blue);
       #local RGB = <RGB.red+W, RGB.green+W, RGB.blue+W>/(1+W);
     #break
    #end // switch
  #end  
  (RGB)
#end 

#macro MapGamutNorm(RGB_Input)
  #local RGB = MapGamut(RGB_Input);
  (RGB/max(RGB.red, RGB.green, RGB.blue))
#end


//=====================================================================
// Wavelength conversion 
// Wavelength to xyz color space aproximation by using the cubic interpolated
// CMF as lookup table where the "border" of the *visible* color space is 
// defined in xy coordinates.
//=====================================================================
#macro Wavelength2xyz(W)
  CMF_xyz(min(825,max(380,W)))
#end 

#macro Wavelength(W)
  #if ((W<380)|(W>825)) // Outside of the spectral data range
    #local RGB = <0,0,0>;
  #else
    #local RGB = MapGamutNorm(xyz2RGB(CMF_xyz(W)));
  #end                                                
  (RGB)
#end

#macro Wavelength2RGB(W)
  #if ((W<380)|(W>825))
    #local RGB = <0,0,0>;
  #else
    #local RGB = xyz2RGB(CMF_xyz(W));
  #end
  (RGB)
#end

#macro Wavelength2Lab(W)
  #if ((W<380)|(W>825))
    #local Lab = <0,0,0>;
  #else
    #local Lab = xyz2Lab(CMF_xyz(W));
  #end
  (Lab)
#end

#macro Wavelength2Lch(W)
  #if ((W<380)|(W>825))
    #local Lch = <0,0,0>;
  #else
    #local Lch = xyz2Lch(CMF_xyz(W));
  #end
  (Lch)
#end


//=====================================================================
// Reflective spectrum calculation
//=====================================================================
// converts a given reflectance spectrum to xyz
// this is the basic function and should not be
// used directly because no chroma match is applied
#macro ReflectanceXYZ(S)
  // calculations based on spreadsheet available at:
  // http://www.digitalcolour.org/understanding/ 
  #local W=380;
  #local XYZ=<0,0,0>;
  #local CMy=0;
  #while (W<=760)
    #local CM = CMF_xyz(W);  
    #local RW  = RefWhiteSP(W).y;
    #local CMy = CMy + CM.y*RW;
    #local XYZ = XYZ + CM*S(W).y*RW;
    #local W=W+CIE_IntegralStep;
  #end
  (XYZ/CMy)
#end

// reflectance spectrum to xyz
#macro Reflective2xyz(S)
  ChromaMatchSource(ReflectanceXYZ(S))
#end

// converts to gamut mapped RGB
#macro ReflectiveSpectrum(S)
  MapGamut(xyz2RGB(ChromaMatchSource(ReflectanceXYZ(S))))
#end 

// converts to RGB
#macro Reflective2RGB(S)
  xyz2RGB(ChromaMatchSource(ReflectanceXYZ(S)))
#end 

// no chroma mapping for Lab and Lch !
// converts to Lab
#macro Reflective2Lab(S)
  #local C = ReflectanceXYZ(S)/ReferenceWhiteXYZ;   
  #local E = 216/24389;
  #local J = 16/116;
  #if (C.y > E)
    #local F = pow(C.y,1/3);
    #local L = 116*F-16; 
  #else
    #local F = 7.787*C.y+J;
    #local L = (24389/27)*C.y;   
  #end     
  < L,
    500*((C.x>E ? pow(C.x,1/3) : 7.787*C.x+J) - F),  // a
    200*(F - (C.z>E ? pow(C.z,1/3) : 7.787*C.z+J)) > // b
#end 

// converts to Lch
#macro Reflective2Lch(S)
  Lab2Lch(Reflective2Lab(S))
#end        


//=====================================================================
// Emissive spectrum calculation
//=====================================================================
#macro Emissive2xyz(S)
  #local W=380;
  #local XYZ=<0,0,0>;
  #while (W<=760)   // integrate over spectrum data
    #local XYZ = XYZ + CMF_xyz(W)*S(W).y;
    #local W=W+CIE_IntegralStep;
  #end
  (XYZ/(XYZ.x+XYZ.y+XYZ.z))
#end

// converts a given spectrum to gamut mapped RGB
#macro EmissiveSpectrum(S)
  MapGamutNorm(xyz2RGB(Emissive2xyz(S)))
#end 

// converts a given spectrum to RGB
#macro Emissive2RGB(S)
  xyz2RGB(Emissive2xyz(S))
#end 

// converts a given spectrum to Lab
#macro Emissive2Lab(S)
  xyz2Lab(Emissive2xyz(S))
#end 

// converts a given spectrum to Lch
#macro Emissive2Lch(S)
  xyz2Lch(Emissive2xyz(S))
#end


//=====================================================================
// Line spectrum calculation
//=====================================================================
#macro LineSpectrum2xyz(LS, Ion)
  #local N = dimension_size(LS,1);
  #local K = 0; 
  #local CMy = 0;
  #local XYZ = <0,0,0>;
  #while (K < N)          
    #if (LS[K][2]<=(Ion)) // check for ionization
      #local CMy = CMy+CMF_xyz(LS[K][0]).y;
      #local XYZ = XYZ + CMF_xyz(LS[K][0])*(LS[K][1]/1000);
    #end  
    #local K=K+1;
  #end    
  (XYZ/CMy)
#end  

#macro LineSpectrum(LS, Ion)
  MapGamutNorm(xyz2RGB(LineSpectrum2xyz(LS,Ion)))
#end

#macro LineSpectrum2RGB(LS, Ion)
  xyz2RGB(LineSpectrum2xyz(LS,Ion))
#end

#macro LineSpectrum2Lab(LS, Ion)
  xyz2Lab(LineSpectrum2xyz(LS,Ion))
#end

#macro LineSpectrum2Lch(LS, Ion)
  xyz2Lch(LineSpectrum2xyz(LS,Ion))
#end           

//=====================================================================
// Combine two or more line spectra
//=====================================================================
/*   
   EXAMPLE: 
   
   #declare Lspd_Mix = array[4] {
     LS_Neon,
     LS_Neon,
     LS_Argon,
     LS_Xenon
   } 
   #declare F_Mix = array[4][2] {
     {0.80, 2}, // 80% Neon   two electrons lost
     {0.15, 1}, // 15% Neon   one electron lost
     {0.04, 2}, //  4% Argon
     {0.01, 2}  //  1% Xenon
   }
   #declare GasColor = LineSpectrumMix(Lspd_Mix, F_Mix);
*/

#macro LineSpectrumMix2xyz(LS,LF)
  #local NE = dimension_size(LS,1);
  #local E = 0; 
  #local XYZ = <0,0,0>;
  #while (E<NE)
    #local XYZ = XYZ+LineSpectrum2xyz(LS[E],LF[E][1])*LF[E][0];
    #local E=E+1;
  #end   
  (XYZ)
#end

#macro LineSpectrumMix(LS,LF)
  MapGamutNorm(xyz2RGB(LineSpectrumMix2xyz(LS,LF)))
#end

#macro LineSpectrumMix2RGB(LS,LF)
  xyz2RGB(LineSpectrumMix2xyz(LS,LF))
#end

#macro LineSpectrumMix2Lab(LS,LF)
  xyz2Lab(LineSpectrumMix2xyz(LS,LF))
#end

#macro LineSpectrumMix2Lch(LS,LF)
  xyz2Lch(LineSpectrumMix2xyz(LS,LF))
#end


//=====================================================================
// Kelvin conversion using Planck's blackbody radiation.
// This calculates the "pressure" of the light by a given wavelength.
//
// 1.4388e-2   = (h*c)/k    (was former 1.4380 but has changed due the
//                           redefinition of the Boltzmann constant)
// 3.74183e-16 = 2*pi*c*(h*c) 
//
// where
//   c = 2.997925e8    - Lightspeed
//   h = 6.626176e-34  - Planck's Constant 
//   k = 1.380650e-23  - Boltzmann Constant
//
// Wlm = wavelength in meter
// K   = temperature of blackbody in kelvin
//
// a google search for "Stefan Boltzmann law", "Wien's law" and 
// "Planck Blackbody" should give you detailed information.
//=====================================================================
#declare PlanckBlackBody = function(Wlm, K) {
  ( (3.74183e-16 * pow(Wlm,-5.0)) / (exp(1.4388e-2 / (Wlm*K))-1.0) )
}; 

#macro Blackbody2xyz(K)
  #local W=380;
  #local XYZ=<0,0,0>;
  #while (W<=760)
    #local XYZ = XYZ + CMF_xyz(W)*PlanckBlackBody(W*1e-9, K);
    #local W=W+5;
  #end
  (XYZ/(XYZ.x+XYZ.y+XYZ.z))
#end      

#macro Blackbody(K)
  MapGamutNorm(xyz2RGB(Blackbody2xyz(K)))
#end

#macro Blackbody2RGB(K)
  xyz2RGB(Blackbody2xyz(K))
#end

#macro Blackbody2Lab(K)
  xyz2Lab(Blackbody2xyz(K))
#end

#macro Blackbody2Lch(K)
  xyz2Lch(Blackbody2xyz(K))
#end 

//=====================================================================
// Illuminant D daylight conversion using S0, S1, S2 (implemented as
// DS012 lookup table) to calculate (ideal) daylight spectra and then
// integrate in the same way as any other emissive spectrum.
//=====================================================================
#macro DaylightM1M2(K,M1,M2)
  #if ((K<4000)|(K>25000))
    #error "CIE.inc: K is limited to the range of 4000-25000 Kelvin (use Blackbody instead)"
  #end
  #if (K<=7000)
    #local X = -4.6070e9/pow(K,3) + 2.9678e6/pow(K,2) + 0.09911e3/K + 0.244063;
  #else
    #local X = -2.0064e9/pow(K,3) + 1.9018e6/pow(K,2) + 0.24748e3/K + 0.237040;
  #end
  #local Y = -3.0*(X*X) + 2.87*X - 0.275;
  #local I = (0.0241 + 0.2562*X - 0.7341*Y);
  #local M1 = (-1.3515 - 1.7703*X + 5.9114*Y) / I;
  #local M2 = ( 0.03   -31.4424*X +30.0717*Y) / I;
#end

#macro Daylight2xyz(K)
  #local M1=0; #local M2=0;   
  DaylightM1M2(K,M1,M2)
  #local I=0;
  #local W=380;
  #local XYZ = <0,0,0>;
  #while (W<=760)
    #local XYZ = XYZ + CMF_xyz(W)*(DS012[I][0]+M1*DS012[I][1]+M2*DS012[I][2]);
    #local W=W+5;
    #local I=I+1;
  #end
  (XYZ/(XYZ.x+XYZ.y+XYZ.z))
#end

#macro Daylight(K)
  MapGamutNorm(xyz2RGB(Daylight2xyz(K)))
#end

#macro Daylight2RGB(K)
  xyz2RGB(Daylight2xyz(K))
#end

#macro Daylight2Lab(K)
  xyz2Lab(Daylight2xyz(K))
#end

#macro Daylight2Lch(K)
  xyz2Lch(Daylight2xyz(K))
#end


//=====================================================================
// Custom whitepoint calculation based on blackbody temperature
//=====================================================================
#if (CIE_MultiObserver)

#macro Blackbody2Whitepoint(K)
//calculate for all observer
  #local W=380; 
  #local I=0;
  #local XYZ31=<0,0,0>;
  #local XYZ64=<0,0,0>;
  #local XYZ78=<0,0,0>;
  #while (W<=760)
    #local BB = PlanckBlackBody(W*1e-9, K);
    #local XYZ31 = XYZ31 + (CMF_table[I][0]*BB);
    #local XYZ64 = XYZ64 + (CMF_table[I][1]*BB);
    #local XYZ78 = XYZ78 + (CMF_table[I][2]*BB);
    #local I=I+1;
    #local W=W+5;
  #end
  #local XYZ31 = xyz2xyY(XYZ31);
  #local XYZ64 = xyz2xyY(XYZ64);
  #local XYZ78 = xyz2xyY(XYZ78);
// the result
  array[3][2] {
    {XYZ31.x, XYZ31.y}
    {XYZ64.x, XYZ64.y}
    {XYZ78.x, XYZ78.y}
  }
#end

#else //no MultiObserver

#macro Blackbody2Whitepoint(K)
  #local W=380; 
  #local XYZ=<0,0,0>;
  #while (W<=760)
    #local XYZ = XYZ + (CMF_xyz(W)*PlanckBlackBody(W*1e-9, K));
    #local W=W+5;
  #end
  #local XYZ = xyz2xyY(XYZ);
// the result
  array[2] {XYZ.x, XYZ.y}
#end

#end //if MultiObserver


//=====================================================================
// Custom whitepoint calculation based on D illuminant recommendation
//=====================================================================

#if (CIE_MultiObserver)

#macro Daylight2Whitepoint(K)  
  #local M1=0; #local M2=0; DaylightM1M2(K,M1,M2)
  #local I=0;
  #local XYZ31=<0,0,0>;
  #local XYZ64=<0,0,0>;
  #local XYZ78=<0,0,0>;
  #while (I<90)  
    #local S = DS012[I][0] + M1*DS012[I][1] + M2*DS012[I][2];
    #local XYZ31 = XYZ31 + (CMF_table[I][0]*S);
    #local XYZ64 = XYZ64 + (CMF_table[I][1]*S);
    #local XYZ78 = XYZ78 + (CMF_table[I][2]*S);
    #local I=I+1;
  #end 
  #local XYZ31 = xyz2xyY(XYZ31);
  #local XYZ64 = xyz2xyY(XYZ64);
  #local XYZ78 = xyz2xyY(XYZ78);
// the result
  array[3][2] {
    {XYZ31.x, XYZ31.y}
    {XYZ64.x, XYZ64.y}
    {XYZ78.x, XYZ78.y}
  }
#end

#else //no MultiObserver

#macro Daylight2Whitepoint(K)  
  #local M1=0; #local M2=0; DaylightM1M2(K,M1,M2)
  #local I=0;
  #local XYZ=<0,0,0>; 
  #while (I<90)  
    #local XYZ = XYZ + (CMF_xyz(380+I*5)*(DS012[I][0] + M1*DS012[I][1] + M2*DS012[I][2]));
    #local I=I+1;
  #end 
  #local XYZ = xyz2xyY(XYZ);
// the result
  array[2] {XYZ.x, XYZ.y}
#end

#end //if MultiObserver 


//=====================================================================
// This one is only added to precalculate the whitepoints from the
// spectra in espd_cie_standard.inc 
//=====================================================================

#if (CIE_MultiObserver)

#macro Illuminant2Whitepoint(ES)
  #local I=0;
  #local XYZ31=<0,0,0>;
  #local XYZ64=<0,0,0>;
  #local XYZ78=<0,0,0>;
  #while (I<77)  
    #local S = ES(380+I*5).y;  
    #local XYZ31 = XYZ31 + CMF_table[I][0]*S;
    #local XYZ64 = XYZ64 + CMF_table[I][1]*S;
    #local XYZ78 = XYZ78 + CMF_table[I][2]*S;
    #local I=I+1;
  #end     
  #local XYZ31 = xyz2xyY(XYZ31);
  #local XYZ64 = xyz2xyY(XYZ64);
  #local XYZ78 = xyz2xyY(XYZ78);
  array[3][2] {
    {XYZ31.x, XYZ31.y}
    {XYZ64.x, XYZ64.y}
    {XYZ78.x, XYZ78.y}
  }
#end

#else

#macro Illuminant2Whitepoint(ES)
  #local I=0;
  #local XYZ=<0,0,0>;
  #while (I<77)  
    #local XYZ = XYZ + CMF_xyz(380+I*5)*ES(380+I*5).y;
    #local I=I+1;
  #end     
  #local XYZ = xyz2xyY(XYZ);
  array[2] {XYZ.x, XYZ.y}
#end

#end   


//=====================================================================
// Find the correlated blackbody color temperature for a given point 
// in CIE XYZ space.
// This implementation uses the method of A.R.Robertson as described in
// "Color Science: Concepts and Methods, Quantitative Data and Formulae"
// by G. Wyszecki and W.S. Stiles. 
//=====================================================================

#macro XYZ2K(XYZ)

  #local LowerLimit = 1.0e-20; // prevent possible division by 0

  // reciprocal temperatures (K) 
  #local ReciprocalTmp = array[31] { 
      LowerLimit,
         10.0e-6,
         20.0e-6,
         30.0e-6,
         40.0e-6,
         50.0e-6,
         60.0e-6,
         70.0e-6,
         80.0e-6,
         90.0e-6,
        100.0e-6,
        125.0e-6,
        150.0e-6,
        175.0e-6,
        200.0e-6,
        225.0e-6,
        250.0e-6,
        275.0e-6,
        300.0e-6,
        325.0e-6,
        350.0e-6,
        375.0e-6,
        400.0e-6,
        425.0e-6,
        450.0e-6,
        475.0e-6,
        500.0e-6,
        525.0e-6,
        550.0e-6,
        575.0e-6,
        600.0e-6
  };

  // correlated UV coordinates
  #local UVT = array[31][3] {  
        {0.18006, 0.26352, -0.24341},
        {0.18066, 0.26589, -0.25479},
        {0.18133, 0.26846, -0.26876},
        {0.18208, 0.27119, -0.28539},
        {0.18293, 0.27407, -0.30470},
        {0.18388, 0.27709, -0.32675},
        {0.18494, 0.28021, -0.35156},
        {0.18611, 0.28342, -0.37915},
        {0.18740, 0.28668, -0.40955},
        {0.18880, 0.28997, -0.44278},
        {0.19032, 0.29326, -0.47888},
        {0.19462, 0.30141, -0.58204},
        {0.19962, 0.30921, -0.70471},
        {0.20525, 0.31647, -0.84901},
        {0.21142, 0.32312, -1.0182},
        {0.21807, 0.32909, -1.2168},
        {0.22511, 0.33439, -1.4512},
        {0.23247, 0.33904, -1.7298},
        {0.24010, 0.34308, -2.0637},
        {0.24792, 0.34655, -2.4681},
        {0.25591, 0.34951, -2.9641},
        {0.26400, 0.35200, -3.5814},
        {0.27218, 0.35407, -4.3633},
        {0.28039, 0.35577, -5.3762},
        {0.28863, 0.35714, -6.7262},
        {0.29685, 0.35823, -8.5955},
        {0.30505, 0.35907, -11.324},
        {0.31320, 0.35968, -15.628},
        {0.32129, 0.36011, -23.325},
        {0.32931, 0.36038, -40.770},
        {0.33724, 0.36051, -116.45}
  };

  // just to make the indexing into the UVT lookup table more readable
  #local U_lut = 0;
  #local V_lut = 1;
  #local T_lut = 2;
  
  #if ( (XYZ.x < LowerLimit) & (XYZ.y < LowerLimit) & (XYZ.z < LowerLimit) )
    #error("CIE.inc: Error in macro XYZ2K() - invalid input with possible devide by zero error as result")
  #end  
  
  // transpose CIE XYZ to UV coordinates
  #local U = (4.0 * XYZ.x) / (XYZ.x + 15.0 * XYZ.y + 3.0 * XYZ.z);
  #local V = (6.0 * XYZ.y) / (XYZ.x + 15.0 * XYZ.y + 3.0 * XYZ.z);

  #local DM = 0.0;

  #local I = 0;
  #while (I < 31)
     #local DI = (V - UVT[I][V_lut]) - UVT[I][T_lut] * (U - UVT[I][U_lut]);
     #if ((I > 0) & (((DI < 0.0) & (DM >= 0.0)) | ((DI >= 0.0) & (DM < 0.0))))
       #break;  // found UV bounding lines  
     #else  
       #local DM = DI;
       #local I = I+1;
       #if (I = 31) 
          #error("CIE.inc: Error in macro XYZ2K() - invalid input, color temperature would be less than minimum of 1666.7° or way too blue")
       #end   
     #end  
  #end     

  #local DI = DI / sqrt(1.0 + UVT[I][T_lut]   * UVT[I][T_lut]);
  #local DM = DM / sqrt(1.0 + UVT[I-1][T_lut] * UVT[I-1][T_lut]);
  #local K = DM / (DM - DI);

  // linear interpolation from the found bounding lines to the blackbody curve
  #local K = 1.0 / ((ReciprocalTmp[I] - ReciprocalTmp[I-1]) * K + ReciprocalTmp[I-1]);
  
  // finally return the color temperature
  (K) 

#end



//=====================================================================
// Mixing of 2 reflective SPD's 
// S1: first SPD
// S2: second SPD
// F1,F2: factor to mix them
//=====================================================================
#macro ReflectiveMixXYZ(S1,S2,F1,F2)
  #local W=380;
  #local XYZ=<0,0,0>;
  #local CMy=0; 
  #local P1 = F1/(F1+F2);
  #local P2 = 1.0-P1;
  #while (W<=760)
    #local CM = CMF_xyz(W);  
    #local RW = RefWhiteSP(W).y;
    #local CMy = CMy + CM.y*RW;
    #local XYZ = XYZ + (S1(W).y*P1+S2(W).y*P2)*CM*RW;
    #local W=W+CIE_IntegralStep;
  #end               
  (XYZ/CMy)
#end      

// reflectance mix spectrum to xyz
#macro ReflectiveMix2xyz(S1,S2,F1,F2)
  ChromaMatchSource(ReflectiveMixXYZ(S1,S2,F1,F2))
#end

// converts to gamut mapped RGB
#macro ReflectiveMix(S1,S2,F1,F2)
  MapGamut(xyz2RGB(ChromaMatchSource(ReflectiveMixXYZ(S1,S2,F1,F2))))
#end 

// converts to RGB
#macro ReflectiveMix2RGB(S1,S2,F1,F2)
  xyz2RGB(ChromaMatchSource(ReflectiveMixXYZ(S1,S2,F1,F2)))
#end 

// converts to Lab
#macro ReflectiveMix2Lab(S1,S2,F1,F2)
  #local C = ReflectiveMixXYZ(S1,S2,F1,F2)/ReferenceWhiteXYZ;   
  #if (C.y > 0.008856)
    #local F = pow(C.y,1/3);
    #local L = 116*F-16; 
  #else
    #local F = 7.787*C.y+16/116;
    #local L = 903.3*C.y;   
  #end  
  < L,
    500*((C.x>0.008856 ? pow(C.x,1/3) : 7.787*C.x+16/116) - F),
    200*(F - (C.z>0.008856 ? pow(C.z,1/3) : 7.787*C.z+16/116)) >
#end 

// converts to Lch  
#macro ReflectiveMix2Lch(S1,S2,F1,F2)
  Lab2Lch(ReflectiveMix2Lab(S1,S2,F1,F2))
#end 


//=====================================================================
// public utilities and helper macros
//=====================================================================

// normalization    
#macro NormRGB(RGB)
  (RGB/max(RGB.red,RGB.green,RGB.blue))
#end

// Kelvin correction due the revised Boltzmann constant
#declare BoltzmannCorrect = function(K) {(1.4388/1.4380)*K}

//=====================================================================
// Free allocated memory
// 
//=====================================================================

#macro CIE_ReleaseMemory() 
  #if (CIE_MultiObserver)
    #undef CMF_table
    #undef CIE_1931
    #undef CIE_1964
    #undef CIE_1978
  #end
 
  #undef Illuminant_A
  #undef Illuminant_B
  #undef Illuminant_C
  #undef Illuminant_D50
  #undef Illuminant_D55
  #undef Illuminant_D65
  #undef Illuminant_D75
  #undef Illuminant_E
  #undef Illuminant_F1
  #undef Illuminant_F2
  #undef Illuminant_F3
  #undef Illuminant_F4
  #undef Illuminant_F5
  #undef Illuminant_F6
  #undef Illuminant_F7
  #undef Illuminant_F8
  #undef Illuminant_F9
  #undef Illuminant_F10
  #undef Illuminant_F11
  #undef Illuminant_F12
  
  #undef sRGB_ColSys
  #undef CIE_ColSys
  #undef ITU_ColSys
  #undef NTSC_ColSys
  #undef SMPTE_ColSys
  #undef Adobe_ColSys
  #undef Match_ColSys
  #undef Apple_ColSys
  #undef Beta_ColSys
  #undef Dell_ColSys
  #undef ShortPersistence_ColSys
  #undef LongPersistence_ColSys
  #undef HOT_ColSys
    
  #undef ColSys2xyz
  #undef ApplyMatchMat
  #undef CreateChromaMat
  #undef ApplyMat
  #undef MultMat33
  #undef InvertMat33
  #undef CIE_ReferenceWhite
  #undef CIE_ChromaticAdaption
  #undef CIE_GamutMapping
  #undef CIE_ColorSystem
  #undef CIE_ColorSystemWhitepoint
  #undef Blackbody2Whitepoint
  #undef Daylight2Whitepoint
  #undef Illuminant2Whitepoint
  
  #undef CMF_xyz
  #undef DS012
  #undef RefWhiteSP
  #undef Scaling_MM
  #undef VonKries_MM
  #undef Bradford_MM
  #undef Scaling_ChromaMatch
  #undef VonKries_ChromaMatch
  #undef Bradford_ChromaMatch
  #undef CIE_Spectrum
  #undef xy_ColSys
  #undef ReferenceWhiteXY
  #undef ReferenceWhiteXYZ
  #undef ChromaMatchFunction
  #undef GamutMapFunction
  #undef csMatXYZ
  #undef csMatRGB
  #undef ChromaMatS
  #undef ChromaMatD
  #undef CIE_IntegralStep
  #undef CIE_MultiObserver
  #undef ChromaMatchSource
  #undef ChromaMatchDest
  
  #undef BoltzmannCorrect
  #undef PlanckBlackBody
    
  #undef xyz2XYZ
  #undef XYZ2xyz
  #undef xyz2xyY
  #undef xyY2xyz
  #undef RGB2xyz
  #undef xyz2RGB
  #undef ReferenceRGB
  #undef ReferenceRGB2RGB
  #undef RGB2xyY
  #undef xyY2RGB
  #undef xyz2Lab
  #undef Lab2xyz
  #undef Lab2Lch
  #undef Lch2Lab
  #undef xyz2Lch
  #undef Lch2xyz
  #undef RGB2Lab
  #undef Lab2RGB
  #undef RGB2Lch
  #undef Lch2RGB
  #undef MapGamut
  #undef MapGamutNorm
  #undef Wavelength2xyz
  #undef Wavelength
  #undef Wavelength2RGB
  #undef Wavelength2Lab
  #undef Wavelength2Lch
  #undef ReflectanceXYZ
  #undef Reflective2xyz
  #undef ReflectiveSpectrum
  #undef Reflective2RGB
  #undef Reflective2Lab
  #undef Reflective2Lch
  #undef Emissive2xyz
  #undef EmissiveSpectrum
  #undef Emissive2RGB
  #undef Emissive2Lab
  #undef Emissive2Lch
  #undef LineSpectrum2xyz
  #undef LineSpectrum
  #undef LineSpectrum2RGB
  #undef LineSpectrum2Lab
  #undef LineSpectrum2Lch
  #undef LineSpectrumMix2xyz
  #undef LineSpectrumMix
  #undef LineSpectrumMix2RGB
  #undef LineSpectrumMix2Lab
  #undef LineSpectrumMix2Lch
  #undef Blackbody2xyz
  #undef Blackbody
  #undef Blackbody2RGB
  #undef Blackbody2Lab
  #undef Blackbody2Lch
  #undef DaylightM1M2
  #undef Daylight2xyz
  #undef Daylight
  #undef Daylight2RGB
  #undef Daylight2Lab
  #undef Daylight2Lch
  #undef ReflectiveMixXYZ
  #undef ReflectiveMix2xyz
  #undef ReflectiveMix
  #undef ReflectiveMix2RGB
  #undef ReflectiveMix2Lab
  #undef ReflectiveMix2Lch
  #undef GammaRGB
  #undef GammaCorrectXYZ2RGB
  #undef GammaXYZ2RGB
  #undef GammaCorrectRGB2XYZ
  #undef GammaRGB2XYZ
  #undef NormRGB
    
  #undef CIE_Observer_String
  #undef CIE_ColorSystem_String
  #undef CIE_Whitepoint_String
  #undef CIE_ReferenceWhite_String
  #undef CIE_ChromaticAdaption_String

  #ifdef (espd_cie_standard_Inc_Temp)     espd_cie_standard_Release()     #end
  #ifdef (espd_lightsys_Inc_Temp)         espd_lightsys_Release()         #end
  #ifdef (espd_sun_Inc_Temp)              espd_sun_Release()              #end
  #ifdef (lspd_elements_Inc_Temp)         lspd_elements_Release()         #end
  #ifdef (rspd_jvp_Inc_Temp)              rspd_jvp_Release()              #end
  #ifdef (rspd_aster_Inc_Temp)            rspd_aster_Release()            #end
  #ifdef (rspd_lunar_Inc_Temp)            rspd_lunar_Release()            #end
  #ifdef (rspd_pantone_coated_Inc_Temp)   rspd_pantone_coated_Release()   #end
  #ifdef (rspd_pantone_uncoated_Inc_Temp) rspd_pantone_uncoated_Release() #end
  #ifdef (rspd_pantone_matte_Inc_Temp)    rspd_pantone_matte_Release()    #end
#end

//=====================================================================
// Some utilities for text information 
//=====================================================================

#macro CIE_csID(CS)
  (CS[0][0]+CS[1][0]+CS[2][0]+CS[2][1])
#end 

#if (CIE_MultiObserver)
 #macro CIE_iID(I)
  I[CIE_Spectrum][0]+I[CIE_Spectrum][1]
 #end
#else
 #macro CIE_iID(I)
  I[0]+I[1]
 #end 
#end

#macro CIE_Observer_String(S)
  concat(S,
 #if (CIE_MultiObserver)
  #switch (CIE_Spectrum)
    #case(CIE_1931) "CIE 1931 2 degree" #break
    #case(CIE_1964) "CIE 1964 10 degree" #break
    #case(CIE_1978) "CIE 1978 2 degree (Judd/Vos)" #break
  #end 
 #else                    
   "CIE 1931 2 degree"
 #end
  )                    
#end  

#macro CIE_ColorSystem_String(S)
  concat(S,
  #switch (CIE_csID(xy_ColSys))  
    #case(CIE_csID(sRGB_ColSys))            "sRGB"                       #break
    #case(CIE_csID(CIE_ColSys))             "CIE RGB"                    #break
    #case(CIE_csID(NTSC_ColSys))            "NTSC"                       #break
    #case(CIE_csID(ITU_ColSys))             "ITU (Pal Secam)"            #break
    #case(CIE_csID(SMPTE_ColSys))           "SMPTE"                      #break
    #case(CIE_csID(Match_ColSys))           "Color Match RGB"            #break
    #case(CIE_csID(Adobe_ColSys))           "Adobe RGB"                  #break
    #case(CIE_csID(Apple_ColSys))           "Apple RGB"                  #break 
    #case(CIE_csID(Beta_ColSys))            "Beta RGB"                   #break 
    #case(CIE_csID(Dell_ColSys))            "Dell Phosphors"             #break 
    #case(CIE_csID(ShortPersistence_ColSys))"Short Persistence Phosphors"#break
    #case(CIE_csID(LongPersistence_ColSys)) "Long Persistence Phosphors" #break
    #case(CIE_csID(HOT_ColSys))             "Hydrogen-Oxygen-Thermal"    #break 
    #else  "unknown"
  #end
  )
#end

#macro CIE_Whitepoint_String(S)
  concat(S,  
  #local X=xy_ColSys[3+CIE_Spectrum][0];
  #local Y=xy_ColSys[3+CIE_Spectrum][1];
  #switch (X+Y)
    #case(CIE_iID(Illuminant_A))   "Illuminant_A",   #break
    #case(CIE_iID(Illuminant_B))   "Illuminant_B",   #break
    #case(CIE_iID(Illuminant_C))   "Illuminant_C",   #break
    #case(CIE_iID(Illuminant_D50)) "Illuminant_D50", #break
    #case(CIE_iID(Illuminant_D55)) "Illuminant_D55", #break
    #case(CIE_iID(Illuminant_D65)) "Illuminant_D65", #break
    #case(CIE_iID(Illuminant_D75)) "Illuminant_D75", #break
    #case(CIE_iID(Illuminant_E))   "Illuminant_E",   #break
    #case(CIE_iID(Illuminant_F1))  "Illuminant_F1",  #break
    #case(CIE_iID(Illuminant_F2))  "Illuminant_F2",  #break
    #case(CIE_iID(Illuminant_F3))  "Illuminant_F3",  #break
    #case(CIE_iID(Illuminant_F4))  "Illuminant_F4",  #break
    #case(CIE_iID(Illuminant_F5))  "Illuminant_F5",  #break
    #case(CIE_iID(Illuminant_F6))  "Illuminant_F6",  #break
    #case(CIE_iID(Illuminant_F7))  "Illuminant_F7",  #break
    #case(CIE_iID(Illuminant_F8))  "Illuminant_F8",  #break
    #case(CIE_iID(Illuminant_F9))  "Illuminant_F9",  #break
    #case(CIE_iID(Illuminant_F10)) "Illuminant_F10", #break
    #case(CIE_iID(Illuminant_F11)) "Illuminant_F11", #break
    #case(CIE_iID(Illuminant_F12)) "Illuminant_F12", #break
    #else "user defined",
  #end
  "  (",str(X,1,4),", ",str(Y,1,4),")"
  )
#end

#macro CIE_ReferenceWhite_String(S)
  concat(S,
  #switch (CIE_iID(ReferenceWhiteXY))
    #case(CIE_iID(Illuminant_A))   "Illuminant_A"   #break
    #case(CIE_iID(Illuminant_B))   "Illuminant_B"   #break
    #case(CIE_iID(Illuminant_C))   "Illuminant_C"   #break
    #case(CIE_iID(Illuminant_D50)) "Illuminant_D50" #break
    #case(CIE_iID(Illuminant_D55)) "Illuminant_D55" #break
    #case(CIE_iID(Illuminant_D65)) "Illuminant_D65" #break
    #case(CIE_iID(Illuminant_D75)) "Illuminant_D75" #break
    #case(CIE_iID(Illuminant_E))   "Illuminant_E"   #break
    #case(CIE_iID(Illuminant_F1))  "Illuminant_F1"  #break
    #case(CIE_iID(Illuminant_F2))  "Illuminant_F2"  #break
    #case(CIE_iID(Illuminant_F3))  "Illuminant_F3"  #break
    #case(CIE_iID(Illuminant_F4))  "Illuminant_F4"  #break
    #case(CIE_iID(Illuminant_F5))  "Illuminant_F5"  #break
    #case(CIE_iID(Illuminant_F6))  "Illuminant_F6"  #break
    #case(CIE_iID(Illuminant_F7))  "Illuminant_F7"  #break
    #case(CIE_iID(Illuminant_F8))  "Illuminant_F8"  #break
    #case(CIE_iID(Illuminant_F9))  "Illuminant_F9"  #break
    #case(CIE_iID(Illuminant_F10)) "Illuminant_F10" #break
    #case(CIE_iID(Illuminant_F11)) "Illuminant_F11" #break
    #case(CIE_iID(Illuminant_F12)) "Illuminant_F12" #break
    #else "unknown"
  #end
  )
#end

#macro CIE_ChromaticAdaption_String(S)
  concat(S,
  #switch (ChromaMatchFunction)
    #case(off) "None" #break
    #case(Scaling_ChromaMatch)  "Scaling"   #break
    #case(VonKries_ChromaMatch) "Von Kries" #break
    #case(Bradford_ChromaMatch) "Bradford"  #break
    #else "unknown"
  #end
  )
#end

//=====================================================================

// EOF
#version CIE_Inc_Temp;
#end