//Parameters real Mu0 = 4.*pi*1.e-7; //Vacuum magnetic permeability real MuC = 1.25e-6; //Copper magnetic permeability real J0 = 500.e4; //Current density //Mesh int nn = 25; //Mesh quality real AirRadius = 2.; //"Infinite" boundary real WireRadius = 8e-3; real CoilYSpacing = 0.018; real CoilRadius = 0.015; // int nWinding = 7; int WindCtrl = nWinding-3; real[int] yWinding(nWinding); for(int i=0; i Laplacian = vLaplacian(Ah, Ah, solver=sparsesolver); real[int] LaplacianBoundary = vLaplacian(0, Ah); Az[] = Laplacian^-1 * LaplacianBoundary; //Magnetic induction Ah Bx, By; Bx = dy(Az); By = -dx(Az); Ah B = sqrt(Bx^2 + By^2); //Magnetic field Ah Hx, Hy; Hx = Nu * Bx; Hy = Nu * By; Ah H = sqrt(Hx^2 + Hy^2); //Current Ah J = dx(Hy) - dy(Hx); //Lorentz force Ah Lx, Ly; Lx = -J * By; Ly = J * Bx; Ah L = sqrt(Lx^2 + Ly^2); //Plot plot(Az, nbiso=30, fill=true, value=true, cmm="A", wait=1); plot(B, nbiso=30, fill=true, value=true, cmm="B - Magnetic induction", wait=1); plot(H, nbiso=30, fill=true, value=true, cmm="H - Magnetic field", wait=1); plot(J, nbiso=30, fill=true, value=true, cmm="J - Current", wait=1); plot(L, nbiso=30, fill=true, value=true, cmm="Lorentz force", wait=1); Ah MuVal = Mu; plot(MuVal, nbiso=30, fill=true, value=true, cmm="Permeability");