Prolate 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-11-08, new results added  2012-07-25,
corrected 2012-10-10

Summary:  A prolate 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 prolate 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 3.1;  Ref.[3], Table 1.06 

Prolate spheroidal coordinates

	x =  ep sinhη sinθ cosψ
	y =  ep sinhη sinθ sinψ
	z =  ep coshη cosθ

Interface = prolate spheroid with semi-axes  ap, ap, cp;  cp > ap

	((x^2)/(Subscript[a, p]^2)) +  ((y^2)/(Subscript[a, p]^2)) + ((z^2)/(Subscript[c, p]^2)) =  1	
	ep = (Sqrt[(cp^2) -  (ap^2)]) 
	
	η = ηp = Arcoth(cp/ap)                            corrected 2012-10-10
	
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)(coshη) 	associated Legendre function of the first kind
	LeQ[n,m, θ]	=  (Subscript[Q, n]^m)(coshη) 	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 Prolate Elliptic Coordinates

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

Transformation from Prolate Ellipsoida to Cartesian Coordinates

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

η for the interface

sep = {η ->  ηp};

Transformation to Particular Solutions

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

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,1,η] -> Sinh[η] ,LeP[1,0,η] ->  Cosh[η],
LeQ[1,0,η]->   Cosh[η] ArcCoth[Cosh[η]] -  1,LeQ[1,1,η]->  Sinh[η] ArcCoth[Cosh[η]] - Coth[η] };

Transformation of the Derivatives w.r.t. η

ToElemD[η_] = Drop[D[{ P[1,0,θ] ->  Cos[θ] , P[1,1,θ]  -> Sin[θ],
LeP[1,1,η] -> Sinh[η] ,LeP[1,0,η] ->  Cosh[η],
LeQ[1,0,η]->   Cosh[η] ArcCoth[Cosh[η]] -  1,LeQ[1,1,η]->  Sinh[η] ArcCoth[Cosh[η]] - Coth[η] },η]//Simplify,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 Potential, [2], eq.(3.10)

Φ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.(3.11)

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

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 /. sep,{A[1,1,0],A[2,1,0]}]//Factor//Flatten

L0,   [2], eq.(3.12), [1], eq.(A21)

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

M0,  eq.(3.13), [1], eq.(A22)

M0 = A[2,1,0]/(ep 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 /. sep,{A[1,1,1],A[2,1,1]}]//Factor

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

L1,  [2] eq.(3.14); [1], eq.(A23)

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

M1,  [2], eq.(3.15); [1], eq.(A24)

M1 = A[2,1,1]/(ep 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 /. sep,{B[1,1,1],B[2,1,1]}]//Factor//Flatten

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

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

L1b == L1

M1b== M1

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

n = 2, m = 0

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

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

dls = μ1 D[ls,η]

drs = μ2 D[rs,η]

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

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

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

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

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

dls = μ1 D[ls,η]

drs = μ2 D[rs,η]

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

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

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

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

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

dls = μ1 D[ls,η]

drs = μ2 D[rs,η]

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

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

All coefficients in  Φ1  and  Φ2  for  n ≠  1  are zero  !

Replacing the Coefficients A[]  and  B[]  by   L0, M0, L1, M1

Substitution Rules Resulting from the Continuity Conditions

sucoff1 = {A[1,1,0]-> L0c(ep H0z), A[2,1,0]-> M0c(ep H0z),
	           A[1,1,1]-> L1c(ep H0x), A[2,1,1]-> M1c(ep H0x),
                    B[1,1,1]-> L1c(ep H0y),  B[2,1,1]-> M1c(ep 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[η_]= ep H0z L0c Cos[θ] Cosh[η](-1/Cosh[η] +ArcCoth[Cosh[η]] )+ep H0x L1c Cos[ψ] Sin[θ] Sinh[η] (-Coth[η]/Sinh[η]+ArcCoth[Cosh[η]] )+ep H0y L1c Sin[θ] Sin[ψ]Sinh[η] (-Coth[η]/Sinh[η]+ArcCoth[Cosh[η]] )

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 f1, f2, f3:

sufi = {ArcCoth[Cosh[η]] -> f1 ,
Coth[η]/Sinh[η] -> f2,  1/Cosh[η] -> f3}

Exterior Reaction Potential

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

Map[Reverse, sufi]

Total Exterior Potential

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

Φ1xyzf[η_]=
(-1 + (f1-f2)L1c) H0x  x+(-1 + (f1-f2)L1c) H0y  y+(-1 + (f1-f3) L0c) H0z z

% - %% //Simplify

Total Interior Potential

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

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

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

L0

L0

hL0 = L0 /.ToElem01[ηp]/. suchi//Simplify

nu = Numerator[hL0]/Sinh[ηp]

de = Denominator[hL0]/Sinh[ηp]//Expand//Simplify

L0f = nu/de

hL0 - L0f //FullSimplify

L1

L1

hL1 = L1 /.ToElem01[ηp]/. suchi//Simplify

nu = Numerator[hL1]/Sinh[2 ηp]

de = Denominator[hL1]/Sinh[2 ηp]//Expand//TrigExpand//Simplify

The above expression, de,  has been simplified by hand giving denew.

denew = (χ1-χ2) (ArcCoth[Cosh[ηp]] - Sech[ηp])  + (2+ χ1+χ2  ) ((Csch[ηp])^2) Sech[ηp]

de - denew //FullSimplify

L1f = nu/denew

hL1 - L1f //FullSimplify

M0

M0

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

nu =  Numerator[hM0](1+χ1)

de = - Denominator[hM0](1+χ1)//Expand//FullSimplify

The above expression, de,  has been simplified by hand giving denew.

denew = 1+χ1+(χ1-χ2) (1- ArcCoth[Cosh[ηp]] Cosh[ηp]) ((Sinh[ηp])^2)

de - denew //FullSimplify

M0f = - nu/denew

hM0 - M0f//FullSimplify

M1

M1

hM1 = M1 /.ToElem01[ηp]/. suchi //FullSimplify

nu =  Numerator[hM1]/Csch[ηp] /2

de = - Denominator[hM1]/Csch[ηp] /2//Expand//FullSimplify

The above expression, de,  has been simplified by hand giving denew.

denew = 2+χ1+χ2-(χ1-χ2) (1-ArcCoth[Cosh[ηp]] Cosh[ηp]) ((Sinh[ηp])^2)

de - denew //FullSimplify

M1f = - nu/denew

hM1 - M1f//FullSimplify

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

η = ηp = Arcoth(ap/cp)

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[η]

Results

Primary Potential

Φ1pxyz

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

Φ1xyzrf[η]

L0c -> L0f
L1c -> L1f

Map[Reverse,sufi], (* ref.[2],eq.(3.26-28) *)

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

Φ2xyzr

M0c -> M0f
M1c -> M1f

Testing Continuity Conditions

Continuity of Potential at interface  η = ηp

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

te1 = Φ1xyzf[ηp]/. Map[Reverse,sufi]/.sep/.sufico ;

te2 = Φ2xyzf[ηp] /. sufico ;

te1 - te2;

Simplify[%]

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

Φ1xyzf[η]/ep /.ToEll/. Map[Reverse,sufi]//Expand ;

D[%,η]//Simplify ;

% /. sep ;

% /. sufico ;

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

D[Φ2xyzf[η]/ep/.ToEll,η]/. sep/.sufico//Expand ;

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

te3 - te4 ;

Simplify[%]

Revised Results

2012-07-25

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

L0f

L1f

M0f

M1f

Substituting semi-axes

suabe ={Sinh[ηp]->(ap/ep),Cosh[ηp]->(cp/ep),Sech[ηp]->(ep/cp), Csch[ηp] ->  ep/ap,Coth[ηp]-> cp/ap,  ArcCoth[(cp/ep)]-> dp}

L0e = L0f //. suabe

L1f //. suabe

L1e =Cancel[ % (ep^3)/{(ap^2) cp}]//FullSimplify

M0e = M0f //. suabe//Simplify

M1e = M1f //. suabe //Simplify

Results for IGTE12 paper

L0i = (((χ1-χ2)cp /ep)/(1+χ2 -(((cp^2) (1+χ1))/(ap^2))+((cp dp (χ1-χ2))/ep)))

L0e - L0i

L1i= 
(((χ1-χ2))/((2+χ1+χ2) - (χ1-χ2)((ap^2) /(ep^2) -  (ap^2) cp dp/(ep^3) )))

L1e - L1i//Simplify

M0e

M0i =- ((1+χ1)/(1+χ1+(((ap^2) (-cp dp+ep) (χ1-χ2))/(ep^3))))

M0e - M0i//Simplify

M1i = -((2 (1+χ1))/(2+χ1- (((ap^2) (-cp dp +ep) (χ1-χ2))/(ep^3))+χ2))

M1e - M1i//Simplify

Numeric Tests

suva = {ap ->  apv, cp ->  cpv, ep -> epv, χ1  -> χ1v,    χ2 -> χ2v,
μ1  ->1+ χ1v,    μ2 -> 1 + χ2v}

sunuv = {apv ->  1.17, cpv ->  2.23, epv ->  1.8984203960134858`, χ1v->  0.23, χ2v -> 0.77}

ep->(Sqrt[-(ap^2)+(cp^2)])/. sunuv

ap/cp/. suva/. sunuv

ηp = ArcCoth[cp/ap]/. suva/. sunuv