function spectrum_plot
[lambda, D65] = illuminant('d65');
figure;
id=lambda>=360;
lambda=lambda(id);
D65= D65(id);
plot(lambda, D65)
hold on
for i=1:(numel(lambda)-1)
    X=[lambda(i) lambda(i) lambda(i+1) lambda(i+1) lambda(i)];
    Y=[0 D65(i) D65(i+1) 0 0];
    RGB=Wavelength2RGB(lambda(i));
    fill(X',Y',[RGB(1)/255 RGB(2)/255 RGB(3)/255],'EdgeColor',[RGB(1)/255 RGB(2)/255 RGB(3)/255])
end
title('D_{65} illuminant')

function RGB=Wavelength2RGB(Wavelength)
Gamma = 0.8;
IntensityMax = 255;

if Wavelength >= 380 & Wavelength <= 439
    Red = -(Wavelength - 440) / (440 - 350);
    Green = 0.0;
    Blue = 1.0;
elseif Wavelength >= 440 & Wavelength <= 489
    Red = 0.0;
    Green = (Wavelength - 440) / (490 - 440);
    Blue = 1.0;
elseif Wavelength >= 490 & Wavelength <= 509
    Red = 0.0;
    Green = 1.0;
    Blue = -(Wavelength - 510) / (510 - 490);
elseif Wavelength >= 510 & Wavelength <= 579
    Red = (Wavelength - 510) / (580 - 510);
    Green = 1.0;
    Blue = 0.0;
elseif Wavelength >= 580 & Wavelength <= 644
    Red = 1.0;
    Green = -(Wavelength - 645) / (645 - 580);
    Blue = 0.0;
elseif Wavelength >= 645 & Wavelength <= 780
    Red = 1.0;
    Green = 0.0;
    Blue = 0.0;
else
    Red = 0.0;
    Green = 0.0;
    Blue = 0.0;
end
  
if Wavelength >= 350 & Wavelength <= 419
    factor = 0.3 + 0.7*(Wavelength - 380)/(420 - 380);
elseif Wavelength >= 420 & Wavelength <= 700
    factor = 1.0;
elseif Wavelength >= 701 & Wavelength <= 780
    factor = 0.3 + 0.7*(780 - Wavelength)/(780 - 700);
else
    factor = 0.0;
end
r=255*(Red*factor)^Gamma;
g=255*(Green*factor)^Gamma;
b=255*(Blue*factor)^Gamma;

R=r;
G=g;
B=b;
RGB=[R G B];