%Chapter 7,  Fig 7.3
%Eric Dubois, updated 2018-12-12
%Plot spectral distribution of a colored light obtained by shining a white
%light through a filter, and of the metameric color displayed on a CRT.
%Data files should be placed in a data directory on the MATLAB path.
%Plot spectral transmission of color filter
clear all, close all;
%read and interpolate the spectral absorption curve
filterName = 'f830_MediumPink.txt';
%filterName = input('Enter filter file name: ','s')
absorp_in = dlmread(filterName,' ');
%interpolate the absorption curve
lami = 400:2:700;
absorp = interp1(absorp_in(:,1),absorp_in(:,2),lami,'spline')';

light = absorp;  % use equal energy white
%compute the XYZ tristimulus values of the light
%read and interpolate xyz color matching functions
xyz=dlmread('ciexyz31.txt',',');
lam = 360:5:830;
lami = 400:2:700;
xyzi = interp1(lam,xyz(:,2:4),lami,'spline');
%compute un-normalized tristimulus values
light_XYZ = light' * xyzi;
%chromaticities
light_xyz = light_XYZ/sum(light_XYZ);
%scale to have desired luminance Y0
y0 = 0.3;
light_XYZ0 = light_xyz *y0/light_xyz(2);
%Scale the light spectral density
light_scaled = light*light_XYZ0(1)/light_XYZ(1);
%Sony RGB to CIE 1931 XYZ  -- matrix from analyze_sony_rgb_primaries.m
A_RGBCRTtoXYZ = [0.464066  0.305612   0.180777; ...
                 0.259640  0.659360   0.081000; ...
                 0.035763  0.142251   0.910739];
%Convert to Sony RGB tristimulus values
light_RGBCRT = A_RGBCRTtoXYZ \ light_XYZ0';
%Gamma corrected
light_RGBCRT_gcorr = (light_RGBCRT).^(0.42);
load Sony_RGB_norm.mat;
%spectral density of displayed color
light_CRT = RGBnorm(16:166,:)*light_RGBCRT;
figure;

plot(lami,light_scaled,'k',lami,light_CRT','k')
xlabel('wavelength (nm)')
axis([400 700 0 .2])
xlabel('Wavelength (nm)');
ylabel('Normalized power');
set(gcf,'Color',[1 1 1]);
set(gca,'fontname','times')          %set the plot font to Times