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