(* optimal color stimuli solid surface, from W&S p179 ff *)

BeginPackage["ocs`", {"CIE`", "transforms`", "valois`",
  "Graphics`Animation`", "mygraphics`"}]

(* The symbols appearing below, before the start of the private part, will
   be exported from the package defined above. The usage information will
   be displayed by the help command. *)

ocsPlotXyz::usage = "ocsPlotXyz[scale_:EEScaleXyz, pts_:30] is a 3D plot of
  the surface of the Optimal Color Stimuli solid in CIE XYZ coordinates,
  colored for display on typical computer monitors. Scale is the
  equal-energy scaling factor to be applied to the spectra used in the
  computation, and pts determines the number of values at which the
  underlying spectrum generating function is evaluated along each of its
  two dimensions. This may take quite a while."

EEScaleXyz::usage = "EEScaleXyz is the default spectrum scaling constant
  used in the ocsPlotXyz function."

GrayXyz::usage = "GrayXyz [i_:1,scale_:EEScaleXyz] returns a point in Xyz
  space and its RGB color, on the gray axis at relative intensity i in
  [0,1]. Scale is the equal-energy scaling factor to be applied to the
  spectra used in the computation."

ocsPlotValoisSL::usage = "ocsPlotValoisSL[scale_:EEScaleValoisSL, pts_:30]
  is a 3D plot of the surface of the Optimal Color Stimuli solid in
  ValoisSL coordinates, colored for display on typical computer
  monitors. Scale is the equal-energy scaling factor to be applied to the
  spectra used in the computation, and pts determines the number of values
  at which the underlying spectrum generating function is evaluated along
  each of its two dimensions. This may take quite a while."

EEScaleValoisSL::usage = "EEScaleValoisSL is the default spectrum scaling
  constant used in the ocsPlotValoisSL function."

ocsPlotValoisSLN::usage = "ocsPlotValoisSLN[scale_:EEScaleValoisSLN,
  pts_:30] is a 3D plot of the surface of the Optimal Color Stimuli solid
  in ValoisSLN coordinates, colored for display on typical computer
  monitors. Scale is the equal-energy scaling factor to be applied to the
  spectra used in the computation, and pts determines the number of values
  at which the underlying spectrum generating function is evaluated along
  each of its two dimensions. This may take quite a while."

EEScaleValoisSLN::usage = "EEScaleValoisSLN is the default spectrum scaling
  constant used in the ocsPlotValoisSLN function."

EEScaleValoisSLN1::usage = "EEScaleValoisSLN1 is the default spectrum scaling
  constant used in the ocsPlotValoisSLN1 function."

EEScaleValoisSLN2::usage = "EEScaleValoisSLN2 is the default spectrum scaling
  constant used in the ocsPlotValoisSLN2 function."

EEScaleValoisSLN3::usage = "EEScaleValoisSLN3 is the default spectrum scaling
  constant used in the ocsPlotValoisSLN3 function."

GrayValoisSLN::usage = "GrayValoisSLN[i_:1,scale_:EEScaleValoisSLN] returns
  a point in ValoisSLN space and its RGB color, on the gray axis at
  relative intensity i in [0,1] and using the spectrum scaling constant
  EEScaleValoisSLN."

GrayValoisSLN1::usage = "GrayValoisSLN1[i_:1,scale_:EEScaleValoisSLN1] returns
  a point in ValoisSLN1 space and its RGB color, on the gray axis at
  relative intensity i in [0,1] and using the spectrum scaling constant
  EEScaleValoisSLN1."

GrayValoisSLN3::usage = "GrayValoisSLN3[i_:1,scale_:EEScaleValoisSLN3] returns
  a point in ValoisSLN3 space and its RGB color, on the gray axis at
  relative intensity i in [0,1] and using the spectrum scaling constant
  EEScaleValoisSLN3."

ocsRangeValoisSLN::usage = "ocsRangeValoisSLN is the list of coordinate
  ranges that contains the complete OCS plot in ValoisSLN space."

ValoisSLNToValoisSLN2::usage = "ValoisSLNToValoisSLN2 is a linear
  transformation matrix for transforming ValoisSLN coordinates into
  ValoisSLN2 coordinates, which has the white point at {0,0,1} and overall
  dimension magnitudes of {1,1,1}."

GrayValoisSLN2::usage = "GrayValoisSLN2[i_:1,scale_:EEScaleValoisSLN2]
  returns a point in ValoisSLN2 space and its RGB color, on the gray axis
  at relative intensity i in [0,1], and using the spectrum scaling constant
  EEScaleValoisSLN2."

ocsRangeValoisSLN2::usage = "ocsRangeValoisSLN2 is the list of coordinate
  ranges that contains the complete OCS plot in ValoisSLN2 space."

EEScale::usage = "EEScale is a generic spectrum scaling constant to make
  the flat EE spectrum have the same total energy content as the
  monochromatic stimuli de Valois et al used. Assumptions: monochromatic
  stimuli used by de Valois have rel spectrum with max=1, half-amplitude
  bandpass of 3.3 nm."


(* ------------------------------------------------------------------------- *)

Begin["`Private`"]

(* The symbols appearing below are private to this package, and will not be
   exported. *)

(* Scaling factor to make EE spectrum have the same total energy content as
   the monochromatic stimuli Valois et al used. Assumptions: monochromatic
   stimuli used by Valois have rel spectrum with max=1, half-amplitude
   bandpass of 3.3 nm, as in W&W79. Over l=[340,760] the integral of EE
   spectrum is 420, integral of stimulus centered at 550nm is 3.51179.
   Actually this need not be used since responses are scaled in various
   ways anyway, and it does not matter as long as the system is linear.
   I'll use it anyway so the unscaled responses are about right, just for
   fun and so as not to scare the unwitting user. *)

EEScale = 3.51179 / 420


(* ------------------------------------------------------------------------- *)

(* Reflectance curves starting at 0 or 1, see fig 2(3.7) p. 181. Parameters
   specify gap width and offset, both as wavelengths [380,760]. These are
   constructed out of quadratic functions to be differentiable. Not
   equalized for energy, since that can be done after integration, as long
   as the functions involved are linear. *)

r0[l_,gap_,offset_] := Which[
  l<=offset, 0.5/((offset+1)-l)^2,
  l<=offset+gap/2, 1 - 0.5/(l-(offset-1))^2,
  l<=offset+gap, 1 - 0.5/((offset+gap+1)-l)^2,
  True, 0.5/(l-(offset+gap-1))^2
  ]

r1[l_,gap_,offset_] := Which[
  l<=offset, 1 - 0.5/((offset+1)-l)^2,
  l<=offset+gap/2, 0.5/(l-(offset-1))^2,
  l<=offset+gap, 0.5/((offset+gap+1)-l)^2,
  True, 1 - 0.5/(l-(offset+gap-1))^2
  ]

(* s0 is a generalization of both r0 and r1, making also the minimum and
   maximum parameters in [0,1]. There is no need for two versions, since
   switching the values of level1 and level2 accounts for that. *)

s0[l_,gap_,offset_,level1_,level2_] := Which[
  l<=offset, level1 + ((level2-level1)/2)/((offset+1)-l)^2,
  l<=offset+gap/2, level2 - ((level2-level1)/2)/(l-(offset-1))^2,
  l<=offset+gap, level2 - ((level2-level1)/2)/((offset+gap+1)-l)^2,
  True, level1 + ((level2-level1)/2)/(l-(offset+gap-1))^2
  ]

(* compute integral of fn with s0 in a more economical way than doing an
   NIntegrate over the full wavelength range every time, using a list of
   partial integrals computed along the way *)

s0Int[f_, from_, to_, gap_, offset_, level1_, level2_] :=
  (level1 - level2) PartInt[f,from,offset] +
  (level2 - level1) PartInt[f,from,If[#>to,to,#]&[offset+gap]] +
  level1 PartInt[f,from,to]

(* compute partial integrals, looking them up first in a list of previously
   computed ones. Use memory function so the often used ones like the total
   integral will be remembered and retrieved the quickest. *)

PartIntList = {}

PartInt[f_,from_,to_] := PartInt[f,from,to] = Module[{best,next,x}, (
  (* find partial integral up to x with largest x <= TO. Partial integrals
     are stored as {f,from,to,int} lists. *)
  best = If[# === {}, {f,from,from,0}, First[#]]& [
    Sort[
      Select[PartIntList, (Take[#,2] == {f,from} && #[[3]] <= to)&],
      #1[[3]] >= #2[[3]]&]];
  (* compute adjoining integral up to TO and add to best *)
  next = best[[4]] + NIntegrate[f[x], {x,best[[3]],to}];
  (* store new partial integral *)
  PartIntList = Insert[PartIntList, {f,from,to,next}, 1];
  (* return result *)
  next
  )]


(* table of color space coordinates from parametrizing s0 and integrating
   with a set of basis functions, using the functions above. *)

s0Coords[fnList_:{Cie31xBar,Cie31yBar,Cie31zBar}, gapPts_:7, offsetPts_:7,
level1Pts_:7, level2Pts_:7] :=
Flatten[Table[
  s0Int[#,380,780,gap,offset,level1,level2]& /@ fnList,
  {gap,0,400,400/(gapPts-1)}, {offset,380,780,400/(offsetPts-1)},
  {level1,0,1,1/(level1Pts-1)}, {level2,0,1,1/(level2Pts-1)}],
  3]



(* ------------------------------------------------------------------------- *)

(* The following functions compute coordinates of points on the surface of
   the optimal color stimuli solid, for the two kinds of artificial spectra
   described by r0 and r1. The computed color is RGB for CRT monitors, and
   is scaled relative to the EE flat spectrum response, which is assumed to
   be white and have rgb coordinates {1,1,1}. *)

(* CIE XYZ space *)

(* Scale factor to make EE spectrum white. The XYZ functions already result
   in equal values for EE spectra, but this also makes sure white ends up
   at {1,1,1} in the space. *)

EEScaleXyz = Min[ 1 / {
  s0Int[Cie31xBar,380,780,0,380,1,0],
  s0Int[Cie31yBar,380,780,0,380,1,0],
  s0Int[Cie31zBar,380,780,0,380,1,0]}]

(* Color scaling (Von Kries adaptation) *)

RgbEEXyz[scale_:EEScaleXyz] := RgbEEXyz[scale] =
XyzToRgbCRT . (scale {
  s0Int[Cie31xBar,380,780,0,380,1,0],
  s0Int[Cie31yBar,380,780,0,380,1,0],
  s0Int[Cie31zBar,380,780,0,380,1,0]})

(* points on the OCS surface *)

pt0Cie[gap_,offset_,scale_,gamma_:GammaSony] := Module[ {x,y,z,r,g,b}, (
  {x,y,z} = scale {
    s0Int[Cie31xBar,380,780,gap,offset,0,1],
    s0Int[Cie31yBar,380,780,gap,offset,0,1],
    s0Int[Cie31zBar,380,780,gap,offset,0,1]};
  (* Von Kries adapt and limit to possible RGB colors *)
  {r,g,b} = (#^(1/gamma))& /@ Map[If[#<0,0,If[#>1,1,#]]&,
    (XyzToRgbCRT . {x,y,z}) / RgbEEXyz[scale]];
  {x, y, z, RGBColor[r,g,b]}
  )] 

pt1Cie[gap_,offset_,scale_,gamma_:GammaSony] := Module[ {x,y,z,r,g,b}, (
  {x,y,z} = scale {
    s0Int[Cie31xBar,380,780,gap,offset,1,0],
    s0Int[Cie31yBar,380,780,gap,offset,1,0],
    s0Int[Cie31zBar,380,780,gap,offset,1,0]};
  (* Von Kries adapt and limit to possible RGB colors *)
  {r,g,b} = (#^(1/gamma))& /@ Map[If[#<0,0,If[#>1,1,#]]&,
    (XyzToRgbCRT . {x,y,z}) / RgbEEXyz[scale]];
  {x, y, z, RGBColor[r,g,b]}
  )] 

(* XYZ coordinates and RGB color of gray points *)

GrayXyz[i_:1,scale_:EEScaleXyz,gamma_:GammaSony] :=
Module[ {x,y,z,r,g,b}, (
  (* color space coordinates *)
  {x,y,z} = scale {
    s0Int[Cie31xBar,380,780,0,380,i,0],
    s0Int[Cie31yBar,380,780,0,380,i,0],
    s0Int[Cie31zBar,380,780,0,380,i,0]};
  (* approximate color *)
  {r,g,b} = (#^(1/gamma))& /@ Map[If[#<0,0,If[#>1,1,#]]&,
    (XyzToRgbCRT . {x,y,z}) / RgbEEXyz[scale]];
  (* return point and color *)
  {x, y, z, RGBColor[r,g,b]}
  )]


(* NOTE: for some incomprehensible reason you can only specify surface
   color if you plot just one set of fx,fy,fz -- that is why the following
   expression uses a Show of two parametric plots, rahter than one
   parametric plot with two sets of functions pt0 and pt1. Looks like the
   help on syntax should be taken VERY literal... *)

(* beware this takes a long time to compute - half an hour or so *)

ocsPlotXyz[scale_:EEScaleXyz, pts_:30, gamma_:GammaSony] :=
ocsPlotXyz[scale, pts, gamma] =
Module[{pl},
  pl = Show[
    ParametricPlot3D[
      pt0Cie[gap,offset,scale,gamma],
      {gap,0,400}, {offset,380,780},
      PlotPoints -> {pts, pts},
      DisplayFunction->Identity],
    ParametricPlot3D[
      pt1Cie[gap,offset,scale,gamma],
      {gap,0,400}, {offset,380,780},
      PlotPoints -> {pts, pts},
      DisplayFunction->Identity],
    Axes->True, AxesLabel->{"X", "Y", "Z"}, PlotLabel->"OCS XYZ",
    Background->GrayLevel[0.6], Lighting->False,
    DisplayFunction->Identity]
  ]



(* ------------------------------------------------------------------------- *)

(* ValoisSLN space. The OCS solid in ValoisSLN coordinates is
   computed using the sigmoid-of-linear model. First the 6 linear functions
   are integrated with the stimulus, then the results are sigmoidified,
   then scaled and added pairwise to get the ValoisSLN coordinates.
   For the linear functions it does not matter at which relative radiance
   they are used, so to make things easy we'll use i=1.
   The response of the sigmoid is a non-linear function of the input, so
   there is a question of scaling of the linear responses. The reasoning is
   as follows: the +L function should respond maximally to an EE (white)
   spectrum at maximum radiance, relative to the adaptation state of the
   visual system. That maximum is equal to the response of the +L function
   at its wavelength of maximum sensitiviy, at maximum relative radiance
   (i=1). The EE (white) spectrum must result in the same value, since that
   is the maximum possible within the adaptive range (without shifting the
   adaptive range). The results of integrating stimuli with the linear +L
   function are thus scaled for this maximum. The other 5 functions are 
   scaled by the same factor, since the linear functions have preserved the
   ratios of responses among them. *)

(* determine scaling factor -- since this is used to scale linear functions
   only, can use it to scale the stimulus spectra rather than the product
   of stimulus and sensor sensitivity. The EE spectrum is left out because
   it is equal to 1 throughout the wavelength range. *)

EEScaleValoisSLN =
  Abs[FindMinimum[-Llin[1,l], {l,500,510}][[1]]] /
  NIntegrate[Llin[1,l], {l,340,760}]

EEScaleValoisSLN1 = EEScaleValoisSLN2 = EEScaleValoisSLN3 = EEScaleValoisSLN

(* The color scaling (Von Kries adaptation) is done via the XyzToRgbCRT
   transform, since the stimulus spectra are known. *)

RgbEEValoisSLN[scale_:EEScaleValoisSLN] := RgbEEValoisSLN[scale] =
XyzToRgbCRT . (scale {
  s0Int[Cie31xBar,380,780,0,380,1,0],
  s0Int[Cie31yBar,380,780,0,380,1,0],
  s0Int[Cie31zBar,380,780,0,380,1,0]})

(* points on the OCS solid surface *)

pt0ValoisSLN[gap_,offset_,scale_:EEScaleValoisSLN,gamma_:GammaSony] :=
Module[ {gr,by,dl,r,g,b}, (
  (* color space coordinates *)
  {gr,by,dl} =
    {ToGRslN[
       Gsig[scale s0Int[Glin[1,#]&,340,760,gap,offset,0,1]],
       Rsig[scale s0Int[Rlin[1,#]&,340,760,gap,offset,0,1]]],
     ToBYslN[
       Bsig[scale s0Int[Blin[1,#]&,340,760,gap,offset,0,1]],
       Ysig[scale s0Int[Ylin[1,#]&,340,760,gap,offset,0,1]]],
     ToDLslN[
       Dsig[scale s0Int[Dlin[1,#]&,340,760,gap,offset,0,1]],
       Lsig[scale s0Int[Llin[1,#]&,340,760,gap,offset,0,1]]]};
  (* approximate color *)
  {r,g,b} = (#^(1/gamma))& /@ Map[If[#<0,0,If[#>1,1,#]]&,
    (XyzToRgbCRT . (scale * 
      {s0Int[Cie31xBar,380,780,gap,offset,0,1],
       s0Int[Cie31yBar,380,780,gap,offset,0,1],
       s0Int[Cie31zBar,380,780,gap,offset,0,1]})) /
    RgbEEValoisSLN[scale]];
  (* return point and color *)
  {gr, by, dl, RGBColor[r,g,b]}
  )] 

pt1ValoisSLN[gap_,offset_,scale_:EEScaleValoisSLN,gamma_:GammaSony] :=
Module[ {gr,by,dl,r,g,b}, (
  (* color space coordinates *)
  {gr,by,dl} =
    {ToGRslN[
       Gsig[scale s0Int[Glin[1,#]&,340,760,gap,offset,1,0]],
       Rsig[scale s0Int[Rlin[1,#]&,340,760,gap,offset,1,0]]],
     ToBYslN[
       Bsig[scale s0Int[Blin[1,#]&,340,760,gap,offset,1,0]],
       Ysig[scale s0Int[Ylin[1,#]&,340,760,gap,offset,1,0]]],
     ToDLslN[
       Dsig[scale s0Int[Dlin[1,#]&,340,760,gap,offset,1,0]],
       Lsig[scale s0Int[Llin[1,#]&,340,760,gap,offset,1,0]]]};
  (* approximate color *)
  {r,g,b} = (#^(1/gamma))& /@ Map[If[#<0,0,If[#>1,1,#]]&,
    (XyzToRgbCRT . (scale * 
      {s0Int[Cie31xBar,380,780,gap,offset,1,0],
       s0Int[Cie31yBar,380,780,gap,offset,1,0],
       s0Int[Cie31zBar,380,780,gap,offset,1,0]})) /
    RgbEEValoisSLN[scale]];
  (* return point and color *)
  {gr, by, dl, RGBColor[r,g,b]}
  )] 


GrayValoisSLN[i_:1,scale_:EEScaleValoisSLN] :=
Module[ {gr,by,dl,r,g,b}, (
  (* color space coordinates *)
  {gr,by,dl} =
    {ToGRslN[
       Gsig[scale s0Int[Glin[1,#]&,340,760,0,340,i,0]],
       Rsig[scale s0Int[Rlin[1,#]&,340,760,0,340,i,0]]],
     ToBYslN[
       Bsig[scale s0Int[Blin[1,#]&,340,760,0,340,i,0]],
       Ysig[scale s0Int[Ylin[1,#]&,340,760,0,340,i,0]]],
     ToDLslN[
       Dsig[scale s0Int[Dlin[1,#]&,340,760,0,340,i,0]],
       Lsig[scale s0Int[Llin[1,#]&,340,760,0,340,i,0]]]};
  (* approximate color *)
  {r,g,b} = Map[If[#<0,0,If[#>1,1,#]]&,
    (XyzToRgbCRT . (scale * 
      {s0Int[Cie31xBar,380,780,0,380,i,0],
       s0Int[Cie31yBar,380,780,0,380,i,0],
       s0Int[Cie31zBar,380,780,0,380,i,0]})) /
    RgbEEValoisSLN[scale]];
  (* return point and color *)
  {gr, by, dl, RGBColor[r,g,b]}
  )] 

GrayValoisSLN1[i_:1,scale_:EEScaleValoisSLN1] :=
Module[ {gr,by,dl,r,g,b}, (
  (* color space coordinates *)
  {gr,by,dl} =
    {ToGRslN1[
       Gsig[scale s0Int[Glin[1,#]&,340,760,0,340,i,0]],
       Rsig[scale s0Int[Rlin[1,#]&,340,760,0,340,i,0]]],
     ToBYslN1[
       Bsig[scale s0Int[Blin[1,#]&,340,760,0,340,i,0]],
       Ysig[scale s0Int[Ylin[1,#]&,340,760,0,340,i,0]]],
     ToDLslN1[
       Dsig[scale s0Int[Dlin[1,#]&,340,760,0,340,i,0]],
       Lsig[scale s0Int[Llin[1,#]&,340,760,0,340,i,0]]]};
  (* approximate color *)
  {r,g,b} = Map[If[#<0,0,If[#>1,1,#]]&,
    (XyzToRgbCRT . (scale * 
      {s0Int[Cie31xBar,380,780,0,380,i,0],
       s0Int[Cie31yBar,380,780,0,380,i,0],
       s0Int[Cie31zBar,380,780,0,380,i,0]})) /
    RgbEEValoisSLN[scale]];
  (* return point and color *)
  {gr, by, dl, RGBColor[r,g,b]}
  )] 

GrayValoisSLN3[i_:1,scale_:EEScaleValoisSLN3] :=
Module[ {gr,by,dl,r,g,b}, (
  (* color space coordinates *)
  {gr,by,dl} =
    {ToGRslN3[
       Gsig[scale s0Int[Glin[1,#]&,340,760,0,340,i,0]],
       Rsig[scale s0Int[Rlin[1,#]&,340,760,0,340,i,0]]],
     ToBYslN3[
       Bsig[scale s0Int[Blin[1,#]&,340,760,0,340,i,0]],
       Ysig[scale s0Int[Ylin[1,#]&,340,760,0,340,i,0]]],
     ToDLslN3[
       Dsig[scale s0Int[Dlin[1,#]&,340,760,0,340,i,0]],
       Lsig[scale s0Int[Llin[1,#]&,340,760,0,340,i,0]]]};
  (* approximate color *)
  {r,g,b} = Map[If[#<0,0,If[#>1,1,#]]&,
    (XyzToRgbCRT . (scale * 
      {s0Int[Cie31xBar,380,780,0,380,i,0],
       s0Int[Cie31yBar,380,780,0,380,i,0],
       s0Int[Cie31zBar,380,780,0,380,i,0]})) /
    RgbEEValoisSLN[scale]];
  (* return point and color *)
  {gr, by, dl, RGBColor[r,g,b]}
  )] 


ocsPlotValoisSLN[scale_:EEScaleValoisSLN, pts_:30, gamma_:GammaSony] :=
ocsPlotValoisSLN[scale, pts, gamma] =
Module[{pl, colors, polygons, EERgb},
  Show[
    ParametricPlot3D[
      pt0ValoisSLN[gap,offset,scale,gamma],
      {gap,0,420}, {offset,340,760},
      PlotPoints -> {pts, pts},
      DisplayFunction->Identity],
    ParametricPlot3D[
      pt1ValoisSLN[gap,offset,scale,gamma],
      {gap,0,420}, {offset,340,760},
      PlotPoints -> {pts, pts},
      DisplayFunction->Identity],
    Axes->True, AxesLabel->{"GR", "BY", "Br"},
    Background->GrayLevel[0.6], Lighting->False,
    DisplayFunction->Identity]
  ]

(* plotting range *)

ocsRangeValoisSLN = {{-1,1},{-1,1},{0,2}}


(* ------------------------------------------------------------------------- *)

(* Linear transform to make gray axis perfectly vertical, i.e. equal to the
   Br axis. determination of the matrix to accomplish the linear transform
   relies on the fact that the color space fits exactly inside the cube
   with magnitude of dimensions {2,2,2}.  *)

ValoisSLNToValoisSLN2 = Inverse[
  {{2,0,#[[1]]}, {0,2,#[[2]]}, {0,0,#[[3]]}}&[GrayValoisSLN[1]] ]

(* gray axis coordinates in ValoisSLN2 *)

GrayValoisSLN2[i_:1,scale_:EEScaleValoisSLN] := Chop[
  Join[ValoisSLNToValoisSLN2 . Take[#,3], {#[[4]]}]&[GrayValoisSLN[i,scale]]
  ]

(* plotting range *)

ocsRangeValoisSLN2 = {{-0.4,0.6},{-0.6,0.4},{0,1}}

(* ------------------------------------------------------------------------- *)

(* ValoisSL space is like ValoisSLN, but not normalized to give a "square"
   color space. *)

(* determine scaling factor -- since this is used to scale linear functions
   only, can use it to scale the stimulus spectra rather than the product
   of stimulus and sensor sensitivity. The EE spectrum is left out because
   it is equal to 1 throughout the wavelength range. *)

EEScaleValoisSL =
  Abs[FindMinimum[-DLlin[1,l], {l,500,510}][[1]]] /
  NIntegrate[DLlin[1,l], {l,340,760}]

(* The color scaling (Von Kries adaptation) is done via the XyzToRgbCRT
   transform, since the stimulus spectra are known. *)

RgbEEValoisSL[scale_:EEScaleValoisSL] := RgbEEValoisSL[scale] =
XyzToRgbCRT . (scale {
  s0Int[Cie31xBar,380,780,0,380,1,0],
  s0Int[Cie31yBar,380,780,0,380,1,0],
  s0Int[Cie31zBar,380,780,0,380,1,0]})

(* points on the OCS solid surface *)

pt0ValoisSL[gap_,offset_,scale_:EEScaleValoisSL,gamma_:GammaSony] :=
Module[ {gr,by,dl,r,g,b}, (
  (* color space coordinates *)
  {gr,by,dl} = {
    GRsig[scale s0Int[GRlin[1,#]&,340,760,gap,offset,0,1]], 
    BYsig[scale s0Int[BYlin[1,#]&,340,760,gap,offset,0,1]], 
    DLsig[scale s0Int[DLlin[1,#]&,340,760,gap,offset,0,1]] }; 
  (* approximate color *)
  {r,g,b} = (#^(1/gamma))& /@ Map[If[#<0,0,If[#>1,1,#]]&,
    (XyzToRgbCRT . (scale * 
      {s0Int[Cie31xBar,380,780,gap,offset,0,1],
       s0Int[Cie31yBar,380,780,gap,offset,0,1],
       s0Int[Cie31zBar,380,780,gap,offset,0,1]})) /
    RgbEEValoisSL[scale]];
  (* return point and color *)
  {gr, by, dl, RGBColor[r,g,b]}
  )] 

pt1ValoisSL[gap_,offset_,scale_:EEScaleValoisSL,gamma_:GammaSony] :=
Module[ {gr,by,dl,r,g,b}, (
  (* color space coordinates *)
  {gr,by,dl} = {
    GRsig[scale s0Int[GRlin[1,#]&,340,760,gap,offset,1,0]], 
    BYsig[scale s0Int[BYlin[1,#]&,340,760,gap,offset,1,0]], 
    DLsig[scale s0Int[DLlin[1,#]&,340,760,gap,offset,1,0]] }; 
  (* approximate color *)
  {r,g,b} = (#^(1/gamma))& /@ Map[If[#<0,0,If[#>1,1,#]]&,
    (XyzToRgbCRT . (scale * 
      {s0Int[Cie31xBar,380,780,gap,offset,1,0],
       s0Int[Cie31yBar,380,780,gap,offset,1,0],
       s0Int[Cie31zBar,380,780,gap,offset,1,0]})) /
    RgbEEValoisSL[scale]];
  (* return point and color *)
  {gr, by, dl, RGBColor[r,g,b]}
  )] 

GrayValoisSL[i_:1,scale_:EEScaleValoisSLN] :=
Module[ {gr,by,dl,r,g,b}, (
  (* color space coordinates *)
  {gr,by,dl} = {
    GRsig[scale s0Int[GRlin[1,#]&,340,760,0,340,i,0]], 
    BYsig[scale s0Int[BYlin[1,#]&,340,760,0,340,i,0]], 
    DLsig[scale s0Int[DLlin[1,#]&,340,760,0,340,i,0]] }; 
  (* approximate color *)
  {r,g,b} = Map[If[#<0,0,If[#>1,1,#]]&,
    (XyzToRgbCRT . (scale * 
      {s0Int[Cie31xBar,380,780,0,380,i,0],
       s0Int[Cie31yBar,380,780,0,380,i,0],
       s0Int[Cie31zBar,380,780,0,380,i,0]})) /
    RgbEEValoisSL[scale]];
  (* return point and color *)
  {gr, by, dl, RGBColor[r,g,b]}
  )] 


ocsPlotValoisSL[scale_:EEScaleValoisSL, pts_:30, gamma_:GammaSony] :=
ocsPlotValoisSL[scale, pts, gamma] =
Module[{pl, colors, polygons, EERgb},
  Show[
    ParametricPlot3D[
      pt0ValoisSL[gap,offset,scale,gamma],
      {gap,0,420}, {offset,340,760},
      PlotPoints -> {pts, pts},
      DisplayFunction->Identity],
    ParametricPlot3D[
      pt1ValoisSL[gap,offset,scale,gamma],
      {gap,0,420}, {offset,340,760},
      PlotPoints -> {pts, pts},
      DisplayFunction->Identity],
    Axes->True, AxesLabel->{"GR", "BY", "Br"},
    Background->GrayLevel[0.6], Lighting->False,
    DisplayFunction->Identity]
  ]


End[]	     		(* private *)
EndPackage[]		(* package *)
