来自 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
