function rho=LHTSS(T,sRGB)
% This is the Least Hyperbolic Tangent Slope Squared (LHTSS) 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 ((tanh(z)+1)/2) = rgb,
% using Lagrangian formulation and Newton's method.
% Reflectance will always be in the open interval (0->1).

% 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->1),

% Written by Scott Allen Burns, May 2019.
% Licensed under a Creative Commons Attribution-ShareAlike 4.0 International
% License (http://creativecommons.org/licenses/by-sa/4.0/).
% For more information, see http://scottburns.us/reflectance-curves-from-srgb/

% 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

% handle special case of (255,255,255)
if all(sRGB==255)
    rho=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 rho=1/2
lambda=zeros(3,1); % starting Lagrange mult
maxit=100; % max number of iterations
ftol=1.0e-8; % function solution tolerance
count=0; % iteration counter

% Newton's method iteration
while count <= maxit
    d0 = (tanh(z) + 1)/2;
    d1 = diag((sech(z).^2)/2);
    d2 = diag(-sech(z).^2.*tanh(z));
    F = [D*z + d1*T'*lambda; T*d0 - rgb]; % 39x1 F vector
    J = [D + diag(d2*T'*lambda), d1*T'; T*d1, zeros(3)]; % 39x39 J 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)
        % solution found
        rho=(tanh(z)+1)/2;
        return
    end
    count=count+1;
end
disp(['No solution found in ',num2str(maxit),' iterations.'])