BeginPackage["ColorFunctions`"]

SpectralColor::usage="SpectalColor[h] will return the RGBColor of sun light with wave length 380+400*h nm."

BlackBodyColor::usage="BlackBodyColor[h] will return the RGBColor of a black body with temperature 1000+10000*h Kelvin."

ChromaDepthColor::usage="ChromaDepthColor[h] will return the RGBColor for ChromaTeck images with black background."

ChromaDepthCMYColor::usage="ChromaDepthColor[h] will return the RGBColor for ChromaTeck images with white background."

LandscapeColor::usage="LandscapeColor[h] wil return a color common for maps changeing from green to brown and than to black."
Begin["`Private`"]

(* This visible color spectrum is adapted from the FORTRAN code 
in http://www.isc.tamu.edu/~astro/color/spectra.html to Mathematica. 
Since it use colors from 0..255  or a wavelength form 380-780 nm we 
calculate the wave length from the intensity value.*)

SpectralColor[h_]:=
   Module[{wl,col=RGBColor[0,0,0],intensity},
     wl=380.+h*400.;
     Which[ 
       wl>380. &&  wl<440.,col=RGBColor[-1.*(wl-440.)/(440.-380.),0,1],
       wl>=440. && wl<490., col=RGBColor  [0,(wl-440.)/50.,1],
       wl>=490. &&    wl<510., col=RGBColor[0,1,-(wl-510.)/20.],
       wl>=510. && wl <580., col=RGBColor[(wl-510.)/70.,1,0],
       wl>=580. && wl<= 645., col=RGBColor[1,-(wl-645.)/65.,0],
       wl>=645. && wl<=780., col=RGBColor[1,0,0]
      ];
     If[wl>700.,	
       intensity=.3+.7* (780.-wl)/80.,
       If[wl<420.,
         intensity=.3+.7*(wl-380.)/40.,
         intensity=1.
	]
      ];
     Map[#*intensity &,col]
   ]			

(* Calculate the color palette due to the black body radiation *)

$bbwave=Partition[{0.0014, 0.0000, 0.0065, 0.0022, 0.0001, 0.0105,
       0.0042, 0.0001, 0.0201, 0.0076, 0.0002, 0.0362,
       0.0143, 0.0004, 0.0679, 0.0232, 0.0006, 0.1102,
       0.0435, 0.0012, 0.2074, 0.0776, 0.0022, 0.3713,
       0.1344, 0.0040, 0.6456, 0.2148, 0.0073, 1.0391,
       0.2839, 0.0116, 1.3856, 0.3285, 0.0168, 1.6230,
       0.3483, 0.0230, 1.7471, 0.3481, 0.0298, 1.7826,
       0.3362, 0.0380, 1.7721, 0.3187, 0.0480, 1.7441,
       0.2908, 0.0600, 1.6692, 0.2511, 0.0739, 1.5281,
       0.1954, 0.0910, 1.2876, 0.1421, 0.1126, 1.0419,
       0.0956, 0.1390, 0.8130, 0.0580, 0.1693, 0.6162,
       0.0320, 0.2080, 0.4652, 0.0147, 0.2586, 0.3533,
       0.0049, 0.3230, 0.2720, 0.0024, 0.4073, 0.2123,
       0.0093, 0.5030, 0.1582, 0.0291, 0.6082, 0.1117,
       0.0633, 0.7100, 0.0782, 0.1096, 0.7932, 0.0573,
       0.1655, 0.8620, 0.0422, 0.2257, 0.9149, 0.0298,
       0.2904, 0.9540, 0.0203, 0.3597, 0.9803, 0.0134,
       0.4334, 0.9950, 0.0087, 0.5121, 1.0000, 0.0057,
       0.5945, 0.9950, 0.0039, 0.6784, 0.9786, 0.0027,
       0.7621, 0.9520, 0.0021, 0.8425, 0.9154, 0.0018,
       0.9163, 0.8700, 0.0017, 0.9786, 0.8163, 0.0014,
       1.0263, 0.7570, 0.0011, 1.0567, 0.6949, 0.0010,
       1.0622, 0.6310, 0.0008, 1.0456, 0.5668, 0.0006,
       1.0026, 0.5030, 0.0003, 0.9384, 0.4412, 0.0002,
       0.8544, 0.3810, 0.0002, 0.7514, 0.3210, 0.0001,
       0.6424, 0.2650, 0.0000, 0.5419, 0.2170, 0.0000,
       0.4479, 0.1750, 0.0000, 0.3608, 0.1382, 0.0000,
       0.2835, 0.1070, 0.0000, 0.2187, 0.0816, 0.0000,
       0.1649, 0.0610, 0.0000, 0.1212, 0.0446, 0.0000,
       0.0874, 0.0320, 0.0000, 0.0636, 0.0232, 0.0000,
       0.0468, 0.0170, 0.0000, 0.0329, 0.0119, 0.0000,
       0.0227, 0.0082, 0.0000, 0.0158, 0.0057, 0.0000,
       0.0114, 0.0041, 0.0000, 0.0081, 0.0029, 0.0000,
       0.0058, 0.0021, 0.0000, 0.0041, 0.0015, 0.0000,
       0.0029, 0.0010, 0.0000, 0.0020, 0.0007, 0.0000,
       0.0014, 0.0005, 0.0000, 0.0010, 0.0004, 0.0000,
       0.0007, 0.0002, 0.0000, 0.0005, 0.0002, 0.0000,
       0.0003, 0.0001, 0.0000, 0.0002, 0.0001, 0.0000,
       0.0002, 0.0001, 0.0000, 0.0001, 0.0000, 0.0000,
       0.0001, 0.0000, 0.0000, 0.0001, 0.0000, 0.0000,
       0.0000, 0.0000, 0.0000},3];

$pconst=1240./8.617*^-5;

planck=Compile[{{wl,_Real},{t,_Real}},
    3.74183*^-16*wl^(-5)/(Exp[$pconst/(wl*t)]-1.)]
			
BlackBody[t_]:=
 Module[{xyz,tt=N[t],i},
   xyz=Inner[Times,Table[planck[i,tt],{i,380,780,5}],$bbwave]; 
   xyz/Max[xyz]
  ]
   

makeChroma[{v1_,v2_}]:={v1,v2,1.-(v1+v2)}
		
$chromaR=makeChroma[{0.64,0.33}];
$chromaG=makeChroma[{0.29,0.6}];
$chromaB=makeChroma[{0.15,0.06}];
$chromaW=makeChroma[{0.3127,0.3291}];
$chromaDet=Det[{$chromaR,$chromaG,$chromaB}];

xyzToRGB[xyzR_,xyzG_,xyzB_,xyzC_]:=
  Module[{rgb},
    rgb={Det[{xyzC,xyzG,xyzB}],Det[{xyzR,xyzC,xyzB}],Det[{xyzR,xyzG,xyzC}]}/$chromaDet;
    Min[Max[#,0.],1.] & /@ rgb
   ]

BlackBodyColor[h_]:=
  Module[{rgb,t},
    t=1000.+10000.*h;
    rgb=BlackBody[t];
    RGBColor @@ xyzToRGB[$chromaR,$chromaG,$chromaB,rgb]
   ]

(* The two color palettes from 
   
   http://www.chromatek.com/color_palettes.html

   for ChromaDepth pictures *)

redFunction=Compile[{{range,_Real}},
     Which[
       range>0.75 || range<0., 0.,
       True, (((-2.13*range-1.07)*range+0.133)*range+0.0667)*range+1
      ]
    ]  

greenFunction=Compile[{{range,_Real}},
       If[
         range<=0.5,(1.6*range+1.2)*range,
	 (3.2*range-6.8)*range+3.6
        ]
      ]                                                                        

blueFunction=
  Compile[{{range,_Real}},If[range>0.5,(-4.8*range+9.2)*range-3.4,0.]]

ChromaDepthColor[h_]:=
  Module[{rgb},
    rgb={redFunction[#],greenFunction[#],blueFunction[#]} &[N[h]];	                                                      
    RGBColor @@ (Max[0.,Min[#,1.]] & /@ rgb)
   ]     

ChromaDepthCMYColor[h_]:=1.-# & /@ ChromaDepthColor[h] 

lscape=Compile[{{v,_Real}},{0.25660999433616455 + v*(-0.15450689601371334 + 
    v*(16.99064278686296 + v*(-49.77485354775084 + 
        (53.311185629847536 - 20.331201294602224*v)*v))), 
 0.5950036100590319 + v*(0.9030709139591158 + 
    v*(-4.7684534226105955 + v*(39.51236857804843 + 
        v*(-117.41187124896338 + (131.37166677177956 - 
            50.03797646732823*v)*v)))), 0.20815128511236033 + 
  v*(0.911878579970309 + v*(-4.479304425265597 + 
      v*(36.027512397958255 + v*(-101.52716099689235 + 
          (110.30681084859079 - 41.31008420034276*v)*v))))}]
          
LandscapeColor[v_]:=RGBColor @@ lscape[v]          

End[] (* Private *)
EndPackage[]