%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%
%% Observer and Sun positions
%% view is also defined elsewhere
%% The longitude is computed in the Ozx plane, beginning at the Oz axis.
\tikzset{%
  view/.style 2 args={%  observer longitude and latitude (y upwards)
                      %  Remark. lomg=0 means x=0
    z={({-sin(#1)}, {-cos(#1)*sin(#2)})},
    x={({cos(#1)}, {-sin(#1)*sin(#2)})},
    y={(0, {cos(#2)})},
    evaluate={%
      \toxx={sin(#1)*cos(#2)};
      \toyy={sin(#2)};
      \tozz={cos(#1)*cos(#2)};
    },
    longitude = #1,
    latitude = #2
  },
  sun/.style n args={3}{% longitude, latitude, light contrast in [0, 1]
    evaluate={%
      real \sunxx, \sunyy, \sunzz;
      \sunxx = sin(#1)*cos(#2);
      \sunyy = sin(#2);
      \sunzz = cos(#1)*cos(#2);
    },
    sun longitude = #1,
    sun latitude = #2,
    contrast = #3
  }
}
\pgfkeys{/tikz/.cd,
  latitude/.store in=\toLatit,  % observer's latitude
  latitude=0
}
\pgfkeys{/tikz/.cd,
  longitude/.store in=\toLongit,  % observer's longitude
  longitude=0  % corresponds to x=0
}
\pgfkeys{/tikz/.cd,
  sun latitude/.store in=\sunLatit,
  sun latitude = 40
}
\pgfkeys{/tikz/.cd,
  sun longitude/.store in=\sunLongit, 
  sun longitude=-60  % on the left of the Oz axis
}
\pgfkeys{/tikz/.cd,
  contrast/.store in=\lightCtr,  % light contrast \in [0, 1]
  % It is introduced (modified) with the style sun.  The default value
  % allows the styles to know if the sun is on or off.
  % It is used under sor in pic objects with su and s variants.
  contrast=-1
}
\pgfkeys{/tikz/.cd,
  % Value of the inner product of the normal vector with the sun
  % direction for which the color is the local color.
  local color threshold/.store in=\lColorThreshold,
  local color threshold=.2
}
\pgfkeys{/tikz/.cd,
  % max value of the constant c in Color!c!LocalColor
  max color cst/.store in=\maxColorCst,
  max color cst=.99
}
\pgfkeys{/tikz/.cd,  % multiplicative coefficients used in ?? to be seen
  seen key/.store in=\seenKey,
  seen key=1
}
\pgfkeys{/tikz/.cd,
  unseen key/.store in=\unseenKey,
  unseen key=.2
}
\pgfkeys{/tikz/.cd,
  near zero/.store in=\nearZero,
  near zero=.01
}


%%%% bounds of the domain
\pgfkeys{/tikz/.cd,
  initialT/.store in=\initialT,
  initialT=0
}
\pgfkeys{/tikz/.cd,
  finalT/.store in=\finalT,
  finalT=1
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% general functions
\newcommand{\setColorsForRSurface}{%
  \colorlet{rgbLeft{0}}{K}  % not sun
  \colorlet{rgbLeft{1}}{W}  % sun
  \colorlet{rgbRight}{.}
}


\tikzmath{%
  function norm3d(\vectorx, \vrctory, \vectorz) {%
    real \res;
    \res = \vectorx*\vectorx +\vrctory*\vrctory +\vectorz*\vectorz;
    return {pow(\res, .5)};
  };
  function ipFaceAndTo(\nox, \noy, \noz) {%
    return {sign(\nFacex*\toxx +\nFacey*\toyy +\nFacez*\tozz)};
  };
  function whatLayer(\ip, \l1, \l2) {%
    if \ip>=0 then {% if seen
      return \l1;
    } else {%
      return \l2;
    };
  };
  %
  %% The next functions verify if the face is seen and then return a
  % color "number" depending on the sun position.  It is used in
  % rgbConstant.
  function faceColorCstI(\t) {%  I = increasing part [\lColorThreshold, 1]
    \tmpFCCst = \lightCtr* \maxColorCst;
    return \tmpFCCst/(pow(1 -\lColorThreshold, 2))*pow(\t -\lColorThreshold, 2);
  };
  function faceColorCstD(\t) {%  D = decreasing part [-1, \lColorThreshold)
    return -\lightCtr* \maxColorCst/(1 +\lColorThreshold)*(\t -\lColorThreshold);
  };
  function rgbConstant(\ipo, \ips) {%
    % \ipo = inner product with observer's direction
    % \ips =  inner product with sun's direction
    \ipsp = abs(\ips);
    if \ipo*\ips >= 0 then {% if seen in sun or unseen and not in sun
      if \ipsp >= \lColorThreshold then {%
        \tmprgbcst = faceColorCstI(\ipsp);
      } else {%
        \tmprgbcst = faceColorCstD(\ipsp);
      };
    } else {%
      \tmprgbcst = faceColorCstD(-\ipsp);
    };
    return min(int(100*\tmprgbcst), 100);
  };
  function rgbLeftIndex(\ipo, \ips) {%
    % \ipo = inner product with observer's direction
    % \ips =  inner product with sun's direction
    \ipsp = abs(\ips);
    \t = 0;
    if \ipo*\ips >= 0 then {% if seen in sun or unseen and not in sun
      if \ipsp >= \lColorThreshold then {%
        \t = 1;
      };
    };
    return \t;
  };
}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% layers and show coordinates
\newcommand{\setlayers}[1]{% layers
  \foreach \k in {#1}{%
    \pgfdeclarelayer{\k}
  }
  \pgfsetlayers{#1}
}


\tikzset{%
  pics/showCoordinates/.style n args={3}{%
    %  #1 = length
    %  #2 = fg layer
    %  #3 = bg layer
    code={%
      \tikzmath{%
        integer \xyLayer, \zLayer;
        real \aLen, \aArrow, \aR;
        \aLen = #1;
        \aArrow = #1/6;
        \aR = #1/3;
        \xyLayer = #3;
        \zLayer = #2;
        if \toLongit>90 then {%
          \xyLayer = #2;
          \zLayer = #3;
        };
        if \toLongit<-90 then {%
          \xyLayer = #2;
          \zLayer = #3;
        };        
        {%
          \begin{pgfonlayer}{\xyLayer}
            \begin{scope}[canvas is xy plane at z=0]
              \draw[B, fill=B!60] (0, \aArrow/4) -- ++(\aLen, 0)
              -- ++(0, 3*\aArrow/4)
              -- ++(\aArrow, -\aArrow) -- ++(-\aArrow, -\aArrow)
              -- ++(0, 3*\aArrow/4) -- ++(-\aLen, 0) -- cycle;
              \draw[B, fill=B!60, rotate=90] (0, \aArrow/4) -- ++(\aLen, 0)
              -- ++(0, 3*\aArrow/4)
              -- ++(\aArrow, -\aArrow) -- ++(-\aArrow, -\aArrow)
              -- ++(0, 3*\aArrow/4) -- ++(-\aLen, 0) -- cycle;
              \draw[B, fill=B!60, even odd rule]
              (0, 0) circle (\aR)
              (0, 0) circle (.8*\aR);
            \end{scope}
          \end{pgfonlayer}
          \begin{pgfonlayer}{\zLayer}
            \begin{scope}[canvas is zx plane at y=0]
              \draw[W, line width=4pt] (0, \aArrow/4) -- ++(\aLen, 0)
              -- ++(0, 3*\aArrow/4)
              -- ++(\aArrow, -\aArrow) -- ++(-\aArrow, -\aArrow)
              -- ++(0, 3*\aArrow/4) -- ++(-\aLen, 0);
              \draw[O, fill=O] (0, \aArrow/4) -- ++(\aLen, 0)
              -- ++(0, 3*\aArrow/4)
              -- ++(\aArrow, -\aArrow) -- ++(-\aArrow, -\aArrow)
              -- ++(0, 3*\aArrow/4) -- ++(-\aLen, 0) -- cycle;
            \end{scope}
          \end{pgfonlayer}
        };
      }
    }
  }
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Code elements
\tikzset{%
  new basis/.code={%
    \tikzmath{%
      \ux = \cay *\caz;
      \uy = \saz;
      \uz = -\say *\caz;
      \vx = -\cay *\saz *\cax +\say *\sax;
      \vy = \caz *\cax;
      \vz = \say *\saz *\cax +\cay *\sax;
      \wx = \cay *\saz *\sax +\say *\cax;
      \wy = -\caz *\sax;
      \wz = -\say *\saz *\sax +\cay *\cax;
    }
  },
  the points/.code={%  It takes one argument, the name of the curve. 
  % It constructs the points \P{\k,\j} with \k defining the meridian
  % and \j the parallel.  There are \NParallels+1 parallels equaly
  % spaced from \initialT to \finalT.  There are \NMeridians (the one
  % with index 0 coincides with the last one, with index \NMeridians).
    \tikzmath{%
      for \j in {0, ..., \NParallels}{% the plane (Ozx) |-> 0
        \s = \initialT +\j/\NParallels*(\finalT -\initialT);
        \yC0 = y#1(\s);
        \zC0 = z#1(\s);
        for \k in {0, ..., \NMeridians}{% (Oz) |-> 0
          \t = \k/\NMeridians*360;
          \wTmp = cos(\t)*\zC0;
          \vTmp = \yC0;
          \uTmp = sin(\t)*\zC0;
          \Px{\k,\j} = \wTmp*\wx +\vTmp*\vx +\uTmp*\ux;
          \Py{\k,\j} = \wTmp*\wy +\vTmp*\vy +\uTmp*\uy;
          \Pz{\k,\j} = \wTmp*\wz +\vTmp*\vz +\uTmp*\uz;
        };
      };
    }
  },
  face normal vector and ips/.code={%
  % It depends on the value of \lightCtr.  If >=0, i.e. if the light
  % source is given, then the inner product with the light direction
  % is computed. 
  % 
  % Elements for each face  \ax, ...
  % \i for meridians (around Oy)
  % \j parallels (from -pi/2 to pi/2 for example)
  % \j = 1 is special
  % 
  % n unit normal vector and the inner products with observer and
  % sun needed for faceColor.
  % 
  % The (a, b, n) vectors  form a direct basis and n points
  % towards the "exterior" of the 3d object to be drawn.
    \tikzmath{%
      integer \i, \tmpInt;
      for \j in {1, ..., \NParallels}{%
        \prevj = \j -1;
        if \j==1 then {%
          \i = \j;
        } else {%
          \i = \prevj;
        };
        for \k in {0, ..., \NMeridians}{%
          \prevk = int(Mod(\k -1, \NMeridians));
          \ax = \Px{\k,\i} -\Px{\prevk,\prevj};
          \ay = \Py{\k,\i} -\Py{\prevk,\prevj};
          \az = \Pz{\k,\i} -\Pz{\prevk,\prevj};
          \aTmpNorm = norm3d(\ax, \ay, \az);
          \ax = \ax/\aTmpNorm;
          \ay = \ay/\aTmpNorm;
          \az = \az/\aTmpNorm;
          \bx = \Px{\prevk,\j} -\Px{\prevk,\prevj};
          \by = \Py{\prevk,\j} -\Py{\prevk,\prevj};
          \bz = \Pz{\prevk,\j} -\Pz{\prevk,\prevj};
          \bTmpNorm = norm3d(\bx, \by, \bz);
          \bx = \bx/\bTmpNorm;
          \by = \by/\bTmpNorm;
          \bz = \bz/\bTmpNorm;        
          \nFacex = \ay*\bz -\az*\by;
          \nFacey = -\ax*\bz +\az*\bx;
          \nFacez = \ax*\by -\ay*\bx;
          \tmp = norm3d(\nFacex, \nFacey, \nFacez);
          \nFacex = \nFacex/\tmp;
          \nFacey = \nFacey/\tmp;
          \nFacez = \nFacez/\tmp; % the normal vector is normalized
          \ipFTo{\k,\j} = \nFacex*\toxx +\nFacey*\toyy +\nFacez*\tozz;
          if \lightCtr>=0 then {%
            \ipFSun{\k,\j} = \nFacex*\sunxx +\nFacey*\sunyy +\nFacez*\sunzz;
          };
        };
      };
    }
  },
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Main style and pic objects
\tikzset{%
  sor/.style args={name=#1,
    meridians=#2, parallels=#3, t_ini=#4, t_end=#5, rAxis=#6}{%
    % It provides the elements for drawing the quadrilaterals
    % forming a surface of revolution, i.e. the points for the
    % rotating surface obtained by rotating the curve given by #1
    % around the Oy axis. 
    % There are #2 meridians and #3+1 parallels. The variables 
    % \NMeridians, \NParallels, ... act as global  variables.  
    % 
    % #6 is the position of the rotational axis.  It is the image of
    % Oy through a rotation of #6_1 around Ox followed by #6_3 around
    % Oz and by #6_2 around Oy. 
    evaluate={%
      real \a, \s, \t, \ipFTo, \ipFTo;
      real \cax, \sax, \cay, \say, \caz, \saz;
      real \ux, \uy, \uz, \vx, \vy, \vz, \wx, \wy, \wz;
      integer \NMeridians, \NParallels;
      integer \j, \k, \i, \prevj, \prevk; % \prev[] needed in codes 
      \NMeridians = #2;
      \NParallels = #3;
      \initialT = #4;
      \finalT = #5;
      \i = 0; 
      for \tmp in {#6}{%  angles that characterise the axis of
                     %  rotation's position; in fact the new basis
                     %  obtained by the compositon R_\ay R_\az R_\ax.
        \i = \i +1;
        \aVar\i = \tmp;
      };
      \cax = cos(\aVar1);
      \sax = sin(\aVar1);
      \cay = cos(\aVar2);
      \say = sin(\aVar2);
      \caz = cos(\aVar3);
      \saz = sin(\aVar3);
    },
    new basis,
    the points={#1},
    face normal vector and ips
  },
}


\tikzset{%
  pics/meridian/.style args={index=#1, layer=#2}{%
    %  index is in 0, ..., \NMeridians -1
    code={%
      \tikzmath{%
        integer \j;
        {%
          \begin{pgfonlayer}{#2}
            \draw[color=.] (\Px{#1,0}, \Py{#1,0}, \Pz{#1,0})
            \foreach \j in {1, ..., \NParallels}{%
              -- (\Px{#1,\j}, \Py{#1,\j}, \Pz{#1,\j})
            };
          \end{pgfonlayer}
        };
      }
    }
  },  
  pics/parallel/.style args={index=#1, layer=#2}{%
    %  index is in 0, ..., \NParallels
    code={%
      \tikzmath{%
        integer \k;
        {%
          \begin{pgfonlayer}{#2}
            \draw[color=.] (\Px{0,#1}, \Py{0,#1}, \Pz{0,#1})
            \foreach \k in {1, ..., \NMeridians}{%
              -- (\Px{\k,#1}, \Py{\k,#1}, \Pz{\k,#1})
            };
          \end{pgfonlayer}
        };
      }
    }
  },
  pics/meridians/.style args={initial=#1, multiple=#2, layer=#3}{%
    % displays meridians #1 + \i \cdot #2
    code={%
      \tikzmath{%
        integer \i, \j, \k;
        \k = #1;
        for \i in {0, ..., int(\NMeridians/#2) -1}{%
          {%
            \begin{pgfonlayer}{#3}
              \draw[color=.] (\Px{\k,0}, \Py{\k,0}, \Pz{\k,0})
              \foreach \j in {1, ..., \NParallels}{%
                -- (\Px{\k,\j}, \Py{\k,\j}, \Pz{\k,\j})
              };
            \end{pgfonlayer}
          };
          \k = \k + #2;
        };
      }
    }
  },
  pics/parallels/.style args={initial=#1, multiple=#2, layer=#3}{%
    % displays parallels #1 + \i \cdot #2
    code={%
      \tikzmath{%
        integer \i, \j, \k, \bound;
        \j = #1;
        if \j == 0 then {%
          \bound = int(\NParallels/#2);
        } else {%
          \bound = int(\NParallels/#2) -1;
        };
        for \i in {0, ..., \bound}{%
          {%
            \begin{pgfonlayer}{#3}
              \draw[color=.] (\Px{0,\j}, \Py{0,\j}, \Pz{0,\j})
              \foreach \k in {1, ..., \NMeridians}{%
                -- (\Px{\k,\j}, \Py{\k,\j}, \Pz{\k,\j})
              };
            \end{pgfonlayer}
          };
          \j = \j + #2;
        };
      }
    }
  }
}


\tikzset{%
  pics/rSurface/.style args={above=#1, below=#2, actions=#3}{%
    % #1 = fg layer
    % #2 = bg layer
    % #3 = actions 
    code={%
      \tikzmath{%  
        integer \j, \k, \prevk, \prevj, \layerI, \rgbCst;
        %% faces
        for \j in {1, ..., \NParallels}{%
          for \k in {1, ..., \NMeridians}{%
            \prevk = \k -1;
            \prevj = \j -1;
            \layerI = whatLayer(\ipFTo{\k,\j}, #1, #2);
            {%
              \begin{pgfonlayer}{\layerI}
                \path[., #3]
                (\Px{\prevk,\prevj}, \Py{\prevk,\prevj}, \Pz{\prevk,\prevj})
                -- (\Px{\prevk,\j}, \Py{\prevk,\j}, \Pz{\prevk,\j})
                -- (\Px{\k,\j}, \Py{\k,\j}, \Pz{\k,\j})
                -- (\Px{\k,\prevj}, \Py{\k,\prevj}, \Pz{\k,\prevj})
                -- cycle;
              \end{pgfonlayer}
            };
          };
        };
      }
    }
  },
  pics/rSurface L/.style args={above=#1, below=#2, actions=#3}{%
    % #1 = fg layer
    % #2 = bg layer
    % #3 = actions (drawn must be provided if full opacity
    code={% !!
      \setColorsForRSurface
      \tikzmath{%  
        integer \j, \k, \prevk, \prevj, \layerI, \rgbCst, \li;
        real \ipToTmp, \ipSunTmp, \ipSunTmpP;
        %% faces
        for \k in {1, ..., \NMeridians}{% 
          for \j in {1, ..., \NParallels}{%
            \prevk = \k -1;
            \prevj = \j -1;
            \layerI = whatLayer(\ipFTo{\k,\j}, #1, #2);
            \li = rgbLeftIndex(\ipFTo{\k,\j}, \ipFSun{\k,\j});
            \rgbCst = rgbConstant(\ipFTo{\k,\j}, \ipFSun{\k,\j});
            {%
              \begin{pgfonlayer}{\layerI}
                \fill[rgbLeft{\li}!\rgbCst!rgbRight,
                line width=.0001pt, #3]
                (\Px{\prevk,\prevj}, \Py{\prevk,\prevj}, \Pz{\prevk,\prevj})
                -- (\Px{\k,\prevj}, \Py{\k,\prevj}, \Pz{\k,\prevj})
                -- (\Px{\k,\j}, \Py{\k,\j}, \Pz{\k,\j})
                -- (\Px{\prevk,\j}, \Py{\prevk,\j}, \Pz{\prevk,\j})
                -- cycle;
              \end{pgfonlayer}
            };
          };
        };
      }
    }
  },
  pics/partial rSurface/.style args={above=#1, below=#2, actions=#3,
    i from=#4, i to=#5}{%
    % #1 = fg layer
    % #2 = bg layer
    % #3 = actions 
    % #4 = index of the first meridian
    % #5 = index of the last meridian 
    code={%
      \tikzmath{%  
        integer \j, \k, \prevk, \prevj, \layerI, \rgbCst, \kini, \kend;
        if #4 < #5 then {%
          \kini = int(#4 +1);
          \kend = #5;
        } else {%
          \kini = #4;
          \kend = int(#5 +1);          
        };
        %% faces
        for \k in {\kini, ..., \kend}{% 
          for \j in {1, ..., \NParallels}{%
            \prevk = int(Mod(\k -1, \NMeridians));
            \k = int(Mod(\k, \NMeridians));  % when the index is
                                             % bigger than the number of meridians
            if \k==0 then {%
              \k = \NMeridians;
            };
            \prevj = \j -1;
            \layerI = whatLayer(\ipFTo{\k,\j}, #1, #2);
            {%
              \begin{pgfonlayer}{\layerI}
                \path[., #3]
                (\Px{\prevk,\prevj}, \Py{\prevk,\prevj}, \Pz{\prevk,\prevj})
                -- (\Px{\prevk,\j}, \Py{\prevk,\j}, \Pz{\prevk,\j})
                -- (\Px{\k,\j}, \Py{\k,\j}, \Pz{\k,\j})
                -- (\Px{\k,\prevj}, \Py{\k,\prevj}, \Pz{\k,\prevj})
                -- cycle;
              \end{pgfonlayer}
            };
          };
        };
      }
    }
  },
  pics/partial rSurface L/.style args={above=#1, below=#2, actions=#3,
    i from=#4, i to=#5}{%
    % #1 = fg layer
    % #2 = bg layer
    % #3 = actions 
    % #4 = index of the first meridian
    % #5 = index of the last meridian 
    code={%
      \setColorsForRSurface
      \tikzmath{%  
        integer \j, \k, \kini, \kend, \prevk, \prevj, \layerI, \rgbCst;
        %% faces
        if #4 < #5 then {%
          \kini = int(#4 +1);
          \kend = #5;
        } else {%
          \kini = #4;
          \kend = int(#5 +1);          
        };
        for \k in {\kini, ..., \kend}{% 
          for \j in {1, ..., \NParallels}{%
            \prevk = int(Mod(\k -1, \NMeridians));
            \k = int(Mod(\k, \NMeridians));  % when the index is
                                % bigger than the number of meridians
            \prevj = \j -1;
            \layerI = whatLayer(\ipFTo{\k,\j}, #1, #2);
            \li = rgbLeftIndex(\ipFTo{\k,\j}, \ipFSun{\k,\j});
            \rgbCst = rgbConstant(\ipFTo{\k,\j}, \ipFSun{\k,\j});
            {%
              \begin{pgfonlayer}{\layerI}
                \fill[rgbLeft{\li}!\rgbCst!rgbRight,
                line width=.0001pt, #3]
                (\Px{\prevk,\prevj}, \Py{\prevk,\prevj}, \Pz{\prevk,\prevj})
                -- (\Px{\prevk,\j}, \Py{\prevk,\j}, \Pz{\prevk,\j})
                -- (\Px{\k,\j}, \Py{\k,\j}, \Pz{\k,\j})
                -- (\Px{\k,\prevj}, \Py{\k,\prevj}, \Pz{\k,\prevj})
                -- cycle;
              \end{pgfonlayer}
              % \path (-5 +.7*\k, \j +\k/3, 0)
              % node[scale=.8] {\k, \j, \ipFTo{\k,\j}, \rgbCst, \li}; 
            };
          };
        };
      }
    }
  },
  pics/rSurface x/.style args={above=#1, below=#2, actions=#3, factor=#4}{%
    % In this version, only the borders of certain quadrilaterals are
    % drawn.  The choice is done through #4 which is a common factor of
    % \NMeridians and \NParallels.
    %
    % Note that the first parallel (for \j = 0) must be drawn and it
    % is nor obtainable using the previous structure (from rSurface).
    %
    % #1 = fg layer
    % #2 = bg layer
    % #3 = actions
    % #4 = factor
    code={%
      \tikzmath{%  
        integer \j, \k, \prevk, \prevj, \layerI, \rgbCst;
        %% faces
        for \k in {1, ..., \NMeridians}{% 
          for \j in {1, ..., \NParallels}{%
            \prevk = int(Mod(\k -1, \NMeridians));
            \prevj = \j -1;
            \layerI = whatLayer(\ipFTo{\k,\j}, #1, #2);
            {%
              \begin{pgfonlayer}{\layerI}
                \path[#3]
                (\Px{\prevk,\prevj}, \Py{\prevk,\prevj}, \Pz{\prevk,\prevj})
                -- (\Px{\prevk,\j}, \Py{\prevk,\j}, \Pz{\prevk,\j})
                -- (\Px{\k,\j}, \Py{\k,\j}, \Pz{\k,\j})
                -- (\Px{\k,\prevj}, \Py{\k,\prevj}, \Pz{\k,\prevj})
                -- cycle;
              \end{pgfonlayer}
            };
            if \k==1 then {%
              {%
                \begin{pgfonlayer}{\layerI}
                  \draw[.]
                  (\Px{0,\j}, \Py{0,\j}, \Pz{0,\j})
                  -- (\Px{0,\prevj}, \Py{0,\prevj}, \Pz{0,\prevj});
                \end{pgfonlayer}
              };
            };
            if int(Mod(\k, #4))==0 then {%
              {%
                \begin{pgfonlayer}{\layerI}
                  \draw[.]
                  (\Px{\k,\j}, \Py{\k,\j}, \Pz{\k,\j})
                  -- (\Px{\k,\prevj}, \Py{\k,\prevj}, \Pz{\k,\prevj});
                \end{pgfonlayer}
              };
            };
            if \j==1 then {%
              {%
                \begin{pgfonlayer}{\layerI}
                  \draw[.]
                  (\Px{\prevk,0}, \Py{\prevk,0}, \Pz{\prevk,0})
                  -- (\Px{\k,0}, \Py{\k,0}, \Pz{\k,0});
                \end{pgfonlayer}
              };
            };      
            if int(Mod(\j, #4))==0 then {%
              {%
                \begin{pgfonlayer}{\layerI}
                  \draw[.]
                  (\Px{\prevk,\j}, \Py{\prevk,\j}, \Pz{\prevk,\j})
                  -- (\Px{\k,\j}, \Py{\k,\j}, \Pz{\k,\j});
                \end{pgfonlayer}
              };
            };      
          };
        };
      }
    }
  },
}