Appendix: Listing of the Mathematica package FoldPlot.m


(*^

::[	Information =

	"This is a Mathematica Notebook file.  It contains ASCII text, and can be
	transferred by email, ftp, or other text-file transfer utility.  It should
	be read or edited using a copy of Mathematica or MathReader.  If you 
	received this as email, use your mail application or copy/paste to save 
	everything from the line containing (*^ down to the line containing ^*)
	into a plain text file.  On some systems you may have to give the file a 
	name ending with ".ma" to allow Mathematica to recognize it as a Notebook.
	The line below identifies what version of Mathematica created this file,
	but it can be opened using any other version as well.";

	FrontEndVersion = "X Window System Mathematica Notebook Front End Version 2.2";

	X11StandardFontEncoding; 
	
	fontset = title, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e8,  24, fontName, "times";
	fontset = subtitle, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, bold, e6,  18, fontName, "times";
	fontset = subsubtitle, inactive, noPageBreakBelow, noPageBreakInGroup, nohscroll, preserveAspect, groupLikeTitle, center, M7, italic, e6,  14, fontName, "times";
	fontset = section, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, grayBox, M22, bold, a20,  18, fontName, "times";
	fontset = subsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, blackBox, M19, bold, a15,  14, fontName, "times";
	fontset = subsubsection, inactive, noPageBreakBelow, nohscroll, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12,  12, fontName, "times";
	fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, fontName, "times";
	fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, fontName, "times";
	fontset = input, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold,  12, fontName, "courier";
	fontset = output, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5,  12, fontName, "courier";
	fontset = message, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23,  12, fontName, "courier";
	fontset = print, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23,  12, fontName, "courier";
	fontset = info, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeOutput, M42, N23,  12, fontName, "courier";
	fontset = postscript, PostScript, formatAsPostScript, output, inactive, noPageBreakInGroup, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287,  12, fontName, "courier";
	fontset = name, inactive, noPageBreakInGroup, nohscroll, preserveAspect, M7, italic, B65535,  10, fontName, "times";
	fontset = header, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic,  12, fontName, "times";
	fontset = leftheader,  12, fontName, "times";
	fontset = footer, inactive, nohscroll, noKeepOnOnePage, preserveAspect, center, M7, italic,  12, fontName, "times";
	fontset = leftfooter,  12, fontName, "times";
	fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, fontName, "times";
	fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, fontName, "times";
	fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, fontName, "courier";
	fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, fontName, "times";
	fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, fontName, "times";
	fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, fontName, "times";
	fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, fontName, "times";
	fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, fontName, "times";currentKernel; 
]
:[font = input; preserveAspect]
(* :Copyright: 1998-1999,   Macquarie University, Sydney*)

(*:Version: Mathematica 2.2 *)

(*:Title: Graphics for Folding Deformations. *)

(*:Author: Ross Moore *)

(*:Keywords:
	Geophysical Deformations, folds, folding, axial, shear, rotation, contours
*)

(*:Requirements: None. *)

(*:Warnings: None. *)

(*:Summary:
	This package provides the functions and coding used to produce
	the examples of simulated geophysical folding deformations presented
	in the paper ``Three-dimensional reconstruction and modelling of complexly
	folded surfaces using \Mathematica'' by Ross Moore and Scott Johnson,
	published in  Computers & Geosciences, (special volume) 1999.	
*)

:[font = input; initialization; preserveAspect]
*)
axial[t_] := Sin[-Pi t];

axial::usage = " axial[t]  generates the deformations used in  axiZY, axiXZ etc.
	The default deformation function is sinusoidal:  axial[t_] := Sin[-Pi t] .
	This should be overridden, when desired."; 
(*
:[font = input; initialization; preserveAspect]
*)
BeginPackage["FoldPlot`",{"Geometry`Rotations`","Graphics`ParametricPlot3D`"}];
Off[General::spell,General::spell1];
(*
:[font = input; initialization; preserveAspect]
*)
FoldPlot::usage = "FoldPlot[fun, Xrange, Yrange] plots the result of a flat surface
	 deformed according to fun, as a parametric function of 2 variables, with ranges
	 given by the range specifiers Xrange and Yrange. Specific options for Plot
	 are automatically enforced, while any extra options are passed to Show .";
(*
:[font = input; initialization; preserveAspect]
*)
FoldPlotEPS::usage = "FoldPlotEPS[name, fun, Xrange, Yrange] does the same
	as FoldPlot[fun, Xrange, Yrange], but saves the graphic into a file 'name.eps'
	using Encapsulated PostScript (EPS) Format.";
(*
:[font = input; initialization; preserveAspect]
*)
FoldContours::usage = "FoldContours[fun, Xrange, Yrange] produces a ContourPlot
	of the deformation function fun of two variables, with ranges given by the range
	specifiers Xrange and Yrange.  FoldContours[fun, Xrange, Yrange, num ]
	presents num contours, instead of the default 7. Any extra options are
	passed to ContourPlot.";
(*
:[font = input; initialization; preserveAspect]
*)
FoldFaceContours::usage = "FoldFaceContours[fun] produces a layout of 3 contour plots
	for the deformation function fun . These plots give the contours in 3 faces of
	a rectangular box, meeting at a single vertex, so sharing an edge pairwise.
	FoldFaceContours[fun, num] specifies that num contours be used, covering all
	of the three faces; not all contours need appear in each face.";
(*
:[font = input; initialization; preserveAspect]
*)
FoldContourMovie::usage = "FoldContourMovie[fun] creates a sequence of Contour plots
	of the deformation function fun, in a series of slices through a rectangular block.
	These may be animated for viewing as a movie.
	 FoldContourMovie[fun, num_frames]  allows the number of frames to be specified;
	the default of 4 is deliberately small, to keep short initial testing time.
	 FoldContourMovie[fun, num_frames, num_levels]  allows the number of contour
	levels to be given for the complete set of slices; not all levels need appear
	in each slice.";
(*
:[font = input; initialization; preserveAspect]
*)
axiZY::usage = axiZX::usage = axiYX::usage = axiYZ::usage = axiXY::usage = axiXZ::usage = 
	"axiZY[][{x,y,z}]  produces a deformation in the z-direction, which varies\n
	along the y-direction according to the function axial[y], to give {x,y,z+axial[y]}.\n
	 axiZY[u][{x,y,z}] magnifies the deformation to give {x,y, z+u*axial[y]}.\n
	 axiZY[u, v][{x,y,z}] rescales to give  {x,y, z+u*axial[v*y]}.\n
	Variants axiZX , axiYX , axiYZ , axiXY , axiXZ  deform along other axis directions.";
(*
:[font = input; initialization; preserveAspect]
*)
shrZY::usage = shrZX::usage = shrYX::usage = shrYZ::usage = shrXY::usage = shrXZ::usage = 
	"shrZY[][{x,y,z}]  produces a shear deformation in the Z-direction, which varies
	along the y-direction to give {x,y,z+y}.\n
	 shrZY[u][{x,y,z}] rescales the shear to give {x,y, z+u*y} .\n
	Variants shrZX , shrYX , shrYZ , shrXY , shrXZ  shear along other axis directions.";
(*
:[font = input; initialization; preserveAspect]
*)
rot::usage = "rot[a,b,c][{x,y,z}] rotates the vector {x,y,z} by Euler-angles a, b, c
	measured in Degrees.";
(*
:[font = input; initialization; preserveAspect]
*)
Begin["`Private`"];
(*
:[font = input; initialization; preserveAspect]
*)
FoldPlot[ fun_ ,Xrange_List,Yrange_List,opts___] := 
Module[{thisplot,allopts} ,
	thisplot = ParametricPlot3D[
		Evaluate[fun],Xrange, Yrange,DisplayFunction-> Identity,opts ];
	Off[FaceGrids::fgstl];
	Show[Graphics3D[{EdgeForm[]
		,FaceForm[SurfaceColor[GrayLevel[1]],SurfaceColor[GrayLevel[.75]]]
		,thisplot//First}]
		, DisplayFunction->$DisplayFunction
		, ViewPoint->{1.96, -2.16, 1.72}
		, LightSources -> {
         {{1., 0., 1.},RGBColor[0,0, 1]}, {{-1., 1.,1.},RGBColor[1, 0, 0]},
         {{-1., -1., 1.},RGBColor[0, 1, 0]},{{.8, 1.2, 0.1},RGBColor[.5,.5,.5]}}
	, AxesLabel->{FontForm["x",{"Helvetica",12}],FontForm["y",{"Helvetica",12}]
		,FontForm["z",{"Helvetica",12}]}
	, BoxRatios -> {1,1,.8}
	, AxesEdge->{{-1,-1},{1,-1},{1,1}}
		,FaceGrids->{{{0,-1,0},{,{-1,0,1}}},{{1,0,0},{,{-1,0,1}}}}
	, Ticks->{{-1,0,1},{-1,0,1},{-1,0,1}}
	, PlotRange->{{-1,1},{-1,1},1.3{-1,1}}
	, thisplot//Last
]]
(*
:[font = input; initialization; preserveAspect]
*)
FoldPlotEPS[name_,fun_ ,Xrange_List,Yrange_List,opts___] := 
Module[{thisplot,allopts,eps,xbm,epsDisplay,xbmDisplay} ,
	eps = StringJoin["!psfix -epsf > ", name,".eps"];
	Print["Output to: ",eps];
	epsDisplay = Display[eps,#]&;
	thisplot = ParametricPlot3D[Evaluate[fun],Xrange,Yrange
			, DisplayFunction-> Identity,opts ];
	Off[FaceGrids::fgstl];
	Show[Graphics3D[{EdgeForm[]
		,FaceForm[SurfaceColor[GrayLevel[1]],SurfaceColor[GrayLevel[.75]]]
		,thisplot//First}]
		(*, DisplayFunction->$DisplayFunction*)
		, DisplayFunction-> epsDisplay
		, ViewPoint->{1.96, -2.16, 1.72}
		, LightSources -> {
		{{1., 0., 1.}, RGBColor[0,0, 1]}, {{-1., 1.,1.}, RGBColor[1, 0, 0]}
	  ,{{-1., -1., 1.}, RGBColor[0, 1, 0]},{{.8, 1.2, 0.1}, RGBColor[.5,.5,.5]}}
	, AxesLabel->{FontForm["x",{"Helvetica",10}],FontForm["y",{"Helvetica",10}]
		,FontForm["z",{"Helvetica",10}]}
	, BoxRatios -> {1,1,.8}
	, AxesEdge->{{-1,-1},{1,-1},{1,1}}
		,FaceGrids->{{{0,-1,0},{,{-1,0,1}}},{{1,0,0},{,{-1,0,1}}}}
	, Ticks->{{-1,0,1},{-1,0,1},{-1,0,1}}
	, PlotRange->{{-1,1},{-1,1},1.3{-1,1}}
	, thisplot//Last
]]
(*
:[font = input; initialization; preserveAspect]
*)
FoldPlot`Private`axial := Global`axial;
axiZY[u_:1,v_:1] := ((#) + {0,0,u*axial[#[[2]] v]})&;
axiZX[u_:1,v_:1] := ((#) + {0,0,u*axial[#[[1]] v]})&;
axiYX[u_:1,v_:1] := ((#) + {0,u*axial[#[[1]] v],0})&;
axiYZ[u_:1,v_:1] := ((#) + {0,u*axial[#[[3]] v],0})&;
axiXY[u_:1,v_:1] := ((#) + {u*axial[#[[2]] v],0,0})&;
axiXZ[u_:1,v_:1] := ((#) + {u*axial[#[[3]] v],0,0})&;
(*
:[font = input; initialization; preserveAspect]
*)
rot[a_:0,b_:0,c_:0] := (
	(#).N[RotationMatrix3D[a Degree, b Degree, c Degree]])&
(*
:[font = input; initialization; preserveAspect]
*)
shrZY[u_:1,v_:1] := ((#) + {0,0,u #[[2]]})&;
shrZX[u_:1,v_:1] := ((#) + {0,0,u #[[1]]})&;
shrYX[u_:1,v_:1] := ((#) + {0,u #[[1]],0})&;
shrYZ[u_:1,v_:1] := ((#) + {0,u #[[3]],0})&;
shrXY[u_:1,v_:1] := ((#) + {u #[[2]],0,0})&;
shrXZ[u_:1,v_:1] := ((#) + {u #[[3]],0,0})&;
(*
:[font = input; initialization; preserveAspect]
*)
FoldContours[fun_,Xrange_List,Yrange_List,levels_Integer:7,opts___]:=
ContourPlot[fun,Xrange,Yrange,opts
	,Axes -> None, Contours->levels, FrameTicks -> None
	,ColorFunction->(GrayLevel[Mod[Round[levels #+.4],2]]&)
	,PlotPoints -> 40
]
(*
:[font = input; initialization; preserveAspect]
*)
FoldFaceContours[fun_,levels_Integer:10,opts___]:= Module[
	{cont1,cont2,cont3,s,t,crange,cfun,defopts,copts, cbase = 2
		,cmin,cmax,cmin1,cmax1,cmin2,cmax2,cmin3,cmax3,legend}
	,
	cont1 = ContourPlot[fun[s,t,1.],{s,-1,1},{t,-1,1}
		,opts,FrameTicks -> None, ContourShading->False
		,PlotPoints -> 40, DisplayFunction->Identity];
	{cmin1,cmax1} = {Min[#],Max[#]}&@(cont1//First//Flatten);

	cont2 = ContourPlot[fun[s,-1.,t],{s,-1,1},{t,-1,1}
		,opts,FrameTicks -> None,ContourShading->False
		,PlotPoints -> 40, DisplayFunction->Identity];
	{cmin2,cmax2} = {Min[#],Max[#]}&@(cont2//First//Flatten);

(*	cont3 = ContourPlot[fun[1.,s,t],{s,-1,1},{t,-1,1}	*)
	cont3 = ContourPlot[fun[1,t,-s],{s,-1,1},{t,-1,1}
		,opts,FrameTicks -> None,ContourShading->False
		,PlotPoints -> 40, DisplayFunction->Identity];
	{cmin3,cmax3} = {Min[#],Max[#]}&@(cont3//First//Flatten);

	{cmin,cmax} = {Min[#],Max[#]}&@(
		{cmin1,cmax1,cmin2,cmax2,cmin3,cmax3});
	crange = Range[cmin,cmax,(cmax-cmin)/levels];
    defopts = {Axes -> None, Contours->crange, FrameTicks -> None
		, ContourShading->True, {ColorFunction-> cfun}
		, DisplayFunction->Identity };

	cfun = GrayLevel[Mod[Round[
		levels*((cmin1-cmin+#(cmax1-cmin1))/(cmax-cmin))+.45],cbase]]&;
    copts = Join[defopts, cont1//Last];
    cont1 = ContourGraphics[cont1//First, Sequence@@copts];

	cfun = GrayLevel[Mod[Round[
		levels*((cmin2-cmin+#(cmax2-cmin2))/(cmax-cmin))+.45],cbase]]&;
    copts = Join[defopts, cont2//Last];
    cont2 = ContourGraphics[cont2//First, Sequence@@copts];

	cfun = GrayLevel[Mod[Round[
		levels*((cmin3-cmin+#(cmax3-cmin3))/(cmax-cmin))+.45],cbase]]&;
    copts = Join[defopts, cont3//Last];
    cont3 = ContourGraphics[cont3//First, Sequence@@copts];

	legend = { GrayLevel[.8]
		, Rectangle[Scaled[{.175,.525}],Scaled[{.475,.825}]]
		, Rectangle[Scaled[{.525,.525}],Scaled[{.825,.825}]]
		, Rectangle[Scaled[{.175,.175}],Scaled[{.475,.475}]]
		, Hue[0]
		, Text[FontForm["Top",{"Helvetica",5}],Scaled[{.325,.675}]]
		, Text[FontForm["Right",{"Helvetica",5}],Scaled[{.675,.675}]]
		, Text[FontForm["Front",{"Helvetica",5}],Scaled[{.325,.325}]]};
(*	Show[GraphicsArray[{{cont1,Graphics[{}]},{cont2,cont3}}] *)
	Show[GraphicsArray[{{cont1,cont3},{cont2,Graphics[legend]}}]
		,DisplayFunction->$DisplayFunction]
]
(*
:[font = input; initialization; preserveAspect]
*)
FoldContourMovie[fun_,frames_Integer:4,levels_Integer:15,opts___] := 
Module[
	{cont,s,t,th,crange,defopts,copts,cbase=2,cmin,cmax,lrange,cfirst,i,j}
	,
	{cmin[0],cmax[0]} = 1000{1,-1}; j=0;
	Do[ ++j;
		cont[j] = ContourPlot[fun[s,t,-(s+t)/Sqrt[2] Tan[th]]
			,{s,-1,1},{t,-1,1}
			, opts, FrameTicks -> None, ContourShading->False
			, PlotPoints -> 40, DisplayFunction->Identity
			];
		{cmin[j],cmax[j]} = {Min[#],Max[#]}&[cont[j]//First//Flatten];
		cmin[0] = Min[cmin[0],cmin[j]];
		cmax[0] = Max[cmax[0],cmax[j]];
	, {th,-2Pi/5,2Pi/5,4Pi/(5*(frames-1))}];
	crange = Range[cmin[0],cmax[0],(cmax[0]-cmin[0])/(levels)];
    defopts = {Axes -> None, FrameTicks -> None, ContourShading->True };
	Do[
		i=0; cfirst[j] = Catch[If[#>=cmin[j],Throw[i],i++]&/@ crange];
		lrange[j] = Select[crange, TrueQ[And[ #>= cmin[j], #<= cmax[j]]]&];
    	copts[j] = Join[
		{ToExpression[StringJoin[{"ColorFunction->(modStartColor["
			, ToString[cbase],",",ToString[cfirst[j]],",#]&)"}]]}
		, defopts, cont[j]//Last];
	, {j,1,frames}]; colIndex = 0;
	Show[GraphicsArray[Evaluate[
		Partition[Table[ ContourGraphics[cont[j]//First
			, Contours-> lrange[j], Sequence@@copts[j] ]
			,{j,frames}]
		, Floor[(frames+1)/2]]
		]],DisplayFunction->(($DisplayFunction[#])&)]
]
(*
:[font = input; initialization; preserveAspect]
*)
modStartColor[base_,start_,j_]:= Module[ {tmp},
	If[Not[TrueQ[colIndex > 0]], colIndex = start];
	If[Chop[j]==1,tmp=++colIndex; colIndex = 0,
		If[Chop[j]==0, tmp=colIndex=start, tmp=++colIndex ]];
	GrayLevel[N[Mod[tmp,base]/(base-1)]]
]
(*
:[font = input; initialization; preserveAspect]
*)
End[];
On[General::spell,General::spell1];
EndPackage[];
(*
^*)


Ross Moore 1999-07-16