(* ::Package:: *)

(* :Name: sccc` *)

(* :Title: sccc Package "sccc" *)
(* :VERSION  1.5, Feb. 2008 *)
(* :Author: Eric Schulz , Mike Pepe
*)

(* :Summary: 
This package provides additional functions for SCCC students.  
*)

(* :Mathematica Version: 6.0 *)


BeginPackage["sccc`"]

DrawVector::usage ="DrawVector[{u,v}] plots the 2D vectors u and v in standard position."
DrawVector3D::usage="DrawVector3D[{u,v}] plots the 3D vectors u and v in standard position."
Axes3D::usage="Axes3D[n] draws three axes each from -n to n."
System3D::usage="System3D[n] draws three coordinate planes from -n to n on each axis."
Point3D::usage="Point3D[u] plots the point u (or list of points).Options: xyPlane->True shows projection onto xy-plane; xyAxes->True shows coordinate lines."
Line3D::usage="Line3D[a,b] draws line from point a to b."
Taylor::usage="Taylor[f[x],{x,a,n}] Taylor poly of degree n for f centered at x=a."
ClearNames::usage="ClearNames[] clears all user defined names."
TF::usage="TF is shortname for the TraditionalForm command."
Mag::usage="Mag[v] computes the magnitude of vector v."
Lim::usage="Lim[f[x],x->number] computes the limit of f[x]. A third option argument of either Left or Right can be specified to compute a one sided limit."
EquationPlot::usage="EquationPlot[eqn,{x,xmin,xmax},{y,ymin,ymax}] graphs implicit plot for equation in two variables"
EquationPlot3D::usage="EquationPlot3D[eqn,{x,xmin,xmax},{y,ymin,ymax},{z,zmin,zmax}] graphs implicit plot for equation in three variables"


(* ::Input:: *)
(**)


Options[DrawVector] = {
   Color -> RGBColor[1, 0, 0], 
   VectorLabel -> None,
   HeadLength -> Large, 
   ShaftThickness -> AbsoluteThickness[2],
   AspectRatio -> Automatic,
   Axes -> True, 
   PlotRange -> All};
Options[DrawVector3D] = { 
	Color -> RGBColor[1,0,0], 
	ShaftColor -> RGBColor[0,0,0], 
	VectorLabel -> None,
	HeadLength -> 0.2,
	HeadAngle -> 0.2, 
	Axes->True,
	PlotRange->All};
Options[System3D]={
	Solid->False,
	Planes->True,
	xyzAxes->True,
	AxesColor->Black,
	AxesTicks->False,
	PlaneColor->Lighter[Gray],
	PlaneGrids->False,
	GridColor->LightGray,
	ViewPoint->{2,2,1.9}};
Options[Point3D]={
	Color->Blue,
	Size->6,
	Axes->False,
	Boxed->False,
	CoordinateLinesColor->Black,
	CoordinateLines->AbsoluteThickness[1.5],
	xyPlane->False,
	xyAxes->False,
	ViewPoint->{2,2,1.9}};


Begin["`Private`"]

Print[Style["sccc package version 1.5.",Italic,11,Gray]]

Print[Style["Axes3D, ClearNames, DrawVector, DrawVector3D, EquationPlot, EquationPlot3D,Lim, Line3D, Mag, Point3D, System3D,Taylor, and TF. 
For command description type ?command.",Italic,11,Gray]]


Axes3D[n_] := Module[
{xMin, xMax, yMin, yMax, zMin, zMax},
If[ IntegerQ[n], 
	(*T*){xMin, xMax, yMin, yMax, zMin, zMax} = {-n, n, -n, n, -n, n}];
If[ VectorQ[n], 
	(*T*){xMin, xMax, yMin, yMax, zMin, zMax} = {-n[[1]], n[[1]], -n[[2]], n[[2]], -n[[3]], n[[3]]}];
DrawVector3D[
	{{xMin, 0, 0}, {xMax, 0, 0}, {0, yMin, 0}, {0, yMax,0}, {0, 0, zMin}, {0, 0, zMax}}, 
	Color -> Blue, HeadLength -> 0, ShaftColor -> Blue, VectorLabel -> {"", "x", "", "y", "", "z"}, 
	ViewPoint -> {2,2,1.9}, Boxed -> False, Axes -> False]]

Line3D[u_,v_]:=
Module[{},
Show[
	Graphics3D[ {Thickness[.005],Line[{u,v}]} ],
	Boxed->False, Axes->False, ViewPoint -> {2,2,1.9}]
]

Point3D[p_,opts___]:=
Module[{pointColor,size,validOptions,lineColor,lines,xyplane,xyaxes},
validOptions=FilterRules[Flatten[{opts,Options[Point3D]}],Options[Graphics3D]];
{pointColor,size,lineColor,lines,xyplane,xyaxes}={Color,Size,CoordinateLinesColor,CoordinateLines,xyPlane,xyAxes}/.{opts}/.Options[Point3D];
Show[Graphics3D[{
	pointColor,
	AbsolutePointSize[size],
	If[	Dimensions[Dimensions[p]]=={1},
		(*T*)Point[p],
		(*F*)Map[(Point[#]&),p]],
	If[	xyplane==True||xyaxes==True,
		(*T*)Flatten[{lineColor,lines,
		If[	Dimensions[Dimensions[p]]=={1},
			(*T*)Line[{{1,1,0}p,{1,1,1}p}],
			(*F*)Map[(Line[{{1,1,0}#,{1,1,1}#}]&),p]]}],
		(*F*)Black],
	If[	xyaxes==True,
		(*T*)Flatten[{lineColor,lines,
		If[	Dimensions[Dimensions[p]]=={1},
			(*T*)Line[{{1,0,0}p,{1,1,0}p,{0,1,0}p}],
			(*F*)Map[(Line[{{1,0,0}#,{1,1,0}#,{0,1,0}#}]&),p]]}],
		(*F*)Black]
	}],
Evaluate[validOptions]]
	]

System3D[n_, opts___] := 
Module[{xMin, xMax, yMin, yMax, zMin, zMax, solid, validOptions, planes, xyzaxes, planecolor, axescolor, axesticks, ticksize, planegrids, gridcolor}, 
validOptions = FilterRules[Flatten[{opts, Options[System3D]}], Options[Graphics3D]];
{solid, planes, xyzaxes, planecolor, axescolor, axesticks, planegrids, gridcolor} = {Solid, Planes, xyzAxes, PlaneColor, AxesColor, AxesTicks, PlaneGrids, GridColor} /. {opts} /. Options[System3D];
If[ IntegerQ[n],
	(*T*){xMin, xMax, yMin, yMax, zMin, zMax} = {-n, n, -n, n, -n, n}];
If[ Dimensions[n] == {3, 2},
	(*T*){{xMin, xMax}, {yMin, yMax}, {zMin, zMax}} = n,
	If[ Dimensions[n] == {3},
		(*T*)
		If[Dimensions[n[[1]]] === {2}, {xMin, xMax} = n[[1]], {xMin, xMax} = {-n[[1]], n[[1]]}];
		If[Dimensions[n[[2]]] === {2}, {yMin, yMax} = n[[2]], {yMin, yMax} = {-n[[2]], n[[2]]}];
		If[Dimensions[n[[3]]] === {2}, {zMin, zMax} = n[[3]], {zMin, zMax} = {-n[[3]], n[[3]]}]]];
ticksize = 0.035 Max[{xMax, yMax, zMax}];
Show[Graphics3D[{
		planecolor, 
		Opacity[If[solid, 1, 0.5]], 
		If[ planes == True,
			(*T*){
			Polygon[{{xMin, yMin, 0}, {xMin, yMax, 0}, {xMax, yMax, 0}, {xMax, yMin, 0}}], 
			Polygon[{{0, yMin, zMin}, {0, yMin, zMax}, {0, yMax, zMax}, {0, yMax, zMin}}], 
			Polygon[{{xMin, 0, zMin}, {xMin, 0, zMax}, {xMax, 0, zMax}, {xMax, 0, zMin}}]},
			(*F*)Black], 
		Opacity[1], 
		If[ xyzaxes == True,
			(*T*){
			axescolor, AbsoluteThickness[2], 
			Line[{{xMin, 0, 0}, {xMax, 0, 0}}], 
			Line[{{0, yMin, 0}, {0, yMax, 0}}], 
			Line[{{0, 0, zMin}, {0, 0, zMax}}], 
			Text[Style["x", Bold], {1.1 xMax, 0, 0}], 
			Text[Style["y", Bold], {0, 1.1 yMax, 0}], 
			Text[Style["z", Bold], {0, 0, 1.1 zMax}]},
			(*F*)Black], 
		If[ axesticks == True,
			(*T*){
			axescolor, AbsoluteThickness[1.5], 
			Map[(Line[{{#1, -ticksize, 0}, {#1, ticksize, 0}}]) &, Range[Round[xMin], Round[xMax], 1]], 
			Map[(Line[{{-ticksize, #1, 0}, {ticksize, #1, 0}}]) &, Range[Round[yMin], Round[yMax], 1]], 
			Map[(Line[{{0, -ticksize, #1}, {0, ticksize, #1}}]) &, Range[Round[zMin], Round[zMax], 1]]},
			(*F*)Black], 
		If[ planegrids == True,
		(*T*){
		gridcolor, AbsoluteThickness[1], 
		Map[(Line[{{#1, yMin, 0}, {#1, yMax, 0}}]) &, Cases[Range[Round[xMin], Round[xMax], 1], Except[0]]], 
		Map[(Line[{{xMin, #1, 0}, {xMax, #1, 0}}]) &, Cases[Range[Round[yMin], Round[yMax], 1], Except[0]]], 
		Map[(Line[{{#1, 0, zMin}, {#1, 0, zMax}}]) &, Cases[Range[Round[xMin], Round[xMax], 1], Except[0]]], 
		Map[(Line[{{0, #1, zMin}, {0, #1, zMax}}]) &, Cases[Range[Round[yMin], Round[yMax], 1], Except[0]]], 
		Map[(Line[{{0, yMin, #1}, {0, yMax, #1}}]) &, Cases[Range[Round[zMin], Round[zMax], 1], Except[0]]], 
		Map[(Line[{{xMin, 0, #1}, {xMax, 0, #1}}]) &, Cases[Range[Round[zMin], Round[zMax], 1], Except[0]]]},
		(*F*)Black]
		}], 
	Evaluate[validOptions], 
	Boxed -> False, 
	Axes -> False,
    BaseStyle -> 14]]

EquationPlot[eq_, {xvar_, xmin_, xmax_}, {yvar_, ymin_, ymax_}, opts___] :=
 Module[
 {validOptions},
 validOptions = FilterRules[{opts}, Options[ContourPlot]];
 ContourPlot[eq, {xvar, xmin, xmax}, {yvar, ymin, ymax}, 
 	Evaluate[validOptions], ContourStyle -> Thick, 
 	AspectRatio -> Automatic, Axes -> True, Frame -> False]
  ]

EquationPlot3D[eq_, {xvar_, xmin_, xmax_}, {yvar_, ymin_, ymax_}, {zvar_, zmin_, zmax_}, opts___] :=
Module[
{validOptions},
validOptions = FilterRules[{opts}, Options[ContourPlot3D]];
ContourPlot3D[eq, {xvar, xmin, xmax}, {yvar, ymin, ymax}, {zvar, zmin, zmax}, 
	Evaluate[validOptions],ViewPoint->{2,2,1.9}]
  ]

Mag[v_] := Sqrt[v.v]
Taylor[a_,b_]:=Normal[Series[a,b]]
ClearNames[]:=Module[{},Map[Remove,Delete[Names["Global`*"],Position[Names["Global`*"],"ClearNames"]]];]

TF = TraditionalForm

SetOptions[Plot,AxesOrigin->{0,0}]
SetOptions[ListPlot, AxesOrigin -> {0, 0}]
SetOptions[ContourPlot, AxesOrigin -> {0, 0}]
SetOptions[ParametricPlot, AxesOrigin -> {0, 0}]
SetOptions[PolarPlot, AxesOrigin -> {0, 0}]
SetOptions[Plot3D,ViewPoint -> {2,2,1.9}]
SetOptions[ContourPlot3D,ViewPoint -> {2,2,1.9}]

(* Start of Lim *)
Lim[a_, b_, dir_: Both] :=
Module[
{LimLeft, LimRight, LimBoth},
LimLeft[aa_, bb_] :=
	If[Head[bb[[2]]] === Symbol,
    	(*T*)"Second argument should be in the form of x \[Rule] n where n is a real number or \[PlusMinus]\[Infinity].",
    	(*F*)
	    If[ bb[[2]] \[Element] Reals ,
    		(*T*)
    		If[ Head[(aa /. (bb[[1]]) -> (bb[[2]] - .00000001))] =!= Complex,
				(*T*)
					If[ Head[Limit[aa, bb, Direction -> 1]] === Interval,						(*T*)"Limit does not exist",
						(*F*)Limit[aa, bb, Direction -> 1]],
	(*F*)"Limit does not exist"]
     	]];
LimRight[aa_, bb_] :=
	If[ Head[bb[[2]]] === Symbol,
		(*T*)"Second argument should be in the form of x \[Rule] n where n is a real number or \[PlusMinus]\[Infinity].",
		(*F*)
		If[ bb[[2]] \[Element] Reals ,
		(*T*)
		If[Head[(aa /. (bb[[1]]) -> (bb[[2]] + .00000001))] =!= Complex,
			(*T*)
				If[ Head[Limit[aa, bb, Direction -> -1]] === Interval,
						(*T*)"Limit does not exist",
						(*F*)Limit[aa, bb, Direction -> -1]],
      		(*F*)"Limit does not exist"]
		]];
LimBoth[aa_, bb_] :=
	If[ Head[bb[[2]]] === Symbol,
		(*T*)"Second argument should be in the form of x \[Rule] n where n is a real number or \[PlusMinus]\[Infinity].",
		(*F*)
		If[bb[[2]] \[Element] Reals ,
			(*T*)
			If[ Head[(aa /. (bb[[1]]) -> (bb[[2]] - .00000001))] =!= Complex && 
				Head[(aa /. (bb[[1]]) -> (bb[[2]] + .00000001))] =!= Complex &&
				LimLeft[aa, bb] == LimRight[aa, bb],
				(*T*)
				If[ Head[Limit[aa, bb]] === Interval,
						(*T*)"Limit does not exist",
						(*F*)Limit[aa, bb]],
				(*F*)"Limit does not exist"],
			(*F*)
				If[ Head[Limit[aa, bb]] === Interval,
						(*T*)"Limit does not exist",
						(*F*)Limit[aa, bb]]
		]];
  Which[
   dir === Both, LimBoth[a, b],
   dir === Left, LimLeft[a, b],
   dir === Right, LimRight[a, b]]]
(* End Lim *)

(* Start DrawVector *)
DrawVector[v_, opts___] :=
Module[
{veclist, colors, vectors, c, label, aspect, plotrange,optPlot,hs,thickness},
{c, label,hs,thickness} = {Color, VectorLabel,HeadLength,ShaftThickness} /. {opts} /. Options[DrawVector];
optPlot = FilterRules[{opts,Options[DrawVector]},Options[Plot]];
If[
	Dimensions[Dimensions[v]] == {1},
	(* Code used for a single vector *)
	Show[
		Graphics[{
			Arrowheads[hs],
			c,
			thickness,
			If[label === None, Text["", 1.1v], Text[label, 1.1v]],
			Arrow[{{0, 0}, v}]}],
		Evaluate[optPlot],Axes -> True, PlotRange -> All,BaseStyle->12],

	(* Code used for a list of vectors *)
	vectors = Dimensions[v][[1]];
	If[Head[c] === RGBColor,colors = Table[c, {i, 1, vectors}],colors = c];
	Show[
		Graphics[ 
			Flatten[
				Table[{
					thickness,
					(* Set color for vector and label *) 
					colors[[i]],
					(* Set label for vector *)
					If[
						Length[Dimensions[v[[i]]]] == 1,
						(* TRUE *)
						If[
							label === None,
							(* TRUE *)
							Text["", 1.1*Flatten[v[[i]]]],
							(* FALSE *)
							Text[label[[i]], 1.1*Flatten[v[[i]]]]
							],
						If[
							label === None,
							(* TRUE *)
							Text["", 1.1*Flatten[v[[i, 2]]]],
							(* FALSE *)
							Text[label[[i]], 1.1*Flatten[v[[i, 2]] - v[[i, 1]]] + v[[i, 1]]]]
						],
					Arrowheads[hs],
					(* Arrow representing vector *)
					Arrow[{
						(* TAIL *) If[Length[Dimensions[v[[i]]]] == 1, {0, 0}, Flatten[v[[i, 1]]]],
						(* HEAD *) If[Length[Dimensions[v[[i]]]] == 1, Flatten[v[[i]]], Flatten[v[[i, 2]]]]}]
					},
					{i, 1, vectors}]
				]
			],
		Evaluate[optPlot],Axes -> True, PlotRange -> All,BaseStyle->12]
	]
]
(* End DrawVector *)

(* Start DrawVector3D *)
DrawVector3D[v_,opts___]:=
Module[
{headlength, headangle, headcolor, shaftcolor, labelfontsize,h11, h12, h21, h22, vectors, vectorList, vectorListLive, vOne,optPlot3D, vectorMagnitudes, avgMagnitude, absHeadLength,  vectorLabel,w},
{headlength, headangle, headcolor, shaftcolor, vectorLabel} = {HeadLength, HeadAngle, Color, ShaftColor,  VectorLabel} /. {opts} /. Options[DrawVector3D];
optPlot3D = FilterRules[{ opts, Options[DrawVector3D]},Options[Plot3D]];

(* Determine how many vectors are to be plotted *)
If[	
	First[ Dimensions[v] ] == 1 || (Length[Dimensions[v]]==1 && First[Dimensions[v]]==3 && Length[Flatten[v]]<=6) || Length[ v ] == 1,
	(* TRUE *)
	vectors = 1; vOne = Flatten[v] (* Support inputs of the form v and {v} *),
	(* FALSE *)
	vectors = Length[v] 
	];

(* Determine the average magnitude of the vectors to be plotted *)
If[
	vectors == 1,
	(* TRUE *)
	avgMagnitude = Sqrt[vOne.vOne],
	(* FALSE *)
	vectorMagnitudes = Table[ 
		If[ 
			Dimensions[Dimensions[v[[i]]]] == {1}, 
			(* TRUE *)
			Sqrt[ v[[i]].v[[i]] ],
			(* FALSE *)
			Sqrt[ (v[[i,2]]-v[[i,1]]).(v[[i,2]]-v[[i,1]]) ]
			],
		{i,1,vectors}
		];
		avgMagnitude = Apply[ Plus, vectorMagnitudes ] / Length[ vectorMagnitudes ]
	];

(* Determine the headlength to be used for all vector heads *)
absHeadLength = headlength * avgMagnitude;

(* Modify absHeadLength if it is greater than the magnitude of the smallest vector *)
If[
	absHeadLength > Min[vectorMagnitudes],
	(* TRUE *)
	absHeadLength = .8 * Min[vectorMagnitudes],
	(* FALSE *)
	Null];

If[
	vectors == 1,
	(* TRUE *)
	If[
		Dimensions[v] == {3} || Dimensions[v] == {1,3},
		(* TRUE *)
		(* Make call to HeadPoints to determine the four points which make up the base of the vector's head *)
		{h11,h12,h21,h22} = HeadPoints[vOne,absHeadLength,headangle];
		vectorList = {shaftcolor,Line[{{0,0,0},vOne}],FaceForm[headcolor],Polygon[{vOne,h11,h12,vOne}],Polygon[{vOne,h21,h22,vOne}]};
		If[
			vectorLabel =!= None, 
			(* TRUE *)
			AppendTo[ vectorList, Text[vectorLabel, 1.1vOne]]
			],
		(* FALSE *)
		w = v[[1,2]]-v[[1,1]];
		{h11,h12,h21,h22} = HeadPoints[w,absHeadLength,headangle];
		vectorList = {shaftcolor, Line[{v[[1,1]],v[[1,2]]}], FaceForm[headcolor],Polygon[{v[[1,2]],h11+v[[1,1]],h12+v[[1,1]],v[[1,2]]}],Polygon[{v[[1,2]],h21+v[[1,1]],h22+v[[1,1]],v[[1,2]]}]};
		If[
			vectorLabel =!= None,
			(* TRUE *)
			AppendTo[	vectorList,Text[vectorLabel, 1.1 w + v[[1,1]]]]
			]
		];
	Show[Graphics3D[vectorList],Evaluate[optPlot3D],AxesLabel -> {"x","y","z"}],
	(* FALSE *)
	If[
		Head[headcolor] === RGBColor,
		(* TRUE *)
		headcolor = Table[headcolor,{i,1,vectors}],
		(* FALSE *)
		Null];
	If[
		Head[shaftcolor] === RGBColor,
		(* TRUE *)
		shaftcolor = Table[shaftcolor,{i,1,vectors}],
		(* FALSE *)
		Null];
	vectorList = {};
	Table[
		If[
			Dimensions[v[[i]]] == {3},
			(* TRUE *)
			(* Make call to HeadPoints to determine the four points which make up the base of the vector's head *)
			{h11,h12,h21,h22} = HeadPoints[v[[i]],absHeadLength,headangle];
			AppendTo[
				vectorList,
				{	shaftcolor[[i]], 
					Line[{{0,0,0},v[[i]]}],
					FaceForm[headcolor[[i]]], 
					Polygon[{v[[i]],h11,h12,v[[i]]}],
					Polygon[{v[[i]],h21,h22,v[[i]]}]
					}
				];
			If[
				vectorLabel =!= None,
				(* TRUE *)
				AppendTo[
					vectorList,
					Text[ vectorLabel[[i]], 1.15 v[[i]]]
					]
				],
			(* FALSE *)
			w = v[[i,2]]-v[[i,1]];
			(* Make call to HeadPoints to determine the four points which make up the base of the vector's head *)
			{h11,h12,h21,h22} = HeadPoints[w,absHeadLength,headangle];
			AppendTo[
				vectorList,
				{	shaftcolor[[i]],
					Line[{v[[i,1]],v[[i,2]]}],
					FaceForm[headcolor[[i]]],
					Polygon[{v[[i,2]],h11+v[[i,1]],h12+v[[i,1]],v[[i,2]]}],
					Polygon[{v[[i,2]],h21+v[[i,1]],h22+v[[i,1]],v[[i,2]]}]}
					];
			If[
				vectorLabel =!= None,
				(* TRUE *)
				AppendTo[
					vectorList,
					Text[vectorLabel[[i]], 1.1 w + v[[i,1]]]
					]
				]
			],
		{i,1,vectors}];
	Show[Graphics3D[vectorList],Evaluate[optPlot3D],AxesLabel -> {"x","y","z"},Axes -> True, PlotRange -> All,BaseStyle->12],
$Failed]
]
(* End DrawVector3D *)

(* Start HeadPoints *)
(* HeadPoints is used by DrawVector3D to build the arrowhead for each vector *)
HeadPoints[v_,headlength_,headangle_]:=
Module[{norm,xangle,zangle,h11,h12,h21,h22},

(* Convert v from rectangular coordinates to spherical coordinates *)
norm=N[Sqrt[v.v]];
xangle = If[v[[1]] != 0, N[ArcTan[ v[[2]]/v[[1]] ]] + If[ v[[1]] < 0, Pi, 0], If[ v[[2]] >= 0, Pi/2, -Pi/2]];
zangle = If[ v[[3]] != 0, N[ArcTan[ Sqrt[ v[[1]]^2+v[[2]]^2 ] / v[[3]] ]] + If[ v[[3]] < 0, Pi, 0], Pi/2];

(* Modify phi and theta in spherical coordinates to compute the base points for the vector head *)
h11= v - headlength * {Sin[zangle+headangle]*Cos[xangle],Sin[zangle+headangle]*Sin[xangle],Cos[zangle+headangle]};
h12= v - headlength *{Sin[zangle-headangle]*Cos[xangle],Sin[zangle-headangle]*Sin[xangle],Cos[zangle-headangle]};
h21= v - headlength * {Sin[zangle]*Cos[xangle+headangle],Sin[zangle]*Sin[xangle+headangle],Cos[zangle]};
h22= v - headlength * {Sin[zangle]*Cos[xangle-headangle],Sin[zangle]*Sin[xangle-headangle],Cos[zangle]};

(* Return the four base points for the vector head *)
{h11,h12,h21,h22}
]
(* End HeadPoints *)
End[]
EndPackage[]