Understanding the MeshShading options

Question

The following code draws the intersection of a surface N: eqNcar=0 with three other surfaces defined by eq1=0, eq2=0, eq3=0.

I would like to shade with dots a specific region on $N$ namely the region defined by $\{eq1>0\}, \{eq2>0\} \cap \{eq3<0\}$. I would like to use the option MeshShading but I don't understand what tensor I should put. The documentation says an array of depth $k$ (here $k=3$ because I have $3$ mesh functions) but such specification raises an error (like "bad tensor dimension"). Anyway I am fighting with what to put in the structure of MeshShading.

Clear["Global`*"];

nn          = 3;
eps         = Table[Subscript[ep,i],{i,nn}];
xs          = Table[Subscript[x,i],{i,nn}];
alp         = 1/2;
bet         = 2;
Subscript[ep,2] = -alp Subscript[ep,1]/bet;
Subscript[ep,1] = 1;
Subscript[ep,3] = 0;

eqNcar      = -(1 + alp + bet)^(-1) - ((1 + alp + bet)^(-1) - Subscript[x, 1])^2 - ((1 + alp + bet)^(-1) - Subscript[x, 2])^2 + Subscript[x, 3];
eqNpar      = {w,s,w^2+s^2} + 1/(1+ alp+bet) {1,1,1};

eq1         = (-((1 + (1 + alp + bet)*s^2 + (1 + alp + bet)*w^2)*\
    (bet*s + s^2 + w*(alp + w))) + 2*s*(1 + (1 + alp + bet)*s)*\
   (s + alp*s^2 + w*(bet + alp*w)) + 2*w*(1 + (1 + alp + bet)*w)*\
   (alp*s + w + bet*(s^2 + w^2)))/(1 + alp + bet); 

eq2         = (alp*(-s + w)*((1 + alp + bet)^(-1) + s^2 + w^2) -
  (2*w*((1 + alp + bet)^(-1) + w)*(-(alp^2*((1 + alp + bet)^(-1) + s)) +
     bet*((1 + alp + bet)^(-1) + w)))/bet -
  (2*s*((1 + alp + bet)^(-1) + s)*(-(alp*((1 + alp + bet)^(-1) + s)) +
     bet^2*((1 + alp + bet)^(-1) + w)))/bet)*Subscript[ep, 1];

eq3 = (-3/(1 + alp + bet) + s + 2*w)*Subscript[ep, 1] + 
 ((-2 + alp + bet)*(3 + (1 + alp + bet)^2*s^2 + (1 + alp + bet)*s*
      (-3 + (1 + alp + bet)*w) + (1 + alp + bet)*w*
      (-3 + (1 + alp + bet)*w)) + (1 + alp + bet)*(-3 + 2*(1 + alp + bet)*s + 
     (1 + alp + bet)*w)*Subscript[ep, 2])/(1 + alp + bet)^2;

pN = ParametricPlot3D[eqNpar,{w,-1,1},{s,-1,1},\
  MeshFunctions -> {Function[{x,y,z,w,s},(eqNcar/. Thread[xs->eqNpar])-eq1],\
    Function[{x,y,z,w,s},(eqNcar/. Thread[xs->eqNpar])-eq2],\
    Function[{x,y,z,w,s},(eqNcar/. Thread[xs->eqNpar])-eq3]},\
  MeshStyle -> {{Blue},{Red},{Black}},\
  MeshShading->{
    {{Automatic,Automatic,Automatic}},\
    {{Automatic,Automatic,Automatic}},\
    {{HalftoneShading[],Automatic,Automatic}}
  },
  Mesh -> {{0}},PlotStyle->Orange,PlotPoints->200,
  BoundaryStyle->{1->None},\
  AxesLabel->xs
 ];

cm              = 72/2.54;
image       = Rasterize[Show[pN,ImageSize->10 cm],"Image",ImageResolution->800];
Export["ex.pdf", Show[image, ImageSize -> 10 cm]]

In my example, MeshShading is an array $M$ of size $(1,3,3)$ and of depth $4$ (note the documentation would expect a array of depth $3$) and I was expecting the coefficient $m_{1ij}$ to be the color of the region between the ith meshcurves and jth mesh curve but it is not the case.

ps: I would like the method to be easily adaptable for a number of mesh functions up to $10$.

Answer 1

I think you are misunderstanding how MeshShading is supposed to work. The way it works is as follows: You give a set of MeshFunctions $\{f_1, f_2, \dots\}$, and a set of mesh divisions $\{\{d_{1,1},d_{1,2},\dots\},\{d_{2,1},d_{2,2},\dots\},\dots\}$ (effectively the setting of Mesh, i.e. the values of $f_1,\dots$ that where mesh contours should be drawn). The mesh divisions effectively define intervals as $$\{\{I_{1,1},I_{1,2},\dots\},\{I_{2,1},I_{2,2},\dots\},\dots\}=\{\{(-\infty,c_{1,1}],(c_{1,1},c_{1,2}],\dots\},\{(-\infty,c_{2,1}],(c_{2,1},c_{2,2}],\dots\},\dots\}.$$

This means that every point can be assigned to a set of intervals $\{I_{1,i},I_{2,j},\dots\}$, whose style is given by the entry $m_{i,j,\dots}$ from the MeshShading option.

To illustrate this, consider the following:

With[{plot = 
   ParametricPlot3D[eqNpar, {w, -1, 1}, {s, -1, 1}, 
    MeshFunctions -> {Function[{x, y, z, w, 
        s}, (eqNcar /. Thread[xs -> eqNpar]) - eq1], 
      Function[{x, y, z, w, s}, (eqNcar /. Thread[xs -> eqNpar]) - 
        eq2], Function[{x, y, z, w, 
        s}, (eqNcar /. Thread[xs -> eqNpar]) - eq3]}, 
    MeshStyle -> {{Blue}, {Red}, {Black}}, 
    MeshShading -> 
     SparseArray[{{i_, j_, k_} :> style[i, j, k]}, {2, 2, 2}, 
      Directive[]], Mesh -> {{0}}, PlotStyle -> Orange, 
    PlotPoints -> 200, BoundaryStyle -> {1 -> None}, AxesLabel -> xs]},
 Manipulate[
  plot /. 
   style[i_, j_, k_] :> 
    If[{i, j, k} == {i0, j0, k0}, Green, Directive[]],
  {i0, 1, 2, 1},
  {j0, 1, 2, 1},
  {k0, 1, 2, 1}
  ]
 ]

As you can see, the first slider selects which "side" of the black mesh line you are on, the second which side of the red mesh line, and the third slider which side of the blue mesh line you are on.