1(* Author: Silvio Levy *)
2
3BeginPackage["BezierPlot`"]
4
5BezierPatch::usage = "\n
6  BezierPatch[ 4 by 4 array ] describes a bicubic Bezier patch\n
7  to be plotted using OOGL.m, for example.  The elements of the\n
8  array are usually 3D vectors, but the function doesn't care."
9
10BezierPlot::usage = "\n
11  BezierPlot[f,{t,tmin,tmax},{u,umin,umax}] \"plots\" the function f\n
12  of two variables t and u, using Bezier patches.  (Usually f returns\n
13  a 3D vector.)  BezierPlot returns an array of BezierPatch'es.\n
14  Options:\n
15    Epsilon is the increment used to compute derivatives (default: 10^-4)\n
16    PlotPoints is a number or pair of numbers indicating the fineness\n
17      of the subdivision (default: 5)";
18
19BezierPlot::baditer = "bad iterator";
20BezierPlot::badpp = "requested number of plot points is < 1";
21
22Options[BezierPlot] = {Epsilon -> .0001, PlotPoints->5}
23
24Begin["`private`"]
25
26BezierPlot[f_,{t_,tmin_,tmax_},{u_,umin_,umax_},options___]:=
27Block[{it,iu,nt,nu,dt,du,loc,dfdt,dfdt0,dfdu,dfdu0,dboth,all,
28      plotPoints=PlotPoints/.{options}/.Options[BezierPlot],
29      eps=Epsilon/.{options}/.Options[BezierPlot]},
30  If[Length[plotPoints]==2,{nt,nu}=plotPoints,nt=plotPoints;nu=plotPoints];
31  If[nx<1 || ny<1,Message[BezierPlot::badpp];Return[{}]];
32  dt=N[(tmax-tmin)/nt]; du=N[(umax-umin)/nu];
33  If[!(NumberQ[dt]&&dt>0&&NumberQ[du]&&du>0),
34    Message[BezierPlot::baditer];Return[{}]];
35  loc=Table[f/.{t->tmin+dt it,u->umin+du iu},{iu,0,nu},{it,0,nt}];
36  dfdt0=-loc+Table[f/.{t->tmin+dt(it+If[it<nt,1,-1]eps),u->umin+du iu},
37                   {iu,0,nu},{it,0,nt}];
38  dfdt=MapAt[-1#&,#,-1]& /@ dfdt0/(3 eps);
39  dfdu0=-loc+Table[f/.{t->tmin+dt it,u->umin+du(iu+If[iu<nu,1,-1]eps)},
40                   {iu,0,nu},{it,0,nt}];
41  dfdu=MapAt[-1#&,#,-1]& @ dfdu0/(3 eps);
42  dboth=MapAt[-1#&,#,-1]& /@ MapAt[-1#&,#,-1]& @
43        (-loc-dfdt0-dfdu0+Table[f/.{t->tmin+dt(it+If[it<nt,1,-1]eps),
44                                     u->umin+du(iu+If[iu<nu,1,-1]eps)},
45                                {iu,0,nu},{it,0,nt}])/(9 eps^2);
46  all=Transpose[{loc,dfdt,dfdu,dboth},{3,1,2}]//N;
47  Map[
48  BezierPatch[
49  {{#[[1,1]],#[[1,1]]+#[[1,2]],#[[2,1]]-#[[2,2]],#[[2,1]]},
50  {#[[1,1]]+#[[1,3]],#[[1,1]]+#[[1,2]]+#[[1,3]]+#[[1,4]],
51    #[[2,1]]-#[[2,2]]+#[[2,3]]-#[[2,4]],#[[2,1]]+#[[2,3]]},
52  {#[[3,1]]-#[[3,3]],#[[3,1]]+#[[3,2]]-#[[3,3]]-#[[3,4]],
53   #[[4,1]]-#[[4,3]]-#[[4,2]]+#[[4,4]], #[[4,1]]-#[[4,3]]},
54  {#[[3,1]],#[[3,1]]+#[[3,2]],#[[4,1]]-#[[4,2]],#[[4,1]]}}]&,
55    Transpose[{Drop[#,-1]&/@ Drop[#,-1]& @ all, Rest/@Drop[#,-1]& @ all,
56                    Drop[#,-1]&/@ Rest @ all, Rest/@Rest @ all}, {3,1,2}],{2}]
57]
58
59Format[x_BezierPatch]:="-BezierPatch-"
60
61End[]
62EndPackage[]
63