function rho=LLSS(T,sRGB)
% This is the Least Log Slope Squared (LLSS) 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.

% Solves min sum(z_i+1 - z_i)^2 s.t. T exp(z) = rgb, where
% z=log(reflectance), using Lagrangian formulation and Newton's method.
% Allows reflectance values >1 to be in solution.

% T is 3x36 matrix converting reflectance to D65-weighted linear rgb,
% sRGB is a 3 element vector of target D65 referenced sRGB values in 0-255 range,
% rho is a 36x1 vector of reconstructed reflectance values, all > 0,

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

% initialize outputs to zeros
rho=zeros(36,1);

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

% 36x36 difference matrix for Jacobian
% having 4 on main diagonal and -2 on off diagonals,
% except first and last main diagonal are 2.
D=full(gallery('tridiag',36,-2,4,-2));
D(1,1)=2;
D(36,36)=2;

% compute target linear rgb values
sRGB=sRGB(:)/255; % convert to 0-1
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

% initialize
z=zeros(36,1); % starting point all zeros
lambda=zeros(3,1); % starting Lagrange mult
maxit=100; % max number of iterations
ftol=1.0e-8; % function solution tolerance
deltatol=1.0e-8; % change in oper pt tolerance
count=0; % iteration counter

% Newton's method iteration
while count <= maxit
    r=exp(z);
    v=-diag(r)*T'*lambda; % 36x1
    m1=-T*r; % 3x1
    m2=-T*diag(r); % 3x36
    F=[D*z+v;m1+rgb]; % 39x1 function vector
    J=[D+diag(v),m2';m2,zeros(3)]; % 39x39 Jacobian matrix
    delta=J\(-F); % solve Newton system of equations J*delta = -F
    z=z+delta(1:36); % update z
    lambda=lambda+delta(37:39); % update lambda
    if all(abs(F)<ftol) % check if functions satisfied
        if all(abs(delta)<deltatol) % check if variables converged
            % solution found
            disp(['Solution found after ',num2str(count),' iterations'])
            rho=exp(z);
            return
        end
    end
    count=count+1;
end
disp(['No solution found in ',num2str(maxit),' iterations.'])