function rgb2cb(args)
*********************************************************************
*********************************************************************
***
*** GrADS Script Name:  "rgb2cb.gs"
***                     by Kyle Ahern
***                     on 2022 August 31
***
*** Purpose:  Convert a given color expressed in RGB to 
***           colorspaces restricted by colorblindness.
***        
*** Usage pattern:
***
***           "run rgb2cb.gs type ID r g b a"
***           
***           Required arguments:
***           type  = type of colorblind colorspace to target.
***                   Must be either p (protanopia), 
***                                  d (deuteranopia), or
***                                  t (tritanopia).
***           ID    = number for the color in GrADS.
***           r     = red value (0 - 255).
***           g     = green value (0 - 255).
***           b     = blue value (0 - 255).
***
***           Optional arguments:
***           a     = alpha value (0 - 255).
***                   
*** Example Usage:
***
***           * Create a new color 30, with R,G,B = 255,177,97
***           * converted to the colorspace for protanopia.
***           "run rgb2cb.gs p 30 255 177 97"
***
*** Products:
***
***           + The user-specified RGB is converted to a colorspace
***             that simulates the a specified type of colorblindness.
***             The converted color is defined in GrADS (color #ID),
***             which can be used in shading. 
***
*********************************************************************
*********************************************************************

type=subwrd(args,1) ;* type of colorblindness
                    ;* p = protanopia (red-blind)
                    ;* d = deuteranopia (green-blind)
                    ;* t = tritanopia (blue-yellow-blind)

ID=subwrd(args,2);  ;* color number
r=subwrd(args,3);   ;* red component (0-255)
g=subwrd(args,4);   ;* green component (0-255)
b=subwrd(args,5);   ;* blue component (0-255)
a=subwrd(args,6);   ;* alpha component (0-255, optional)
                    ;* NOTE:  The mixing of translucent colors
                    ;*        that have been converted to simulate
                    ;*        colorblindness may yield results that
                    ;*        do not exist within the target
                    ;*        color space.

if (type = "" | ID = "" | r = "" | g = "" | b = "")
  say "Error: Not enough arguments specified (5 required)."
  say "arguments received: "
  say "   type: "type
  say "   ID  : "ID
  say "   r   : "r
  say "   g   : "g
  say "   b   : "b
  return;
endif

if (type = "p")
  srgb=convprotan(r, g, b);
else; if (type = "d")
  srgb=convdeuter(r, g, b);
else; if (type = "t")
  srgb=convtritan(r, g, b);
else
  say "Error: Unrecognized colorspace type ["type"]"
  say "Valid types are...
  say "   type = p (protanopia, red colorblindness)"       
  say "   type = d (deuteranopia, green colorblindness)"       
  say "   type = t (tritanopia, blue-yellow colorblindness)"
  return;
endif; endif; endif;

newr=subwrd(srgb,1);
newg=subwrd(srgb,2);
newb=subwrd(srgb,3);

'set rgb 'ID' 'newr' 'newg' 'newb' 'a

return

*************************************************

function convprotan(sr,sg,sb)

  ;* convert sRGB to raw RGB
  rgb=srgb2rgb(sr,sg,sb); r=subwrd(rgb,1); g=subwrd(rgb,2); b=subwrd(rgb,3);

  ;* convert raw RGB to XYZ
  xyz=rgb2xyz(r,g,b); x=subwrd(xyz,1); y=subwrd(xyz,2); z=subwrd(xyz,3);

  ;* convert XYZ to protanopia-perceived raw RGB
  prgb=xyz2prorgb(x,y,z); pr=subwrd(prgb,1); pg=subwrd(prgb,2); pb=subwrd(prgb,3);

  ;* convert protanopia-perceived raw RGB to sRGB
  srgb=rgb2srgb(pr,pg,pb);

return srgb

function convdeuter(sr,sg,sb)

  ;* convert sRGB to raw RGB
  rgb=srgb2rgb(sr,sg,sb); r=subwrd(rgb,1); g=subwrd(rgb,2); b=subwrd(rgb,3);

  ;* convert raw RGB to XYZ
  xyz=rgb2xyz(r,g,b); x=subwrd(xyz,1); y=subwrd(xyz,2); z=subwrd(xyz,3);

  ;* convert XYZ to deuteranopia-perceived raw RGB
  prgb=xyz2deurgb(x,y,z); pr=subwrd(prgb,1); pg=subwrd(prgb,2); pb=subwrd(prgb,3);

  ;* convert deuteranopia-perceived raw RGB to sRGB
  srgb=rgb2srgb(pr,pg,pb);

return srgb

function convtritan(sr,sg,sb)

  ;* convert sRGB to raw RGB
  rgb=srgb2rgb(sr,sg,sb); r=subwrd(rgb,1); g=subwrd(rgb,2); b=subwrd(rgb,3);

  ;* convert raw RGB to XYZ
  xyz=rgb2xyz(r,g,b); x=subwrd(xyz,1); y=subwrd(xyz,2); z=subwrd(xyz,3);

  ;* convert XYZ to tritanopia-perceived raw RGB
  prgb=xyz2trirgb(x,y,z); pr=subwrd(prgb,1); pg=subwrd(prgb,2); pb=subwrd(prgb,3);

  ;* convert tritanopia-perceived raw RGB to sRGB
  srgb=rgb2srgb(pr,pg,pb);

return srgb


**************************************************

function srgb2rgb(sr,sg,sb)

if (sr>255); sr=255; else; if (sr<0); sr=0; endif; endif;
if (sg>255); sg=255; else; if (sg<0); sg=0; endif; endif;
if (sb>255); sb=255; else; if (sb<0); sb=0; endif; endif;

i=sr/255.0;
j=sg/255.0;
k=sb/255.0;

if (i>0.04045)
  r=math_pow((i+0.055) / 1.055, 2.4);
else
  r=i / 12.92;
endif

if (j>0.04045)
  g=math_pow((j+0.055) / 1.055, 2.4);
else
  g=j / 12.92;
endif

if (k>0.04045)
  b=math_pow((k+0.055) / 1.055, 2.4);
else
  b=k / 12.92;
endif

retval=r" "g" "b

return retval

function rgb2srgb(r,g,b)

if (r>0.003486)
  sr=(1.055 * math_pow(r, 1.0/2.4)) - 0.055;
else
  sr=12.92 * r;
endif

if (g>0.003486)
  sg=(1.055 * math_pow(g, 1.0/2.4)) - 0.055;
else
  sg=12.92 * g;
endif

if (b>0.003486)
  sb=(1.055 * math_pow(b, 1.0/2.4)) - 0.055;
else
  sb=12.92 * b;
endif

sr=sr*255;
sg=sg*255;
sb=sb*255;

if (sr>255); sr=255; else; if (sr<0); sr=0; endif; endif;
if (sg>255); sg=255; else; if (sg<0); sg=0; endif; endif;
if (sb>255); sb=255; else; if (sb<0); sb=0; endif; endif;

retval=sr" "sg" "sb

return retval

function rgb2xyz(r,g,b)

X = (0.4124 * r) + (0.3576 * g) + (0.1805 * b);
Y = (0.2126 * r) + (0.7152 * g) + (0.0722 * b);
Z = (0.0193 * r) + (0.1192 * g) + (0.9505 * b);

retval=X" "Y" "Z

return retval

function xyz2prorgb(X,Y,Z)

* inverse of C-matrix (constant):
Ci.1.1 =  3.24062547732005;
Ci.1.2 = -1.53720797221032;
Ci.1.3 = -0.498628598698248;
Ci.2.1 = -0.968930714729319;
Ci.2.2 =  1.87575606088524;
Ci.2.3 =  0.0415175238429539;
Ci.3.1 =  0.0557101204455106;
Ci.3.2 = -0.204021050598487;
Ci.3.3 =  1.05699594225439;

* D-matrix for protanopia colorblindness:
D.1.1 = -0.381227242724273;
D.1.2 =  1.12278919714022;
D.1.3 =  0.0668758459598379;
D.2.1 = -0.468976237623762;
D.2.2 =  1.38122724272427;
D.2.3 =  0.0227067506750672;
D.3.1 =  0;
D.3.2 =  0;
D.3.3 =  1;

* multiply C-1 by D's first column to get M"s first column
M.1.1=(Ci.1.1 * D.1.1) + (Ci.1.2 * D.2.1) + (Ci.1.3 * D.3.1);
M.2.1=(Ci.2.1 * D.1.1) + (Ci.2.2 * D.2.1) + (Ci.2.3 * D.3.1);
M.3.1=(Ci.3.1 * D.1.1) + (Ci.3.2 * D.2.1) + (Ci.3.3 * D.3.1);

* get second column...
M.1.2=(Ci.1.1 * D.1.2) + (Ci.1.2 * D.2.2) + (Ci.1.3 * D.3.2);
M.2.2=(Ci.2.1 * D.1.2) + (Ci.2.2 * D.2.2) + (Ci.2.3 * D.3.2);
M.3.2=(Ci.3.1 * D.1.2) + (Ci.3.2 * D.2.2) + (Ci.3.3 * D.3.2);

* get third column...
M.1.3=(Ci.1.1 * D.1.3) + (Ci.1.2 * D.2.3) + (Ci.1.3 * D.3.3);
M.2.3=(Ci.2.1 * D.1.3) + (Ci.2.2 * D.2.3) + (Ci.2.3 * D.3.3);
M.3.3=(Ci.3.1 * D.1.3) + (Ci.3.2 * D.2.3) + (Ci.3.3 * D.3.3);

* protanopia raw RGB
r = (M.1.1 * X) + (M.1.2 * Y) + (M.1.3 * Z);
g = (M.2.1 * X) + (M.2.2 * Y) + (M.2.3 * Z);
b = (M.3.1 * X) + (M.3.2 * Y) + (M.3.3 * Z);

retval = r" "g" "b

return retval

function xyz2deurgb(X,Y,Z)

* inverse of C-matrix (constant):
Ci.1.1 =  3.24062547732005;
Ci.1.2 = -1.53720797221032;
Ci.1.3 = -0.498628598698248;
Ci.2.1 = -0.968930714729319;
Ci.2.2 =  1.87575606088524;
Ci.2.3 =  0.0415175238429539;
Ci.3.1 =  0.0557101204455106;
Ci.3.2 = -0.204021050598487;
Ci.3.3 =  1.05699594225439;

* D-matrix for deuteranopia colorblindness:
D.1.1 =  0.188418821882189;
D.1.2 =  0.659769581974319;
D.1.3 =  0.0392952711325292;
D.2.1 =  0.231773597359735;
D.2.2 =  0.811581178117809;
D.2.3 = -0.0112220522052205;
D.3.1 =  0;
D.3.2 =  0;
D.3.3 =  1;

* multiply C-1 by D's first column to get M"s first column
M.1.1=(Ci.1.1 * D.1.1) + (Ci.1.2 * D.2.1) + (Ci.1.3 * D.3.1);
M.2.1=(Ci.2.1 * D.1.1) + (Ci.2.2 * D.2.1) + (Ci.2.3 * D.3.1);
M.3.1=(Ci.3.1 * D.1.1) + (Ci.3.2 * D.2.1) + (Ci.3.3 * D.3.1);

* get second column...
M.1.2=(Ci.1.1 * D.1.2) + (Ci.1.2 * D.2.2) + (Ci.1.3 * D.3.2);
M.2.2=(Ci.2.1 * D.1.2) + (Ci.2.2 * D.2.2) + (Ci.2.3 * D.3.2);
M.3.2=(Ci.3.1 * D.1.2) + (Ci.3.2 * D.2.2) + (Ci.3.3 * D.3.2);

* get third column...
M.1.3=(Ci.1.1 * D.1.3) + (Ci.1.2 * D.2.3) + (Ci.1.3 * D.3.3);
M.2.3=(Ci.2.1 * D.1.3) + (Ci.2.2 * D.2.3) + (Ci.2.3 * D.3.3);
M.3.3=(Ci.3.1 * D.1.3) + (Ci.3.2 * D.2.3) + (Ci.3.3 * D.3.3);

* deuteranopia raw RGB
r = (M.1.1 * X) + (M.1.2 * Y) + (M.1.3 * Z);
g = (M.2.1 * X) + (M.2.2 * Y) + (M.2.3 * Z);
b = (M.3.1 * X) + (M.3.2 * Y) + (M.3.3 * Z);

retval = r" "g" "b

return retval

function xyz2trirgb(X,Y,Z)

* inverse of C-matrix (constant):
Ci.1.1 =  3.24062547732005;
Ci.1.2 = -1.53720797221032;
Ci.1.3 = -0.498628598698248;
Ci.2.1 = -0.968930714729319;
Ci.2.2 =  1.87575606088524;
Ci.2.3 =  0.0415175238429539;
Ci.3.1 =  0.0557101204455106;
Ci.3.2 = -0.204021050598487;
Ci.3.3 =  1.05699594225439;

* D-matrix for tritanopia colorblindness:
D.1.1 =  0.914109117508657
D.1.2 =  0.17251541488986
D.1.3 = -0.0904742417651181
D.2.1 =  0.029312797581535
D.2.2 =  0.941124025162149
D.2.3 =  0.0308770041508312
D.3.1 = -0.811907653788819
D.3.2 =  1.63075033906832
D.3.3 =  0.144766857329192

* multiply C-1 by D's first column to get M"s first column
M.1.1=(Ci.1.1 * D.1.1) + (Ci.1.2 * D.2.1) + (Ci.1.3 * D.3.1);
M.2.1=(Ci.2.1 * D.1.1) + (Ci.2.2 * D.2.1) + (Ci.2.3 * D.3.1);
M.3.1=(Ci.3.1 * D.1.1) + (Ci.3.2 * D.2.1) + (Ci.3.3 * D.3.1);

* get second column...
M.1.2=(Ci.1.1 * D.1.2) + (Ci.1.2 * D.2.2) + (Ci.1.3 * D.3.2);
M.2.2=(Ci.2.1 * D.1.2) + (Ci.2.2 * D.2.2) + (Ci.2.3 * D.3.2);
M.3.2=(Ci.3.1 * D.1.2) + (Ci.3.2 * D.2.2) + (Ci.3.3 * D.3.2);

* get third column...
M.1.3=(Ci.1.1 * D.1.3) + (Ci.1.2 * D.2.3) + (Ci.1.3 * D.3.3);
M.2.3=(Ci.2.1 * D.1.3) + (Ci.2.2 * D.2.3) + (Ci.2.3 * D.3.3);
M.3.3=(Ci.3.1 * D.1.3) + (Ci.3.2 * D.2.3) + (Ci.3.3 * D.3.3);

* tritanopia raw RGB
r = (M.1.1 * X) + (M.1.2 * Y) + (M.1.3 * Z);
g = (M.2.1 * X) + (M.2.2 * Y) + (M.2.3 * Z);
b = (M.3.1 * X) + (M.3.2 * Y) + (M.3.3 * Z);

retval = r" "g" "b

return retval

;* \^o^/