function [L_D65,a_D65,b_D65]=C2D65O(L,a,b)
A=[0.8951 0.2664 -0.1614
    -0.7502 1.7135 0.0367
    0.0389 -0.0685 1.0296];
iA=inv(A);
i_C=[98.074 100 118.232]';
XW=i_C(1);
YW=i_C(2);
ZW=i_C(3);
R_C=A*i_C;

i_D65=[95.047 100 108.883]';
XW_D65=i_C(1);
YW_D65=i_C(2);
ZW_D65=i_C(3);

R_D65=A*i_D65;

M=diag(R_D65./R_C);
B_C_D65=iA*M*A;


corte=6/29;
coefi=1/3*(29/6)^2;
arg=(L+16)/116;
q=(arg.^3).*(arg>corte)+...
    ((arg-4/29)/coefi).*(arg<=corte);
Y=YW*q;

arg=a/500+(L+16)/116;
q=(arg.^3).*(arg>corte)+...
    ((arg-4/29)/coefi).*(arg<=corte);
X=XW*q;

arg=(L+16)/116-b/200;
q=(arg.^3).*(arg>corte)+...
    ((arg-4/29)/coefi).*(arg<=corte);
Z=ZW*q;

XYZ_C=[X;Y;Z];
XYZ_D65=B_C_D65*[X;Y;Z];

temp=XYZ_D65(1)/XW;
escala=1/3*(29/6)^2;
corte=(6/29)^3;
hhX=(escala*temp+4/29).*(temp<=corte) + (temp.^(1/3)).*(temp>corte);
temp=XYZ_D65(2)/YW;
hhY=(escala*temp+4/29).*(temp<=corte) + (temp.^(1/3)).*(temp>corte);
temp=XYZ_D65(3)/ZW;
hhZ=(escala*temp+4/29).*(temp<=corte) + (temp.^(1/3)).*(temp>corte);


L_D65=116*hhY-16;
a_D65=500*(hhX-hhY);
b_D65=200*(hhY-hhZ);

