BeginPackage["FacetSurface`"];
Clear["FacetSurface`*"];
 
Options[VisualizeFacetSurface]={Ticks      ->Automatic,
                                Axes       ->True, 
                                ImageSize  ->Large,
                                PlotLabel  ->None,
                                PlotRange   ->All};
VisualizeFacetSurface::usage =
"VisualizeFacetSurface[data, (options...}]\t creates a triangulated faceted surface from a rectangular array of data values";
TriangulatedInterpolation::usage =
"TriangulatedInterpolation[{x,y},data]\t interpolates data grid in a triangulated faceted fashion consistent with VisualizeFacetSurface[]\n"<>
"\t\t\t\t\t\t\tRange of values for {x,y} is: 1\[LessEqual]x\[LessEqual](#columns in data),  1\[LessEqual]y\[LessEqual](#rows in data)";  

Options[FacetSurfaceGridData]={Global`IC          ->Automatic, 
                               Global`ICDim       ->{21,21},
                               Global`SemiRandomIC->True,
                               RandomSeed         ->5};
FacetSurfaceGridData::usage =
"FacetSurfaceGridData[rule,(nsteps, ...options)]      generates the grid data for the triangular faceted-surface model.\n"<>
"Options:\t\tIC\t\t->Automatic  uses combined random/periodic initial condition."<>
"(IC->2d array will overide every other initial condition option.)\n"<>
"\t\t\tSemiRandomIC-> False    (False will create purely random IC)\n"<>
"\t\t\tICDim\t\t-> {21,21}    (dimensions for automatically generated IC {ny,nx})\n"<>
"\t\t\tSeedRandom\t-> 5       (random-number generator seed)."

Begin["`Private`"]

(*===================== FacetSurface routines =================================*)
(*                                                                             *)
(*- June/2008 NKS Summer School 2008 project nugget/Vallorie J. Peridier       *)
(*                                                                             *)
(*- Objective:  Simulate/visualize irregular surfaces using cellular automata  *)
(*=============================================================================*)   


(*==============================================================================
  | Visualize the data[[y-indices, x-indices]] as triangular facets            |
  ==============================================================================*)
VisualizeFacetSurface[data_, opts___Rule]:= Module[{nx, ny, testdata}, 
              {ny,nx}=Dimensions[data];
              testdata=Table[{x,y, data[[y,x]]},{y,1,ny},{x,1,nx}];
              Graphics3D[Flatten[Join[
                   Table[{Polygon[Extract[testdata,{{jj,ii},{jj,ii+1},{jj+1,ii+1}}]],
                          Polygon[Extract[testdata,{{jj,ii},{jj+1,ii},{jj+1,ii+1}}]]},
                        {jj,1,ny-1,1}, {ii,Mod[jj,2,1],nx-1,2}],
           
                   Table[{Polygon[Extract[testdata,{{jj,ii},{jj,ii+1},{jj+1,ii}}]],
                          Polygon[Extract[testdata,{{jj+1,ii},{jj+1,ii+1},{jj,ii+1}}]]},
                       {ii,1,nx-1,1}, {jj,Mod[ii-1,2,1],ny-1,2}]
                    ]],
                     Ticks    -> (Ticks     /.{opts}/.Options[VisualizeFacetSurface]), 
                     Axes     -> (Axes      /.{opts}/.Options[VisualizeFacetSurface]),
                     PlotLabel-> (PlotLabel /.{opts}/.Options[VisualizeFacetSurface]),
                     PlotRange-> (PlotRange /.{opts}/.Options[VisualizeFacetSurface]),
                     ImageSize-> (ImageSize /.{opts}/.Options[VisualizeFacetSurface])
                     ]];                 


(*==============================================================================
  | Data interpolation consistent with the FacetSurface visualization          | 
  ==============================================================================*)
TriangulatedInterpolation[{x_,y_},data_]:=
      If[!MatrixQ[data],
            Print["interpolation data is not a matrix"];                Abort[],                                           
            Print["field point ",{x,y}, " out of data range: [1<=x<=",
                   #[[2]],"], [1<=y<=",#[[1]],"]"]&[Dimensions[data]];  Abort[]
         ]/;(!(MatrixQ[data]&&InRange[x,y,data]));
  
TriangulatedInterpolation[{x_,y_},data_]:=data[[y,x]]/;MatrixQ[data]&& 
                                                       InRange[x,y,data]&&
                                                       IntegerQ[x]&&
                                                       IntegerQ[y];
TriangulatedInterpolation[{x_,y_},data_]:=Module[
         {xf=FractionalPart[x],ii=Floor[x], jj=y,  v1,v2, nx, ny},
         {ny,nx} = Dimensions[data];
         {v1,v2}=Extract[data,{{jj,ii},{jj,Min[ii+1,nx]}}];
         ((v1-v2)(1-xf) +  v2)//N]/;MatrixQ[data]&&
                                    InRange[x,y,data]&&
                                    IntegerQ[y];

TriangulatedInterpolation[{x_,y_},data_]:=Module[
         {ii=x, yf=FractionalPart[y],jj=Floor[y],v1,v4, nx, ny},
         {ny,nx} = Dimensions[data];
         {v1,v4}=Extract[data,{{jj,ii},{Min[ny,jj+1],ii}}];
         ((v1-v4)(1-yf)  + v4)//N]/;MatrixQ[data]&&
                                    InRange[x,y,data]&&
                                    IntegerQ[x];

TriangulatedInterpolation[{x_,y_},data_]:=Module[
        {xf=FractionalPart[x],ii=Floor[x], 
         yf=FractionalPart[y],jj=Floor[y],nx,ny,v1,v2,v3,v4},

        {ny,nx}=Dimensions[data];
        If[ii==nx && jj==ny,Return[TriangulatedInterpolation[{ii,jj},data]]];
        If[ii==nx,          Return[TriangulatedInterpolation[{ii,y },data]]];
        If[          jj==ny,Return[TriangulatedInterpolation[{ x,jj},data]]];

        {v1,v2,v3,v4}=Extract[data,{{jj,ii},{jj,ii+1},{jj+1,ii+1},{jj+1,ii}}];
        If[EvenQ[ii+jj], 
                Return[EvenInterpolation[xf,yf,{v1,v2,v3,v4}]],
                Return[OddInterpolation[xf,yf,{v1,v2,v3,v4}]]]]/;MatrixQ[data]&&InRange[x,y,data];


EvenInterpolation[xf_,yf_,{v1_,v2_,v3_,v4_}]:=
         If[(xf==yf), xf v3 +(1-xf)v1,
                      If[yf>xf,  (v1-v4)(1-yf) + (v3-v4)(xf) + v4,
                                 (v1-v2)(1-xf) + (v3-v2)(yf) + v2]]//N;
OddInterpolation[xf_,yf_,{v1_,v2_,v3_,v4_}]:=
        If[(xf+yf==1), xf*v2+(1-xf)v4,
                       If[yf>(1-xf), (v2-v3)(1-yf) + (v4-v3)(1-xf)+v3,
                                     (v4-v1)(yf) + (v2-v1)(xf) + v1]]//N;
InRange[x_,y_,data_]:=If[x<1 || y<1 || x>Dimensions[data][[2]] || y>Length[data], False,True];


 
(*==========================================================
  | The irregular-surface grid-data generator              |
  ==========================================================*)

FacetSurfaceGridData[rn_Integer, nsteps_Integer:5, opts___Rule]:=Module[
     {ic,cum},
     
     (*--- Automatic IC? -----------------------------------------------------*)
     ic= (Global`IC     /. {opts} /. Options[FacetSurfaceGridData]); 
     ic=If[  ic===Automatic, FacetSurfaceGridAutomaticIC[opts],     ic];

     (*-- cum is the result; set it to zero and accumulate the answer --------*)
     cum = ic/.{_Integer->0};
     Nest[(cum = cum+#; CellularAutomaton[{rn,{2,1},{1,1}},#,{{{1}}}])&,ic,nsteps];
     cum]/;(0<=rn && rn<2^10);
     
FacetSurfaceGridAutomaticIC[opts___Rule]:=Module[
      {ny,nx,ic},
      {ny,nx} = Global`ICDim   /.{opts} /. Options[FacetSurfaceGridData]; 
      {ny,nx} = If[OddQ[#],#,#+1]&/@{ny,nx};   (* make them odd numbers for an even number of intervals*)
      SeedRandom[ RandomSeed /. {opts} /. Options[FacetSurfaceGridData]];
      ic = RandomInteger[{0,1},{ny,nx}];
      If[(Global`SemiRandomIC /. {opts} /. Options[FacetSurfaceGridData]),
             ic = ((ic + Table[If[EvenQ[j],0,1],{j,-Floor[ny/2],Floor[ny/2]},{nx}])/.{2->1})];
      ic];

End[];
EndPackage[];
