function rho=ILSS(B11,B12,sRGB)
% This is the Iterative Least Slope Squared (ILSS) algorithm for generating
% a "reasonable" reflectance curve from a given sRGB color triplet.
% The reflectance spans the wavelength range 380-730 nm in 10 nm increments.

% It solves min sum(rho_i+1 - rho_i)^2 s.t. T rho = rgb, K1 rho = 1, K0 rho = 0,
% using Lagrangian formulation and iteration to keep all rho (0-1].

% B11 is upper-left 36x36 part of inv([D,T';T,zeros(3)])
% B12 is upper-right 36x3 part of inv([D,T';T,zeros(3)])
% sRGB is a 3-element vector of target D65-referenced sRGB values in 0-255 range,
% rho is a 36x1 vector of reflectance values (0->1] over wavelengths 380-730 nm,

% Written by Scott Allen Burns, 4/26/15.
% Licensed under a Creative Commons Attribution-ShareAlike 4.0 International
% License (http://creativecommons.org/licenses/by-sa/4.0/).
% For more information, see http://www.scottburns.us/subtractive-color-mixture/

rho=ones(36,1)/2; % initialize output to 0.5
rhomin=0.00001; % smallest refl value

% handle special case of (255,255,255)
if all(sRGB==255)
    rho=ones(36,1);
    return
end

% handle special case of (0,0,0)
if all(sRGB==0)
    rho=rhomin*ones(36,1);
    return
end

% compute target linear rgb values
sRGB=sRGB(:)/255; % convert to 0-1 column vector
rgb=zeros(3,1);
% remove gamma correction to get linear rgb
for i=1:3
    if sRGB(i)<0.04045
        rgb(i)=sRGB(i)/12.92;
    else
        rgb(i)=((sRGB(i)+0.055)/1.055)^2.4;
    end
end

R=B12*rgb;

% iteration to get all refl 0-1
maxit=10; % max iterations
count=0; % counter for iteration
while ( (any(rho>1) || any(rho<rhomin)) && count<=maxit ) || count==0
    
    % create K1 matrix for fixed refl at 1
    fixed_upper_logical = rho>=1;
    fixed_upper=find(fixed_upper_logical);
    num_upper=length(fixed_upper);
    K1=zeros(num_upper,36);
    for i=1:num_upper
        K1(i,fixed_upper(i))=1;
    end
    
    % create K0 matrix for fixed refl at rhomin
    fixed_lower_logical = rho<=rhomin;
    fixed_lower=find(fixed_lower_logical);
    num_lower=length(fixed_lower);
    K0=zeros(num_lower,36);
    for i=1:num_lower
        K0(i,fixed_lower(i))=1;
    end
    
    % set up linear system
    K=[K1;K0];
    C=B11*K'/(K*B11*K'); % M*K'*inv(K*M*K')
    rho=R-C*(K*R-[ones(num_upper,1);rhomin*ones(num_lower,1)]);
    rho(fixed_upper_logical)=1; % eliminate FP noise
    rho(fixed_lower_logical)=rhomin; % eliminate FP noise
    
    count=count+1;
end
if count>=maxit
    disp(['No solution found after ',num2str(maxit),' iterations.'])
end