I'm writing an application to create and manipulate spectral power distributions (the physical basis of color). What I have thus far is kinda cool, you can click on the either of the spectrums on the left and it'll give you an spd (a gaussian distribution for the color spectrum and plank's formula for the blackbody spectrum).
The problem is the results are inaccurate, and I don't have any papers on the subject I can mimick...
The steps I'm taking to generate these colors are as follows
1) create a SPD
2) normalize it (divide by the max)
3) create an XYZ color by taking the dot product of the samples of the SPD with the color matching function from 
xyz.X = 0.0f;
xyz.Y = 0.0f;
xyz.Z = 0.0f;
for (u_int i=0; i<n_cmf_samples; ++i)
float I = spd.intensity(cmf_wl_samples[i]);
xyz.X += cmf_x_samples[i]*I;
xyz.Y += cmf_y_samples[i]*I;
xyz.Z += cmf_z_samples[i]*I;
4) convert the XYZ color to RGB using the matrix form  or maybe 
5) tone map the RGB color (C = C/(1 + C))
6) clamp the RGB color
The result I get isn't quite accurate, especially the black body spectrum. Any references or advice on how to improve my steps would be greatly appreciated!
Also if anyone understands how renderers use spds in their calculations let me know. I understand that luxrender uses sampling, but from their uncommented source code I can't figure it out...