p21.html

> q:=(R,rho,lambda)->sqrt(R^2+rho^2+2*rho*R*cos(lambda));

(1)

> k:=(q,z)->2*sqrt(q)/sqrt((1+q)^2+z^2);

(2)

> H:=(k)->((2-k^2)*EllipticK(k)-2*EllipticE(k))/k;

(3)

> Kern:=(R,rho,lambda,z)->(1+(R^2-rho^2)/q(R,rho,lambda)^2)*H(k(q(R,rho,lambda),z))/sqrt(q(R,rho,lambda));

proc (R, rho, lambda, z) options operator, arrow; `/`(`*`(`+`(1, `/`(`*`(`+`(`*`(`^`(R, 2)), `-`(`*`(`^`(rho, 2))))), `*`(`^`(q(R, rho, lambda), 2)))), `*`(H(k(q(R, rho, lambda), z)))), `*`(sqrt(q(R, ... (4)

Phase shift:

> Phi:=(R,rho,z)->evalf(Int(Kern(R,rho,lambda,z),lambda=0.0..2*Pi));

(5)

Hand-crafted numerical guidance

Unfortunately Maple cannot handle this integration - even numerically. Hence we must help out with pre-evaluation of H(k). This will have the additional advantage of speeding up the integration in the case that the value k is repeated during evaluation thereof. First find maximum value of k:

> diff(k(q,0),q);

(6)

> solve(,q);

(7)

> k(,0);

(8)

So k is in [0,1]. Now tabulate H over that range. First get limiting values of H:

> limit(H(k),k=0,right);

(9)

> limit(H(k),k=1,left);

(10)

Hence q=0 => k=0 => H=0 and q=1 => k=1 => H=infinity. These extreme values cause a problem in the integration. Ideally, the integral should really be re-formulated so that integration through these values does not contribute an infinite part that presumably, due to sign change on passing through these points ultimately gives a finite value. They correspond to physical coincidence of detector ring path and source. For now I will assume that no damage is done if they are simply excluded. Justification for this approach will come from plots which do not show any pathology as the geometry is varied so that the detector passes through the source.

Now pre-compute H(k):

> with(LinearAlgebra):
Num_k_pts:=20000; # to precompute H
Num_lambda_pts:=5000; # to do integration by hand

(11)

> make_H:=proc()
local i, res;
res:=Vector(Num_k_pts);
for i from 1 to Num_k_pts do
res[i]:=evalf(H(i/(Num_k_pts+1))); #do not pass zero or 1 as arguments - see above comments
end do;
return res;
end proc:

> H_lut:=make_H();

Typesetting:-mrow(Typesetting:-mverbatim( (12)

Now make a pseudo-function that is a nicely-behaved version of H.

> evalf_H:=proc(k)
local kt;
kt:=round(k*Num_k_pts);
if kt<1 then kt:=1 elif kt>Num_k_pts then kt:=Num_k_pts end if;
return H_lut[kt];
end proc:

And redefine Kern to call evalf_H(k) instead of H(k)

> Kern:=proc(R,rho,lambda,z)
local qt,kt;
qt:=evalf(q(R,rho,lambda));
if (qt<1/Num_lambda_pts) then qt:=1/Num_lambda_pts end if; #This deals with the pathology at q = 0.
return (1+(R^2-rho^2)/qt^2)*evalf_H(k(qt,z))/sqrt(qt);
kt:=k(qt,z);
end proc:

Check it is working:

> evalf(Kern(1,1,1,0));

(13)

And do the integration numerically 'by hand':

> Phi:=proc(R,rho,z)
local lambda,s,d_lambda,i;
lambda:=0;
s:=0;
d_lambda:=2*Pi/Num_lambda_pts;
for i from 1 to Num_lambda_pts do
lambda:=lambda+d_lambda;
s:=s+evalf(Kern(R,rho,lambda,z));
end do;
return s*d_lambda;
end proc:

Plots

Plot of detector reading (flux through loop) when:
R = 1 (=> FOG radius = source ring radius) and
z = 0 (in plane with source)
for various rho (rho=0 => no displacement, rho=2 => detector path is wholly exterior to source, just grazing source ring at a tangent).

First define 1D function:

> f0:=(rho)->Phi(1,rho,0);

(14)

> plot(f0,-3..3,
labels=["horizontal displacement of sensor","flux"],font=[TIMES,NORMAL,10],
title="In-plane flux versus offset when FOG radius = source radius");

Plot_2d

> f1:=(rho)->Phi(1,rho,1);
f2:=(rho)->Phi(1,rho,2);
f3:=(rho)->Phi(1,rho,3);

(15)

> plot([f1,f2,f3],-3..3,
labels=["horizontal displacement of sensor","flux"],font=[TIMES,NORMAL,10],
title="Flux versus offset for various heights when FOG radius = source radius",
legend=["z=1","z=2","z=3"]);

Plot_2d

> f0:=(z)->Phi(1,0,z,Num_lambda_pts);

(16)

> plot(f0,-3..3,
labels=["height of sensor","flux"],font=[TIMES,NORMAL,10],
title="Flux versus height when no offset and FOG radius = source radius");

Plot_2d

> f1:=(z)->Phi(1,1,z);
f2:=(z)->Phi(1,2,z);
f3:=(z)->Phi(1,3,z);

(17)

> plot([f1,f2,f3],-3..3,
labels=["height of sensor","flux"],font=[TIMES,NORMAL,10],
title="Flux versus height for varius offsets when FOG radius = source radius",
legend=["rho=1","rho=2","rho=3"]);

Plot_2d

> frz:=(rho,z)->Phi(1,rho,z);

(18)

> plot3d(frz,-3..3,-1..1,axes=boxed,labels=["offset","height","flux"]);

Plot

> with(plots):

> Num_frames:=50;

(19)

> P:=animate(plot3d,
[frz,-3..3,-1..1,
axes=boxed,
labels=["offset","height","flux"],
orientation=[t*360,t*360/2]],
t=0..1-1/Num_frames,frames=Num_frames):

> P;

Plot