Setup:

Clear["Global`*"];
Off[ General::spell ];
Off[ General::spell1];
Off[Clear::wrsym];

Load Packages:

Needs["LinearAlgebra`Orthogonalization`"]
?LinearAlgebra`Orthogonalization`*

Needs["Calculus`VectorAnalysis`"]
?Calculus`VectorAnalysis`*

Needs["Graphics`PlotField`"]
?Graphics`PlotField`*

Needs["Graphics`PlotField3D`"]
?Graphics`PlotField3D`*

Initialization: 
(Don't read this; just execute it)
(Click on the outer bracket--far right--and hit the ENTER key.

Unprotect[{VectorPlot,VectorPlot3D,VectorLength}];

VectorPlot[list__List?VectorQ]:= VectorPlot[{list}];

VectorPlot[list_List?MatrixQ]:= 
Module[{length,totList},
 length=Length[list];
 totList=Table[{{0,0},list[[i]]} ,{i,1,length}];
 Return[
  Show[Graphics[{Thickness[0.015],Line /@ totList }
          ,{GridLines->Automatic
           ,Axes->True
           ,AspectRatio->1}]]
 ]
];

VectorPlot3D[list__List?VectorQ]:= VectorPlot3D[{list}];

VectorPlot3D[list_List?MatrixQ]:= 
Module[{length,totList},
 length=Length[list];
 totList=Table[{{0,0,0},list[[i]]} ,{i,1,length}];
 Return[
  Show[Graphics3D[{Thickness[0.015],Line /@ totList}
          ,{ FaceGrids->{{0,0,-1},{0,1,0},{-1,0,0}}
            ,BoxRatios->{1,1,1}
            ,Axes->True}]]
 ]
];

VectorLength[vector_List?VectorQ]:=Sqrt[vector.vector]

Protect[{VectorPlot,VectorPlot3D,VectorLength}];

2-D Example:

v1={1,0} //Normalize
v2={1,1} //Normalize
v3={0,1} //Normalize
VectorPlot[v1,v2,v3];

VectorLength /@ {v1,v2,v3,{1,1}}

ArcCos[v1.v2]

ArcCos[v1.v2] /Degree //N

3-D Example:

v1={1,0,0} //Normalize
v2={1,1,0} //Normalize
v3={0,1,0} //Normalize
VectorPlot3D[v1,v2,v3];

VectorLength /@ {v1,v2,v3,{1,1,1}}

CrossProduct[v1,v2] //Normalize

VectorPlot3D[v1,v2,CrossProduct[v1,v2]];

{CrossProduct[v1,v2].v1,CrossProduct[v1,v2].v2}

vx={1,0,0};
vy={0,1,0};
vz={0,0,1};
VectorPlot3D[vx,vy,vz];
{CrossProduct[vx,vy].vz,CrossProduct[vz,vy].vx}

Gradient:

?Grad

Plot3D[x^2+y^2+z^2 /.{z->0} 
      ,{x,0,1}
      ,{y,0,1}
      ,BoxRatios->{1,1,1}
      ];

ContourPlot[- (x^2+y^2+z^2) /.{z->0} //Evaluate
      ,{x,0,1}
      ,{y,0,1}
      ];

Grad[x^2+y^2+z^2]

PlotGradientField[x^2+y^2+z^2 /.{z->0} //Evaluate
   ,{x,0,1}
   ,{y,0,1}
   ,VectorHeads->True
   ,PlotPoints->4
   ,Axes->True
   ];

PlotGradientField3D[x^2+y^2+z^2
   ,{x,0,1}
   ,{y,0,1}
   ,{z,0,1}
   ,VectorHeads->True
   ,PlotPoints->4
   ,Axes->True
   ];

Div:

Curl:

Multiple Operations:

Grad[x^2+y^2+z^2]

Div[{x,0,0}]

Curl[{y,-x,0}]

Dirac Delta:

f[1]= n/Sqrt[Pi] Exp[-n^2 x^2];
f[2]=n/Sqrt[Pi] 1/(1+n^2 x^2);
f[3]=Sin[n x]/(Pi x);

Do[

Print["================================"];

ftemp=f[i]; Print["f[",i,"] =  ",ftemp];

Plot[ Table[ ftemp ,{n,2,10,2}] //Evaluate 
  ,{x,-1,1}
  ,PlotRange->All];
  
iftemp= Integrate[ftemp,x]; Print["int-f[",i,"] =  ",iftemp];

Plot[ Table[ iftemp ,{n,2,10,2}] //Evaluate 
  ,{x,-1,1}
  ,PlotRange->All];
  
,{i,1,3}]