function [D,As,Ad]=adaptation_prep_xy(A,xy_wps,xy_wpd,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 chromaticity
% coordinates, xy.

% A is nx3 matrix of color matching functions (CMFs) by columns.
% xy_wps is two-element vector of source illuminant white point
%        chromaticity coordinates, (x,y).
% xy_wpd is a two-element vector of destination illuminant white
%         point chromaticity coordinates, (x,y).
% 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
%
% A list of xy chromaticity coordinates for various standard
% illuminants can be found at Wikipedia:
% https://en.wikipedia.org/wiki/Standard_illuminant
% 
% 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;
% compute illuminant tristimulus values from chromaticity coords
XYZwps = zeros(3,1);
XYZwps(1) = xy_wps(1)/xy_wps(2);
XYZwps(2) = 1;
XYZwps(3) = (1-xy_wps(1)-xy_wps(2))/xy_wps(2);
XYZwpd = zeros(3,1);
XYZwpd(1) = xy_wpd(1)/xy_wpd(2);
XYZwpd(2) = 1;
XYZwpd(3) = (1-xy_wpd(1)-xy_wpd(2))/xy_wpd(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;
