A theoretical study of the electromagnetic fields in a limb, excited by artificial sources

A theoretical study of the electromagnetic fields in a limb,

excited by artificial sources

Citation for published version (APA):

Amelsfort, van, A. M. J., & Scharten, T. (1986). A theoretical study of the electromagnetic fields in a limb, excited by artificial sources. (EUT report. E, Fac. of Electrical Engineering; Vol. 86-E-156). Eindhoven University of Technology.

Document status and date: Published: 01/01/1986 Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers) Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

providing details and we will investigate your claim.

j

I

A Theoretical Study of the

Electromagnetic Fields in a

Limb, Excited by artificial

. Sources

by

AM.J. van Amelsfort ar,d Th. Scharten

EUT Report 86-E-156 ISBN 90-6144-156-0 ISSN 0167-9708 February 1986

Department of Electrical Engineering Eindhoven The Netherlands

A THEORETICAL STUDY OF THE ELECTROMAGNETIC FIELD IN A LIMB, EXCITED BY ARTIFICIAL SOURCES

by

A.M.J. van Arnelsfort

and Th. Scharten

EUT Report 86-E-156 ISBN 90-6144-156-0 ISSN 0167-9708 Coden: TEUEDE Eindhoven February 1986

8'-'

'; i~t~

IOTH'--"

I t:..~:-".(\ "_,.A _ _ ,,,-,_ 0

ilO3391

0 ----~

T.

~

.. ;

Ef~JDHC!\-'r:N

I

I

i

I

1 '"-":,-: .~~ I~-' , -d-'~'

'-1

" .. , '. ! . .' . :~ -.-~ J

i

f.

l ! ' ',' , " '1',

-=,-,J,

i "I .... , l ' \ .... 1 ~

I

,

CIP-GEGEVENS KONINKLIJKE BIBLIOTHEEK, DEN HAAG

Amelsfort, A.M.J. van

A theoretical study of the electromagnetic field in a limb,

excited by artificial sources I by A.M.J. van Amelsfort and

Th. Scharten. - Eindhoven: University of Technology. - Fig.

(Eindhoven University of TeChnology research reports I

Department of Electrical Engineering, IS5N 0167-9708, 86-E-156) Met lit. opg., reg.

ISBN 90-6144-156-0

5ISO 539.1 UDC 537.8:615.84 UGI 650

ABSTRACT

I:SL!lg integral transfurm techniques, a procedure ,is developed to computt: the

electromagnetic field both inside and outside a r.ircularly cYlindrical limb,

pxclt.ed by any electric or electromagnetic source. The ~omponent layers of the

limb are assumed to b", homogeneous, anisotropic and dissipative. As an illus"

~ratj on, the time'-harmonic field of an axial pair of ring-shaped electrodes,

located symmetrically inside the limb, is calculated. The influences of

anisotropy and boundaries on field distl~ibution and surface charges are inves-tigated.

Amelsfort, A.M.J. van and

Tbl

Scharten .

A THEORETICAL STUDY OF THE

E~CTROMAGNETIC

FIELD IN A LIMB, EXCITED

BY ARTIFICIAL SOURCES.

Department of Electrical Engineering, Eindhoven University of Technology (The Netherlands), 1986.

EUT Report 86-E-156

Address of the authors:

Group Electromagnetism and Circuit Theory, Department of Electrical Engineering,

Eindhoven University of ~echnology,

P.O. Box 513,

5600 MB EINDHOVEN,

CONTENTS

List of principal symbols Part 1: Theoretical analysis.

By A.M.J. van Amelsfort and Th. Scharten 1. Introduction

2. Statement of the problem

3. Analysis 3.1. Theory

3.2. Excitation by electrodes

3.3. Generalization

Part 2: Numerical algorithm, analysis and results concerning two intrameddulary placed

v 1 1 3 8 8 12 18

electrodes. By A.M.J. van Ameslfort 24

1. Introduction

2. Statement of the problem

3. Numerical procedures

3.1. Numerical algorithm 3.2. Numerical analysis 4. Numerical results

4.1. Tissue of infinite extent

4.2. Composition of tissues surrounded by lossy air

5. Discussion and conclusions 5.1. Discussion 5.2. Conclusions 6. References 24 26 33 33 39 48 50 57 64 64 68 70

LIST OF PRINCIPAL SYMBOLS

E electric field phasor

f frequency

H magnetic field phasor

I identity matrix

=

J electric current density

J Bessel function, order n

n

j

1-1

K magnetic current density

K modified Bessel function, order n

n

N number of tissue layers

Q electric charge

r radial coordinate

r position vector

t time coordinate

u unit vector

U unit step function

X,Y transverse cartesian coordinates

Y Neumann function, order n

n

z axial coordinate

S

axial phase-change vector, axial wavenumber

o

impulse function

~ complex tensor permittivity

£ permittivity of free space

o

K transverse wavenumber

~o permeability of free space

v normal coordinate

p electric charge density

a tensor conductivity

=

T tangential coordinate

w radian" frequency

a

_i partial derivation with respect to some coordinate i

V del operator

Subscripts

i layer number t boundary number

s source

S surface density

t transverse component

Superscripts

e belonging to the axial electric field

h belonging to the axial magnetic field

J belonging to the Bessel amplitude coefficients

Part I: Theoretical Analysis 1. Introduction

As it became known that time-dynamic stress imposed upon bone influences not only the (re)modeling of bone itself [39J, but also generates an electromagnetic field [40,4lJ, numerous experiments has been carried out confirming bone

Ire)modeling can be influenced by electromagnetic fields [2,3,6,8,15,19,21]. In most of these experiments these electromagnetic fields where being generated by electrodes or coils whichwere, on their part, being driven by static or time dependent currents or voltages. A number of experiments, methods and models have already been summarized in recent reviews [16,37J. In spite of all these

activities as yet, no satisfying explanation for the observed phenomena has been

found.

.

"

From an electrophysic point of' view we can distinguish three research areas: - determination of the biological electromagnetic sources [22,34]

determination of the macroscopic electromagnetic properties of the tissues

involved [10,11,23,24,35,42J

- determination of the macroscopic electromagnetic field distribution generated by artificial electromagnetic sources [13,17,27,28,29,30J. In this study we restrict ourselves to the third area only.

The difficulty in designing a model to determine the electromagnetic field distribution in a limb, lies in the combination of a non-uniform geometry, the existence of different layers and the anisotropic property of bone and muscle

t issue (see fig. 1.1).

W~III'I--- Fat/skin

tfll"Ht--- Bone

"l1IC-_- Muscle _ _ Air

2

-Another aspect is the iarg~ variety of electromagnetic sources, biological and

art ificial, for which the electromagnetic field distribut ion should be deter-mined.

The models, recently proposed, consider a very simplified configuration or give

only a part of the solution of the boundary value problem. Dierikx et.al. have calculated and measured th" elf'c! ric field distribution of two Helmholtz coils and two circular met.al plates placed in a conductor [29]. in the same article

he describes a possible change in the field distribution when a dielectric rod

is placed between the coils or plates. Lunt has calculated the electromagnetir

field distribution of two square coj Is placed in a conductor [27·!. In order to

investigate the influence of bone on the held distribution he also placed a

dielectric rod between the coils. At first he used a rod with a square cross

section, and later Oil a circular cross section [28J. An essential aspect in his

model is that the p<'rmittivity of both tissues must be identical. Mcloed has

calculated th .. field distribution around a bone fracture [:30 i and around a

culture dish lTL ~ using Helmholtz coils or conducting plates. In the papers

mentioned thus far, the proposed methods to determine the field distribution

at'" all based upon a static or quasi-static approach to the boundary value

problem. Hagmann has chosen a model, where the various tissues are modelled as

concentrlc circularly cylindrical layers inside an infinitely long concentric helix. Departing from a time-harmonic formulation of Maxwell's equations, he calculates the dispersion characteristics of this configuration. Although his method could be used to determine the Specific Absorbtion Rate (SAR) of the various tissues, we have some doubt about the reliability of his numerical results. The first reason to question his results is that he states to hav" an

overdetermined eigen value problem, which 1s not. Secondly, he considers only

one mode. EH om ' and he has not explained in his paper for which m [1,2,3 •... J

thE' results apply. Because the solution of the eigenvalue problem form a discrete set of spectral components, this approach cannot be used to find the complete solution of an active' houndary value problem. which will lead to a continuous spectrum. In all the above mentioned models it is supposed that t'he various tissues are isotropic. Moreover the methods cannot easily be generalized 1n or'der to determine the field distribution generated by other

stimulators.

With our study we want to make the next logical step in the development of a model to determine the electric field distribution in a limb. In our proposal the various tissues, the anisotropy of bone and muscle will be taken into

consideration. Furthermore, the results of the derivation can be used for diff~­

rent stimulators.

Using a layer'ed circularly cylindrical geometry and assuming the electrical propet'ties to be anisotropic with respect to the cyl inder axis, we will show how the complete solution to :'-1a'{well· s phasor equations can be found for an arbi tral'y electromagnet ic source distribution (Part I). Wi th the proposed method, based upon the theory of guided waves, the electric field distribution

will then be numerical determined for the case of two current driven

intramed-du I ar'y placed cir{'ular electrodes (Part II). We have chosen thi.s particular

configuration in order to compare the numerical results with the experimental

resuHs obtained from the Utl'echt \;niv"rsity Hospital [19J. However, most of

the attention will be focussed upon the influence of source frequency and a.nisotropy on the field distribution.

2. Statement of the problem

Fi~. 2.1 shows the cross-section and longitudinal section of a circularly

cYlindrical limb surrounded by air. The various tissue layers are assumed to bl~

coaxially arranged, where the boundary surface between two adjacent layers

marks an abrupt change in the electromagnetic properties of the media.

Further, we suppose the tissues to be independent of time and signal strength,

reciprocal, homogeneous, disspative, and in general, anisotropic with respect

to the axial direction. For a more correct description of muscle and bone tissue the last mentioned property is an essential one [10,35,42]. The

permeability of each layer e.quals that of free space, since a tissue exhibits no substantial magnetic effects.'

So, in cribed • - - - - a i r

--:::::g~~=

faUskin f' muscle -:~I=: bone

h'%'b%>",*,',,*sS§§%'§

I- marrow - - z ax is (a) (b)

Fig. :!.l Schematic rlrawing of the mod,,1 of a limb (a) cross-, (bl longitudinal section. the frequency domain,

by the permeability "

the electromagneti.c behaviour of a and a complex tensor permittivity

( ' t fi - : o ' t o

~)

_z

.

(2.1 ) tissue is des~· § given by where e:

t and e:oz; denote the transverse and axial permittivity, respecti\'ely. In

view of the above suppositions, ' t and " are layer-dependent complex constants.

i ,. ~. - ji"

(2.2)

we notice that the imaginary parts are representing both electric conduction and dielectric hysteresis of a tissue. For, as is well-known,

(2.3 )

where 2 is the tensor conductivity and l; is the imaginary part of electric tensor susceptibility, representing dielectric hysteresis. The real part of permittivity is related to the real part of susceptibility, ,~ , as

(2.4)

The two before mentioned loss mechanisms cause an exponential decay of the electromagnetic field in axial directions off the source. It is therefore that the physical length of the limb does not enter the problem as long as the decay is sufficiently steep.

The stimulator is represented by a bounded region inside or outside the limb, where electric and magnetic current densities are prescribed functions of posi-tion and time. If coils are used as a stimulator

[2,21),

then the electric current density ~ is to be specified accordingly and the magnetic current densi ty ~ is zero everywhere. In that case the stimulating electric field is generated by the time-varying magnetic field, emanuating from the electric current. On the other hand, employing electrodes [3,15,19), the electric field is excited directly by a distribution of static or time-varying electric charge. In the latter case the charge density p is related to the electric current by

the equation of continuity, which in the frequency domain reads as follows:

(2.5 )

This indicates that the currents which are supporting the time-varying charges must be taken into account, unless the configuration allows discarding these

currents.

Apart from this, it is possible in some cases to replace a given distribution of static or time-varying charge by a distribution of magnetic current .. As an

example, let us consider two circular disk electrodes with opposite charge, located closely to one another (Fig.

2.2).

z-axis

-

"

Fig. 2.2 Stimulator consisting of two circularly disk electrodes with opposite charge.

The local magnetic current density inside the cylinder <. [0.<.1 , ' • [-L.L] :'is related to the electric field between the two disks according to ~ -, § '~<'

Supposing this field to be axially directed we have

< • [0.<,1. , • [-L.LI.

(2.6)

elsewhere,

where Ps is the local surface charge density.

Finally we notice that surface distributions of:! and ~ form a useful means

in representing irradiated sources, such as electromagnetic horn antennas.

We now may formulate the problem as to determine the field (§.~) from any source distribution (:!.~) located inside or outside the limb structure as specified in Fig. 2.1. In other words, the problem is to find the solution (§.~) to Maxwell' s phasor equations for a given distribution of :! and ~

v x ~ + jwu~ - - ~

v x tl - jw~'~ • d (2.7)

satisfying the conditions at the boundary surface between two adjacent tissue layers, namely the continuous transition of the tangential components of § and ~

through that surface.

In eqns. 2.7, :! and ~ are the spectral functions of the real, time-dependent current densities

I

and

g ;

so

"

-,

f Jexp(-jl.llt) dt (2.Ba)

-~

-,

f !eltp(-jwt) dt (2.8b)

'"

,~ ,;!

i

/,.

from the inverse Fourier transforms

~ os:

k

L

§ exp(j",t) dill (2.9a)

E

a:

4w

I

~ exp(jwt) dw . (2.9b)

If sources with simple harmonic time-depence (w-w

o) are concerned, we have

instead of eqn. 2.8:

~:- Re(::!. exp(jwot)] (2.10a)

As a consequence the time-harmonic electromagnetic field (g'.!!) is found from

(2.11a)

(2.llb)

The same arguments yield for the electric conduction current density ~.' given by

J - O'~E

-& - - (2.12 )

Finally, something has to be said about the presence of electric surface charges, being of interest from biophysical and bio-electromechanical points of view. A discontinuity of the normal electric field components across a tissue boundary is supported by electric surface charges in that boundary. We name these charges secundary as distinct from the primary ones which, may the case arise, are

representing the stimulating source.

Disgarding dielectric hysteresis (being of minor importance comparedwita electric conduction), the secundary surface charge density due to conduction currents in the tissues is given by (Fig. 2.3)

P -: - __ 1__ 0 E -0 E

Sc j w ( t2 v2 tl vl)

(2.13 )

Besides, we have a secondary surface charge density due to electric polariza-tion, given by

PSp·: - EO( X~t2EY2- X~tlEV1)

,. (£~l-£O)<l- (£~2-to)EY2

(2.14 )

The limb structure, as shawn in Fig. 2.1, is taken as a composite, circularly cylindrical open boundary waveguide of infinite length. For our purpose the theory of uniform waveguides has to be extended so as to comprise anisotropic matter, as specified by eqns. 2.1 and 2.2. This will be done in Sec. 3.1, paying special attention to boundary equations. The theory will then be applied to the limb, where an axial pair of ring-shaped electrodes is acting as a source (Sec. 3.2). In Sec. 3.3 the procedure for solving the boundary value problem will be generalized in order to include other source configurations.

3.1 Theory

In connection with the uniformity of the structure, the boundary value problem can be reduced to a problem in terms of the coaxial electric and axial magnetic field components. First, all phasor fields and the del operator will be sepa--rated into their transverse and axial components:

F :- F + F u - -c z-z 'V :- V + u a t -z z' where ~ is ~,~,:! ~ and f_t " (Fx,Fy'O) v_t " (ax'CIy'O). (3.1 ) (3.2 )

Secondly, we carry out a Fourier transformation with respect to the axial coordinate according to the transformed pair

f -: f f exp(-jBZ) dz (3.3a)

1

-~.. z. _~ E exp(jBz) dB . (3.3b)

Bearing in mind the exponential tions off the source, eqns. 3.1

decay of the field components in axial direc-then give

(3.4.)

V·F. - (V + jBu )o(r +

F

u )

t -z -t z-z (3.4b)

Field Equations

Using prescriptions 3.4, the differential equations 2.7 become

-

- -VtX ~t+ _{jwuoKz2z - - Kz2 z} v x

H -

jwe'E u - J u t -t ~ z-z z-z ja~,' ~t- ~zx V E ~ jwu H - -K t Z o - t - t jau _-z x _{- t}

H -

_-z u x v

H -

jw£ E - J t Z t - t - t (3.58) (3.5b) (3.Se) (3.5d)

In order to derive differential equations for the axial field components, §t and tit

are successively eliminated from eqns. 3.5, using the vector identities

VtX <2_zx ~) • (V

t" E)2z

v_t" <~zx ~) - -(VtX E)"2_z

The differential equations are then found to be

where the transverse wave numbers ICe and < h are defined by

(3.6)

(3.7 a)

(3.7b)

(3.8b)

If the stimulator is represented by electric current and electric charge only, then ~=~ in the above equations. Using eqn. 2.5, we have

v .J _{t -t} + JaJ - -_z jwp

with which eqns. 3.7 change into

(3.lOa)

(3.10b)

where p is the Fourier transformed volume charge density. For given sources,

the axial field components are to be determined from the above equations. When these components have been found, the transverse fields are obtained from eqns. 3.5 c/d as

(3.Ua)

(3.1.lb)

Boundary equations

In order to complete the theory, the boundary conditions have to be expressed in terms of the axial field components. These conditions read (see Fig. 3.1)

u x E - u x E

-\I -1 -\I -2

u )( H. .. u x H -v -1 -v -2

Fig. 3.1 Boundary surface

(.1.12a)

(3.12b)

where ~v is the unit vector along the normal on the uniform boundary C. the subscripts 1 and 2 referring to the regions on either side of C.

There are no surface currents in the boundary because the permeabilities and permittivities of the tissues involved are finite-valued.

Applying the Fourier transformation 3.3a, separating into transverse and axial phasors and subsequently substituting expressions 3.11. boundary conditions 3.12 become

(3.13a)

(3.13b)

! t C 3.13c)

! Eo C (3.13d)

where an orthonormal set (~v'~"~z) on C has been adopted, as shown in Fig.' 3.1, d_r

ffild

d,

denoting partial derivatives with respect to the transverse tangential and the normal coordinate off the surface, respectively. In virtue of Ez and ii_z

being continuous through C, the following equalities also hold:

a E: • a _{. z l 1z2 E : - a .E:z• .} ! f C (3.14a)

a _{• zl} K - a _{• 2:2} K :- a _{• z'} K ! f C . (3.14b)

Therefore. the buwldary equations 3.13 can be rewritt",n to give

From eqns. 3.15 it is seen that, irrespective of the source, the axial field components E_z and ii_z are coupled at a boundary if the field is depending on the , .. coordinate along C. If not, these components are not being coupled. As a

consequence, solutions of the homogeneous differential equations {3.,

1< ,

2; ~ , 2) are representing hybrid wave modes if the field components depend on the , coordinate, whereas the solutions represent pure transverse electric and transverse magnet ic wave modes, i f the fields are independent of, .

Source-boundary equati 0115

I f charge and currents are confined to an infinitely thin, axially directed uniform sheet, the source may give rise to boundary equations as well. We assume that the source sheet is not coinciding with a tissue boundary. The boundary equations may then be derived directly from the field equations 3.7/10, namely

by integration from v - v- to v - v+, yielding

-+ -- -E - E .. K z Z ST -+- - -Ii - H .. -J z Z 5"[

at E ... E_{s '} and, if a charge distribution is specified:

r - r - -s (J .lba) (3.16b) (3.16c) (3.16d) (3.l6e)

The index S points to a surface density of the quantity involved.

3.2 Excitation by electrodes

The theory, developed in the previous paragraph, will now be applied to the limb structure, where an axial pair of identical ring-shaped electrodes is acting as the stimulating source. This configuration, as shown in Fig. 3.2, is fairly a good model for theoretical verifying some experiments carried out at the Utrecht Uni versi ty Hospital [19].

_ _ _ _ _ air

~

~~S~==:fat/slcin

muscle

bone _{W$ , !!:>sss'} ~> <g.$' S

"-. ~ - .. source .") ~,

+A--H _ _ marrow - - - .! .. .!!~~

Fig. 3.2 Cylindrical limb, excited by ring-shaped electrodes.

FOI7J1uiation of the boundary value probleJJI

We wish to determine the time-harmonic electromagnetic field. at. extreme low source frequencies as an approximation for the static case (freq~encies~ 1 HzI

To allow for comparison also higher frequencies, say 1

MHz,

are of interest. The electrodes are at any moment oppositely charged. In view of the relatively low frequencies, the charge will be distributed·homogeneously along the ring. The current wires, supplying the charges, are left aside: the· electric field

between the electrodes (the region where the fracture is to be lor.atedl will be predominantly determined by the charged electrodes. If, in the future, it will be found necessary to include the field caused by the current wires, these contributions can easily be added to the field without the current wires. Consequently, the source is represented by the spectral functions

P - ~ _Z'II'r 6(r-r ) [6(z-L) - '(Z+L)l. _s

,

r, denoting the radius of the 1lectrode.

Thus, after Fourier transformation,

implying a surface density Ps given by

p - -_S ~ sin (~L).

>r

,

On substituting, eqns. 3.10 give

(3.17)

(3.18a)

(3.18b)

(3.19a)

(3.19b)

This means the electrodes are not generating an axial magnetic field. Because of symmetry, the fields are independent of the • coordinate (Fig. 3.2). As has been pointed out in Sec. 3.1, this_ independence implies that ii_z is neither generated at the tissue boundaries. So Hz will be zero everywhere. Hence, using circular

cylindrical coordinates, the boundary value problem is written as

[ d 2 _r + r -1 d _r + (K,) e 2 -

1

E _z, " 0 (3.20a) i o(r) = i 0+1(') Z,l 8 Z,l 8 (3.20b) d E , ( c ) .. d E '+l(C) - _Q_- a siu(aL) c z, s c Z, S WE: r t , l s (3.Z0c) (3.20d) (3.20e)

where use has been made of eqns. 3.19a, 3.16b, 3.16e, 3.13a and 3.14b,

succes-sively. The numbering (il of the various layers and boundaries is illustrated in Fig. 3.3. The layer containing the source has been split into two regions (of the same electrical properties) and the numbering goes from the z-axis outwards. Let N denote the number of layers, then the infinite region surrounding the limb may be considered as the (N+I)th layer.

4

3 3

2 2

1 1

boundary layer

Fig. 3.3 Cross-section. Numbering of layers and boundaries.

Derivation of the solution for E,

The solution of differential equation 3.20a is a linear combination of Bessel and Neumann functions of order 0:

where

r • 0

o r . . l

-where the C's and D's are complex amplitude functions of a and., to be determined from the boundary equations 3.20 ble, r being the radius of observation, We

name the complex amplitude coefficients the Bessel amplitude coefficients. Since the solutions must be bounded and, for the exterior region (i=N+Il, must

represent outward traveling waves, we have the following additional constraints:

().2ib)

().2ic)

Since the determination of the Bessel amplitude coefficients in case of more than some two or three layers is lengthy without resulting in' sureyable' formulas, it might be convenient to call in a computer making use of a

and Reduce are examples of such a language [9,18]. Depending on the complexity of the procedures it may also be possible to evaluate the Bessel amplitude

coefficients numerically. For both purposes the boundary equations, substituting solutions 3.21, are written in matrix form. We then obtain

where

e e T li,' (C"D , )

,,[l,N],".

In these equations use has been made of

and (3.22) (3.23) (3.24) (3.25a) (3.25b) (3.26a) (3.26b) (3.27) (3.28 )

We name ~,the Bessel amplitude vector. The set of equations 3.22 - 3.23 is solved in two steps. First, using eqn. 3.22 the Bessel amplitude vector for the k-th layer is expressed in terms of either the Bessel amplitude vector for the regionQ,,)or that for the exterior region (~N+l)' namely by repeated eliminations. It then follows that

i f • C [2,.] with. l [2.NJ

(3.29) i f k. [O+l,N] with. £ [l,'HJ

where the matrices ~k and ~k are defined as and 1 -1 p~: n (B • A ) -k l~k-l -1 ~i' ~l" ! . k f [2,.] k - 1 k E [o+"N] B - N • k - Ni-l. (3.30a) (3.30b)

Secondly, substituting expressions 3.29 into eqn. 3.23 yields

(3.31)

from which the vectors '1 and 'N+l

of the constraints 3.21, we have

can be fOWld in the following way. Making use

(3.32 ) e T £N+l- CN+l·(l.-j) Next we introduce p

:_(Pll

PL2 ) _{i f • E ["N]} -s P" Pn (3.33a) P .. I s-l

-

- '

t'"

'u)

i f • E [l,N-l] ~O+l'-'21 '22 (3.33b) gN+l- ~ •

.-ft.'

-where the elements Pi) and ' i j are known by comparison with definitions 3.30. By

means of eqns. 3.32 and 3.3~, eq. 3.31 is reduced to

{c.,. ( ) T .

2 • Pll,P12 - cN+l

or, written in a more appropriate form:

2.

(p

11

P"

.

(c~

)_

S . • -p c .... , (3.34)

Thus, the solution reads

(3.35)

where the transfer matrix

.p

is defined by

T -_'p

12 •

(3.36)

At several stages of the derivation, matrices have to be inverted. In this connections it may be noticed that

2t;t,l+1 It l dec R .. - - - - (K ) 51 'II r 1 1 2 det 0 ..

-,

'.

(3.378) (3.J7b) (3.37c)

In the derivation use has been made of the Wronskian of Bessel functions [1)

2

- - -

..

. (3.38 )

FieldS in spatial domain

Summarizing the results, obtained so far, the transformed axial electric field

~component in the k-th t issue layer reads,

18

-(3.39a)

the Bessel amplitude vector ~. being known from eqns, 3,29 and 3,35. From eqns. 3.10 the transformed transverse field components are found as

and _ WEtk._ H - - _ _ E ,k B rk (3.39b) (3.39c)

Using the inversion formula (3.3b), the spatial distribution of the field components is obtained from

l "J

-E - - E cos(sz) dB zit 'II 0 zit

E ." _r .1 _{11 0}

j

e

_r .sln(Sz) dB l • -H - _

,It

f H ,.cos(sz) dS. 11' 0 , .. (3.40a) (3.40b) (3.40c)

since ii_z• and iii. are even functions of a and Erk is an odd function of a. These

integrals are to be evaluated numerically. A numerical algorithm will be given in Part II, as well as some numerical results.

3.3 Generalization

In view of further investigations, we will now generalize the procedure, out-lined in the previous paragraph, so as to include sources which are more complex than the ring-shaped electrodes. We restrict ourselves to surface distributions of electric and magnetic currents over a circular cylindrical surface of finite axial extent. In general the surface currents will be -dependent. As was shown in Sec. 3.1 this dependency results in hybrid field solutions.

As a first step the • -dependency is removed by applying a finite Fourier trans-form according to

,

F -: f F exp(-jn,)d, -. i' • ...l _2,

L

...

-(3 ... 41a) (3.41b)

19

-where F denotes any component of ~, ~

Electric current sources Let

,

E_z·: J Ezexp(-jn+)d+ -,

,

K1-: J H1exp(-jn+)d+ -'II • -;.

then .. and

a.

are solutions of the transformed equations

respectively, yielding

where, just as in the previous case,

D~- 0 De _{N+l- -}J Ce _N+l' h h DN+l--j CN+l' (3.42a) (3.42b) (3.43a) (3.43b) (3.44a) (3.44b) (3.45.) (3.45b) (3.45c) (3.45')

It will be understood that the Bessel amplitude coefficients

depend on the value of n. Next the boundary equations

3.13 alb, 3.15 alb and 3.16 are transformed making use of the transformations

3.43. Again the results are written in matrix form, where the following

assignments for~, ,2, and 2 hold:

"

J

We then have where A" "I a," a₂ -a₃" B • -I b •

,

b," lot ,.w 2nd row 3rd row 4th row C " 0 -s "

•

,

;;

_,

d

, ,

E

d _{r ,} ii

•

J_a("1(1) 0 alJ~(K~ri)

•

-a₂J n(1(1 r 1) h 2 • ("1+1) Et,1"1

-

20

-lot column C· 2nd column 0" 3rd column Ce 4th column De l£['.N].l •• 1 " •

•

Yo("1 r 1' 0 0 In(Kir

•

_l) alY~(K~ri) h &2Jo(K_{1 (1)}

•

-a_{21 n(1(1(1)} a3J~(IC~ri) .

-,

Jar aW~o(€t.i+l-Et.l) h 2 h ).10("1+1) lei

•

J o(I(1+1 r1) YU(K!+l

•

r 1) 0 0 0 _{I n(''1+1}h r 1) blJ~(lCi+l (i) bIY'(lCi+l r_{1 )} 0

0 0 b2J~(IC~+l _{r 1 )} h 2 • (rei) Et.i+11<:1+1

•

h 11 0(1(1) "1+1 (3.46a) (3.46b) (3.47) 0

•

Y n(K1(1) (3.48a) h &2Yn("1 r 1) a3Y~(K~ri) (3.48b) (3.48c) (3.48d) 0 h Yn("1+1 (1) (3.49a) 0 b2Y~(IC~+lrl) (3.49b) (3.49c)

21

-J (Ker )

•

0 0 n " Yn(KSrS) J (Khr ) h 0

"

n s s Yn("sr s ) 2 • K~' (Ker ) (3.50) KeJ'(Ker ) 0 0 s n s s 8 n S S 0 0 "hJ • (Khr ) s n s s "hy'(lChr s n s s

and where the source vector for electric current sources is given by

5 --J _ 2 . . Sn T ( 0 _{• s+'} J -J(w. - _a_ ) J ... _ _ _ 1 ,0) 0 WE Sz WE r Sf t,s t,s S (3.51)

js. and jsz being the transformed azimuthal and axial electric surface current densities.

In case of magnetic current sources all formulas are the same as in case of electric current sources, with exception of the source vector ~J which has to be replaced by

(3.52)

where Ks. and Ksz are the transformed azimuthal and axial magnetic surface current densities.

Solutions of equations 3.46

Referring to the derivation described in Sec. 3.2, formulas 3.29 - 3.30, the Bessel amplitude vector for the k-th layer is written as

c --k where i f Of [2 •• ] i f • f [ .... , •• ] (3.53)

,

,

I

/

(3.5.)

~,and ~, being given by eqns. 3.48 and 3.49 respectively. Introducing

P11 P12 P13

")

1:.s:- Pn P22 P23 P2. Pll P32 P33 P3• P4l P42 P43 P44 (3.55) and '11 '12 '13 '14 '21 '22 '23 '2. 9.+1 ,- _{q31 '32 '33 '3.} (3.56 )

'.,

'.2 '.3 '44

the Bessel amplitude vectors for the central and exterior region are found as

(3.51)

where the transfer matrix .Jis defined by

P11 P13 -qll+ jq 12 -q13+ jq 14 -1 P21 P23 -q21 +jq22 -qn + jq₂₄ IJ -: R Pll P33 -q41+ jq 42 -q43+ jq44 (3.58)

p.,

P43 -q41+jq42 -q4J+ jq44

With view to inversion of matrices, the following determinants are easily derived from 3.47, 3.48 and 3.50:

dec n .. _ ( __

.

..

2 __ )2

•

(3.59a)

(3.59b)

23

-Inverse transformation

Summarizing the results obtained so far, he transformed axial field components in the k-th layer reads

E 0

','

(Ck,Dk)o e e T (J n (Kkr). Y n ("kc» (3.60a) ii h h h h T (3.60b) z,k 0 (c. ,D.), (In(~r)'Yn(Kkr)) I

the Bessel amplitude coefficients being known from eqn, 3.53. From eqn. 3.10 c/d the transformed transverse field components are found to be

- h -z[ e[ e e • e , e T] Er,k- (lC k) -j~lCk (Ck,D_k) • (In(lCk,c),Yn(lCkr)) + ,1 -1 h h It h T + wlJoor [(Ck,D k) ' (J~(lCkr).Y~(lCkr})]} (3.61a) (3.Gle) ().6ld)

Finally, using the inversion formula 3.41b and 3.3b, the complex field in the k-th layer reads

-E -_-k ....!

2lf

1 [

2; _....-

1

!k

up(jn,)] exp(jSz) dB (3.62a)

-

-H _ 1

-k, "2'i

1

h~

1

tlk

exp(jn,> J exp(JSz) de (3.62b)

n

-I

/ .I

intrameddulary placed electrodes.

1. Introduction

The numerical evaluation of the electric field distribution involves the compu-tation of the inverse Fourier transform given by

§ •

k

J § exp(jBz)dS. (1.1)

In the development of numerical procedures, the integral can be evaluated separately from the computation of the integrand itself. Once the integrand is known, the integration can be carried out.

The integrand is the complete solution of the active boundary value problem given by eqns. 3.45 - 3.58 Part I. However, straightforward numerical imple-mentation of this procedure gives no accurate results. Especially, if the model

consists of three or more tissue layers, numerical evaluation will give rise to

combinations of Bessel and Neumann functions {In''n}' which are equal to, or

numerical approach to Wronskians

(1.2)

or Hankel functions of the second kind

(1.3)

where

~ 2 2 €'zi , - r₁ "£ 1" 6

-Z 0 EeL (1.4)

If standard computer library routines are being used to compute _{I n(,)} and _Yn(,) the loss of significant digits in the determination of a combination can be as large as the total number of significant digits, used by the computer for

internal representation of reals. So, an alternative numerical algorithm should be developed with which reliable results can be ontained. Essential for the

integrand ~ is the computation of the Bessel amplitude coefficients. These coefficients depend on the number of tissue layers (N) and the location of the source.(s). Because, apriori it is not known which tissues are dominant with respect to the electric field distribution, it should be possible to use numerical procedures, concerning the Bessel amplitude coefficients, for an arbitrary number of tissue layers and source location. It is then possible, on the one hand, to investigate the influence of each tissue on the field

distribution separately. On the other hand, during the development of the program itself, numerical instabilities can be recognized and abolished. Therefore, the algorithm for the coefficients must be a compromise between direct numerical implementation of the original procedure and, following the

same procedure, a straightforward analytical derivation of the coeffecients in terms of functions and parameters for any combination of N and s.

The first method leads to inaccurate results. The second method can only be used practically together with computer programs suitable for processing mathematical expression [9,18J. However, these programs are not designed to find numerically stable algorithms for numerical problems. We therefore must develop a new

algorithm that should include a method to detect and replace those combinations of In('land Yn(,l, introducing a loss of significant digits, by more accurate alternatives.

At first sight, the use of excisting F.F.T. algorithms for the calculation of the inverse transform seems to be favourable. In this case however, this is not a very efficient method. A correct use of F.F.T. demands the number of samples to be such that. the samp 1 ed- ibtegrand is a good representation of the continuous integrand. So, if the integrand is not a smooth function of the axial wavenumberN,. the number of samples must be very large, thus increasing the required

processing time. The processing time, predominantly determined by the computation of the integrand, can be limited by deviding the integral into separate parts, using a minimum number of samples for each part.

In this part of the report, I shall propose an algorithm to compute the electric field distribution generated by circularly symmetric charge distributions, which answers the above consideration (Sec. 3). The results obtained with this algo-rithm will be shown in Sec. 4 and concern the field distribution generated by two intremeddulary placed electrodes. These include not only the electric field distribution itself, but also the secondary current density and secondary sur-face charge density. The results will be given in the transformed domain (ampli-tude spectrum, 3D plot) as well as in the space-time domain (distribution func-tions-, 3D plot). In order to understand more about the quality of the field distribution, the results will also be represented in a directional field pattern.

In Sec. 5 the results will be discussed and conclusions will be formulated.

"

Bessel amplitude functions

The set of equations (Part I, eqns. 3.21-3.28) from which the Bessel amplitude coefficients are to be determined, can also be written as

where we define

b~).

i

₁

so, from comparison we learn that

where J __ -_''''' ',;:-,-ill'" h 2 2£ti(K_i+1 ) 1"'1 (2.1a) !E[2.S] (2.1b) 1-s (Z.ic) l'[S+l.N] (2.1d) (Lie) (2.1£ ) (2.tg) (2.ih) (2.11) (2.1j) (2.11<) (2.1m) (2.ln) (2.tp) (2.1q)

In order to explain where the numerically unstable combinations stem from, an example will be examined in more detail. Let us consider a configuration of two tissues, the source being located in the central tissue (N=2, s=l). Straight-forward analysis of eqn. 2.1 gives

(2.2a)

where

(2.2b)

The attention will now be focussed upon an auxiliary function, which appears simultaneously in various forms and at different places in the expressions of the Bessel amplitude coefficients. This auxiliary function is given by

(2.3 )

which can also be written as

(2.4)

For large arguments <" i t follows that

1<,1»1

(2.5)

Substituting in eqn. 2.4 gives

If the absolute value of <, is large, numerical evaluation of the intermediate

result~(a)according to eqn. 2.3 introduces a loss of significant digits. However, using eqn. 2.6 more accurate results are found.

I

auxiliary function can be recognized. Even worse, increasing the axial wave-·

number (a) , some combinations numerically approach an auxiliary function, thus

increasing the loss of significant digits. Due to the variety in the values of the tissue parameters and boundary radii, the range of', can be very large for each a. Therefore, simultaneous replacement of the auxiliary functions by their counterparts will in general not lead to more reliable results. The substitution must only be made when the absolute value of the argument involved is large enough.

Hankel amplitude coefficients

If, for large absolute values of', ' the combinations of I_n and Y_n are nearly

equal to Hankel functions, an obvious alternative, is to write the solution of the Bessel differential equation (Part I, eqn. 3.20a) as a linear combination of Hankel functions of the first and second kind "~'I,),"~21,) according to

Using similar set of equations (eqn. 2.1), the Hankel amplitude coefficients are defined by H. H T H T (C"D,) • (C"O) , where we define with

b~),

d"

,

1"1 (2.8&) if[2,.] (2.8b) 1"s (2.Se) '<['+',NJ (2.8d) (lo8e) (2.Sf) (2.8g) (L8h) (2.8') (2.8l)

where (2.8k) (2.8m) (2.8.) and (2.8.) (2.8q)

Equation 2.8a meets the constraint of a bounded solution in the interior region. Equation 2.8e describes an outward travelling wave in the exterior region. Using only the last set of equations to compute the field solution, again numerical instabilities can be expected. To clarify the cause of these inaccuracies we take the same example as before N=2, s=l. It follows that

(2.9a)

where

(2.9b)

Now the auxiliary function can be recognized to be

(2.11 )

Both auxiliary functions bJ(6) and bH(6) are related by

(2.12 )

So, it can be expected that for small absolute values of " expression 2.10 will introduce a loss of significant digits and 2.11 will be more accurate.

Summarizing, there are several methods to our disposal to obtain a numerical result. The accuracy of these methods depends on the numerical values of the arguments " involved. For instance, if these arguments are small, say

(2.[3 )

the best approach is to use the Bessel amplitude coefficients. If, however, these arguments are large, say

(2.14 )

it is better to use the Hankel amplitude coefficients. Clearly, there is a methodical gap between small and large arguments. If an argument lies in the

interval between small and large values, additional measures have to be taken to replace the numerical unstable combinations by their more accurate counterparts. Both, Bessel and Hankel amplitude coefficients originate from the boundary conditions. The occurring numerical unstable combinations stem also from the boundaries. This leads to the suggestion that the algorithm to be developed for the determination of the coefficients has to be constructed starting from the boundary, or boundaries, causing the numerical instabilities.

Inverse Fourier Transfo171l

As soon as the transformed field solution is available, the inverse transform is computed according to

(2.15 )

Separating!l in its even and odd functions (!le'~o) of the axial wavenumber,

(2.16) gives 8 L aL E - 11m { .!. f E c08(6z)d8) + 11m { 1 f il'oSlQ(8z)d8) - lI' -e 'I\' 8L- 0 _8L- 0 (2.17)

Because the computation of the integrand (or better the coefficients) consumes most of the required processor time, it was suggested to divide the integration interval into separate parts. Defining

and yields where B"" I -! f E coa(az)d6 -em.... -e B. B ... , I - 1 f g s1n(Sz)dS -om 11 -0 B. H-l IH § - r fem + E 1_0m + 1Il-0 mao B - 0 _o (2.18a) (2.18b) (2.19) (2.20)

It is now possible to determine the first M definite integrals with standard numerical algorithms. The various values '. and corresponding number of the samples (N.) can be chosen independently and are only determined by the

behaviour of the function !(S)

For the determination of the tails of the integrals formed by

BL . I - lim { ~

f

E

co.(B.)dS} -te B +- 11' 8 -e L H BL I - lim { 1 f

E

'in(S.)dS} -to '11'-0 6L- 8_t (2.21a) (2.21b)

another method is being used. We suppose that we can choose SM so large that the numerical value of these two integrals is small compared with the sum of the first M parts. I f we then neglect the contributions of these tails to the final result, this will be deteriorated by a sine function.

Therefore, as a compromise we determine an estimation of these integrals.

Numerical analysis of! shows that for large axial wavenumbers !. and!o can be approximated by

E • E •

A exp(G S) -e -ea -e -e (2.22a) and (2.22b) i.

i

/

respectively, where ~e J ~e ,Ao and Q_o are functions of B for which a A .. 0 e > a-e - eM (2.23a) a G .. 0 B-e - e > eM (2.23b) a A .. 0 6-0 - e > eM (2.23c) a G .. 0 e- - 8 > 8M (2.23d)

Substitution of these functions into eqns. 2.21 gives rise to an explicit formulation of the infinite integrals in terms of Ae,Ge,Ao,Go and eM • What is

left is how the auxiliary variables Ae,Ge,Ao,Go can be determined and how the

3. Numerical procedures

In the first paragraph of this section we will clearify the algorithm for the determination of the amplitude coefficients and the one for the evaluation of the inverse Fourier transform. The construction of the algorithm will not be described in detail. Instead, the method leading to the algorithm will be

sketched roughly and is illustrated by an example. The algorithm itself will be' given as a set of equations. In the second paragraph the numerical analysis of the most important functions and procedures will be described.

3.1. Numerical algorithm

Amp] i tude coefficients "

If the active boundary is not taken into consideration, the amplitude coeffi-cients of the i-th tissue layer can, in general, be expressed in terms of the coefficients of the (i+2)-th tissue layer by

where use has been made of eqn. 2.1 bid. Numerical evaluation of expression 3.1 involves a sum of products. However, it can also be written as a product of sums following

C1+2

c x - - + d

1+1 01+2 H1

(J

Defining a recurrence relation by

and substituting into eqn. 3.2 gives

where C1+₂ R -1+2 01+2 ci+1Ri+2+ di+1 )D i+2) c₁+_1R1+₂+ _d1+_1)Di+₂ (3.2 ) (3.3) (3.4) (3.5) , " I , , . !

The advantage of reorganizing the original equation 3.1 is that the intermediate results are now formed by expressions and functions related to the same boundary i or i+l. The influence of the other boundaries is reduced to only one complex number "HI. In order to illustrate the reformulation of eqn. 2.1, the expression

involvIng the smplitude coefficient

ot

is given for the case the model consists of five tissue layers, the source beIng located in the exterior tissue

(N=5,s=5). Equation 3.6 and 3.7 show

ot

as a sum of products and a product of sums, respectively.

--\ d IC2- )_ .. Y O(O:~fs)-d IC2A) ~ JO( O:;fS) .. d Id 2c l-"Y:<O:SfS)-d ,d:zC lb .. JO( .. ;f S)+d Id2d,c .. Y O(o:;rS)-d ,d2dld .. JO(O::fS~ lIus_

Jdlc~J." ... dtC2.,cI.. + C!dld 2c,... ... d 1d 2cJb.. + Jd,d 2d,c" ... dldzd,d.. } • ~ Ps (l"e ) ,,( .. ;rs) -', _ _ . _ + b, JO(O:SfS) " YO<O:;fS) +b, -0, _ _ . _ + d, Jo(O:sts) " 'o(o:;rs) + b, - ' , _ _ . _ + b, Jo(O:sts) 0, 'O(O:;fS) + d,

ut-

_{j ... b ..} . , jc .. + d~ + b,

.,

_{j ... + b ..} + b, Cl~+d, -0, _ _ . _ + d, 0 _ _ ' _ 0 ( ° JO(O:SfS) + d l ) ja .. + b .. jc .. + d.. ' _Yu(o:~ts) c,~+d3 -" _ _ ,_+ b, Jo(o:srs) " Yo{o:;rs) + b,

0,

_{j ... + b ..} + d, -, jc .. + d .. + b,

0,

_{J ... + b ..} + d, Cl~+dJ ~,--,-

...

J O(O:5 t S) 0, 'o(o:;rs) + d, - ' , _ _ . _ + b. Jo{ests) 0, 'o(e;rs) + d, ~. _ _ . _ + d, Ja(ests) (3.7) 0

Substitution of eqn. 3.3 into 3.7 gives ( q l ) + dz} y,(o:~rs)

.,----Jo(o:;rs) Yo(o:~rs)

.,----Jo(4 c s) yo(o:~rs)

',-

--.-JO(0:5 CS) y,(0:5 C5) 1.5 . -J,(0:5 C5' (J.a)

The same can be done for the Hankel amplitude coefficients. Using this method, the set of equations defining the Bessel and Hankel amplitude coefficients can be written as 8-1 IT (

.-1

.-1 IT (

'-1

"

C~Rk+;1

+

d~)-l]

RN+l-j + d~)-l) R -R S 8+l RN+l-j x [ R +11 x J,«e r ) + y,«e r ») , _{8 S S 5 8} '1<+I'j J i-2 J + d~») di-d IT ( ckRk+l k-I J b1-l R_{i -l J} -d • if[2 •• ) i-I J i-2 J +dJ ») , i-I[ IT ( _{ckR k+} t

,-1

k J 8 i_1 K₁__l" J

-,

i-I

(

:~)

-8-1 J + d~l) 1-1 J +d~)-l) [ IT ( _ckRk+l

,

[ IT (c_k".+! k-1 k-I i YO(Ker ) kh R 'R s+1 R . - • 8

,

•

_{JO(K:rs )} RN+l_j -S1I'r s _

(~il~I_J

) • 1<[s+1 ... 1) x-zEPs

_,.

, (3.9a) (3.9b) , (3.ge) I! ! ,I ,,'

. i

and

(:;)_(R1

_ 1 ) -lx R_s·RS+1 R_N+_{1 ..(}} (3.10a) R +11 •• _s ~11. _{s s} , ) + "~21. _{s s} ,

)J •

"+1=0 • x( R1

•

·i-l [ H c_i_₁[ .-1 [ n (

'-1

i-2

•

n ( ckRk+l '-1 i-2 _H n ( ckRk+l '=1 i-I _H + d~)j. • b i_1 Ri _1 "" - . -_d i_1 +d~)] x

•

a i_1 R i_1'" H -C i-I -1 - 1) x [ n ( '-1 c_k

,,+!

+ .~) -1] R -R _{• 0+1 k>.} Rs·Rs+1 ~l-o RN+₁-o ( Rl ( _Rl - 1) H b i_1 Ri -1 "" . -d i_1 - 1) H a_i_l Ri _ 1'" . -C i-I -1 i.[2,.J (3.10b) x (R1 - 1) x

-F!.", )

R _ S 8 s _~1!

•• ,)

• •

• (:i!RN+1-O ) , i.[0+I,N+1J (3.l0c)

It can be seen that the algorithms defined by eqns. 3.9 and 3.10 is being constructed around the recurrence relation given by eqn. 3.3 This relation is independent of both the number of tissue layers and the locations of the source. Only the stop-terms, with which the recurrency is stopped, are dependent of the tissue and source configuration. From the numerical analysis, given in the next paragraph, it will become clear that this particular relation can be used to minimize the loss of significant digits. Furthermore, it preserves the

flexibi-lity of the algorithm with respect to the number of tissues. In order to incor-porate one extra tissue in a typical model the recurrent relation needs to be nested only one level deeper.

Inverse Fourier Transfol71l

For the computation of the first M finite integral parts (eqn.2.1B) a five

points integration method is applied [25]. As has been suggested in the previous section, the infinite integral parts will be approximated by an analytic expres-sion which will be derived here. Because the transformed field distribution is proportional to the transformed source distribution, the knowledge about the source distribution itself will be used in the derivation of the analytic expression. In case the two ringshaped charge distribution at distance 2L it follows that ~L f A exp(G a)cos(az)sln(8L)dSt = ~'

,

H (3.lla) A_z GzSln[8M(L~)] - (L+z)coa[aH(L+z>] - u j - exp(G a ) l( { 2 2 + -z 11' Z H [ G + (L+Z) J +

,

GzSin(SM(L-z>] - (L-z)cos[BH(L-z>] [ 02

_,

+ (L_.)2] (3.11b) l!

,

I'

II

I

I

I

.1 / ,1

111m { • 8 -_L

'L

f ~osin(az)deJ 8. 8L

'L

I E_rs1n(Sz)dS J '" 8. .. u _-, 1.

_,

lim I f Arexp(GrB)sin(ez)sin(SL)dst '" ~M ~ ... L A '" ~r ~ _exp(GrBM) x I GrcOslSH(L+z)j + (L+z)8in[8H(L~)] [ G; + (L+,)2₁ Grcos[BL(L-Z)] + (L-z)8in[B M(L-z)] [ G; + (L-Z)2J 0·12a) + (3.12b)

What is left, is the determination of the functions A"G"A_r and G,. Because we

have already used the information about the source distribution in the axial

direction, in order to determine Az.Gz.A r and G

r we consider the field is being

generated by only one ring-shaped electrode in the

z=O

plane. Then the asympto-tic functions are found by calculating the corresponding field component at two adjacent. 8 samples. The asymptotic behaviour of " and " can be described

as and and, also as and leading to and E (6)· A exp(G 5)

,

,

,

E (6)· A exp(G 5)

,

,

,

1 E,CB+M) Gr "

n

In J:;

_,

(B) B > 8

•

8 > 8

•

8 > 8

•

(3.13a) (3.Ub) (3.14a) (3.14b) (3.15a) (J.15b)

On substituting the numerical result of 3.15 into 3.13 it follows that A .. E (t3)exp(-G 6) "

,

and A ~ E (a)exp(-G B) ~

"

,

a > , m a > , m (3.16a) (3.16b)

Substitution of 3.15 and 3.16 into 3.11 and 3.12 results in an estimation of the remainder of the inverse Fourier transform.

3.2 Numerical analysis

The numerical analysis will be restricted here to the cylindrical functionsJ

n, Yn,

and "~2). the recurrence relation", and the stop terms. Because the computation of the electric field distribution in the space-time domain. involves

a numerical integration with respect to the axial wavenumber a • the numerical analysis concerns only the numerical behaviour of the various expressions as a function of the axial wavenumber. with system parameters fixed.

ArgW1Jents

The arguments for which the cylindrical functions must be completed satisfy the following equation:

(3.17a)

which can also be written as

(3.17b)

From eqn. 3.17b i t is seen that there exists an axial wavenumber aa such that

, .2£,,"0

I

- - - «1, a2 a > a a (3.la) H(l) n ~ :

! '

~

I

I

,

j

I t: ; I I,

'I

(3.19)

showing that for large axial wavenumbers the arguments are proportional to the product of the axial wavenumber and the radius, with system parameters fixed. It

must be noted, that the solution of the complex root has to be chosen such that

~vlindrical functions

For the determination of the Bessel and Neumann functions of complex arguments, standard computer library procedures are being used [l,43J. Because the compu-tation of the Hankel functions, involves exponential functions we do not use standard routines. It is known [38J that, for large arguments

(ll l-J nwJ H 1;)" - exp(- - ) W (-2jt;) n 1\'t; 2 on • larg(,) I <, (3.2l) (2) l+; .q..t' H tt;) - ~ exp( ) ~ (2)1;), n TI; on larg(')1 < w (3.22)

where W (±2J') are the Whittaker functions [38J given by on N ED 2 2 1 + E (-J)m[ n ( 4n - (2k-l) + Bkl; m. ... l k-l

lar

g

(!;)I

<" (3.23)

In order to extract the exponential function, two auxiliary functions ~n and are introduced such that

larg(1;)1 < '/I" (3.24)

larg(')1 <, (3.25)

WP

"

This separation of functions is advantageous for the clarification of both the

numerical analysis and the numerical evaluation of the recurrence relation. From