Oblate Spheroid, c < a, in Homogeneous Magnetic Field.
(Evaluation of Coefficients in Kuchel and Bulliman, and New Expressions)

Markus Kraiger

mkraiger@student.tugraz.at

Bernhard Schnizer

schnizer@itp.tu-graz.ac.at

2011-12-07, new results added  2012-07-25

Summary:  An oblate spheroid (symmetry axis  ||  z-axis, permeability  μ2)  is in 
a medium of permeabilty  μ1 and in an external homogeneous magnetic field  

				H  =  (H0x, H0y, H0z).
				
The solution of this problem in oblate spheroidal coordinates given in ref.[1] is repeated 
below at first in particular solutions consisting of Legendre Functions and Polynomials.
Thereafter these functions are replaced by the corresponding elementary functions.

[1]  Philip W. Kuchel, Brian T. Birman: Perturbation of Homogeneous Magnetic Fields by 
       Isolated Single and  Confocal Spheroids. Implications for NMR Spectroscopy of Cells. 
       NMR in Biomedicine, 2(4) (1989) 151-160.
       
[2]  M. Kraiger, B. Schnizer:  Reaction Fields of Homogeneous Magnetic Spheroids of 
       Arbitrary Direction in a Homogeneous Magnetic Field. A Toolbox for MRI and MRS 
       of Heterogeneous Tissue.
       Report ITPR-2011-21, Graz University of Technology, Institute of Theoretical and 
       Computational Physics, November 2011.

[3] Perry Moon and D.E.Spencer : Field Theory Handbook. Including coordinate systems
      differential equations and their solutions.Springer 1988.

$Version

"7.0 for Mac OS X x86 (64-bit) (February 19, 2009)"

Descriptions

Ref.[2], Paragraph 4.1;  Ref.[3], Table 1.07 

Oblate spheroidal coordinates

	x =  eo coshη sinθ cosψ
	y =  eo coshη sinθ sinψ
	z =  eo sinhη cosθ
	
Interface = oblate spheroid with semi-axes  ao, ao, co;  co < ao

	((x^2)/(Subscript[a, o]^2)) +  ((y^2)/(Subscript[a, o]^2)) + ((z^2)/(Subscript[c, o]^2)) =  1	
	eo = (Sqrt[(ao^2) -  (co^2)]) 
	
	η = ηo = Artanh(co/ao)
	
Exterior permeability:		μ1 = μ0 (1 + χ1)
Interior permeability:		μ2 = μ0 (1 + χ2)	
	
Functions occuring in the particular solutions

	P[n, m, θ ] 	=  (Subscript[P, n]^m)(cosθ)   	associated Legendre polynomial
	LeP[n,m, θ]	=  (Subscript[P, n]^m)(i sinhη) 	associated Legendre function of the first kind
	LeQ[n,m, θ]	=  (Subscript[Q, n]^m)(i sinhη) 	associated Legendre function of the second kind
	
Expansion coefficients

	A[d,n,m]	=  (Subscript[, d])(Subscript[A, n]^m)  ,  d = 1 (Exterior),  d = 2 (Interior of spheroid)
	B[d,n,m]	=  (Subscript[, d])(Subscript[B, n]^m)

Transformations and Substitutions

Transformation fromCartesian to Oblate Elliptic Coordinates

ToEll = {x -> eo Cosh[η] Sin[θ] Cos[ψ],
             y -> eo Cosh[η] Sin[θ] Sin[ψ],
	    z -> eo Sinh[η] Cos[θ] };

Transformation from Oblate Ellipsoida to Cartesian Coordinates

Toxyz =  Map[Reverse,ToEll]

η for the interface

seo = {η ->  ηo};

Transformation to Particular Solutions

ToPart[η_] = { Cos[θ] -> P[1,0,θ], Sin[θ] -> P[1,1,θ],
                          Cosh[η] -> - I LeP[1,1,η], Sinh[η] -> -I LeP[1,0,η]}

Transformation from Particular Solutions to Elementary Functions.
Legendre Polynomials and Functions Expressed by Elementary Functions

Transformation of the Functions

ToElem[η_] ={ P[1,0,θ] ->  Cos[θ] , P[1,1,θ]  -> Sin[θ],
LeP[1,0,η] -> I Sinh[η] ,LeP[1,1,η] -> I Cosh[η],
LeQ[1,0,η]->   Sinh[η] ArcCot[Sinh[η]] -  1,
LeQ[1,1,η]->   Cosh[η] ArcCot[Sinh[η]] - Tanh[η]}

Transformation of the Derivatives w.r.t. η

ToElemD[η_] = Drop[D[%,η]//FullSimplify,2]

Transformation of Functions and Derivatives

ToElem01[η_] = Join[ToElem[η],ToElemD[η]] ;

Substitutiing magnetic susceptibilities in place of permeabilities

suchi = {μ1 ->  μ0 (1 + χ1), μ2 ->  μ0 (1 + χ2)}

Solution of the Potential Equation in Spheroidal Coordinates

Primary Potential

Primary Potential in Cartesian Coordinates

Φ0C = - (H0x x + H0y y + H0z z )

Primary Potential in Elliptic Coordinates

Φ0E[η_] =   Φ0C/. ToEll

Primary Potential in Particular Functions

Φ0P[η_] =   Φ0E[η] /. ToPart[η]

Total External Potentiall, [2], eq.(4.6)

Φ1[η_]  = Φ0P[η]  +
Sum[Sum[A[1,n,m]LeQ[n,m,η]  P[n,m,θ] Cos[m ψ ],{n,1,2}],{m ,0,2}] +
Sum[Sum[B[1,n,m]LeQ[n,m,η]  P[n,m,θ] Sin[m ψ ],{n,1,2}],{m ,1,2}]

Internal Potential, eq.(4.7)

Φ2[η_]  = I Sum[Sum[ A[2,n,m]LeP[n,m,η]  P[n,m,θ] Cos[m ψ ],{n,1,1}],{m ,0,1}] +
I Sum[Sum[B[2,n,m]LeP[n,m,η]  P[n,m,θ] Sin[m ψ ],{n,1,1}],{m ,1,1}]

Solution of the Particular Systems

n = 1, m = 0

ls = Coefficient[Φ1[η] ,P[1,0,θ]]

rs = Coefficient[Φ2[η] ,P[1,0,θ]]

dls = μ1 D[ls,η]

drs = μ2 D[rs,η]

eqs = {ls == rs, dls ==drs}

so10 = Solve[eqs /. seo,{A[1,1,0],A[2,1,0]}]//Factor//Flatten

L0,   [2], eq.(4.9)

L0 = A[1,1,0]/(eo H0z)/. so10

M0,  eq.(4.10)

M0 = A[2,1,0]/(eo H0z)/. so10

n = 1, m = 1, Cos[ψ]

ls = Coefficient[Φ1[η] ,P[1,1,θ]Cos[ψ]]

rs = Coefficient[Φ2[η] ,P[1,1,θ]Cos[ψ]]

dls = μ1 D[ls,η]

drs = μ2 D[rs,η]

eqs = {ls == rs, dls ==drs}

so11c = Solve[eqs /. seo,{A[1,1,1],A[2,1,1]}]//Factor//Flatten

A[1,1,1]/.so11c

L1,  [2] eq.(4.11)

L1 = A[1,1,1]/(eo H0x)/. so11c

M1,  [2], eq.(4.13)

M1 = A[2,1,1]/(eo H0x)/. so11c

n = 1, m = 1. Sin[ψ]

ls = Coefficient[Φ1[η] ,P[1,1,θ]Sin[ψ]]

rs = Coefficient[Φ2[η] ,P[1,1,θ]Sin[ψ]]

dls = μ1 D[ls,η]

drs = μ2 D[rs,η]

eqs = {ls == rs, dls ==drs}

so11c = Solve[eqs /. seo,{B[1,1,1],B[2,1,1]}]//Factor//Flatten

L1b = B[1,1,1]/(eo H0y)/. so11c

M1b = B[2,1,1]/(eo H0y)/. so11c

L1b == L1

M1b== M1

The coefficients   L1  and   M1  are the same  for  H0x  and  H0y  !

Replacing The coefficients A  and  B    by   L0, M0, L1, M1

Substitution rules resulting from the continuity conditions

sucoff1 = {A[1,1,0]-> L0c(eo H0z), A[2,1,0]-> M0c(eo H0z),
	           A[1,1,1]-> L1c(eo H0x), A[2,1,1]-> M1c(eo H0x),
                    B[1,1,1]-> L1c(eo H0y),  B[2,1,1]-> M1c(eo H0y),
                    Table[{A[dom,1,2] ->  0 ,B[dom,1,2] ->  0 },{dom,1,2}]}//Flatten

sucoff2 = Table[{A[dom, 2, mm] -> 0,B[dom, 2, mm] -> 0},{dom,1,2},{mm,0,2}]//Flatten

sucoff = Join[sucoff1, sucoff2]

Coefficients   for   μ2 = μ1

In that case the reaction is zero so that:

NoReact1 = {L0c -> 0 , L1c -> 0}

The Total Outer Potential

Φ1P[η]= Φ1[η]/. sucoff

The Total Outer Potential  for   μ2 = μ1

Φ1P[η]/.NoReact1

The Total Inner Potential

Φ2P[η] = Φ2[η]/. sucoff

Expressing Potentials in Elementary Functions

Exterior Potentials

Total Exterior Potential

Φ1Ef[η_] = Φ1P[η] /. ToElem01[η]

Part proportional to z

hh = Map[Factor,Coefficient[Φ1Ef[η] , Cos[θ]]] Cos[θ] //Factor

Primary potential

Pz = hh  /. NoReact1

Reaction potential

Rz = hh - Pz // Simplify

Part proportional to x

hh = Map[Factor,Coefficient[Φ1Ef[η] , Sin[θ] Cos[ψ]]] Sin[θ] Cos[ψ]//Factor

Primary potential

Px = hh  /.NoReact1

Reaction potential

Rx = hh - Px // Simplify

Part proportional to y

hh = Map[Factor,Coefficient[Φ1Ef[η] , Sin[θ] Sin[ψ]]] Sin[θ] Sin[ψ]//Factor

Primary potential

Py = hh  /. NoReact1

Reaction potential

Ry = hh - Py // Simplify

Total Primary Potential

Φ1p = Px + Py + Pz

Φ1pxyz =  %/. Toxyz

Total Exterior Reaction Potential

Φ1Er[η_] = Rx + Ry + Rz

The above expression is transformed  by hand such that it is linear in x, y, z .

Final Form of the Exterior Reaction Potential in Elliptical Coordinates

Φ1Erf[η_]= eo H0z L0c Cos[θ] Sinh[η](-1/ Sinh[η] + ArcCot[Sinh[η]] )+eo H0x L1c Cos[ψ] Sin[θ] Cosh[η] (-Tanh[η]/Cosh[η] +ArcCot[Sinh[η]] )+eo H0y L1c Sin[θ] Sin[ψ]Cosh[η] (- Tanh[η]/Cosh[η]+ArcCot[Sinh[η]] )

The equivalence of  Φ1Erf[η_]  with  Φ1Er[η_]  is checked:

Φ1Er[η] -Φ1Erf[η]//Simplify

Final Form of the Exterior Reaction Potential in CartesianCoordinates Apart of f1, f2, f3

Φ1xyzr[η_] = Φ1Erf[η]/. Toxyz

Defining the functions g1, g2, g3:

sugi = {ArcCot[Sinh[η]]-> g1 ,
Tanh[η]/Cosh[η] -> g2,  1/Sinh[η] -> g3}

Φ1xyzrf[η_] =Φ1xyzr[η]/. sugi

Total Exterior Potential

Map[Factor, Φ1pxyz + Φ1xyzrf[η],{4}]

Φ1xyzf[η_]=
(-1 + (g1-g2)L1c) H0x  x+(-1 + (g1-g2)L1c) H0y  y+(-1 + (g1-g3) L0c) H0z z

% - %% //Simplify

Total Interior Potential

Φ2Ef[η_] = Φ2P[η] /. ToElem01[η]//Expand

Φ2xyzf[η_] =  Φ2Ef[η] /. Toxyz

The coefficients L0, L1, M0, M1 expressed in Elementary Functions

L0

L0

hL0 = L0 /.ToElem01[ηo]/. suchi//Cancel

nu = Numerator[hL0]/Cosh[ηo]

de =  Denominator[hL0]/Cosh[ηo]//Simplify

L0f = nu/de

hL0 - L0f //FullSimplify

L1

L1

hL1 = L1 /.ToElem01[ηo]/. suchi//FullSimplify

nu = Numerator[hL1]Cosh[ηo]

de = Denominator[hL1]Cosh[ηo]//Simplify

den =  -(2 + χ1+χ2  +  (χ1-χ2) ((Cosh[ηo])^2)(1 - ArcCot[Sinh[ηo]] Sinh[ηo]))

de - den //Simplify

L1f = nu/den

hL1 - L1f //Simplify

M0

M0

hM0 = M0 /.ToElem01[ηo]/. suchi //FullSimplify

nu = Numerator[hM0]/Cosh[ηo] //Simplify

de = Denominator[hM0]/Cosh[ηo] //FullSimplify

den = 1+χ2+(χ1-χ2) ArcCot[Sinh[ηo]] Sinh[ηo]-(1+χ1) ((Tanh[ηo])^2)

M0f = nu/den

hM0 - M0f //Simplify

M1

M1

hM1 = M1 /.ToElem01[ηo]/. suchi//Cancel //Simplify

nu = Numerator[hM1]Cosh[ηo]//Expand//FullSimplify

de = Denominator[hM1]Cosh[ηo] //FullSimplify

M1f = nu/de

nu = Numerator[%]Cosh[ηo]

de =  Denominator[%%]Cosh[ηo]

den = ((2+χ1+χ2) -(χ1-χ2) Cosh[ηo]^2 (-1+ArcCot[Sinh[ηo]] Sinh[ηo]))

de - den //Simplify

M1f = nu/den

hM1 - M1f//FullSimplify

Final formulas for the coefficients  L0, L1, M0, M1

η = ηo = Artanh(co/ao)

L0f

L1f

M0f

M1f

Checks for Φ1 and Φ2

Interior coefficients for interior medium = exterior medium

M0f /. χ2  -> χ1//Simplify

M1f /. χ2  -> χ1//Simplify

NoReact2 = {M0c -> 1, M1c -> 1}

Exterior coefficients for interior medium = exterior medium

L0f /. χ2  -> χ1//Simplify

L1f /. χ2  -> χ1//Simplify

NoReact1

NoReact =Join[NoReact1,NoReact2];

Total potential in the  interior for interior medium = exterior medium

Φ2xyz0 = Φ2xyzf[η] /. NoReact

prinary exterior potential

Φ0C

Φ2xyz0 ==Φ0C

Final form of  total interior  potential

Φ2xyzf[η]

Final form of  interior reaction potential

Φ2xyzr = Φ2xyzf[η] -Φ2xyz0//Simplify

Final form of  total exterior  potential

Φ1xyzf[η]

Final form of  exterior reaction potential

Φ1xyzrf[η]

Map[Reverse,sugi]

Results

Primary Potential

Φ1pxyz

Exterior Reaction Potential, ref.[2], eq.(4.14, 20)

Φ1xyzrf[η]

L0c -> L0f
L1c -> L1f

Map[Reverse,sugi]          (* ref.[2],eq.(4.22 - 24)*)

Interior Reaction Potential, ref.[2], eq.(4.21)

Φ2xyzr

M0c -> M0f
M1c -> M1f

Testing Continuity Conditions

Continuity of Potential at interface  η = ηo

sufico ={L0c -> L0f,L1c -> L1f,M0c -> M0f,M1c -> M1f};

te1 = Φ1xyzf[ηo]/. Map[Reverse,sugi]/.seo/.sufico ;

te2 = Φ2xyzf[ηo] /. sufico

te1 - te2;
Simplify[%]

Continuity of Normal Component of Induction, B1η = B2η, 
at  interface   η = ηo

Φ1xyzf[η]/eo /.ToEll/. Map[Reverse,sugi]//Expand

D[%,η]//Simplify ;

% /. seo ;

% /. sufico  ;

te3 = (1+χ1) % //Expand ;

D[Φ2xyzf[η]/eo/.ToEll,η]/. seo/.sufico//Expand ;

te4 = (1 +χ2)% //Expand ;

te3 - te4 ; 
Simplify[%]

Revised Results

2012-07-25

The expressions depending on the parameter  ηo   are replaced by equivalent expressions depending
on (Subscript[a, o]) , (Subscript[c, o]) , (Subscript[e, o])  , the semi-axes and the excentricity.

L0f

L1f

M0f

M1f

Substituting semi-axes

suabe ={Sinh[ηo]->(co/eo),Cosh[ηo]->(ao/eo),Sech[ηo]->(eo/ao), Tanh[ηo]-> co/ao,  ArcCot[(co/eo)]-> do}

L0e = L0f //. suabe

L1e = L1f //. suabe

M0e = M0f //. suabe

M1e = M1f //. suabe