dansmath-ematica notebooks

dansmath > mathematica > notebooks

Some of Dan's Mathematica Notebooks

I give you now a few of my favorite cool Mathematica notebooks, some of which I showed as a speaker at the CMC Math Conferences at Lake Tahoe and in Monterey in 1998. These emphasize nice pictures to explain concepts. Pick the one that interests you by clicking on the title.
Note: You have to have Mathematica for these to run. You can Copy the text of a notebook here in your Web browser, Open Mathematica (Quit browser if you have limited RAM), make a New document, Paste it in, Enter the cells, and Watch the show!

1. Rotation of Axes with Conic Sections . . . . Enter your own conic section and this will explain steps and rotate to get rid of the xy-term, finally graphing the original and rotated curves.

2. Divisor Plots and Supercomposite Numbers . . . . Some numbers, like 12, have lots of divisors: 1, 2, 3, 4, 6, and 12; while others, like 13, have few: 1, 13. Wouldn't you like to see a ListPlot of this and some cool animation?

3. Graphic Display of Primes . . . . Listing primes was never easier! Create graphic pictures of the arrangement of prime numbers, where you can pick and animate the length and width of the table.

4. Disco Ball . . . . Teaches spherical coordinates and Hue and RGB options by a series of colored, mirrored balls. But how do they stay up there?

5. Color Me In(creasing) . . . . Any curve of your choosing can be displayed with colors depending on slope. The concave up parts are made thicker because that's where the water collects!

6. Curves in 3D . . . . Any 3D parametric curve of your choosing can be displayed, along with tangent, normal, and binormal vectors. Choose osculating circle and three-plane options!

In the notebooks below, import the text into Mathematica and then run it.

The stuff like this is Mathematica input commands;

(* The stuff like this is comments and won't interfere. *)

back to mathematica

(* 1. Rotation of Axes - Conic Sections (c) 1995-2000 *)
(* by Dan Bach, Diablo Valley College *) (* back to list of notebooks *)

(* A general conic section has equation: a x^2 + b x y + c y^2 + d x + e y + f = 0. *)

Clear[a,b,c,d,e,f,g]
g[x_,y_] := a x^2 + b x y + c y^2 + d x + e y + f
g[x,y]

(* Here's where you can put in your own coefficients or just enter this cell and use these! *)

a = 1;   b = 3;   c = -2;
d = 2;   e = 0;   f = -6;
g[x,y]

(* We can use the discriminant to predict the shape: *)

disc = b^2 - 4 a c
Print["The discriminant is ",disc,"."]

If[(disc < 0), shape := "n ellipse"];
If[(disc > 0), shape := " hyperbola"];
If[(disc = 0), shape := " parabola"];
Print["Your curve is a",shape,"."]

(* Let's see about the graph: we can solve for y and graph two functions... *)

func := Solve[g[x,y] == 0, y];
y1 = func[[1,1,2]]
y2 = func[[2,1,2]]
r = 5;

Off[Plot::plnr]
Plot[{y1,y2},{x,-r,r}, PlotRange->{{-r,r},{-r,r}},
			PlotStyle->{{Hue[.0],Thickness[.015]}, 
			            {Hue[.4],Thickness[.015]}},
			AspectRatio->Automatic];

(* How about not separating it into two functions? *)

ContourPlot[g[x,y],{x,-5,5},{y,-5,5},
	PlotPoints->40, Contours->{0}, ContourShading->False];

(* If the equation has an xy-term in it, the axes of the shape won't line up with the x and y axes. *)

(* The angle theta that does this is given by cot(2 theta) = (a - c) / b. *)

theta = If[c == a, Pi/2, ArcTan[b/(a - c)]]/2
ct = Cos[theta];  st = Sin[theta];

(* The new coords (u,v) are given by simple trig... *)

g[x,y] /. {x -> ct u - st v, y -> st u + ct v}

(* Now clean up the expression! *)

gr[u_,v_] := Evaluate[Expand[%]]
gr[u,v]
Simplify[gr[u,v]]
N[gr[u,v]]

(* We can solve for one variable as before, plot two more functions, and compare results! *)

funcr := Solve[gr[x,y] == 0, y];
r = 5;

yr1 = funcr[[1,1,2]];  yr2 = funcr[[2,1,2]];
Simplify[yr1]
Simplify[yr2]

Off[Plot::plnr]
Plot[{y1,y2,yr1,yr2},{x,-r,r}, 
			PlotRange->{{-r,r},{-r,r}},
			PlotStyle-> {{GrayLevel[.7],Thickness[.01],Dashing[{.04}]},
			             {GrayLevel[.7],Thickness[.01],Dashing[{.04}]},
		 	     {Hue[.0],Thickness[.015]}, {Hue[.4],Thickness[.015]}},
			Epilog->{GrayLevel[.6],Dashing[{.04}],
			     	     Line[{{-r ct,-r st},{r ct, r st}}],
					     Line[{{-r st, r ct},{r st,-r ct}}]},
			AspectRatio->Automatic];

(* Now there's a cool picture. Change your a, b, c, d, e, f and go through it again! *)

(* 2. Divisor Plots, Supercomposite Numbers (c) 1996-2000 *)
(* by Dan Bach, Diablo Valley College *) (* back to list of notebooks *)

(* A "supercomposite number" is one with more divisors than any smaller number. *)
(* Here are a couple of tables of how many divisors the first few numbers have. *)

Clear[d,k,m,n]
Div[n_] := Divisors[n]
d[n_] := Length[Div[n]]
p = Table[{n,d[n]},{n,1,10}];
	TableForm[p, TableSpacing->{0,5},
	TableHeadings->{None,{"n","d[n]"}}]

dees = Table[{n,d[n]}, {n,1,100}];
TableForm[dees,TableSpacing->{0,5},
	TableHeadings->{None,{"n","d[n]"}}]

(* Why not make a plot of d(n) vs n and see this? Here's a dot plot: *)

ds[n_] := Table[{k,d[k]}, {k,1,n}]
dp[n_] := ListPlot[ds[n], 
			PlotStyle->{PointSize[.015],Hue[0]},
			DisplayFunction->Identity];

Show[dp[50], DisplayFunction->$DisplayFunction];

(* Here's a joined plot: *)

jp[n_] := ListPlot[ds[n], PlotJoined->True,
			   PlotStyle->{GrayLevel[.5],AbsoluteThickness[.5]},
		   	DisplayFunction->Identity];
			
Show[jp[50], DisplayFunction->$DisplayFunction];

(* This highlights the supercomposites: *)

sc[n_] := ListPlot[
			  Table[If[d[k]>Max[Table[d[m],{m,1,k-1}]],
			  {k,d[k]},{1,1}], {k,1,n}], 
		     PlotStyle->{Hue[.4],PointSize[.02]},
		     AxesOrigin->{0,0}, DisplayFunction->Identity];

Show[sc[50], DisplayFunction->$DisplayFunction];

(* And this puts them all together! *)

Show[jp[50],dp[50],sc[50],
	DisplayFunction->$DisplayFunction];

(* Here's the animation I know you were hoping for! *)

Do[Show[jp[m],dp[m],sc[m], AxesOrigin->{0,0},
	PlotRange->{{-.1 m,1.1 m},{-.2 Log[m], 4.2 Log[m]}},
	 Epilog->{Hue[0],Text["n's up to "<>ToString[m],
	 {.1m,3.5 Log[m]},{-1,0}]}, DisplayFunction->$DisplayFunction], 
{m,10,400,10}];

(* 3. Graphic Display of Primes (c) 1997-2000 *)
(* by Dan Bach, Diablo Valley College *) (* back to list of notebooks *)

(* This is a way to see the distribution of primes. The natural numbers are put in *)
(* rows of k (you pick the k), and the primes are highlighted and then labeled! *)
(* First we use black dots to show the primes... *)

k = 10;
m[n_] := Mod[n,k]
a[n_] := {m[n],(m[n] - n)/k}
g[t_] := Table[Graphics[Disk[a[Prime[n]],.5]],{n,1,t}];

(* g[t] shows the first t primes *)

Show[g[26],AspectRatio->Automatic];

(* Now for the numbers themselves, but only in "prime locations." *)

k = 20;
m[n_] := Mod[n,k]
a[n_] := {m[n],(m[n] - n)/k}
txt[t_]:=Table[Graphics[Text[Prime[n],a[Prime[n]]]],{n,1,t}];
Show[txt[100], AspectRatio->Automatic, PlotRange->All];

(* Here's a nice gray combination: *)

k = 10;
m[n_]   := Mod[n,k]
a[n_]   := {m[n], (m[n] - n)/k}
gr[t_]  := Table[Graphics[
			{GrayLevel[.8], Disk[a[Prime[n]],.5]}], {n,1,t}];
txt[t_] := Table[Graphics[
			Text[Prime[n], a[Prime[n]]]], {n,1,t}];
			
Show[gr[25],txt[25], AspectRatio->Automatic, PlotRange->All];

(* This puts the primes under 1000 into neat columns! *)

k = 30;
m[n_]   := Mod[n,k]
a[n_]   := {m[n], (m[n] - n)/k}
gr[t_]  := Table[Graphics[
			{Hue[.2], Disk[a[Prime[n]],.6]}], {n,1,t}];
txt[t_] := Table[Graphics[
			Text[Prime[n], a[Prime[n]]]], {n,1,t}];
			
Show[gr[168],txt[168], AspectRatio->Automatic, PlotRange->All];

(* Animation showing patterns of dots with various numbers of columns: *)

Clear[m,n,k,a,gr,t]
m[n_]   := Mod[n,k]
a[n_]   := {m[n], (m[n] - n)/k}
gr[t_]  := Table[Graphics[
			{Hue[.2+.02k], Disk[a[Prime[n]],.6]}], {n,1,t}];
Do[Show[gr[168], Frame->True, FrameTicks->{Automatic,None},
   AspectRatio->Automatic, PlotRange->All],{k,25,40}]

(* I suppose now you want to color them in by, say, congruence mod d? *)
(* ¡No hay problema! *)

k = 12;
d = 10;
m[n_] := Mod[n,k]
a[n_] := {m[n],(m[n] - n)/k}
c[t_] := Table[Graphics[
	{Hue[Prime[n]/d],Disk[a[Prime[n]],.5]}],{n,1,t}];
txt[t_] := Table[Graphics[Text[Prime[n],a[Prime[n]]]],
			{n,1,t}];
Show[c[46],txt[46],AspectRatio->Automatic,PlotRange->All];

(* 4. Disco Ball (c) 1992-2000 *)
(* by Dan Bach · Diablo Valley College *) (* back to list of notebooks *)

(* There's a package for you! *)

Needs["Graphics`ParametricPlot3D`"]

(* Here's a generic Mathematica sphere. Note the spherical coordinates (hence the name).*)

ParametricPlot3D[{Cos[u] Sin[v], Sin[u] Sin[v], Cos[v]},
	{u,0,2Pi,Pi/15},{v,0,Pi,Pi/15}];

(* We free the ball from its go-go cage. *)

ParametricPlot3D[{Cos[u] Sin[v], Sin[u] Sin[v], Cos[v]},
	{u,0,2Pi,Pi/20},{v,0,Pi,Pi/12}, Axes->None, Boxed->False];

(* Let's put more and more mirrors on the ball. *)

Do[ParametricPlot3D[
	{Cos[u] Sin[v], Sin[u] Sin[v], Cos[v]},
	{u,0,2 Pi, Pi/n},{v,0,Pi,Pi/n}, Axes->None, Boxed->False],
{n,4,12,1}];

(* But these colors aren't funky enough. *)

Do[ParametricPlot3D[
	{Cos[u] Sin[v], Sin[u] Sin[v], Cos[v], Hue[Random[]]},
	{u,0,2 Pi, Pi/n},{v,0,Pi,Pi/n},
	 Lighting->False, Axes->None, Boxed->False],
{n,8,20,1}];

(* Maybe too random and bright. *)

Do[ParametricPlot3D[
	{Cos[u] Sin[v], Sin[u] Sin[v], Cos[v],
	 RGBColor[t=Random[],t (1+Sin[2v])/2, n/20]},
	{u,0,2 Pi, Pi/n},{v,0,Pi,Pi/n},
	 Lighting->False, Axes->None, Boxed->False],
{n,8,20,1}];

(* Let's open the ball up and look inside! *)

Do[ParametricPlot3D[
	{Cos[u] Sin[v], Sin[u] Sin[v], Cos[v],
	 RGBColor[Random[],(1+Sin[2v])/3,n/20]},
	{u, -Pi/3 + (n-8) Pi/90, 5 Pi/3 - (n-8) Pi/90, Pi/n},
	{v, (n-8) Pi/50, Pi, Pi/n},
	 Lighting->False, Axes->None, Boxed->False],
{n,8,20,1}];

(* 5. Color Me In(creasing)! (c) 1997-2000 *)
(* by Dan Bach, Diablo Valley College *) (* back to list of notebooks *)
(* For fast downloading of this (text) document, this notebook contains links to GIF graphs; just click! *)

(* Let's plot a function that has its ups and downs. *)

f[x_] := x^3 - 3x
a = Plot[f[x],{x,-5,5}];

click here for figure 1

(* Here's a crude way of showing the decreasing parts of a function...*)

b = Plot[f[x],{x,-1,1}, PlotStyle-> 
	{{Hue[.6], Thickness[.01], Dashing[{.02}]}},
	DisplayFunction->Identity];
Show[a,b];

click here for figure 2

(* Run this module in v2.2 any time you want your polynomials in descending order. *)

(* Thanks to Alan Deguzman and Bill Davis for this. Version 3 has this built in. *)

Unprotect[Plus];
Format[Literal[Plus[p__]]] := 
	Module[{s1,s2}, s1 = Hold[p]; s2 = Reverse[s1];
		ReplacePart[HoldForm[Evaluate[s2]],Plus,{1,0}]/;
			OrderedQ[s1] && s1 =!= s2]

(* We abandon the Plot command for the Table setup *)

f[x_] := Product[x-c,{c,-2,2}]; f[x]
h = .1;
curve =  Table[Graphics[
			 {Hue[If[f'[x]>0,.3,.0]], Thickness[.007],
			 Line[{{x,f[x]},{x+h,f[x+h]}}]}], {x,-2,2,h}];
Show[curve, Axes->True, Frame->True];

click here for figure 3

(* Yay! Now let's look at concavity: *)

f[x_] := Product[x-c,{c,-2,2}]; f[x]
h = .1;
curve =  Table[Graphics[{Thickness[If[f''[x]>0,.015,.005]],
			Line[{{x,f[x]},{x+h,f[x+h]}}]}], {x,-2,2,h}];
Show[curve, Axes->True, Frame->True];

click here for figure 4

(* The concave-up parts are thicker; they collect more water. Now show both features:*)

f[x_] := Product[x-c,{c,-2,2}]; f[x]
h = .1;
curve =  Table[Graphics[{Hue[If[f'[x]>0,.3,.0]],
			Thickness[If[f''[x]>0,.015,.007]],
			Line[{{x,f[x]},{x+h,f[x+h]}}]}], {x,-2,2,h}];
Show[curve, Axes->True, Frame->True];

click here for figure 5

(* Try this next function, then play with different f[x]'s and x-ranges. *)

(* If you're red/green color-blind, then change the green to blue (change .3 in Hue to .6).*)

f[x_] := Sin[x^2]; f[x]
h = .01;
curve =  Table[Graphics[{Hue[If[f'[x] > 0, .3, 0, 0]],
			Thickness[If[f''[x] > 0,.015,.007,.007]],
			Line[{{x, f[x]}, {x+h, f[x+h]}}]}], {x, -3, 3-h, h}];
Show[curve, Axes->True, Frame->True];

click here for figure 6

(* The right color scheme is greener upper, redder steeper; concave up parts are thicker. *)

f[x_] := Sin[x^2]
h = .02; {a,b} = {-3,3}; f[x]
maxm = Max[Table[f'[x], {x,N[a],b,h}]];
minm = Min[Table[f'[x], {x,N[a],b,h}]];
redgreen =  Table[Graphics[{Hue[.3(f'[x]-minm)/(maxm-minm)],
			Thickness[If[f''[x] > 0,.018,.006,.006]],
			Line[{{x, f[x]}, {x+h, f[x+h]}}]}],{x,N[a],b-h,h}];
Show[redgreen, Axes->True, Frame->True];

click here for figure 7

(* I'm hoping you use this in your calculus class, whether you're a teacher or a student. *)

(* Please improve on this, and add features. E-mail me with your embellishments. *)

(* 6. Curves in 3D! (c) 1995-2000 *)
(* by Dan Bach, Diablo Valley College *) (* back to list of notebooks *)

(* Comments are enclosed like this *)

(* Let's define a curve in 3D space.  *)
(* This one's a helix; please feel free to *)
(* change the coordinate functions x, y, and z.*)

x[t_] := Cos[t]
y[t_] := Sin[t]
z[t_] := t/3

r[t_] := {x[t],y[t],z[t]}

curve = ParametricPlot3D[Evaluate[r[t]],{t,0,4Pi}]  (* Helix!*)

(* This is a "Graphics3D" object, then. What things can we control? *)

Options[Graphics3D]

(* Gluing on tangent, normal, and binormal vectors *)

(* Let's free it from its cage and add a tangent vector. *)

curve = Show[curve, AspectRatio->Automatic,
			Axes->None, Boxed->False];
len[{a_,b_,c_}] := Sqrt[a^2+b^2+c^2]
T[t_] := r'[t]/len[r'[t]]
tvec[t_] := Graphics3D[{Hue[0],Thickness[.02],
			Line[{r[t], r[t] + T[t]}]}]
Show[curve,tvec[1]];

(* Here's a movie shot while the tangent 
vector was travelling along the curve. *)

Do[Show[curve,tvec[t]], {t,0,2Pi,Pi/8}]

(* A little too animated, isn't it?  Let's steady it 
down by making the "invisible box" a fixed size. *)

Do[Show[curve,tvec[t],
	PlotRange->{{-1.2,1.2},{-1.2,1.2},{-.5,Pi}}], 
{t,0,2Pi,Pi/8}]

(* We've got the red tangent vector. Now we want the normal 
vector (blue) pointing towards the center of curvature... *)

T[t_] := r'[t]/len[r'[t]]
Nor[t_] := T'[t]/len[T'[t]]
tvec[t_] := Graphics3D[{Hue[0],Thickness[.02],
			Line[{r[t], r[t] + T[t]}]}]
nvec[t_] := Graphics3D[{Hue[.6],Thickness[.02],
			Line[{r[t], r[t] + Nor[t]}]}]
Show[curve,tvec[1],nvec[1]];

(* ... and another movie shot while the two vectors climbed up the 
curve together.  The plane they determine is called the "osculating 
plane"; it's the plane that best holds the curve at that point. *)

Do[Show[curve,tvec[t],nvec[t],AspectRatio->Automatic,
	PlotRange->{{-1.2,1.2},{-1.2,1.2},{-.5,Pi}}], 
{t,0,2Pi,Pi/8}];

(* The third vector in the "moving trihedron" (like a varying 
set of x-y-z-coordinate axes) is called the "binormal" vector.  
It's perpendicular to the osculating plane. *)

B[t_] := Cross[T[t],Nor[t]]/len[Cross[T[t],Nor[t]]]
bvec[t_] := Graphics3D[{Hue[.8],Thickness[.02],
			Line[{r[t], r[t] + B[t]}]}]

Show[curve,tvec[1],nvec[1],bvec[1]];

(* Here's how the three unit vectors create a "moving frame" along the curve. *)

Do[Show[curve,tvec[t],nvec[t],bvec[t],AspectRatio->Automatic,
	PlotRange->{{-1.2,1.2},{-1.2,1.2},{-1,Pi}}], 
{t,0,2Pi,Pi/8}];

(* Adding the osculating circle and three planes *)

(* Let's review the functions we've defined. *)

x[t_] := Cos[t]
y[t_] := Sin[t]
z[t_] := t/3

r[t_] := {x[t],y[t],z[t]}

curve = ParametricPlot3D[Evaluate[Flatten[
{r[t],{{RGBColor[.8,.2,.2],Thickness[.01]}}},1]],{t,0,4Pi}];  (* Helix!*)

len[{a_,b_,c_}] := Sqrt[a^2+b^2+c^2]
T[t_] := r'[t]/len[r'[t]]
tvec[t_] := Graphics3D[{Hue[0],Thickness[.02],
			Line[{r[t], r[t] + T[t]}]}]
Show[curve,tvec[1], AspectRatio->Automatic,
			Axes->None, Boxed->False];

Nor[t_] := T'[t]/len[T'[t]]
nvec[t_] := Graphics3D[{Hue[.6],Thickness[.02],
			Line[{r[t], r[t] + Nor[t]}]}]
Show[curve,tvec[1],nvec[1],Axes->None];

B[t_] := Cross[T[t],Nor[t]]/len[Cross[T[t],Nor[t]]]
bvec[t_] := Graphics3D[{Hue[.8],Thickness[.02],
			Line[{r[t], r[t] + B[t]}]}]

Show[curve,tvec[1],nvec[1],bvec[1],Axes->None];

(* Here's the osculating circle along the curve, with the three 
unit vectors. The bigger the curvature, the smaller the circle.*)

k[t_] := len[r''[t]]/(len[r'[t]])^1.5
rad[t_] := 1/k[t]
cent[t_] := r[t] + rad[t] Nor[t]
oscul[t_] := ParametricPlot3D[Evaluate[Append[
		r[t]+ rad[t]*(Cos[u] T[t] + (Sin[u]+1) Nor[t]),
			Hue[.8]]],{u,0,2Pi},DisplayFunction->Identity];
cenpt[t_] := Graphics3D[{PointSize[.03],Point[cent[t]]}]
Show[curve,tvec[1],nvec[1],oscul[1],cenpt[1],
	DisplayFunction->$DisplayFunction,Axes->None];


(* Time for a rolling osculating circle! *)

Do[Show[curve,tvec[t],nvec[t],bvec[t],oscul[t],cenpt[t],
		PlotRange->{{-1.2,1.2},{-1.2,1.2},{-.5,5}},
	DisplayFunction->$DisplayFunction,Boxed->True ],
	{t,0,2Pi,Sqrt[Pi]/8}]

(* This one accelerates down the hill. *)

top = 4Pi;
Do[Show[curve,
		tvec[top-t^2],nvec[top-t^2],bvec[top-t^2],oscul[top-t^2],
		PlotRange->{{-1.3,1.3},{-1.3,1.3},{-1,5.5}},
		DisplayFunction->$DisplayFunction,Axes->None ],
	{t,0,c=Sqrt[top],c/20}]

(* Now we put in the three planes:    *)
(* The osculating plane with T and N  *)
(* The normal plane with N and B      *)
(* The rectifying plane with B and T  *)

oscupl[t_] := Graphics3D[{Hue[.8],
			Line[{r[t],r[t]+Nor[t]}],Line[{r[t]+T[t],r[t]+ T[t]+ Nor[t]}],Table[
			Line[{r[t]+a Nor[t],r[t]+a Nor[t]+T[t]}],{a,0,1,.2}]}]
normpl[t_] := Graphics3D[{Hue[0],
			Line[{r[t],r[t]+Nor[t]}],Line[{r[t]+B[t],r[t]+ B[t]+ Nor[t]}],Table[
			Line[{r[t]+a Nor[t],r[t]+a Nor[t]+B[t]}],{a,0,1,.2}]}]
rectpl[t_] := Graphics3D[{Hue[.6],
			Line[{r[t],r[t]+T[t]}],Line[{r[t]+B[t],r[t]+ B[t]+T[t]}],Table[
			Line[{r[t]+a T[t],r[t]+a T[t]+B[t]}],{a,0,1,.2}]}]

Show[curve,oscupl[1],rectpl[1],normpl[1],oscul[1],cenpt[1],
	bvec[1],tvec[1],nvec[1], Axes->None]

(* And now the whole cast of characters. *)
(* Choose which things to leave in the "Show" command.*)
(* Change the x[t], y[t], and z[t] if you wish. *)

Do[Show[curve,oscupl[t],rectpl[t],normpl[t],oscul[t],cenpt[t],
			bvec[t],tvec[t],nvec[t],
			PlotRange->{{-1.5,1.5},{-1.5,1.5},{-.5,4Pi/3}}],
	{t,0,2Pi,Pi/12}]

(* Now look at Bach's Car on your own curve. *)
(* Change the x[t], y[t], and z[t], and re-enter the cell. *)
(* Adjust the PlotRange, and the {t, a, b} if you need to! *)

x[t_] := Cos[t]
y[t_] := Sin[t]
z[t_] := Cos[2t]

r[t_] := {x[t],y[t],z[t]}

curve = ParametricPlot3D[Evaluate[Flatten[
{r[t],{{RGBColor[.8,.2,.2],Thickness[.01]}}},1]],{t,0,4Pi}];  (* Helix!*)

Do[Show[curve,oscupl[t],rectpl[t],normpl[t],oscul[t],cenpt[t],
			bvec[t],tvec[t],nvec[t],
			PlotRange->{{-2,2},{-2,2},{-2,2}}, Axes->None],
	{t,Pi,5Pi/4,Pi/20}]

dansmath > mathematica > top of page