Clear["Global`*"]
(* © Marko Horbatsch, York University, Toronto, Canada *)

(* We solve for the vector potential of a spinning spherical shell with surface charge σ. The colume current density is j_phi = σ δ(r-R) 
It can be used to verify that in the x-z plane the y-component of A satisfies a Poisson equation (5.64). A_phi is the only component given that j=j_phi. To use (5.64) one has to convert it into A_y and j_y by multiplying on Cos[φ], and then use the Laplacian in SPC. Ultimately one gets an equation for A_phi(r,θ), which can be solved by separation of variables. The simplest solution is (5.69), where j_phi and A_phi are proportional to Sin[θ], as a result of geometry. *)
(* we have generalized this result by finding another form that allows to modulate the surface charge density into σ(θ) by multiplying with Cos[θ], and the appropriate radial function associated with the δ(r-R) profile. 
Two other cases are separable, but they involve a singularity at θ=0,π: 1/Sin[θ], and Cos[θ]/Sin[θ]; they go with a radial function that is a constant=R^2 (for r<R) matched to R/r for r>R *)

c0=1;

(* one way to interpret the above plot: how many Amps/m surface current density do we have in the shell as a function of polar angle θ. To express it as a volume current density we multiply with a Dirac δ(r-R) which has the needed extra dimension of /per unit length, i.e., 1/m in SI. We can then think of the spherical shell as a conductor with different resisitivity, such that it passes (much) more current at the equator, when c2 and c4 are turned on to diminish the current density closer to the poles of the sphere. This is how we can connect the spinning sphere with K(θ)=σ(θ)*ω*R*Sin[θ] to a model for an electromagnet. *)

(* three cases are pre-set:

case 1:  j_ϕ~Sin[θ] Griffiths, Example 5.11, dipole
case 2:  j_ϕ~Cos[θ]Sin[θ] (no dipole, but quadrupole)
case 3:  j_ϕ~|Cos[θ]|Sin[θ] (quadrupole)
*)
iCase=2;

Clear[A2];
A2[r_,θ_]:=Which[iCase==1,Sin[θ]Piecewise[{{c0/3 r,r<R},{c0/3 R^3/r^2,r>R}}]/.R->1,iCase==2,Sin[θ]Cos[θ]Piecewise[{{c0/5 r^2,r<R},{c0/5 R^5/r^3,r>R}}]/.R->1,iCase==3,Piecewise[{{Sin[θ]Cos[θ],-Pi<θ<-Pi/2},{Sin[θ]Cos[θ],-Pi/2<θ<Pi/2},{-Sin[θ]Cos[θ],Pi/2<θ<3/2Pi},{0,θ==0}}]Piecewise[{{c0/5 r^2,r<R},{c0/5 R^5/r^3,r>R}}]/.R->1]

(* Clear[A2];
A2[r_,θ_]:=Cos[θ]/Sin[θ]Piecewise[{{c0 R^2,r<R},{c0 R^3/r,r>R}}]/.R->1 *)
(* Clear[A2];
A2[r_,θ_]:=1/Sin[θ]Piecewise[{{c0 R^2,r<R},{c0 R^3/r,r>R}}]/.R->1 *)

Clear[Br,Bθ];
Br[r_,θ_]:=D[Sin[θ1]A2[r,θ1],θ1]/(r Sin[θ1])/.θ1->θ

Bθ[r_,θ_]:=-D[r1 A2[r1,θ],r1]/r1/.r1->r

Plot[{Br[r,Pi/4],Br[r,Pi/3],Br[r,7Pi/16]},{r,0,3},PlotStyle->{Red,Blue,Green},PlotRange->All]

(* the radial component of the magnetic field is continuous at the charge surface (R=1), since it is the normal component. The Bθ-field component - this is B-parallel which is a 1d vector = scalar should have a jump discontinuity at R=1 due to the surface current *)

Plot[{Bθ[r,Pi/4],Bθ[r,Pi/3],Bθ[r,7Pi/16]},{r,0,3},PlotStyle->{Red,Blue,Green},PlotRange->All]

(* we picked three angles in the first quadrant of the x-z plane; do we understand the magnitude of the jumps? *)
(* We should form the cross product of K_phi with the normal vector (6.27), to see how (5.75) comes about. ϕ-hat cross r-hat gives θ-hat, that is why Bθ has the jump discontinuity, and it is the sole component of the B-parallel vector, since B_ϕ vanishes by symmetry here. *)

η=10^-8;
Discontinuity=N[{Bθ[r+η,Pi/4]-Bθ[r-η,Pi/4],Bθ[r+η,Pi/3]-Bθ[r-η,Pi/3],Bθ[r+η,7Pi/16]-Bθ[r-η,7Pi/16]}/.r->1]

(* if we are not including mu_0, then we are working with the auxiliary field H (given in Amps/m), which (IMHO) we should always, and only convert to B, when using the magnetic field for force calculations, or for measurement comparisons in Tesla. The we will never be confused about which field is which when using the presence of medium other than vacuum or air. A free current produces H, the real field follows when magnetization is involved by tacking on the right factor, or in the case of ferromagnetic material by usually considering the magnetization to be saturated. IMHO, the focus on B-fields as following from currents directly (with the mu_0 there right away - assuming air=vacuum), is more confusing than helpful. *)

Clear[Kϕ];
Kϕ[θ_]:=Which[iCase==1,c0 Sin[θ],iCase==2,c0 Sin[θ] Cos[θ],iCase==3,c0 Sin[θ]Abs[Cos[θ]]];

DiscFun[θ_]:=(Bθ[r+η,θ]-Bθ[r-η,θ])/Kϕ[θ]/.r->1

N[{DiscFun[Pi/4],DiscFun[Pi/3],DiscFun[7Pi/16]}]

(* - we used c0 to incorporate all other physical constants such as R, σ, ω, and μ0. We do need the latter since B=curl(A),and not H. We set c0=1. 
The H-discontinutity would have been consistent with surface charge density given by K_ϕ=σωRSin[θ] for Ex.5.11 *)

(* now that we verified the physics principles at work here, we are allowed to go on with visualization tools provided my Mma *)

(* we need two steps: step 1 = conversion of {Br,Bθ} into {Bx,By}, which will still be in terms of SPC;
step 2 = define {r,θ} in terms of {x,z} to get a cut in a Cartesian coordinate plane, which our brain is used to for visualization. *)

Bx[r_,θ_]:=Br[r,θ]Sin[θ]+Bθ[r,θ]Cos[θ]
Bz[r_,θ_]:=Br[r,θ]Cos[θ]-Bθ[r,θ]Sin[θ]

(* the back page in Griffiths gives the unit vector conversion = step 1; in step 2 we define new functions - this works better than doing conversions 'on the fly'. *)

Clear[BX,BZ];
BX[x_,z_]:=Bx[r,θ]/.{r->Sqrt[x^2+z^2],θ->ArcTan[z,x]}
BZ[x_,z_]:=Bz[r,θ]/.{r->Sqrt[x^2+z^2],θ->ArcTan[z,x]}

(* DensityPlot[{BX[x,z],BZ[x,z]},{x,0.01,2},{z,0.01,2},PlotRange->All,MaxRecursion->5,FrameLabel->{x,z}] *)
(* looks nice, but is useless! *)

(* the density plot is supposed to reveal where the magnetic field is strong, but BEWARE! it does show uniformity inside the spinning shell - it is supposed to be a field strength plot *)
(* You may get fooled by the glow: Apparently, the two magnetic field components conspire to give biggest magnetic field strength just outside the spinning shell, not at the equator, and not at the poles!??? *)
(* can this be right? The 2d plots of the magnetic field components are strongest inside the shell. So what is going on? *)
(* A_ϕ is a 1d vector = scalar. It can be plotted as a contour plot with values attached to the contours. We convert it to Cartesians: *)

Ac[x_,y_]:=A2[r,θ]/.{r->Sqrt[x^2+z^2],θ->ArcTan[z,x]}

ContourPlot[Ac[x,z],{x,0,3},{z,-1.5,1.5},Contours->15,ContourLabels->All,FrameLabel->{x,z}]

(* inside the shell, A_ϕ grows with x (uniform {H,B}-field);
in the region where the DensityPlot for {B_x,B_z} showed a glow we have strong contour discontinuity along the shell surface. *)

Clear[B];
B[x_,z_]:=Sqrt[BX[x,z]^2+BZ[x,z]^2]

(*ContourPlot[B[x,z],{x,η,3},{z,-1.5,1.5},Contours->15,ContourLabels->All,PlotRange->All,FrameLabel->{x,z},MaxRecursion->5] (* the last switch makes this command  run forever, but the plot comes out nice! *)

(* For c2=0: since B is constant inside the shell, we don't get a contour value for the field there. *)

{B[0.25,1.35],B[-0.25,-1.35]}

(* the field is the strongest inside the sphere!
Message here: beware of density plot! *)

Plot3D[B[x,z],{x,-1.5,1.5},{z,-1.5,1.5},PlotRange->{{-1.5,1.5},{-1.5,1.5},Automatic},MaxRecursion->6,AxesLabel->Automatic]

(* there is a very steep drop in field strength at the equator, more so than in the direction towards the poles. *)
(* this observation is perhaps related to the question of the B field outside a solenoid in the radial direction: it drops off so fast that the Ampere-law result for an infinitely long solenoid says that the field is zero there - at the equator our spinning charged sphere shows a dramatic drop in strength just outside the surface *)

(* Can we get a field line plot? *)

(* StreamPlot[{BX[x,z],BZ[x,z]},{x,0.01,2},{z,0.01,2},PlotRange->All] *) (* beats me why this fails! *)

PL1=ListStreamPlot[Table[{{x,z},{BX[x,z],BZ[x,z]}},{x,-2.0001,2,0.02},{z,-2.0001,2,0.02}],StreamPoints->150,FrameLabel->{x,z}];
PL2=ParametricPlot[{Cos[θ],Sin[θ]},{θ,0,2Pi},PlotStyle->Red];
Show[PL1,PL2]

(* this gives a general idea of the field direction, but it is not yet what we consider a field line plot in physics;
also note: inaccuracies at the shell surface were prevented by using a large table - the spacing dx, dz were selected as 0.02, i.e., quite small.

given the uniformity of the field inside the shell, we should choose equidistant seeds for stream lines there. *)

K0=FindMaximum[Kϕ[θ],{θ,0.9Pi/2}][[1]]
(* where the current density is below 50% of max we will show the surface in black vs red *)

PL1=ListStreamPlot[Table[{{x,z},{BX[x,z],BZ[x,z]}},{x,0.0001,1.5,0.02},{z,0.0001,1.5,0.02}](*,StreamPoints->Table[{{(j-1/2)/10,0.05},Blue},{j,1,10}]*),FrameLabel->{x,z},PerformanceGoal->"Quality",StreamStyle->LightGray,VectorPoints->9,VectorColorFunction->Hue,VectorScale->{Medium,1}];
PL2=ParametricPlot[{Sin[θ],Cos[θ]},{θ,0,Pi/2},ColorFunction->Function[{x,y,θ},If[Kϕ[θ]/K0<0.5,Black,Red]],PlotStyle->Thick,ColorFunctionScaling->False];
Show[PL1,PL2]

(* what are field lines? *)

ODE={x'[t]==BX[x[t],z[t]],z'[t]==BZ[x[t],z[t]]};

Nr=15;
tmax=560;

sol={};
Do[IC={x[0]==(j-0.5)/Nr*Sqrt[1-0.5^2],z[0]==0.5};
sol=AppendTo[sol,NDSolve[Flatten[Join[ODE,IC]],{x,z},{t,0,tmax}][[1]]];
,{j,1,Nr}];//Quiet

{x[20],z[20]}/.sol[[6]]

PL3=ParametricPlot[Table[{x[t],z[t]}/.sol[[j]],{j,1,Nr}],{t,0,tmax},PlotRange->{{0,1.6},{0,1.6}},PlotStyle->Blue];
Show[PL3,PL2,AxesLabel->{x,z}]

(* note: the 2d fieldlines construction works only properly if the field is uniform over the region where the lines are started. *)

(* Now we have a representation that let's us feel the field strength from the separation of the parametric curves. We based the field line density on the x-z plane, and used the fact that equispaced seed values in the constant-field region are appropriate. *)
(* Note how strong the discontinuity is in the tangential field component in the regions of appreciable current density. *)

(* verify the calculated magnetic field is curl-free and divergence-free: *)

curlBϕ[r_,θ_]:=(D[r1 Bθ[r1,θ],r1]-D[Br[r,θ1],θ1])/r /.{θ1->θ,r1->r}

curlBϕ[0.5,0.5]

divB[r_,θ_]:=D[r1^2 Br[r1,θ],r1]/r1^2+D[Sin[θ1]Bθ[r,θ1],θ1]/r /Sin[θ1]/.{θ1->θ,r1->r}

divB[0.5,0.5]

(* ASSIGNMENTS: 
1) Compare the three possible scenarios:
a) Griffiths textbook example 5.11, j_ϕ ~ Sin[θ] -> uniformly charged shell σ=const.
b) modulated surface charge: j_ϕ ~ Sin[θ]Cos[θ]
c) modulated surface charge: j_ϕ ~ Sin[θ]Abs[Cos[θ]]

the three cases can be studied in the large-r limit using the multipole expansion. Is it true that case (b) leads to a vanishing dipole moment?
 *)

(* Follow-up: how did we find the examples (b,c) - i.e., when does separation of variables work? *)

Clear[LHS];
LHS[g_]:=(D[f[r],{r,2}]Sin[θ]g[θ]+f[r]/r^2(D[Sin[θ]D[Sin[θ]g[θ],θ],θ]/Sin[θ]-g[θ]/Sin[θ]))/(Sin[θ]g[θ])

LHS[Cos]//FullSimplify

(* there are other cases, which are somewhat pathological in that the surface charge blows up at the poles: (j=-1 does not separate, but j=0,1 works)*)

Myf[x_]:=Table[Cos[x]^j Sin[x]^-2,{j,-1,1}]

LHS[Myf]//FullSimplify

Plot[1/Sin[x],{x,0,Pi}]