
function [xyz] = inv_ciecam_jc(JCh, xyzw, LA, c, Nc, FLL, F, Yb, D);

%
%  INV_CIECAM_JC	function [xyz] = inv_ciecam_jc(JCh, xyzw, LA, 
%                                                      c, Nc, FLL, F, 
%                                                      Yb, D)
%
%       calculates XYZ from CIECAM97s JCh coordinates
%
%  inputs:
%       
%       JCh     JCh color vector
%       xyzw    referencey white in XYZ colorspace
%       LA      absolute adapting luminance
%       c       impact of surround              (see table)
%       Nc      chromatic induction factor      (see table)
%       FLL     lightness contrast factor       (see table)
%       F       degree of adaptation            (see table)
%       Yb      background relative luminance (in the 
%                 same units of xyzw)
%       D       =1 to discount the illuminant
%
%  outputs:
%
%       xyz     color in XYZ colorspace
%
%
%  ref: Condition               c       Nc      FLL     F
%       Avg., Stim.>4deg.       0.69    1.0     0.0     1.0
%       Avg.                    0.69    1.0     1.0     1.0
%       Dim                     0.59    1.1     1.0     0.9
%       Dark                    0.525   0.8     1.0     0.9
%       Cut Sheet               0.41    0.8     1.0     0.9
%
%  (adapted from IDL code by Mark Fairchild)
%  Krishna S. Nayak 3/13/98
%  Leland Stanford Junior University
%

MB =    [       0.8951,         0.2664,        -0.1614 ;  ...
               -0.7502,         1.7135,         0.0367 ;  ...
                0.0389,        -0.0685,         1.0296 ];

MH =    [       0.38971,        0.68898,       -0.07868 ; ...
               -0.22981,        1.18340,        0.04641 ; ...
                0.00,           0.00,           1.00    ];

rgbw = MB * (xyzw/xyzw(2));

if (D == 1.0) 
  D = F - (F/(1.+2.*LA^0.25+(LA^2.)/300.));
end
p = rgbw(3)^0.0834;

rgbcw(1) = (D*(1.0/rgbw(1))+1.0-D)*rgbw(1);
rgbcw(2) = (D*(1.0/rgbw(2))+1.0-D)*rgbw(2);
rgbcw(3) = (rgbw(3)/abs(rgbw(3)))*((D*(1.0/(rgbw(3)^p))+1.0-D)*abs(rgbw(3))^p);

MBinv = inv(MB);
MHinv = inv(MH);

rgbprimew = MH * (MBinv * (rgbcw'*xyzw(2)));

k = 1.0/(5.0*LA +1);
FL = 0.2*(k^4.0)*5.0*LA+0.1*((1-k^4.0)^2.0)*(5.0*LA)^(1./3.);
n = Yb/xyzw(2);
Nbb = 0.725*(1.0/n)^0.2;
Ncb = 0.725*(1.0/n)^0.2;
z = 1.0+FLL*n^0.5;

rgbaw = (40.0*(FL*rgbprimew/100.).^0.73) ./     ...
        ((FL*rgbprimew/100).^0.73+2.0)   + 1.0;

Awhite = (2*rgbaw(1)+rgbaw(2)+(1.0/20.0)*rgbaw(3)-2.05)*Nbb;

h = JCh(3);

if (h <= 20.14) 
  e1 = 0.8565; e2 = 0.8; h1 = 0.0; h2 = 20.14;
elseif (h <= 90.0) 
  e1 = 0.8; e2 = 0.7; h1 = 20.14; h2 = 90.0;
elseif (h <= 164.25) 
  e1 = 0.7; e2 = 1.0; h1 = 90.0; h2 = 164.25;
elseif (h <= 237.53) 
  e1 = 1.0; e2 = 1.2; h1 = 164.25; h2 = 237.53;
else
  e1 = 1.2; e2 = 0.8565; h1 = 237.53; h2 = 360.0;
end

e = e1 +(e2-e1)*(h-h1)/(h2-h1);

Astim = ((JCh(1)/100.0)^(1.0/(c*z)))*Awhite;

s = (JCh(2)^(1.0/0.69)) / ((2.44*((JCh(1)/100.)^(0.67*n))*(1.64-0.29^n))^(1.0/0.69));

if ((h >= 90.0) & (h < 270.0)) 
  sign = -1.0;
else 
  sign = 1.0;
end

a = (s*((Astim/Nbb)+2.05))/(sign*sqrt(1.0+tan(pi*h/180.)^2.0)*(50000.0*e*Nc*Ncb/13.0)+ (s*((11.0/23.0)+(108.0/23.0)*tan(pi*h/180.))));
		
b = a*tan(pi*h/180.);

rgba(1) = (20.0/61.0)*((Astim/Nbb)+2.05)+(41.0/61.0)*(11.0/23.0)*a+(288.0/61.0)*(1.0/23.0)*b;
rgba(2) = (20.0/61.0)*((Astim/Nbb)+2.05)-(81.0/61.0)*(11.0/23.0)*a-(288.0/61.0)*(1.0/23.0)*b;
rgba(3) = (20.0/61.0)*((Astim/Nbb)+2.05)-(20.0/61.0)*(11.0/23.0)*a-(20.0/61.0)*(315.0/23.0)*b;

if ((rgba(1)-1.0) < 0.0) 
  rgbprime(1) = -100.0*((2.0-2.0*rgba(1))/(39.0+rgba(1)))^(1.0/0.73);
else
  rgbprime(1) = 100.0*((2.0*rgba(1)-2.0)/(41.0-rgba(1)))^(1.0/0.73);
end
	
if ((rgba(2)-1.0) < 0.0) 
  rgbprime(2) = -100.0*((2.0-2.0*rgba(2))/(39.0+rgba(2)))^(1.0/0.73);
else
  rgbprime(2) = 100.0*((2.0*rgba(2)-2.0)/(41.0-rgba(2)))^(1.0/0.73);
end
	
if ((rgba(3)-1.0) < 0.0) 
  rgbprime(3) = -100.0*((2.0-2.0*rgba(3))/(39.0+rgba(3)))^(1.0/0.73);
else
  rgbprime(3) = 100.0*((2.0*rgba(3)-2.0)/(41.0-rgba(3)))^(1.0/0.73);
end

rgbcY = MB * (MHinv * (rgbprime'/FL));

Yc = 0.43231*rgbcY(1)+0.51836*rgbcY(2)+0.04929*rgbcY(3);

messR = (rgbcY(1)/Yc)/(D*(1.0/rgbw(1))+1.0-D);
messG = (rgbcY(2)/Yc)/(D*(1.0/rgbw(2))+1.0-D);
messB = ((rgbcY(3)/Yc)/(abs(rgbcY(3)/Yc)))*((abs(rgbcY(3)/Yc))^(1.0/p))/((D*(1.0/rgbw(3)^p)+1.0-D)^(1.0/p));

Yprime = 0.43231*messR*Yc+0.51836*messG*Yc+0.04929*messB*Yc;

temparray = [messR*Yc, messG*Yc, (messB*Yc)/((Yprime/Yc)^((1.0/p)-1.0))];

xyz = MBinv * temparray';

