MAGNETIC FIELD DUE TO A CIRCULAR CURRENT

Thesis for ”'12 Degree OI M. S.
MICHIGAN STATE UNIVERSITY

Bartin T. Smith
3.960

mmmmmInwnuummu '

3 1293 01764 0073

 

L I B R A R Y
Michigan State
University

\CEINHET’IRN 8000* 4-H» hutfron y-u
' "OID Fl" ' manor

MAGNETIC FIELD DUE TO A CIRCULAR CURRENT

By

Bartin T. Smith

AN ABSTRACT

Submitted to the College of Science and Arts of
Michigan State University of Agriculture and
Applied Science in partial fulfillment of
the requirements for the degree of

MASTER OF SCIENCE
Department of Physics and Astronomy

1960

%/m j MW;
fl

Bartin T. Smith

ABSTRACT

Mathematical expressions for the magnitude of the magnetic field
due to a circular current are derived. The difficulty of evaluating elliptic
integrals of modulus near 1 is obviated by successive applications of
Landen's Transformation, resulting in expressions readily calculable
on computers. The expressions for complete elliptic integrals of the
first and second kinds are programmed, with explanations of order
pairs. Directions for change of parameters are included. Programs
for the magnitudes of the components of the magnetic field parallel and
perpendicular to the axis are presented, with explanations as to scaling

and use as a subroutine.

MAGNETIC FIELD DUE TO A CIRCULAR CURRENT

BY

Bartin T. Smith

A THESIS

Submitted to the College of Science and Arts of
Michigan State University of Agriculture and
Applied Science in partial fulfillment of
the requirements for the degree of

MASTER OF SCIENCE

Department of Physics and Astronomy

1960

The work on this thesis was partially supported by the
Atomic Energy Commission under Contract AT(11-l) - 872

ii

ACKNOWLEDGMENTS

I wish to express my appreciation to Dr. H. G. Blosser
for suggesting the problem and for his ready assistance and
advice. I would also like to thank Dr. M. M. Gordon for
making available his knowledge of computer programming

through his many helpful suggestions.

II.

III.

IV.

TABLE OF CONTENTS

INTRODUCTION

DERIVATION OF MAGNETIC FIELD MAGNITUDES

TRANSFORMATION OF ELLIPTIC INTEGRALS .

PROGRAMS

A. Subroutine for K, E .
B. Routine for BP . . . .
C. Routine for B

Z

iii

18

18

22

26

I. INTRODUCTION

In many physical experiments it is desirable to know the magnitude
at any location of a magnetic field due to a circular current. Though
the mathematical expressions for the magnitudes are well known, the
numerical difficulties involved in the evaluation of these expressions in
any given case may be quite large.

It is the purpose of this thesis to consider such a case, and to
show whereby the numerical difficulties may be alleviated. The mathe-
matical operations necessary are readily adaptable to computer pro-
gramming, and use has been made of the Michigan State University
digital computer, Mistic.

The need for the ready calculation of magnetic fields due to
circular currents arose in connection with the problem of current
control in the trimming coils of the proposed Michigan State University

cyclotron.

II. DERIVATION OF MAGNETIC FIELD MAGNITUDES

N

 

 

 

Fig. 1

We wish first to write an expression for the vector potential at
any point P not on a coil of radius L. Noting that there is no loss of
generality if we let P be in the x-z plane, we may write for the distance

of P from any line element ds (see fig. 1):

 

 

1/ 2 2 ‘ 2 2 2 2
r: L +P -2LPcose+z = «(IO-Loos 9) +(Lsin9) +z
Using the expression for r, the vector potential is then:

ds

r
A :Eif
_ 4n CW/(IO - L cos 9)2+(L sin9)2 + z2

 

 

When we consider opposing elements on opposite sides of the x axis,
we see that there is no contribution from x components of the current,

since the resultant current (21 cos 9) is perpendicular to the/0 - 2 plane.

Consequently, A]O = A = 0. So since ds = adcp,

1T

A :PoI LcosGdG
o (P-LcosG)2+(Lsin9)2+z

 

 

9 211' 2

2
In the manner of Smythe (5), let 9 : 1T + 24>, so cos 9 = 2 sin cp-l, and

d9 = qu) so we have

 

 

 

 

 

 

1r
/— 2
A z’JoLI 2 (2 sin (buds
9 11 2 .1.
{(L+)O)Z+zz} 1— 41")“; (I); 2
o (Li-P) +2 I
2
Defining k E 4L; 2 , and substituting, we get
(L +p) + z
i g 1
Ik Z , 2 2
A _po (L) 2 81H (19qu _ dL [1]
9‘ 211’ p 2 2 1/2 2 2 .1/2

 

 

Examining the first integral, we note that

 

 

 

 

2 i
511143 1; £2 1 .. [1-k sin ¢]
_ k '—
2 2 2 2 2
[l-k sin .¢]2 [1-k sin ¢l
so that
1 1 1
m ¢ ¢ ; _ 2 - — [l-k sin ¢] dqa =—(K-E)
1 2 1 2 2
_ k — k k
0 0

Z
O [l-k sinch]2 [l-kzsinch]Z

where K and E are complete elliptic integrals of the first and second

kinds, respectively.

4
The second integral of equation 1 is K, giving us a final expression

for the vector potential:

1
pl 2. 2
o L k
Ae'nkIfi) 1'7K'E

 

The magnetic induction vector may be calculated from B = VXA.

Since there is no component of A in either the Z or/O direction, this

 

 

becomes
B : -éaAe +3 3(pAe)
— az )0 W—
aAe 1 8(PA9)
orB 2- 82 “1desz [2]
Rewrite A

G in the form:

 

 

 

 

 

1
)1 I E
A9: 0 [(L+/O)2+z] 1- 2142/9 K-E
)0 [(L+)O) +z ]
.1
3A p I 2 \
Tie—‘2; z[(L+p)2+zz] 1-[ :5 2]) K-E
(L )+2
1
F01 L 2]E 1 ZLP \BK§_15_+K 4apz [@513
J'Efi [I +10) +2 -[(L+/O) +2 1} 8k 62 '

 

 

—-_ .1. .1. 5
L 2 5k 2(L )‘2
Also, since k = 2 ’0 , we have —- = - [O z
(L+p)2+z 57‘ 2 2 3
-- ‘ [(L+p) +2 ]2
and £21: _ _k_ 1.5: .15
:30 " 2p 4p ' 4L
Substituting these expressions in Eq. [2] and simplifying, we
obtain:

}1 Iz
13/0: 0 -K

2 2 2
L
1 + +p +2 E [3]

(L-P)2+ z2

 

2 2 2
2P1T[(L+P) + z ]
For the z component, it can be shown similarly that

)11 2 2 2
Bz= O 1 K+L P2 22E [4]
3:- (L-P)+z

 

 

III. TRANSFORMATION OF ELLIPTIC INTEGRALS

Using the preceding formulae for BP and Bz, one may calculate
the modulus k and refer to tables for the corresponding K and E.
However, since in a practical problem, k is not likely to be an even
number with its complete elliptic integrals listed, it is necessary to
interpolate between known values or to calculate them directly. With
the availability of automatic calculating machinery, this latter course
is preferable since it provides greater accuracy with a negligible
sacrifice of speed.

In the regions where k is small, that is, whenf) is not near L
or when z is large, the power series representation would be adequate,
since its repetitive form lends itself well to computer calculation.
When k is large, however, a situation that occurs at points near the
coil, the series representation converges very slowly.

The power series is obtained by expanding the radical in a

binomial series:

NI=I
NI:

{'

d4) ,= Z .
————+—— = dcp 14-1-<—sin2q>+-}—%k4sin4¢+. . .
.1] 2 , 2 2 2'4
0 l-k Sln q) -

K

Hg

1 3 . . . (Zn-1) 2n _ 2n
+.
2-4. . . (2n) 1‘ 5m ‘1’

 

+

NI=I

. (Zn-l) _1_r_
(Zn) 2

 

2 1- 3 .
Since by Wallis' Theorem: sin nedcp + 2- 4
O .

we have:

2 2
1 2 .
b-ttg) k +(l—3) k4+... +

 

1T
K-‘z'

1° 3' 5. . .(Zn-l)‘ 2 2n+
2.4

2°4'6...(2n) k

 

. . . . . . th .
As an illustration of 121115, con51der the coeff1c1ent of the n— term in
the series for K:
2 2

'1 3-5 ...(Zn-1I _ (2n)!
2 4 6... (2n) " 22n(n')2

 

. . . . . - ‘— n -n
Usmg Stirling's approx1mation, n! : ‘_/21Tnn e , we may show
that this fraction may be written l/nTr. It can be seen that if this is
. . 2n . ' . .
the coeff1c1ent of k , the result 15 a very slowly converging series
when k is nearly 1.
That k can never be greater than 1 is true by definition. The
modulus k is equal to the numerical eccentricity of an ellipse. For
example, in calculating the arc length of an ellipse, the situation from

which, according to Hancock (4), the "elliptic integral" derives its

 

 

 

 

2 2
name, we have: ifE-z- +y-2- = 1 where a and b are the major and minor
a b
axes respectively,
x X a2 a2- b2 2
— 2 I X
’d 2 a
: - 1 ‘1) Z —.
s O / -+[dx dx 2 2 dx
o a - x

a2 b2 b2

Upon introducing the numerical eccentricity e = —;—2—— = 1 - -—7:,
a a

we have, with x = a sinq>,

 

{4)
S = aJ W/l - ezsinch do

0

which, it can be seen, is an elliptic integral of the second kind.
2

Since 1 - —2- can never be greater than 1, the series will converge,

a
b .
however slowly. The quotient g- is called the complementary modulus

k‘, and setting e = k, we have k2 + k'2 = l.

A method of changing a slowly convergent elliptic integral into
a highly convergent one is the method of Landen's transformation.
Writing now F(k, (p) for the incomplete elliptic integral of the first
kind, we wish to find a F(kl’ cpl) such that it will be more rapidly
convergent than F(k, (1:) and can be related to it by F(k, W = C F(kl' cpl),

where C is a constant of proportionality to be determined.

 

2 2
Instead of using the radical-\[l -k sin ¢, substitute for the
2 . . 2 2 2
modulus k the eccentr1c1ty (a - b )/a , where, as before, a > b.

This give 5:

2 2 1 2 2 2 2
“/l-ksin ¢=;_\/a cos ¢+b sin¢

 

 

and now: F(k, <13) 2 aF(a, b, (p) and E(k, (p) = ’a-1 E(a, b, (1))

With this notation, we shall consider a geometrical derivation
of Landen's transformation as given by Cayley (1).

On a circle of radius AOB, let P be any point, and Q be any
point on the diameter other than the center, as in fig. 2. Considering

a, b, and c1 as noted in the figure, we may write
1 . .
a1 = % (a+b), b1 = \/ab, and c1 = E (a - b) where a1 is the radius,

1
andOQ—al-b---2-(a-b)—c1

 

 

 

 

 

Fig. 2

We see that QP sinq>1= a1 sin Zcp and QP cos cpl = c1+ a1 cos 24>.

Upon substituting the preceding expressions for a1 and c , we

1

 

 

 

 

2 2 2 2 2
obtain: OP = a cos 4) + b sin 4) so that
a sianp c +a cos 24>
. l l l
Sincpl = , cos (pl [5]
_\/ 2 2 2 , 2 V 2 2 2 ,
a cos<p+b51n¢ acos¢+bsm¢
and: 2

2( 2+b'2)
2.2_a1acos¢ Sinq)
15mq" ’5 2 2,2
acos¢+b81n¢

 

2 2
a1 cos q>+b

We may write:

sin(2q>-¢1)= sin 24>cos ¢1_ c0524) sin 4,1 : (a-b) sin a)

 

 

ZV—z 2 2_2
acos¢+bsm¢

and

cos(2q'>-¢1)= cos 24> cos (121+ sin2<p sin 491

2 ,2
acos ¢+bsm ¢

2 2 2
a cos ¢+b2sin ¢

 

1 2 2 2 Z
= a cos Cb + b sin CP
al 1 1 l 1

Consider the point P' which is on the circle next to P. If ¢1
is incremented by d¢1i then:

Pqu)1 = PP' sin P'PQ = PP' cos(2¢ - cpl)

and sinc e

PP' = a1d(2¢) = Zaldcp

2 alcos( 243-491) (14; = PQdcp1

 

 

10

11
Substituting for PO and cos (2¢-¢l) we get:

 

 

 

 

chp _ d¢
2 Z 2. Z 2 Z 2 . 2
a cos q>+b Sincp a1 cos ¢1+b151n cpl
Upon integrating we have
(I)
<I> d
If d<1> = 1 4’1

 

 

 

 

2 Z 2 2, 2 Z 2 2.2
acos¢+bsm¢ acos¢+b Sinq>
0 o l l l l

01'

l
F(a, b, (p) = EF(a1, b1, cpl)
Mk. 4» = ama. b. s) =§F<a1b1¢1) = i?) F(k,-pl)

which is the desired direction.

1 - k' a - b 2a
W ' ' z = ' 1 =
riting k1 1 + k' a +b , we have, Since + k1 a +b

 

 

 

 

 

 

(1+k1)
F k, = F k,
( <1) 2 < 1 4.1)
When (I) = %, (bl = 17, and the complete integral is:
w l+k
K(k.-2-) = 2 F(k,-r)
(1 +k1) (1 +k2)
- 2 2 F(kl, 211)
F(kn, 2n'ln
=(1+k1)(1+k2). . . (1+kn) n
2
l-k' 1
— n-
where kr1 — ———1 +k'
n-l

It will be recalled that al, b1, and c1 were derived from a, b,

and c. Similarly a2, b2, and c2 can be derived from-a1, b1, and c1.

12

This is equivalent to calculating the succeeding kn.

1 1

a1 -— E(a-I-b) b1 —- Vab C1 — E(a-b)

a -l( +b) b = ab --1-( b)
2-2311 2\I11 C2'2‘aL1'1
—i( +b) b — b —l( b)

a3‘2a2 2 3’\Ia22 C3‘2"“‘2'2

It can be seen that as n increases, a and b approach the same
n n

limit. Hancock (3) shows that

a_b _(\/Z-\/?)2
11‘ 2

and

a1+b1

a1'b1
az'b2: 2 'fi:—2_'(Val_vgl—)VEI

 

SO

 

 

OI'

a-b 2
a-b<11 a_b<(\/2—2‘\/E)

 

 

2 9
2 2 2 2 2
Also we see that
-b 2
a2 2 (VS-VB)
a-b <———<
3 3 2 3
2
and in general
(‘2 V132
a-b< "' ,orlim(a-b)=0
n n n n n
2 n—aoo
2
2 b n-1 n-1
Whena =b,k =1-—n—2-=o, sothatlimF(k,2 n),=2 11'
n n n a n90<> n
n

1T
andK-(l+kl)(1+k2) ..... (1+kn)2

13
As a further illustration of the convergence of this series, we

examine the log of the infinite product:

an =ln(l+kl)+1n(1+k2)+1n(1+k3)+ ....... +ln(l+k )+ . . . .
n

When kn is very small, ln(1+kn):::' kn’ and the ratio of successive

terms will be:

 

 

 

 

1
Letting k = —, we see that
n m

 

lim = 0

Thus we see that continued application of this transformation will
reduce an elliptic integral of any modulus k to a continuous product
representation. Since the terms (1+ki) are calculated the same way,
the method is readily adaptable to such machines as the Michigan State
University digital computer, Mistic. In fact, the rapidity of convergence
of the series for any k quickly exceeds the capacity of this machine, as
it would any other. An example of the fast convergence is given in

the following table:

O.

.—|

WWWWW‘WF‘WW
00000000

999999

.999990
.991391
.768452
.219581
.012353
.000038
.000000
.000000

999990
655984
303175
352529
346615
653210
156098
000365
000000

14

CDNO‘thU-JN

In order to derive an expression for E in a more rapidly converging

form, we begin with equation 5:

al sin 24>

 

sin (bl =

 

.\[ 2 2 2 _ 2
a cos ¢+b Sln q:
Squaring.we get:

2 , 2
2 4&1 Sln ¢cos <1)
sin cpl =

 

2 2 2 _ 2
a cos ¢+b Sin q>
. 2 2 2 _ 2 . .
Letting D = a cos (1) + b Sin q> throughout the follow1ng, we write:
2 2 2 2 2 2
D-a = -a sin ¢+b sin cp=(b2-a )sin q»

2 2 2
(a -b )cos cp

U
0"
11

(D - a2)(D - b2) 2 -(a2- b2)(a2- b2) sinch cosch

2 2 , 2 2.
-4:a1 (a-b) 8111 (1) cos q)

2 .
Since this expression is [-(a-b) ] times the numerator of the quotient

2 2 2 2
above, we get: (D-a )(D-b )+ D(a-b) sin 431: 0.

15
, 1 Z 2 2 . 2 .
Adding [E(a +b )cos 421+ ab Sln cpl] to each Slde, we get. after

much simplification:

l 2 2 Z , Z 2 2 2 2 2 2
D - [-2-(a +b )cos ¢1+abSin cpl] = 4 C1 cos (131(a1 cos cpl‘l' b1 sin (1)1)

Replacing D by its value, and rearranging, gives:

 

2 2 2 2 2 2 2 2 2 2 2 2 2
+b ' = 2 + ' — +2 + ' '
a cos q: 8111 ¢ (a1 cos $1 b1 Sln cpl) b1 Clcoscp1 a1 cos cpl b181n q>1
Multiplying this equation by

d¢ _

 

 

 

1

2d¢l

T2 2 2 2 \/ 2 2 2 2
+ ° + '

i/al cos cpl b1 Sln cpl a1 cos cpl b1 5m ((11

we get:

 

 

 

 

l 2
—b
2 2 2 2 2 2 2 2 2.
do a cos cp+b sin cpzdcp Ja cos ¢+b sin q: - 1 +C coscp
l l l 1 l 2 2 2 2 1 1
9L1 cos q>1+b1 Sln cpl

Upon integrating, we obtain:
1 2
: 9 "' — 1 b D .
E(a, b, (p) E(al, b1 cpl) 2bl F(a1 1 cpl) + Clsmcp1

From before, we had:
1
Ma. b, 41) = ;F(k. <1), E(a, b, <1) = aE(k. <1»)
Using these expressions, we can write the transformed integral as:

a b2 C

___1 1 _1 .
E(k’¢) ’ a E(k1’¢l) ' zen1 F(k1,q>1) + a sm‘h

Continuing, however, with the integral in the form above:

b2

1 .
E(a, b, (p) — E (a1, b1, (1)1) - —2— F (al, 131, (pl) + C151n (pl

16

we rewrite it as:
E(b)2F(b)-IE(b)2F(b)] F(b+'
6" "1’21 a' "I —. aL1' 1'¢1’a1 31' 1'¢1 'aici 2:11’ 1’45) Cis’m‘I’i
2

2
b 1

2
which, since a - 92—- - —;—- = -—(a2- b2) = -a c

1 4 1 1, is equivalent to:

E(ab¢)- a2F(ab¢)=[E(a1b1¢l) - a2

1 F(albl¢1)] - a ch(alb1q>l) + C1 sin (pl

1

As n increases, we see that

. 2
11m [E(an, bn, in) - an F(an, bn, en” = o

n—>oo
2
2 b
sincea =b meansk =1-—r-1—=0
n n n a2
n

_ fl _.
and F(an, anPn) " an 1 E(an! bn! (151,1) _ an¢n

Now we can write, substituting for the expression in brackets its

value obtained from the same equation with the subscripts increased by l:

1C1F(albl¢)

E(a, b, ¢)-azF(a, b, ‘1’) = [E(aZquDZ) - a22F(a2b2¢2)] - a
+ c1 sinc|>1 - a2c2F(a2b2¢2) + casinq)2

= [E(az, b2, cpl) - a22F(a2. b2, ¢Z)]-(2a1c1+4a2c2)F(a2b2¢2)+c1 sin ¢1+c2 sin cpz

As n —> 00, we have:

E(a,b,¢)-a2F(a,b,¢) = -(2a1c +4a2c2+8a c +. . . )F(a,b, <t>)+c1 sin¢1+. . . .

1 3 3

Since in our problem, 4) 2172., ((11 = 11', (132 = 2n, . . . . , (p = 2 it,

17 2 TI'
_ = - - 4 -0 o o p ,—
E(a, b, 2) (a 2a1c1 aZCZ )F(a b 2)

17

but
iF<k¢>=F<abm=iF<a b ¢)= =—1—F(a b «M
a ' ’ ’ 2 1'1'1 2n n’n'n
n-l TT

. 1
and Since (1) = 2 1T, we have: —K =—
n a Zan

With this expression, and E (k, 4)) = aL E(a, b, (p)

 

 

 

 

 

 

 

 

n
we may restate E(k,g—) as:
2ac 4a c
‘n' _ 11 2 2 1r
E(k.§)- 1- 2 - 2 K(k,-2-)
a a
2 2
a1C1_ aZ-bz k2 aL2C2 aL1'b1 k1
Weseethat 2— 2=—4and c =__2_.Z:.:1_
a 4a a1 1 a -b
2
ac2_a1c1 a22_k k1 a3c3 _ 1
so that - -—— and also - —k
2 2 a c 4 4 a c 4 2
a a 1 1 2 2

Succeeding terms are calculated in the same manner, allowing

the final expression to be:

k
12 1 1 l
= -— — — _kkk 0.0
E [1 2k(1+2+4k1k2+8123+ )]K

Inasmuch as the k must be calculated in order to determine K,
n
this form of E is easily evaluated without the necessity of complex

programming, since each additional term of E consists of a number

obtained from one multiplication and one right shift.

18

IV . PROGRAMS

A. Subroutine for E, K

The following is a program for the computing of K and E. This
is a subroutine which is entered by a standard entry after the modulus
k has been placed in the accumulator. K and E are then calculated in
less than 40 milliseconds and placed in the accumulator and quotient
register, respectively. A total of 59'Vmemory locations are used to
store this subroutine, and locations 0, 1, and 2 are used for temporary
storage.

In the state in which the subroutine is presented here, the
subroutine gives K and E scaled by 20, 10, or 5, but the scaling can
be easily changed to a multiple of two by replacing the (. 2w) by 1r/ 4,
11/8, etc., depending upon the range of the k to be used.

The square root subroutine from the Mistic library has been
incorporated in this subroutine, and can be entered at (p+45), where
the K, E subroutine begins at p. The standard entry is used and
exit addresses will be computed normally.

As stated before, the representations for K and E to be pro-

grammed are:

 

1T

: 1 O... —

K (1+k1)( +k2)(1+k3) (1+kn)2

2 k kk kk k

_ k 1 12 12 3
E _ 1- 2 (1 +7 + + 8 + ) K

19

The K and E corresponding to sixteen randomly spaced k's were

calculated on Mistic, the results agreeing to 10 places when k < 0. 999

and to 8 places when kN O. 9999 with the values listed in Complete

Elliptic Integrals, C. W. K. Gleerup, Lund University Press, 1957.

 

0

10

11

12

13

14

15

16

17

402F +5F
4236L L5 2F
4041L 41 38L
4139L 4140L
502F 7 JZF
102F 40F
191F LOF
40F 507L
2645L 40F

LJF 401F

L9F 661F
-5F 402F
101F 40F
5038L 7JF
4038L L439L
4039L LJF
401F 5040L

7J1F 4040L

 

2 2
k (or k )
n

1 kg

4 ' 4
k1 (or k 1)

N

k k
_1_ _2

2 2 k k k 5 This will be done on

.1. +__1_ __2
2 2 2 second passage. First
k
.1. , __1_
2 2 passage multiplies by
zero.
1 .51 1,112
2 2 2 2

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

*Scaling counter n must be picked such that

L537L 3222L

L5F 4038L

4039L LJF

4040L L137L

4037L L542L

L043L 4042L

3626L L540L

001F 4040L

L5F LO37L

364L 5041L

7J41L 101F

4041L 5039L

7J41L L441L

4041L 2659L

LJ41L 4041L

5040L7J44L

40F 5041L

7JF 401F

26 54L 22L

LL409 5FLL409 4F

00F 00F

00F 00F

-3

20

(.-2 8) < O? For one passage only

Change sign of (-2-

Count for two jumps’k

 

 

 

1,51 (1,52 (.1,
2 2 2 2 2
k

n -38

—<

2 2

k_2

K

E

38

k
3

7

VI
N

k

. n
Test constant for Size of ——(used as

po sitive number)

2

 

K
n-2

2

<1, eg.

TT k>0.99999, use 3
0.99999>,k>0..9841 , use 2
0.9841 > k>0.802 , usel

4O

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

00F 00F

00F 00F

00F OOZF

00F 001F

21

Scaling counter n

40F 00 1283 1853 0718.1 .ZTr

401F+5F

4253L511F

lOlF-JF

402F 50F

, L51F662F

-5FL02F

101F3653L

L42F2648L

L52F22F

L557L 4037L

L558L 4042L

2661L 2661L

LL409 5F LL409 4F

00F 002F

L941L 4039L

LJ39L 2232L

L5F 501F

2236L 2236L

7

 

R1 square root subroutine from
computer tape library

Restore counters

Counter restoring constants

22

B. Routine for BP

The routine for BPis straightforward, and provisions have been
made for any normalization the user may require.

The formula

)1 Iz 2 2 2
BJO = O -K +£23.22 E
(L-P) +z

NIH

2P11[(L+p)2+ 22]

k 1

ZLP=
I

 

has been changed with the aid of and»O =(4Tr)10-7

.1.
2 2]2
(L4?) +z

to read:

2 2 2
BF kz K+L_+£+z E

110’7 VLp (L-P)2+z2

 

 

Care must be exercised when this code is used because of the
coefficient of E. WhenPR‘eL and Zis small, a large number results.

For example, when Z = J- L andP = L,

 

25
2 2 2 2+—1—
L_+p+z 625
= 1251
L 2+2 1
I 7p) “ 625

Clearly, the numerator must be scaled before the division takes
place. A scaling factor of 10-m has been used, with m determined by
the user as appropriate for the particular location where the field is
to be calculated. pr never approaches L, or if Z is large, then a

smaller m would be sufficient, but if Z is small relative top when/O

23

is near L, then m will be large. The expression is then:

L2+Pz+z2 : 2L2+zz

(L-P)2+z2 22
Using this expression, an m can be chosen to fit the case at hand.

The first three words of this and the following program consist
of carriage returns and delays to facilitate a change to a subroutine,
should the user so desire.

Five random values of Bf) were calculated on a desk calculator,
the results agreeing with the Mistic calculated results to 9 places.

The routine for Bp will stop if the scaling factor is too large

by hanging up at the division order at the 19th word.

0 921311? 9231‘
1 921311“ 923]?
2 921311? 9231-“
3 L542L L443L (L+p)

4 40F 50F

2
5 NF 40F (L+p)
6 50441.. 7.1441. 22

2 2
7 L4F 40F (L+p) +2

8 5042L 7J43L LP

9 66F -5F Lp
2 2
(L430) +z

 

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

40F 5010L

26 KE + 45F 2650L

- 5F 4049L

L542L L043L

40F 50F

7JF 40F

5044L 7J44L

L4F 40F

5049L 7J45L

66F -5F

40F 5042L

7J42L 404F

5043L 7J43L

L44F 404F

5044L 7J44L

L44F 404F

504F 7JF

40F 5048L

7J45L LOF

405F 5042L

7J43L 5030L

Store E

<L-p)

2
(L-p)
2

Z

(L-JD)2+z2

E x (scale factor)

 

LIE}.

(L-p)‘2 + z2

2
L

2

2
L+p

2 2 2
L+P +z
Z 2 2

2 (E) (scale factor)

(L-p) +z

[( )K -( )EI
LP

26 KB + 45F 40F VLfD

24

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

50 54L 7J5F
66F 7041L
6643L 7J44L
52115F 5035L
26P1F L547L
L043L 3639L
OFFOFF
L543L L446L
4043L262L

B (norm)

L (radius)

p

scaling factor
AP
p max

00F 00F

00F 00F

00 1 F 4054L
4054L 5051L
26K, EF 4048L
2612L 2612L

OOFOOF

25

k[( )K-( )E] ZanormU )E'( )K]

 

 

i
O

II

69

Via

Print

Comparepmax-p < 0?

Increase}?

k
Transfer to K, E

Store K

26

C. Routine for Bz

As before with BP’ the routine for B2 was programmed in a
straightforward manner, the only trouble being in calculating the
coefficient of E.

The expression to be evaluated is, since )10 = (4w)10-

Z 2 2
z = 2 K+L 'Pz'zz’ E
[ 2 2] (L—p)+z
(L+p)+z

Unlike the coefficient of E in BP' P = L presents little difficulty,
1

 

NIH

as long as Z is not zero. When ZFSJ‘(1/25)L, a scale factor of 10-
is sufficient, since E is already scaled by 10 in the range 0. 99999
>k>0. 9841, and a smaller k would mean even less difficulty, since
the coil would not be approached so closely.

When Z -> 0, however,

2 2 2 2 2

L -B -z L -P L+p

2 2 —2 2 : L-P
(L-P) +z (L-p)

Keeping gg < 1E0 (scaling factor) will suffice when z << L.

-l -2

For all but regions extremely near the coil, 10 or 10 would be
ample scaling.

As with g), allowance has been made in the program for normal-
ization with respect to any location the user might choose.

The values of B as calculated here agree to 9 places with those
z

calculated by M. M. Mosely of Texas Christian University.

The routine for Bz will hang up if the proper scaling factor is

not picked by failing to execute the divide order at the 28th word.

10

11

12

l3

14

15

16

17

18

19

92131F 923F

92131F 923F

92131F 923F

L542L L443L

40F 50F

7JF 40F

5044L 7J44L

L4F 409F (L-I:F))Z+z2
5042L 7J43L

669F -5F LP
2 2
(L+p) +z

 

40F 5010L

26KB +45F001F k1(or kn)
40F 5012L

26KEF 4049L Store K
-5F 4050L Store E
L542L L043L

40F 50F

7JF 40F

5044L 7J44L

2 2
L4F 402F (L-P) +2

27

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

5045L 7J50L
401F 50 43L
7J43L 40F
5044L 7J44L
2 2 2

L4F 40F L -(p +z )
5042L 7J42L
LOF 40F
50F 7511?

2 2 2
L -(P +2 1

(L-P)2+ 22

(E) (scaling factor)

 

662F -5F

408F 5049L

7J45L L48F [( )K +( )E]
408F L59F

L59F 5032L

 

26KE+45F 409F '\((L+p)2+ 22
L58}? 001F

669F 7J48L Bz

521151“ 5036L Print
26P1F L543L

LO46L 3641L

L543L L447L

4043L 262L
OFF OFF

L (radius)

28

43

44

45

46

47

48

49

50

p

Z
scaling constant
pmax
AP
B (0)
z
00F 00F

00F 00F

Temporary storage for K and E

29

30
BIBLIOGRAPHY

Arthur Cayley, An Elementary Treatise on Elliptic Functions,
Deighton, Bell, and Company: London, 1876.

 

Herbert Bristol Dwight, Tables of Integrals and Other Mathe-
matical Data, The Macmillan Company: New York, 1949.

 

 

Harris Hancock, Elliptic Integrals, Dover Publications, Inc. :
New York, 1958.

 

Harris Hancock, Theory of Elliptic Functions, John Wiley &
Sons: New York, 1910.

 

William R. Smythe, Static and Dynamic Electricity, McGraw-
Hill Book Company: New York, 1939.

 

MICHIGAN STRTE UNIV. LIBRQRIES

 

31293017