function [D,As,Ad]=adaptation_prep_XYZ(A,XYZwps,XYZwpd,degree)
% Creates preparatory arrays to pass to adaptation.m, which
% predicts how tristimulus values will change when moving from one
% illuminant to another, given the white points of the two
% illuminants. Illuminant white points are specified as tristimulus
% values, XYZ.
% A is nx3 matrix of color matching functions (CMFs) by columns.
% XYZwps is 3x1 vector of source illuminant white point tristimulus values.
% XYZwpd is 3x1 vector of destination illuminant white point tristimulus
% values. Note: XYZwps and XYZwpd should both have Y=1.
% degree is a value in the range [0-1] that indicates the degree
% of adaptation. Use degree=1 for full adaptation.
% D is nxn matrix of finite-differencing constants.
% As is nx3 matrix of Ws-referenced CMFs.
% Ad is nx3 matrix of Wd-referenced CMFs.
% Notes
% -----
% Version 2020_01_18 Scott Allen Burns
%
% This chromatic adaptation transform (CAT) is based on smoothest
% reflectance reconstruction (References [1] and [2]) and is
% developed and compared to existing CATs (Bradford, CAT02, CAT16)
% in Reference [3].
%
% References
% ----------
% [1] Burns SA. Numerical methods for smoothest reflectance
% reconstruction. Color Research & Application, Vol 45,
% No 1, 2020, pp 8-21.
% [2] Generating Reflectance Curves from sRGB Triplets, 2015
% http://scottburns.us/reflectance-curves-from-srgb/#LLSS
% [3] Burns SA. Chromatic adaptation transform by spectral
% reconstruction. Color Research & Application, Vol 44,
% No 5, 2019, pp 682-693.
n=size(A,1);
% form tri-diagonal matrix of finite differencing constants
D=full(gallery('tridiag',n,-2,4,-2));
D(1,1)=2;
D(n,n)=2;
% adjust destination illum for degree
XYZwpd = degree * XYZwpd + (1 - degree) * XYZwps;
% reconstruct illuminant spectra
Ws = method_2(D,A,XYZwps);
Wd = method_2(D,A,XYZwpd);
% form illuminant-weighted CMFs
As=diag(Ws)*A;
Ad=diag(Wd)*A;