Vector Initializaton

Off[General::spell];
Off[General::spell1];
Off[ParametricPlot3D::ppcom];
Needs["Graphics`Colors`"];

SetOptions[Graphics,Axes->False,AxesOrigin->{0,0},AspectRatio->Automatic];
SetOptions[Graphics3D,ViewPoint->{5,2,1},Boxed->False,AspectRatio->Automatic];
SetOptions[ParametricPlot3D,ViewPoint->{5,2,1},AspectRatio->Automatic];
SetOptions[Plot3D,ViewPoint->{5,2,1},AspectRatio->Automatic];
SetOptions[ParametricPlot,AspectRatio->Automatic];
SetOptions[Plot,AspectRatio->Automatic];

colorQ[a_]:=(Head[a]===RGBColor);
grayQ[a_]:=(Head[a]===GrayLevel);
thickQ[a_]:=(Head[a]===Thickness);
absthickQ[a_]:=(Head[a]===AbsoluteThickness);
notruleQ[u_] := !(Head[u]===Rule);
ruleQ[u_] := (Head[u]===Rule);

Norm::usage = "Norm[expr] returns the Euclidean length of the
		vector, 'expr'.";

Norm[a_?VectorQ] := Sqrt[a.a];

Angle[a_?VectorQ,b_?VectorQ] := ArcCos[Simplify[(a.b)/(Norm[a] Norm[b])]];

Proj[a_?VectorQ,b_?VectorQ] := ((a.b)/(b.b))b;

tol = .5*10^-6;

OrthUp[a_?VectorQ] :=
	Module[{b,c,d,f},
		d = N[a/Norm[a]];
		If[Chop[d[[3]]^2 - 1,tol] == 0, b={1,0,0} ,
		b = {0,0,1} - d[[3]] d; b = b/Norm[b];];
	Return[b]];
	
OrthSet[a_?VectorQ] :=
	Module[{b,c,d},
		d = N[a/Norm[a]];
		If[Chop[d[[1]],tol]==0,b={1,0,0},
			 If[Chop[d[[2]],tol]==0,b={0,1,0},
			  b={-d[[2]],d[[1]],0}]];
			b=b/Norm[b];c = Cross[d,b];
	Return[{b,c/Norm[c],d}]];

				
Rotate[a_?VectorQ,theta_] :=
	Module[{u,v,r},
		u = OrthSet[a];
		v = Transpose[u];
		r = N[{{Cos[theta],-Sin[theta],0},{Sin[theta],Cos[theta],0},{0,0,1}}];
	Return[v.r.u]];

Vector2D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,arrow},
	ntail = N[tail];
	ntip = N[tip];
	
	sopts = Join[Select[{opts},colorQ],Select[{opts},
				thickQ],Select[{opts},absthickQ],Select[{opts},grayQ]];

If[Chop[ntail.ntail+ntip.ntip,tol] == 0,
			(*too small to plot or not number*)
			Return[Graphics[Point[{0,0}]]];
	,						
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = 0.075 (ntip-ntail) ;
	up = 0.15 {-Punit[[2]],Punit[[1]]};
	arrow = Line[{ntail,ntip,ntip-Punit+up,ntip,ntip-Punit-up}];
	]];
Return[Graphics[{Join[{AbsoluteThickness[2]},Join[sopts,{arrow}]]}]];
];

Vector3D[tail_,tip_,opts___]:=
	Module[{ntail,ntip,sopts,Punit,up,u2,u3,u4,rotmat,arrow},
	ntail = N[tail];
	ntip = N[tip];
	
	sopts = Join[Select[{opts},colorQ],Select[{opts},
				thickQ],Select[{opts},absthickQ],Select[{opts},grayQ]];

If[Chop[ntail.ntail+ntip.ntip,tol] == 0,
			(*too small to plot or not number*)
			Return[Graphics3D[Point[{0,0,0}]]];					
	,		
If[Chop[(ntip-ntail).(ntip-ntail),tol] == 0,
		Return[Graphics3D[{Join[sopts,{Line[{ntail,ntip}]}]}]];
	,
	Punit = (ntip-ntail);
	rotmat = Rotate[Punit,N[Pi/2]];
	up = OrthUp[Punit];
	Punit = .075 Punit;
	up = .15 Norm[Punit] up;
	u2 = rotmat.up;
	u3 = rotmat.u2;
	u4 = rotmat.u3;
	arrow = Line[{ntail,ntip,
				ntip - Punit + up,ntip,
				ntip - Punit + u2,ntip,
				ntip - Punit + u3,ntip,
				ntip - Punit + u4}];
	]];
Return[Graphics3D[{Join[{AbsoluteThickness[2]},Join[sopts,{arrow}]]}]];
];

aRCtAN[y_,x_] := Module[{v},
	If[ Chop[N[x],tol] == 0 ,If[ N[y] < 0, v = -N[Pi/2], v = N[Pi/2]],
		v = N[ArcTan[y/x]];
		If[ x < 0 && y < 0, v = N[Pi -  v],v]];
	Return[v]];

Arc2D[X_,Y_,opts___] := Graphics[{opts,
		Disk[{0,0},1,{Min[
							aRCtAN[
						N[X[[2]]],N[X[[1]]]
								]
							,
							aRCtAN[
							N[Y[[2]]],N[Y[[1]]]
								]
						]
						,
					Max[
							aRCtAN[
						N[X[[2]]],N[X[[1]]]
								]
							,
							aRCtAN[
							N[Y[[2]]],N[Y[[1]]]
								]
						]}
				]}
				];
Arc3D[X_,Y_] := 
	Module[{t,U,V},
		U[V_]:= N[(1/Sqrt[V.V]) V];
		Return[Graphics3D[
		Polygon[Append[Table[U[t X + (1-t) Y],{t,0,1,0.1}],{0,0,0}]]
				]
		]
	];

(*Zero Arrow causes trouble. Put in Vector?*)

Vector[tail_,tip_,opts___]:=
	If[Length[N[tail]]==2,Return[Vector2D[tail,tip,opts]],
		Return[Vector3D[tail,tip,opts]]];
		
Arrow[vec_,opts___]:= 
	Module[{a,b},
				If[Length[N[vec]]==2,
				a=ToExpression["Tail"]/.Select[{opts},ruleQ]/.
					{ToExpression["Tail"]->{0,0}};
				Return[Vector2D[a,a+vec,opts]],
				a=ToExpression["Tail"]/.Select[{opts},ruleQ]/.
					{ToExpression["Tail"]->{0,0,0}};
				Return[Vector3D[a,a+vec,opts]]]
				];

Frame2D[delta_] := 
      Graphics[{
     Green,Line[{{-delta/4,0},{delta,0}}],
                    Text["x",{1.01 delta,0}],Line[{{0,-delta/4},{0,delta}}],Text["y",{0, 1.01 delta}]}];

Frame3D[delta_] := 
      Graphics3D[{Green,Line[{{-delta/4,0,0},{delta,0,0}}],Text["x",{1.01 delta,0,0}],Line[{{0,-delta/4,0},{0,delta,0}}],Text["y",{0, 1.01 delta,0}],Line[{{0,0,-delta/4},{0,0,delta}}],Text["z",{0,0,1.01 delta}]}];

View2D[stuff_] := Show[stuff,Frame2D[1]];
View3D[stuff_] := Show[stuff,Frame3D[1]];

Slice Initialization

Needs["Graphics`Colors`"];

SetOptions[Plot3D,ViewPoint->{5,3,1},
	BoxRatios ->Automatic,PlotPoints->30,
	AxesLabel->{"x","y","z"}];
SetOptions[ContourPlot,PlotPoints->30,
	AxesLabel->{"x","y"}];

spaceCurve[ fun:{_, _, _}, ul:{_, u0_, u1_, du_}, opts___] :=
	Show[ Graphics3D[Line[Table[N[fun], ul]]], opts,ViewPoint->{5,3,1} ] /;
			NumberQ[N[u0]] && NumberQ[N[u1]] && NumberQ[N[du]]

spaceCurve[fun_, {u_, u0_, u1_}, opts___] :=
Module[{plotpoints},
    plotpoints = PlotPoints /. {opts} /. Options[Plot3D];
	spaceCurve[ fun, {u, u0, u1, (u1-u0)/(plotpoints-1)}, opts ]
    ];

\(\(sliceSurface[z_, xLo_, xHi_, yLo_, yHi_, zLo_, zHi_, opts___] := 
        Module[{web, box, 
            slicePlane}, 
box = {{xLo -  .01, 
                xHi +  .01}, {yLo -  .01, yHi +  .01}, {zLo -  .01, 
                zHi +  .01}}; 
web = {}; 

          Do[AppendTo[web, 
              spaceCurve[{x, y, z}, {y, yLo, yHi}, PlotRange -> box, 
                DisplayFunction -> 
                  Identity]]; 
slicePlane = 
              Graphics3D[
                Polygon[{{x, yLo, zLo}, {x, yHi, zLo}, {x, yHi, zHi}, {x, 
                      yLo, zHi}, {x, yLo, zLo}}]]; 
            Show[slicePlane, web, opts, ViewPoint -> {5, 3, 1}, 
              PlotRange -> box, 
              DisplayFunction -> $DisplayFunction];
	, {x, xLo, 
              xHi, \((xHi - xLo)\)/5}]; 

          Do[AppendTo[web, 
              spaceCurve[{x, y, z}, {x, xLo, xHi}, PlotRange -> box, 
                DisplayFunction -> 
                  Identity]]; 
slicePlane = 
              Graphics3D[
                Polygon[{{xLo, y, zLo}, {xHi, y, zLo}, {xHi, y, zHi}, {xLo, 
                      y, zHi}, {xLo, y, zLo}}]]; 
            Show[web, slicePlane, opts, ViewPoint -> {5, 3, 1}, 
              PlotRange -> box, 
              DisplayFunction -> $DisplayFunction];
	, {y, yLo, 
              yHi, \((yHi - yLo)\)/5}]; 
Show[web, opts, 
            ViewPoint -> {5, 3, 1}, PlotRange -> box, 
            DisplayFunction -> $DisplayFunction];];\)\)

plot3DContour[f_,xCnt_,yCnt_,scale_,opts___]:= 
Module[{zSs,surface,contour,cont,baseContour,combined,crvs3Dplus,fPlus,crvs3Dminus,fMinus},

zSs = Table[N[f[x,y]],{x,xCnt-scale,xCnt+scale,scale/30},
			{y,yCnt-scale,yCnt+scale,scale/30}];
zMIN = Min[zSs];

contour = ContourPlot[ f[x,y] ,{x,xCnt-scale,xCnt+scale},
			{y,yCnt-scale,yCnt+scale},PlotPoints->30,
			DisplayFunction -> Identity];
topo=Show[contour,ContourShading->False,
      ContourStyle->AbsoluteThickness[2],DisplayFunction->Identity];
      			
surface = Plot3D[f[x,y] ,{x,xCnt-scale,xCnt+scale},{y,yCnt-scale,yCnt+scale},opts,ViewPoint->{5,3,1},Mesh->False,BoxRatios->Automatic,DisplayFunction->Identity];

fPlus[x_,y_]=f[x,y]+.02;
crvs3Dplus = Show[Graphics3D[First[Graphics[topo]]/.Line[pts_]:>(val=Apply[fPlus,First[pts]];
Line[Map[Append[#,val]&,pts]]
)],DisplayFunction->Identity];

fMinus[x_,y_]=f[x,y]-.02;
crvs3Dminus = Show[Graphics3D[First[Graphics[topo]]/.Line[pts_]:>(val=Apply[fMinus,First[pts]];
Line[Map[Append[#,val]&,pts]]
)],DisplayFunction->Identity];

cont = Show[contour,ContourShading->False,
				DisplayFunction -> Identity];

baseContour = First[Graphics[cont]];
					
baseContour = Graphics3D[baseContour /. Line[x_]:>
			Line[Map[Append[#,zMIN]&,x]]];

combined = Show[surface,baseContour,crvs3Dplus,crvs3Dminus,DisplayFunction->$DisplayFunction];
Return[combined]
		];

sliceAt[f_,xCnt_,yCnt_,scale_,z_,zMIN_] :=
	Module[{surface,plane,contour,pLz,contLz,w},

	surface = Plot3D[ f[x,y] , {x,xCnt-scale,xCnt+scale},
						{y,yCnt-scale,yCnt+scale},
				PlotPoints->30,BoxRatios->Automatic,Mesh->False,
				AxesLabel->{"x","y","z"},
				ViewPoint->{5,3,1},
				DisplayFunction->Identity];


	plane[w_] := Graphics3D[Polygon[{{xCnt-scale,yCnt-scale,w},{xCnt+scale,yCnt-scale,w},{xCnt+scale,yCnt+scale,w},{xCnt-scale,yCnt+scale,w},{xCnt-scale,yCnt-scale,w}}]];
						
	contour[w_] := ContourPlot[ f[x,y] ,
		{x,xCnt-scale,xCnt+scale},{y,yCnt-scale,yCnt+scale},
		PlotPoints->30, Contours->{w},DisplayFunction->Identity];

				
	pLz = plane[z];

	contLz = Show[Evaluate[contour[z]],
			ContourShading->False,
			DisplayFunction -> Identity];

	contLz = First[Graphics[contLz]];
					
	contLz = Graphics3D[contLz /. Line[x_]:>
			Line[Map[Append[#,zMIN]&,x]]];

Show[surface,pLz,contLz,
		DisplayFunction->$DisplayFunction];]

Cell[TextData[{
  "Section ",
  
  CounterBox["Chapter"],
  ".1: Coordinates in 3D"
}], "Subsection"]

Section 0.4   Plots of explicitly specified functions in 3-D

Definition .4.1:  Graph of an explicitly defined function in 3D

Mathematica can draw plots of of surfaces, so that we can easily visualize them.

A circular paraboloid

<< Default3D`
f[x_, y_] := (x^2) + (y^2)
a =  -2;
b = 2;
c =  -2;
d   =  2;
Plot3D[f[x, y], {x, a, b}, {y, c, d},  Boxed  ->  True,  
        Axes ->  True];

The same, but with real time viewpoint selection

<< RealTime3D`
f[x_, y_] := (x^2) + (y^2)
a =  -2;
b = 2;
c =  -2;
d   =  2;
Plot3D[f[x, y], {x, a, b}, {y, c, d},  Boxed  ->  True,  
        Axes ->  True];

The same, but with  circular domain

It's a little unnatural to have the "handkerchief corners" on the graph of the circular parabola.
Here's a trick to draw the circular parabola over a circular domain

<< Default3D`
<< RealTime3D`

CircularPlot3D[f_,  R_,  
        options___]  :=  
ParametricPlot3D[  ({x,  
              y,  f[x,  y]})  /.  {x  ->  r Cos[θ],  
            y  ->  r Sin[θ]},  
{r,  0,  
          R},  {θ,  0,  2 π},  
        options]
f[x_, y_] := (x^2) + (y^2)

CircularPlot3D[f, 3,  Boxed  ->  True,  
      Axes ->  True]

A hyperboloid

<< Default3D`
f[x_, y_] := (x^2) -  (y^2)
a =  -2;
b = 2;
c =  -2;
d   =  2;
Plot3D[f[x, y], {x, a, b}, {y, c, d},  Boxed  ->  True,  
        Axes ->  True];

The same, but with real time viewpoint selection

<< RealTime3D`
f[x_, y_] := (x^2) -  (y^2)
a =  -2;
b = 2;
c =  -2;
d   =  2;
Plot3D[f[x, y], {x, a, b}, {y, c, d},  Boxed  ->  True,  
        Axes ->  True];

The same, but with circular domain

<< RealTime3D`
f[x_, y_] := (x^2) -  (y^2)
CircularPlot3D[f,  1,  Boxed ->  False]

Something more complicated

<< RealTime3D`
f[x_, y_] = 
        Sin[x y] (E^(-(x^2) - (y^2)));
a =  -3;
b = 3;
c =  -3;
d   =  3;
Plot3D[f[x, y], {x, a, b}, {y, c, d},  Boxed  ->  True,  
        Axes ->  True];

Section 0.5  Contour Plots in 3-D by Horizontal Slices

Another way to represent the graph of a function  f[x, y]  is to draw the curves  c  =  f[x, y], where  c  is constant, for various values of c. 

When  z  represents height, the curve  c  =  f[x, y]  is the set of (x, y)-points where the  surface  z  =  f[x, y] has constant height  z  =  c. 

This is the method used to make topographical maps.  The graph of the geographic elevation function is the landscape itself!  A topographic map is the contour plot of the elevation function.

Definition ..: Contour Graph

A Circular Paraboloid

The function  f[x, y] = (x^2) + (y^2)  is graphed explicitly below with curves at various  z-heights shown on the surface and in a base plane below the surface.  The collection of curves in the base plane is the contour graph of the function.

Clear[f, x, y];
f[x_, y_] = (x^2) + (y^2);
xCnt = 0;
yCnt = 0;
scale = 2;
plot3DContour[f, xCnt, yCnt, scale, Ticks -> None, 
        ViewPoint -> {5, 2, 2}];

We illustrate the construction of a countour plot on the function  f[x, y] = (x^2) + (y^2).  We begin by fixing the "level" at  z  =  1, so  1  =  f[x, y]  is the equation  1 = (x^2) + (y^2).  This is an implicit circle of radius  1  centered at the origin.

Clear[f, x, y];
f[x_, y_] = (x^2) + (y^2);
xCnt = 0;
yCnt = 0;
scale = 2;
z = 1;
sliceAt[f, xCnt, yCnt, scale, z, -1]

Next, we set  z  =  2, so the equation is  2 = (x^2) + (y^2).  This is an implicit circle of radius  (Sqrt[2]).  Continue with various constants, say z = 1, 2, 3, 4, 5.

z = 2;
sliceAt[f, xCnt, yCnt, scale, z, -1]

The contour graph of  f[x, y] = (x^2) + (y^2)  is shown next.

ContourPlot[x^2 + y^2, {x, -2, 2}, {y, -2, 2}, 
        ContourShading -> False];

Something more complicated

Here is the contour plot of f[x, y] = x y (e^(-(x^2) - (y^2))).  The plot on the right shows the graph of the function in 3D, with several curves of constant height, as well as the projections of these curves on the x--y plain.  The graph on the left shows the (2D) contour plot of the same function, shaded according to the height; it conveys, or encodes the same information as the graph plot on the right.

Clear[f, x, y];
f[x_, y_] = 
        9  x y Exp[-((x^2) + (y^2))];
xCnt = 0;
yCnt = 0;
scale = 2;
fig1 = 
        plot3DContour[f, xCnt, yCnt, scale, Ticks -> None, 
          ViewPoint -> {5, 2, 2}];
fig2 = 
        ContourPlot[
          f[x, y], {x, xCnt - scale, xCnt + scale}, {y, yCnt - scale, 
            yCnt + scale}, FrameLabel -> {"x", "y"}, 
          FrameTicks -> {{-1, 1}, {-1, 1}}, 
          RotateLabel -> False];
Show[GraphicsArray[{fig2, fig1}]];