Mathematica处理若干个点电荷的电位电场的一个程序包

来自 Mathematica for Theoretical Physics Electrodynamics, Quantum Mechanics, General Relativity and Fractals by Gerd Baumann这本书

程序包

PointCharge.wl(修改版，适用Mathematica12.0)

做的修改有：
删除了Needs["Graphics`PlotField`"];，12.0版本已经内置这些
弃用GradientFieldPlot这一过时的函数，改为使用StreamPlot函数

BeginPackage["PointCharge`"];
(* --- load additional standard packages --- *)
Clear[Potential,Field,EnergyDensity,FieldPlot];
(* --- export functions --- *)
Potential::usage = "Potential[coordinates_List]
creates the potential of
an assembly of point charges. The cartesian
coordinates of the locations of
the charges are given in the form of 
{{x,y,z,charge},{x,y,z,charge},...}.";
Field::usage = "Field[coordinates_List] calculates
the electric field for
an ensemble of point charges. The cartesian
coordinates are
lists in the form of {{x,y,z,charge},{...},...}.";
EnergyDensity::usage =
"EnergyDensity[coordinates_List] calculates the
density of the energy for an ensemble of point
charges. The cartesian
coordinates are lists in the form of
{{x,y,z,charge},{...},...}.";
FieldPlot::usage =
"FieldPlot[coordinates_List,typ_,options___] creates
a
ContourPlot for an ensemble of point charges. The
plot type (Potential,
Field, or Density) is specified as string in the
second input variable. The
third argument allows a change of the Options of
ContourPlot and
PlotGradientField.";
(* --- define the global variables x,y,z --- *)
x::usage;
y::usage;
z::usage;
Begin["`Private`"];
(* --- determine the potential --- *)
Potential[coordinates_List]:=
 Block[{x,y,z},
 Fold[Plus,0,Map[(#[[4]]/Sqrt[(x-#[[1]])^2 +
(y-#[[2]])^2 +
(z-#[[3]])^2])&, coordinates]]];
(* --- calculate the field ---*)
Field[coordinates_List]:=
 Block[{field,x,y,z},
 field = -
Fold[Plus,0,Map[(#[[4]]*({x,y,z}-Take[#,3])/
 (Sqrt[(x-#[[1]])^2 +
(y-#[[2]])^2 +
(z-#[[3]])^2
])^3)&,coordinates]];
 Simplify[field]
 ];
(* --- calculate the energy --- *)
EnergyDensity[coordinates_List]:=
 Block[{density,x,y,z,field},
 field = Field[coordinates];
 density = field.field/(8*Pi)
 ];
(* --- create plots --- *)
FieldPlot[coordinates_List,typ_,options___]:=
 Block[
 {pot, ncharges, xmin, xmax, zmin, zmax, xcoord
= {}, zcoord = {},
 pl1, pl2},
 ncharges = Length[coordinates];
(* --- determine limits for the plot --- *)
 Do[
 AppendTo[xcoord,coordinates[[i,1]]];
 AppendTo[zcoord,coordinates[[i,3]]],
 {i,1,ncharges}];
 xmax = Max[xcoord]*1.5;
 zmax = Max[zcoord]*1.5;
 xmax = Max[{xmax,zmax}];
 zmax = xmax;
 xmin = -xmax;
 zmin = xmin;
 Clear[xcoord,zcoord];
(* --- fix the type of the plot ---*)
 If[typ == "Potential",pot =
Potential[coordinates] /. y -> 0,
 If[typ == "Field",pot =
-Potential[coordinates] /. y -> 0,
 If[typ == "EnergyDensity",pot =
EnergyDensity[coordinates] /. y -> 0,
 Print[" "];
 Print[" wrong key word! Choose "];
 Print[" Potential, Field or EnergyDensity "];
 Print[" to create a plot "];
 Return[]
]]];
(* --- plot the pictures --- *)
 If[typ == "Field",
 pl1 =
StreamPlot[
 Evaluate@{-D[pot, {x}], -D[pot, {z}]}, {x, xmin, xmax}, {z, zmin, 
  zmax}],
 pl1=
ContourPlot[pot,{x,xmin,xmax},{z,zmin,zmax},
 options,
 PlotPoints->50,
 ColorFunction->Hue,
 Contours->15]
 ]
 ];
End[];
EndPackage[];

测试

所用mma版本号12.0