(*Вычисление поля параллелепипеда через сумму полей магнитных \
моментов.*)
\[Mu] = 4 \[Pi] 10^-7;
(*Ток момента*) Ic = 10000000;
(*Ширина основания по горизонтали*) s1 = 0.8;
(*Ширина основания поперечная*) s2 = 0.8;
(*Высота параллелепипеда*) h = 0.8;
(*Радиус магнитного момента*) Rm = (s1 + s2)/40;
(*Диапазон рисунка по вертикали*) a = 0.9;
(*Диапазон изображения по горизонтали*) b = 0.8;
(*Поперечный размер области.*) c = b;
(*Аспект*) v1 = a/b;
(*Аспект*) v = c/b;
(*Относительный шаг таблицы*) \[Lambda] = 0.02;
(*Относительный шаг таблицы*) \[Lambda]1 = 0.05;
(*Количество магнитных моментов*) nm = 6000;
(*Линейная концентрация*) (nm/(s1 s2 h))^(1/3)
(*Концентрация моментов*) nm/(s1 s2 h)
(*Магнитный момент*) Ic \[Pi] Rm^2
F8[x_, y_, z_, \[Beta]_, R_] =
1/(x^2 + y^2 + z^2 + R^2 - 2 R (x Cos[\[Beta]] + y Sin[\[Beta]]))^(
1/2);
(*Векторный потенциал магнитного момента по оси y.*)
Aym[x_, y_, z_,
R_] := -(( \[Mu] Ic R)/(4 \[Pi])) NIntegrate[
F8[x, y, z, \[Beta], R] Cos[\[Beta]], {\[Beta], 0, 2 \[Pi]}];
LAy = Table[{{x, z}, Aym[x, 0, z, Rm]}, {x, -b/2,
b + s1/2, \[Lambda] b}, {z, -a - h/2, a + h/2, \[Lambda] a}];
(*Интерполяция векторного потенциала магнитного момента.*)
Ay = Interpolation[Flatten[LAy, 1]];
Ayo[x_, y_, z_] := Ay[Sqrt[x^2 + y^2], z];
Ayov[x_, y_, z_] := Ayo[x, y, z]/Sqrt[x^2 + y^2] {y, -x, 0};
(*Случайные координаты момента x,y,z. Загружаем область магнита \
моментами.*)
Nz = Table[{s1 RandomReal[] -s1/2, s2 RandomReal[] -s2/2,
h RandomReal[] - h/2}, {n, 1, nm, 1}];
(*Суммарное поле векторного потенциала.*)
Ays[x_, y_, z_] :=
Sum[Ayov[x -Nz[[i]][[1]], y -Nz[[i]][[2]], z -Nz[[i]][[3]]], {i, 1,
nm}];
LAys = Table[{{x, z}, Ays[x, 0, z]}, {x, -b, b, \[Lambda]1 b}, {z, -a,
a, \[Lambda]1 a}];
(*Интерполяция суммарного поля векторного потенциала.*)
Ayss = Interpolation[Flatten[LAys, 1]];
Aysv[x_, z_] := {Ayss[x, z][[1]], Ayss[x, z][[2]]};
(*Вычисление индукции.*)
Bx[x_, y_, z_, R_] := ( \[Mu] Ic R)/(4 \[Pi])
NIntegrate[(
z Cos[\[Beta]])/(x^2 + y^2 + z^2 + R^2 -
2 R (x Cos[\[Beta]] + y Sin[\[Beta]]))^(
3/2), {\[Beta], 0, 2 \[Pi]}];
Bz[x_, y_, z_, R_] := ( \[Mu] Ic R)/(4 \[Pi])
NIntegrate[(
R - y Sin[\[Beta]] -
x Cos[\[Beta]])/(x^2 + y^2 + z^2 + R^2 -
2 R (x Cos[\[Beta]] + y Sin[\[Beta]]))^(
3/2), {\[Beta], 0, 2 \[Pi]}];
LBy = Table[{{x, z}, {Bx[x, 0, z, Rm], Bz[x, 0, z, Rm]}}, {x, -b/2,
b + s1/2, \[Lambda] b}, {z, -a - h/2, a + h/2, \[Lambda] a}];
Byg = Interpolation[Flatten[LBy, 1]];
Byo[x_, y_, z_] := {Byg[Sqrt[x^2 + y^2], z][[1]] x/Sqrt[x^2 + y^2],
Byg[Sqrt[x^2 + y^2], z][[1]] y/Sqrt[x^2 + y^2],
Byg[Sqrt[x^2 + y^2], z][[2]]};
(*Суммарное поле индукции.*)
Bys[x_, y_, z_] :=
Sum[Byo[x -Nz[[i]][[1]], y -Nz[[i]][[2]], z -Nz[[i]][[3]]], {i, 1,
nm}];
LBys = Table[{{x, z}, Bys[x, 0, z]}, {x, -b, b, \[Lambda]1 b}, {z, -a,
a, \[Lambda]1 a}];
Bysv = Interpolation[Flatten[LBys, 1]];
DensityPlot[Ayss[x, z][[2]], {x, -b, b}, {z, -a, a}, ImageSize -> 600,
PlotLegends -> Automatic, ColorFunction -> "SunsetColors",
Epilog -> {{RGBColor[0.5, 0.4, 0], Thick,
Line[{{-s1/2, -h/2}, {-s1/2, h/2}, {s1/2,
h/2}, {s1/2, -h/2}, {-s1/2, -h/2}}]}}]
DensityPlot[Ayss[x, z][[1]], {x, -b, b}, {z, -a, a}, ImageSize -> 600,
PlotLegends -> Automatic, ColorFunction -> "SunsetColors",
PlotRange -> {-17, 7},
Epilog -> {{RGBColor[0.5, 0.4, 0], Thick,
Line[{{-s1/2, -h/2}, {-s1/2, h/2}, {s1/2,
h/2}, {s1/2, -h/2}, {-s1/2, -h/2}}]}}]
DensityPlot[Ayss[x, z][[3]], {x, -b, b}, {z, -a, a}, ImageSize -> 600,
PlotLegends -> Automatic, ColorFunction -> "SunsetColors",
PlotRange -> {-17, 7},
Epilog -> {{RGBColor[0.5, 0.4, 0], Thick,
Line[{{-s1/2, -h/2}, {-s1/2, h/2}, {s1/2,
h/2}, {s1/2, -h/2}, {-s1/2, -h/2}}]}}]
StreamPlot[{Bysv[x, z][[1]], Bysv[x, z][[3]]}, {x, -b, b}, {z, -a, a},
ImageSize -> 600, StreamColorFunction -> "Rainbow",
AspectRatio -> v1, StreamPoints -> Fine,
Epilog -> {{RGBColor[0.5, 0.4, 0], Thick,
Line[{{-s1/2, -h/2}, {-s1/2, h/2}, {s1/2,
h/2}, {s1/2, -h/2}, {-s1/2, -h/2}}]}}]
