////?

DEPARTMENT OF ELECTRICAL

AND COMPUTER ENGINEERING

-'/$%%

(NASA-CR-177236) CN THE ELECIECMAGNETIC N86-27515

SCATTERING FfiCM INFINITE RECIANGU1&B

CCNDUCTING GRIDS (North Carolina State

Univ.) 171 p HC A08/MI A01 CSCL 20N Unclas

G3/32 43245

V

SCHOOL OF ENGINEERING
ON THE ELECTROMAGNETIC SCATTERING

FROM INFINITE RECTANGULAR CONDUCTING GRIDS

by

CHRISTOS CHRISTODOULOU

DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING

NORTH CAROLINA STATE UNIVERSITY

Raleigh, North Carolina

May 1985

This research was supported by the National Aeronautics and Space

Administration through grant NSG 1588.
ABSTRACT

This report describes the study and development

of two numerical techniques for the analysis of

electromagnetic scattering from a rectangular wire mesh. Both

techniques follow from one basic formulation and they are

both solved in the spectral domain. These techniques were

developed as a result of an investigation towards more

efficient numerical computation for mesh scattering. These

techniques are efficient for the following reasons:

a) They make use of the Fast Fourier Transform.

b) They avoid any convolution problems by converting

integrodifferential equations into algebraic

equations.

c) They do not require inversions of any matrices.

The first method, the "SIT" or Spectral Iteration Technique,

is applied for regions where the spacing between wires is

not less than two wavelenghs. The second method, the

"SDGC" or Spectral Domain Conjugate Gradient approach, can

be used for any spacing between adjacent wires. A study of

electromagnetic wave properties, such as reflection

coefficient, induced currents and aperture fields, as

functions of frequency, angle of incidence, polarization and

thickness of wires is presented. Examples and comparisons

of results with other methods are also included to support

the validity of the new algorithms.
Ill

ACKNOWLEDGEMENTS

I would like to thank Dr. J. F. Kauffman, the chairman

of my advisory committee, for his guidance, support and

invaluable assistance with this work. Many helpful

discussions with Jerry Brand and Dr. J. Dunn are also

gratefully acknowledged. Finally, I would like to thank Dr.

M. C. Bailey of NASA Research Center at Langley for his

constructive criticism and comments.
IV

TABLE OF CONTENTS

1. INTRODUCTION ........ 1

2. THE SPECTRAL ITERATION APPROACH 4

2.1 Formulation 4

2.2 Convergence of iterative scheme 11

2.2.1 One Dimensional case .13

2.2.2 Two Dimensional case 14

a) First method 15

b) Second method .17

c) Third method 21

d) Fourth method 23

2.3 Comments. 26

3. THE SPECTRAL DOMAIN CONJUGATE GRADIENT APPROACH. . . ..27

3.1 Review of the conjugate gradient theory 27

3.2 Current density formulation 31

3.3 Numerical properties of the conjugate gradient

method ' .37

3.3.1 Singular operator. 37

3.3.2 Convergence rate 38

3.3.3 Stability 39

3.3.4 Global versus local minimum 39

3.4 Stopping procedures and Initial guess 40

3.4.1 Stopping procedures 40

3.4.2 Initial guess 41
4. FORMULATION OF THE S.D.C.G. METHOD FOR THIN WIRES WITH

FINITE CONDUCTIVITY AND FOR APERTURE FIELDS. .... .42

4.1 Equivalent radius concept and internal

impedance 42

4.2 Solution of aperture fields 46

4.3 Reflection coefficients 49

5. RESULTS AND COMMENTS . .52

5.1 One dimensional case 52

5.1.1 Current densities 52

5.1.2 Aperture fields 58

5.2 Two dimensional case . .62

5.2.1 Current densities and reflection coefficients .62

5.2.2 Aperture fields 91

5.3 CPU time and storage requirements 102

6. CONCLUSIONS AND SUGGESTIONS FOR FUTURE RESEARCH. . . 103

7. LIST OF REFERENCES 110

8. APPENDICES 114

8.1 Derivation of the equation for H as a function

of the electric vector potential F 114

8.2 Derivation of the equation for the scattered

magnetic field Hs 116

8.3 Minimization in the domain 119

8.4 Code 120

8.5 Flow-Chart for the S.D.C.G. method 124

8.6 Listing of the S.I.T. method 125
VI

8.7 Listing of the S.D.C.G. method for thin strips

with the square-shaped unit cell 132

8.8 Listing of the S.D.C.G. method for thin strips

with the cross-shaped unit cell 139

8.9 Listing of the S.D.C.G. method for Aperture

Fields 145

8.10 Listing of One, Two and Three dimensional

complex fast Fourier transforms 151
1. INTRODUCTION

A new technology for large space-based systems requires

antennas with 100 meters or larger in diameter for radio

frequency operation, communication, earth observation and

radio astronomy applications.

A new type of antenna, the MESH DEPLOYABLE ANTENNA,

which appears to be more cost-effective and easier to ,

transport into space compared to a solid reflector of 100

meters in diameter, was the motivation for the study

reported herein. The mesh used to construct large space

reflector antennas is usually made of gold-plated molybdenum

wire about one mill in diameter. The wires run and cross in

a weave that is periodic in nature, forming a reflecting

surface that behaves differently depending on the number of

openings per wavelength and polarization of the incident

energy. The undesirable effects resulting from such a

surface include transmission loss, resistive loss, and cross

polarization loss.

Here a wire mesh structure is used as a simplified model

of the knitted (woven) material. A rectangular mesh

structure is a periodic structure, and scattering from

periodic structures is a subject that has a long and

illustrious history dating back to Lamb and Rayleigh in the

last century [1-5] .
Constructing solutions to the problem of mesh scattering

can be achieved using a variety of methods. One possible

method which has been widely used is the METHOD OF MOMENTS

(MOM) [6-8] . This method, when applied to periodic

surfaces, has the disadvantage of requiring the inversion of

a very large matrix, a fact that renders the method

unwieldy. Other methods involve COUPLED INTEGRAL EQUATIONS.

These methods will not usually yield a solution due to the

complexity of inverting the integrals for a periodic mesh.

Another popular technique used for estimating the reflection

coefficient from a wire mesh is based on AVERAGED BOUNDARY

CONDITIONS [9-10]. This method offers good results when the

number of mesh openings per wavelength is large [11].

However, even this method fails for certain applications

when the number of openings per wavelength becomes small.

This dissertation includes the analysis and formulation

of two new models for studying scattering from wire meshes

that are more efficient and simpler to apply than the

previous methods. The first method is based on the SPECTRAL

ITERATION APPROACH (SIT) [12-18] which is valid for cases

where the spacing between adjacent wires is larger than two

wavelengths. This limitation on the size of spacing between

wires for the SIT method led to the development of the

second model which is valid for all spacings. This new

model is the SPECTRAL DOMAIN CONJUGATE GRADIENT method

(SDCG) [19-22] and is a combination of the SIT and the
Conjugate Gradient method. Both methods utilize the fast

Fourier transform and avoid any convolution problems and any

inversion of matrices.

These two techniques offer new accurate models which can

be extended and applied to the more difficult problems of

knitted mesh surfaces. A number of examples are computed

and compared with other methods. Also, comments and

suggestions are made for possible extension of the SDCG

method to the more complicated problem of the knitted

structure.
2. THE SPECTRAL ITERATION APPROACH

2.1 FORMULATION

Any scattering problem could be expressed in the form of

the integral equation:

«(x)= /K(x,x') <?(x«) dx + «inc(x) (2.1)

with the constitutive equation W(x)=K(x) <3>(x) (2.2)

where K(x,x') is the kernel of the integral transform

i nc

* (x) is the externally applied field

<D(x) is the field quantity, and

W(x) is the source quantity

The S.I.T. method is a frequency domain (Spectral

Domain) solution, and consists of casting the general basic

global equations (i.e. the second order partial differential

equation or its integral representation, such as equation

(2.1)) as a local algebraic equation in the Fourier

transform space, and leaving the local constitutive equation

as a local algebraic equation in real space. That is,

taking the Fourier transform of equation (2.1) and keeping

(2.2) the same, one arrives at:

*(k) = K(k) W(k) +«inc(k) (2.3)

V (x) = K(x) 4> (x) (2.4)

Equations (2.3) and (2.4) show how the original set of

equations are converted into a set of two simultaneous

algebraic equations in two unknowns (the fields and the

induced currents) in two different domains connected by the
Fourier transform which is given by:

C°° A A

F(k) = / f(x) exp(jk.x) dx (2.5)

-Lao

The operation in equation (2.5) from now on will be denoted

by the transform pair:

**

F(k)-f->f(x) (2.6)

By virtue of the numerical Fast Fourier transform and the

local algebraic representation, the number of required

complex multiply and add operations and the number of

required storage locations are of the order of Nlog2N and N

respectively (where N represents the number of Floquet modes

or cells into which the problem is discretized).

For periodic structures, the Floquet theorem [23] is

used to account for the periodicity of the wire mesh and the

coupling between adjacent wires. The specific equations for

a wire mesh (See Figure 2.1) are formulated as follows:

Fig. 2.1. Wire mesh geometry
The electric field E due to a magnetic current M is

given by:

"E(x,y)=-l/e VxF(x,y,z) (2.7)

where F is the associated electric vector potential of the

source and e is the permitivity of the medium in which the

wire mesh is placed. F" and ft are related by the free space

Green's function G=exp(-jit.r)/4n r by

"£("?) = /G(rMT') M(f') dr1 (2.8)

From this the, magnetic field intensity, H, can be derived

(See Appendix 8.1) from Maxwell's equations and expressed

as:

-K -»» p-Tf—i ""*•

H(x,y,z)= -j to F(x,~y,z) + \/\/.F(xfy,z)/jcoeu (2.9)

whereu is the permeability of the medium. Since we have a

planar structure F is set equal to zero. Now expanding

equation (2.9) in terms of its Cartesian coordinates x and y

yields:

2 2

;

)Fx* +(ko2+-^+-r:}F

yy ' Oxdy dy2

(2.10)

A planar periodic structure such as that shown in Figure

(2.1) could be considered to be the source distribution for

the magnetic field of the equation (2.10). Substituting

equation (2.8) into equation (2.10) and taking Fourier

transform of equation (2.10) yields the transformed

scattered tangential fields at z=0 in the following form:

1

ii ( x , y i —

JOT|1C

r

 2 2

IV L

 a .,

° 1 A

2

d

dx
2 2

o ' mn mn mn

-a B k 2- 6 2

mn mn o mn

M(

(2.11)

where

= 2 Tt m/a -ko sind coscp

3 =2 TI n/c - 2 rt m/a cotn -ko sinO sincp are the Floquet

modes [24] and

amn' 3mn}

transform of Green's function.

is the Fourier

Fig. 2.2 FSS Geometry

Taking the inverse Fourier transform of equation (2.11)

yields:
H

s(x,y)=

mn

2 _ 2 B

o mn mn

 p mn

a 0 k 2- 0 2

mn mn o mn

G( amn'

amn'

.exp[j( amnx

(2.12-a)

Now, by using the equivalence theorem and applying the

appropriate boundary conditions on Iis(x,y) at z=0 (See

Appendix 8.2) leads to:

H

inc_ -2

mn

a

 mn ^ mn o ~

 a mn

2 - • 2 '

mn * o ~

 a mn ^ mn

.exp[j(amnx

(2.12-b)

where E represents the transformed electric aperture field

and Hlnc is the incident tangential magnetic field. To

extend the formulation over the full range (i.e. to include

conducting regions), the current densities have to be added

to equation (2.12) to give:

e[J(x,y)]= Ht1110^,

a

mn mn

? 2 —I

k

o-a mn^

2 2

mn ~

ko ~

 a

mn mn

•

exp[J( amnx

 *Pmny )l

(2.13)
where 6 is the complement of the truncation operator

defined as:

9[X(?)]=X(r) for r in the aperture (2.14)

and Q[X(r)]=0 for r in the conducting regions

and 9[X(r)]=X(r)-9[X(r)] (2.15)

Note that in equation (2.13) J and E are both the unknowns

to be solved for.

Equation (2.13) can be recognized as the discrete

Fourier series for a periodic sequence [25]. Note that

there is a direct duality between the (x,y,z) domain and the

(k rk ,k ) domain. Since all the functions involved here

x y z

have a 2tt /m and a 2ti /ri periodicity in their exponents, one

period (i.e. one cell) of the structure is sufficient to

completely specify the transform. That leads to the use of

the discrete Fourier transform which can be evaluated very

efficiently by the Fast Fourier transform. It should be

noted here that because of the exactness of the duality

between the two domains, no aliasing effects will appear

when the FFT is performed. By aliasing we mean overlapping

of spectral components.

Besides equation (2.13), the boundary condition that

governs the behavior of the tangential components of the

electric field, E, over the conducting regions has to be

satisfied. Equation (2.13) can now be rewritten as :

+ 3(J)] (2.16)
10

where E. is the Fourier transform of E.

F is the Fourier transform and p ~ is its inverse

I, is the product of the Floquet expanding modes

and Green's dyadic in the spectral domain.

If the induced currents were available, the solution of If.

could be immediately obtained from (2.16). In practice,

however, J is an unknown to be solved together with "£. and

equation (2.16) cannot be solved directly. Instead, using

equation (2.16) a recursive relationship between the (p+l)th

approximate value of E* and the pth approximation of E is

now derived and both E. and J are computed simultaneously,

via the following iterative procedure:

a) Start with a guess for IT. in the (x,y) domain and

apply the truncation operator (i.e. apply the

boundary condition that E =0 over any perfectly

conducting surfaces).

b) Take the Fourier transform of E.

c) Solve for J(P

)

-F"1[iF<8it

(P))l+HtinC

(2.17)

d) Set currents J equal to zero everywhere except over

the conducting surfaces, that is find:

^IFo St

(P)-)]+5tlnc| (2.18)

Substituting equation (2.18) into (2.16) yields:
11

(2.19-a)

And finally taking the inverse Fourier transform of (2.19-a)

yields:

^(P+lJfr-il-iprSfr-if ek'P' + F'1?"1 Ft-Ht

inc+0(Htinc)]

(2.19-b)

Note that once E is evaluated J can also be computed.

Equation (2.19-b) could be cast in a more convenient form

(operator form) as:

Et

(P+1>-I. VP>+* (2'20)

where L = F~1Z~1F^F~1Z

 8 is an operator

and C - F'1^'1 F[-Ht

inC+®<Htinc)] is a constant that

depends on the initial conditions and the incident wave.

The two most attractive features of this method are the

following:

a) No extreme amount of computer memory storage is

required.

b) No explicit knowledge of appropriate basis functions

is needed.

However, like most iterative techniques, the basic iterative

scheme suffers from convergence problems. These problems

and the attempts to alleviate them is the subject of the

next section.

2.2 CONVERGENCE OF ITERATIVE SCHEME

To achieve convergence the important condition that has
12

to be satisfied is that p (L) < 1 or that the spectral

radius of the operator L has to be less than one. As it

turns out, for two dimensional cases where the wire spacing

is greater than two wavelengths, p (L) < 1 and hence

equation (2.20) converges very quickly for any type of

incident polarization, angle of incidence and wire

thickness. However, for spacings less than two wavelengths

the method fails miserably. To achieve convergence in those

cases the Successive Relaxation method could be employed to

"relax" the process and force p (L) < 1 for some relaxation

factor ©. The choice, not only of the optimum relaxation

factor, but simply of a relaxation factor that would produce

a convergent scheme is a difficult task indeed.

In the one dimensional problem (parallel grid) the

Contraction Mapping Theory was used very successfully to

obtain the optimum relaxation factor 9 which forces the

spectral radius to be less than 1 [26]. To show how this

theory was used, equation (2.20) is rewritten as:

g(xn) = xn+1 = L xn + C (2.21)

Define a new mapping G(xn) so that: ~~

G(xn) =6x + (1-0) g(xn) (2.22)

According to the contraction mapping theory [27-31] a

transformation G of a metric space X onto itself is Lipshitz

continuous if there exists a p, independent of x and y such

that

d(G(x),G(y))<Pd(x,y) for all x,y,e X where d(x,y) is a
13

proper metric in X. For strictly contractive mappings

less than one.

2.2.1 One Dimensional Case

For the one dimensional case the simplest possible

metric d that can be used to obtain the optimum 9 is

chosen as follows:

is

|G(y)-G(x y-x for P< 1

Let y=x +

G(V 6 }

"

G(xo

then

< P |6| or G(x )-G(x < P

(2.23)

(2.24)

So the necessary and sufficient condition for contraction

mapping becomes:

d (G(x)) <p (2.25)

'dx

Now substitute (2.22) in (2.24) to obtain:

|8(x, g(xQ-i-6)-ex0-Kl-e) g(xQ; p|6'| or

9+ (1-6) dg(x)

1 dx

<P

Setting p =0 in the above equation and solving for 6

yields:

8 = (dg(x)/dx)/(dg(x)/dx-l) (2.26)

This value of 6 is called the "contraction" factor since it

will yield a convergent scheme even in those cases where the

basic iterative scheme of equation (2.20) fails to converge.

It should be noted here that in the above analysis © is

treated as a constant when in fact it is a function of x.

The reason for that treatment is that 9 is solved in the
14

neighborhood of a solution (root) x where the values that 9

acquires are approximately equal. Therefore 9 can be

assumed to be constant within that particular neighborhood.

(For examples and results on the one dimensional problem see

[26]).

2.2.2 Two Dimensional Case

In two dimensions, the basic iterative scheme of

equation (2.21) is given by:

(2.27)

or

x"*1 = Lll xn + L12 .yn +.C1

(2.28)

Equation (2.28).can be expressed in the more convenient

form:

x

n+1

 = g(xn,yn)

y

n

+

l

 = h(xn,yn) (2.29)

To achieve convergence in the two dimensional problem, the

four partial derivatives QX'QV'^X

 and h

v

 must satisfy the

following condition [32-35] :

|gy|<

" xn+

f

-

 xn+1-

=

Lll L12

L21 L22

~xn~

y

n_

Cl

C2

x

n+1

 = Lll xn + L12 .yn + -C1

y

n+1

 = L21 xn + L22 yn + C2

hy|< k2 (2.30)

for kl and k2 less than one and for all points (x,y) in the

neighborhood R of the root (xo,yo), where R consists of all

(x,y) with|x-xoj< e , |y-yo|< e , for some positive e. For

wire spacings less than two wavelengths condition (2.30) is
15

not satisfied. Thus, one has to construct new mappings

(functions) for the system in (2.27) to obtain convergence

in a manner similar to the one dimensional case. Now, to

apply the method of contraction mappings the system (2.29)

is rewritten as:

G(xn,yn) = 9vx

n

 + (1-6 ) g(xn,yn)

X A

H(xn,yn) = 6yy

n

 + (1-6 ) h(xn,yn) (2.31)

where 6 and 6 are relaxation factors.

x y

Unlike the one dimensional problem, Gx, Gy, Hx and Hy

cannot be separately set equal to zero since they would

produce a system of equations that are impossible to solve

for gx=0, gy=°/ hx=0 and "h -0, i.e.

8

x

(1-6) g = 0 for 9 (2.32)

A y A

and

8y * (1'V hy

 = °

(1-6 ) hv = 0 for 6 (2.33)

y x y

One way to avoid this difficulty is to set kl and k2 to

nonzero values but their absolute value must always be less

than one.

a) First Method

Let kl and k2 less than one in equations (2.32) and

(2.33) to obtain:

l

e + (1

- 8

(1-6) gv I < 1/2 (2.34)

A y
16

and

I 9 + M-8 ) h

1(1-9) hv I < 1/2 (2.35)

y xi

Since hx, hy, gx and gy are complex numbers that implies

that 9 and 9 can acquire complex values and hence & and

A y A

9 are expressed as:

9x = a + jb (2.36-a)

9 = c + jd (2.36-b)

Moreover, let

real (gx) = a imaginary (gx> = 0

real (g ) -Y imaginary (g ) =6 (2.37)

Upon substituting equations (2.36) and (2.37) into (2.34)

one obtains:

a + jb + (1-a-jb) ( a+ jb) | < 0.5 (2.38-a)

|(l-a-jb) (Y+J6) | < 0.5 (2.38-b)

Taking absolute values yields:

[(a +a -a a + b3 )2 + (b +p - a£- ba) 2]

 1/2

 < 0.5

or

[(Y - aY+ b6) 2 + (6 - a6- b Y )

2

] < 0.5

(a +a-a a +• b3 ) 2 + (b +3- a 3- ba) 2 < (0.5) 2

(2.39-a)

(Y-aY + b6 )2 + ( 6- a6- bY )2 < (0.5) 2

(2.39-b)

Equation (2.39) can be expanded to yield two nonlinear

equations in two unknowns a and b of the form:
17

Al a2 + A2 b2 + A3 a + A4 b + A5 =0.402

a

2 -I- b2 -2 a +1 = .42 /( Y2 + 62) (2.40)

where Al, A2, A3, A4 and A5 are constants that depend on g

and g . Similarly, to solve for 9 =c+jb another set of

nonlinear equations is to be solved:

Bl c2 + B2 d2 + B3 c + B4 d + B5 = 0.402

c

2

 + d2 -2 c + 1 . = .42 /(e2 +n2) (2.41)

where e=real (h ), n=imaginary (h ) and Bl, B2, B3, B4

and B5 are constants that depend on h and h .

x y

The solution of these nonlinear equations give 9 and 9

that would be expected to yield a convergent scheme but,

unfortunately, they fail to do so for a wire mesh. This

failure is attributed to the fact that the chosen metric is

not the appropriate one for this type of geometry, whereas

it could be a good choice for other geometries of frequency

selective surfaces. This fact leads to another choice of a

metric space,

b) Second Method

This time the Euclidean norm is chosen as follows:

2 . i «„ i 2 . I u_, i 2 % 1/2 M 2

(2.42 )

It is desired to solve for Gx, Gy, Hx and Hy that are

functions of 9x and 9y with the hope to yield ||M||2 <1. So
18

the basic minimization scheme for solving for x and y in

this case is the following:

an 2 .„ and

3

9

 <2.43)

It was found previously that Gx, Gy, Hx and Hy can be

written in terms of x, y as:

Gx = Q + (1-6 )g

x x/yx

Gy =

 (1

~

e

x

)

9y

Hx = (1-61 )h

y

 x

Hy =9 + (1-6 )h (2.44)

and hence

= (Gx Gx + Gy Gy + Hx Hx + Hy Hy ) (2.45)

Substituting equations (2.44) into (2.45) and after a long

and tedious manipulation one obtains:

l|M|2=<al

2 A1+a

1BH31Al+j31d + a2

2A2+ a2B2+022A2+j02C2+d)1/2

, , * <2.46)

where Al= |g I +1+ |gxf -9x-gx

Bl=gx - |9x|

2

 + gx* - [gj2

 -2 |

A2=l+ |h/-hy-hy* + |hj2

B2

=V ihy|2 +hy* - lhy|2 ^

C2=h * -h

y y

Now to solve for 9x= a.^+j 3^ and 9 = <*2 +j 32 one

needs to solve the following system of equations:
19

ML 1 1

= o = — - - [2aA,+B ]

d ||M|| i

F,(a ,0 ,a ,0 ) = - - = 0 = - - [2 a0A9+B,]

* i i * * ,5 „ o iiui i * * * d a 2

d ML 1

F (a, ,0,,0-re,)- ^-= 0 = — [2 01A,+jC1]

•3 1 ••• * * ^ 0 1 XI 1 d 0 2

ML 1 1

= o - - — [2

d 02 2

(2.47)

By using Newton's method or any other minimization method

one can solve for a,, a2 0, and 02 which will give

the values for Q and 6". Unfortunately, once more the

x y

values of Q and Q obtained by this method yield values

x y

I|M|L>1 for some points inside the cell. It should be noted

here that the condition that ||M|L<1 should be satisfied at

all sampled points in the cell, and the violation of this

condition at one point is enough to affect all the other

points since they are all related together via the two

dimensional Fourier Transform.

^Figure (2.3) shows a 16 by 16 array of sampled cell

points and the value of ||M|L at each point. It can be seen

that the condition the |JM L<1 is violated at numerous

points, which implies that a contraction mapping cannot be

achieved by this method. It was observed that the smaller

the wire spacing the larger the values of ||M L become,

especially near the edges of the wires.
20

A= 0.25610000 B= 0.25010002 C= 0.25610000 D= 0.25010002

FREQ = 0.2998E+09

PHI = 0.0 IHETA = 0.0 PSI = 90.0

NX = 14 NX1 = 2 NX2 = 15 NY = 14 NY1 = 2 NY2 = 15

0.3 2.8 6.9 4.3 6.8 4.8 6.5 5.2 6.1 5.7 5.6 6.1 3.6 0.5

2.8 0.1 1.8 2.2 1.3 2.0 1.3 1.9 1.1 1.4 1.3 2.3 0.1 7.8

6.9 1.8 0.0 0.8 1.6 0.6 1.6 0.6 1.5 0.5 1.1 0.0 2.3 4.0

4.3 2.2 0.8 0.0 0.1 0.6 0.2 0.5 0.4 0.4 0.0 0.9.2.1 6.8

6.8 1.3 1.6 0.1 0.0 0.1 0.5 0.3 0.5 0.0 0.4 1.2 1.1 5.0

4.8 2.0 0.6 0.6 0.1 0.0 0.2 0.1 0.0 0.4 0.6 0.6 1.7 6.3

6.5 1.3 1.6 0.2 0.5 0.2 0.0 0.0 0.1 0.5 0.1 1.4 0.9 5.5

5.2 1.9 0.6 0.5 0.3 0.1 0.0 0.0 0.2 0.2 0.7 0.1 1.9 5.8

6.1 1.1 1.5 0.4 0.5 0.0 0.1 0.2 0.0 0.6 0.2 1.6 1.2 6.0

5.7 1.4 0.5 0.4 0.0 0.4 0.5 0.2 0.6 0.0 0.4 0.6 1.9 5.3

5.6 1.3 1.1 0.0 0.4 0.6 0.1 0.7 0.2 0.4 0.0 1.7 1.2 6.7

6.1 2.3 0.1 0.9 1.2 0.6 1.4 0.1 1.6 0.6 1.7 0.0 1.5 3.9

3.6 0.0 2.4 2.1 1.1 1.7 0.9 1.9 1.2 1.9 1.2 1.5 0.0 8.2

0.2 7.8 4.0 6.8 5.0 6.3 5.5 5.8 6.0 5.3 6.7 3.9 8.2 0.5

Fig. 2.3 The values for ||MJL at each sample point inside

aperture
21

Since neither one of the previous chosen metric spaces

appear very promising for this particular geometry of

frequency selective surfaces (i.e. a planar mesh) the trial

of different metric spaces is put to an end and a different

line of thought is followed in the next method,

c) Third Method

Instead of using 9 and 9 , four different relaxation

x y

factors 911, 912, 921 and 922 could be utilized to offer

more degrees of freedom in satisfying condition (2.30).

Thus, the new modified system of equations becomes:

.n+1

n+1

911 912

921 922

n

n

(l-9'Il) - 912

-921 (1 - 922)

n ..n,

g(x",y")

n n,

h(x",y")

n n,

=G(xll,y")

=H(xn,yn)

(2.48)

Now it is easy to set all four partial derivatives Gx, Gy,

Hx and Hy equal to zero to obtain:

Gx=911+(1-911) g -912 h =0

A *V

Gy=912+(1-911) g -912 h =0

Hx=921-921 gx +(1-922) hx=0

Hy=922-921 g +(1-922) h =0 (2.49)

Solving this system of equations for 911, 912, 921, 922

yields:

911 = -

h x g y - (2.50)
22

912=

gy - d-gx) d-hy) (2.51)

621=

 X

hx gy - d-gx) d-hy) (2.52)

h

x

 g

922= — —

hx g - U-gx) d-h ) (2.53)

Again, this choice of 8's works very well for the one

dimensional problem but it does not lead to convergence for

the two dimensional wire mesh problem.

To explain why this'method does not work for the two

dimensional problem the theory for constructing convergent

iterations for a pair of trancendental equations is invoked,

According to this theory the original system of equation

x

n+1=g(xn

fy

n)

y

n+1=h(xn

,yn)

can be written as:

x

n+1=xn

+a[g(xn,yn)-xn]+e[h(xnfy

n)-yn]=G(xn,yn)

y

n+1-yn

+Y [g(xn,yn)-xn]+6 [h(xn,yn)-yn]=H(xn,yn)

(2.54)

Note the similarity of the above equation with equation

(2.48). The parameters a, 0 , Y and 6 play the same role in
23

equation (2.54) as the relaxation factors ©11, 912, 821 and

922 in equation (2.48). To find the root of equation (2.54)

it is desired to determine a ,0 , Y an<3 6 , by the four

conditions that the first partial derivatives of G and H are

zero at some point (x,y) that hopefully is near the root.

Note that the unknown parameters enter linearly in the same

way as 9's do in equation (2.48), so the calculation of the

partial derivatives G ,G ,H and H posses no problem. For

x y x y

the case of trancendental equations, it is known that this

method of constructing convergent schemes works provided

that the partial derivatives gv,g_,,h and h DO NOT vary x y x y

very rapidly in the neighborhood of the root (x

0fY0)« Thus,

although it is easy to produce a G and an H that are well

behaved at the root (x ,y ) they may behave quite badly a

small distance away. If this strategy is to be successful G

and H must not only have small partial derivatives in some

region, but this region must also include the desired root.

For the two dimensional wire mesh it was found that the

derivatives g ,g ,h and h vary very rapidly, especially at

x y x y

points close to the edges of the wire. So this fact, and

the lack of knowledge of the region within which a root

exists, causes this method to fail,

d) Fourth Method

Finally, another method that could be tried to solve for

x and y is Newton's method. In this case, we start with the

basic iterative scheme:
24

.n+1

n+1

Lll L12

L21 L22

x

n

y

n

+

+

Cl

C2

which gives:

c

n+1

 = Lll xn -»- L12 yn + Cl

.n+1 L21 xn + L22 yn + C2 (2.55)

Since convergence means that for large n x-—>. xn equation

(2.25) can be rewritten as:

x-Lll x -L12 y - Cl =0

y-L21 x -L22 y - C2 =0 (2.56)

Note that Lll, L12, L21 and L22 are operators so one can

solve the above equation for a root (X

0»y0) by employing

Newton's method.

The convergence of this formulation though suffers since

the derivatives, gv,g.,,h and h are much larger than one

x y x y

for wire spacings less than two wavelengths or so. This

fact by itself causes this method to fail.

Figure (2.4) shows how the relaxation factors 9V and 9

x y

contract the basic iterative scheme, but still not enough to

push the iteration into the region, of convergence.
25

10 - &

p_.,T(r+6)-T(r)|

R " for 6=(o.ol,o.ol)

r+6 -r

R <1>-*CONTRACTION

R> 1-* DIVERGENCE

No relaxation factor

r n+1 n n

x =o.5x +o.5 Lx

n+1 n n

y =o.5y +o.5 Ly

n+1 n

x =9 +(1-9 )Lx

X X

n+1 n

y =9 +(1-6 )Ly

y y

CONVERGENCE REGION

3X

Figure 2.4. Contraction effect of different relaxation

factors.
26

2.3 COMMENTS

It is believed that, unlike the one dimensional problem,

the two dimensional problem has functions and partial

derivatives that are very steep, so any method that depends

in a critical way on magnitudes of derivatives will have

difficulty to converge. It is also believed that all the

above mentioned methods for obtaining a convergent scheme

can be very effectively applied to other geometries of

frequency selective surfaces, such as an array of metallic

patches, cross dipoles, circular apertures, etc.

In conclusion, this method works very well for large

spacings between adjacent wires without making use of any

relaxation, contraction or variational factors, but it fails

miserably to converge for two dimensional problems where the

mesh spacing is less than two wavelengths or so.
27

3. THE SPECTRAL DOMAIN CONJUGATE GRADIENT APPROACH

Unlike the previous (S.I.T.) method, in this method the

induced currents and the aperture fields are solved

separately. The common features of the S.D.C.G. method and

the S.I.T. method are that they are both solved in the

spectral domain and that both make use of the fast Fourier

transform. In the S.D.C.G. approach, the conjugate gradient

method is employed to improve upon each previous iterate.

Hence., the method is basically an iterative technique.

This part of the the dissertation includes the analysis

and formulation of the problem in the spectral domain for

both current densities and electric fields and their

solution via the conjugate gradient technique. Moreover, a

number of numerical properties for the conjugate gradient

method are discussed, and ways of terminating the iterative

process are suggested.

3.1 REVIEW OP THE CONJUGATE GRADIENT THEORY

Suppose that the system that is to be solved is given

by:

A "x = "y (3.1)

Let x be and initial guess for x and the residual error

vector be:

r

(0) =?- Ax< 0 ) (3.2)

If A is symmetric positive definite then A is also
28

symmetric positive definite. Now define the quadratic error

functionals as:

ERRF1 = r A"1 ? =<r,A~1r>

ERRF2 = r* r = <r,r> =11 rll2 (3.3)

ERRF3 = r* (A A*)"1? = <r, (A A*)~1r>

where the asterisk * means the conjugate transpose.

All error functionals in equation (3.3) are positive for all

values of x( except for x* ' =x . where "x. is the exact

6 6

solution of Ax=y. In the case where x* is equal to the

exact solution x all the error functionals in (3.3) would

be equal to zero.

Now, substitute equation (3.2) in the first error

functional of equation (3.3) to obtain:

ERRF1 =<(y -At(0)) , A'1 (y - A x(0)) > or

ERRF1 =<(?, A"1^> - 2<y-,*(0)> +<x(0), A1c(0) >

(3.4)

ERRF is now a quadratic equation function in x . Let

-»»( n \

x

v be a point in N-dimensional space. Then the equation

^

(0)

 =t(n) ^04, ^(n) (3.5)

is the equation of the line through point x ' in the

direction of ~p , called the direction vector. For a two

dimensional interpretation see Figure (3.1). The parameter

is proportional to the distance pc — ^ . Substituting

equation (3.5) in equation (3.4) leads to:
29

ERRF1-

terms

> -2 an <p"(n), t(n) > + Other

(3.6)

,1 - ERRF is now a quadratic function with respect to a and

has a local minimum which is found by differentiating

equation (3.6) with respect to a , i.e.

1 O( ERRF )

A p<n) =0

line x(0)-x(n)+oi,p(n>

correct solution x

minimum

ntoura of

constant error

t solution

Fig. 3.1. The error functional and the conjugate gradient

method in two dimensions
30

from which one can solve for a t«? obtain:

<p(n),r(n>> a =

n —_,./nx ;—r—

 <p^n ; ,Apinj > (3.7)

Once the position of the local minimum has been found, the

next trial vector can be defined as:

x

 = x + <x p (3.8)

From Figure (3.1) one can see now that at each iteration a

new local minimum is found until the global minimum is

reached.

There are two basic methods that can be used here to

obtain the next trial vector. The first is the steepest

descent method and the other one is the conjugate gradient

method. These methods differ only in the choice of their

> / — \

direction vector p . Sarkar showed how the steepest

descent method method can be applied for electrostatic

problems [36] . In the conjugate gradient method, the

direction vectors, ~p , must mutually orthogonal with

respect to the the matrix A. That is,

<p(i) , A "p(j) > =0 for i * j (3.9)

The iterative scheme of the conjugate gradient method, which

can now be used to yield successive approximations towards

the correct solution, is given by Hestenes and Stiefel [37] ,

and A. Jennings [38] as:

First, let the initial vector (i.e. for n=0) be:

= y - A t(0) . (3.10)
31

The equations for the nth iteration are:

an = <s(n)

 , 3<n)

n

r ' = r(n) - a s(n)

n

n

It can be shown [38] that the following othogonal

relationships are also satisfied:

<r(l) , p(^> =0 for i > j (3.12)

<r(l) , ~r(i

}> =0 for i ^ j (3.13)

3.2 CURRENT DENSITY FORMULATION

The magnetic field H due to an electric current density

-*•

J is given by:

A (x,y,z)

H(x,y)= - (3.14)

U

•*• -».

where A is the associated magnetic vector potential. A and

J are related by the free space Green's function
32

G(r) =

exp(-jk.r)

4 TI r

as follows;

A (r) =u / G(r,r'). J(r' ) (3.15)

From this the electric field intensity Es can be derived

from Maxwell's equations and expressed as:

E (x,y,z) = -jw A(x,y,z) +

jcoue (3.16)

For a planar structure we set the z-component of the

magnetic vector A equal to zero. Now, upon expanding

equation (3.16) in cartesian- coordinates we obtain, for z=0,

~ES(x,y)

jcoe

/•G.Jx

d*dy • Jy

(3.17)

Considering the periodicity of the two dimensional structure

shown in Figure (3.2) (planar structure), and taking the

Fourier transform of equation (3.17) leads to:

E

s(a ,3 ) =

mn' mn'

k

 2- a2

o mn

-a 3 mn mn

— o. S k ~ B

mnp mn o p mn

(3.18)

where the sign (~) denotes the Fourier transformed quantity.
33

a and 0 represent the Floquet coefficients which were

mn mn

defined in the previous chapter as:

a =2n m/a -ko sindcoscp mn

and

mn

=2n n/c -2n m/a cot Q -ko sin.d

G< <W &mn

)=

 -J/2 (k°

2 - "'mn mn'

is the

Fourier transform of Green's function, and J , J are the

x y

unknown current densities.

Fig. 3.2 Frequency selective surface geometry
34

Notice that the spectrum of E*s is discrete. That is, it

exists for discrete values of a and &„,„. Note, also, that

mn mn

the convolution problem is avoided and instead of dealing

with an integrodifferential equation we have to consider

algebraic equations.

Taking the inverse Fourier' transform of equation (3.18)

yields:

I2s(x,y) =-

joe

mn

2 • 2

t - a - a B o mn mn

 K mn

mn

k

 2

k

-

o mn

G . J

(3.19)

To enforce the boundary condition over the surface of all

metallic regions we require that the total tangential

electric field should satisfy the condition:

E

1 (x,y) + ES(x,y) =0 (3.20)

where E is the incident electric field and

E is the scattered electric field

Substituting for the value of Es from equation (3.20) into

equation (3.19) yields:

-E1-

, ^ — '

J_\

j coe L i

mn

k — ex — a £

o mn mn

-

amn &mn "o^0

mn

2

mn

.exp[+j(ax+6 mn m n y ) ]

(3.21)
35

Equation (3.21) can be recognized as the inverse

discrete Fourier transform which can be performed via the

fast Fourier transform (FFT). Equation (3.21) could now be

written in an operator form as:

- "E1 = Zmn "Jmn (3.22)

where Z is the product of G, the Floquet modes and the

mn

inverse Fourier transform.

A solution of the above equation will yield the unknown

current densities J and J from which the reflected and

x y

transmitted fields can be obtained and hence the reflection

and transmission coefficients could be calculated, it

should be mentioned here that like the spectral domain

iteration approach, the spectral domain conjugate gradient

method is independent of basis functions.

Now, one way to solve for J and J is to use the

x y

conjugate gradient method [19,20,37,38]. The conjugate

gradient method in the spectral domain was used by other

investigators [21,22] on different geometries. In our case,

the algorithm of equation (3.11) cannot be directly applied

on equation (3.22) since Z is symmetric but not self

mn

adjoint or positive definite. To over come this difficulty

and guarantee a convergent scheme, equation (3.22) has to be

properly modified. To do that, multiply both sides of

*

equation (3.22) by Z (i.e. the conjugate transpose of

Z ) to obtain:

mn
36

-z mnmn

= Z mn Zmn J inn mn

(3.23)

where the product Z Z is a Hermitian matrix and

mn mn

therefore positive definite. That also means that the

algorithm (3.11) can now be applied to the transformed

equation (3.23). In fact, one can apply the previous

algorithm on equation (3.23) without actually forming

*

Z Z explicitly via the following algorithm [37,38]:

Let J* ' be the initial guess and let the initial residual

vector r^ be:

*<°> = z

mn

k

mn

ERRF =

th , The equations for the n iteration are:

a

n

 =

mn

mn

n

ERRF

(n+1)

(3.24)

= :n)

3

mn

n

mn
37

In the above algorithm the root mean square error

|r r|J was chosen as the quadratic functional to be

minimized. This minimization is also called minimization in

the range [39]. For a minimization of the functional

I f * * - 1 111/?

r (A A ) r / one could refer to the work done by

Hestenes and Stiefel, T. Sarkar, J. W. Daniel, T. Cwik and

Appendix [8.3].

3.3 NUMERICAL PROPERTIES OF THE CONJUGATE GRADIENT METHOD

3.3.1 Singular Operators.

*

Although the transformation Z Z appears to render

the conjugate gradient method universally applicable for the

solution of linear operator equations, one must be careful

of the condition number of Z . If Z is almost singular,

mn ran

*

Z Z will be even more ill-conditioned than Zm . For

mn mn . mn

example, let a matrix A be:

1 1

A =

.99 1

This matrix has a condition number approximately equal to

400, i.e.

A.2 / A.]. = 400

*

whereas A A has a condition number of ^2/^-1

 = 160,000. That

means that there is a strong risk of facing poor convergence

rates.

One way to check whether or not Z is nearly singular

mn

 3

is to slightly perturb the coefficients of Z and apply the
38

conjugate gradient method again. If the results of the

perturbed system are very different from those obtained from

the original system, then the matrix Z could be considered

mn

to be nearly singular and hence poor convergence rates

should be anticipated.

3.3.2 Convergence rate

J. W. Daniel, T. Sarkar and Westlake [40] have shown that

the convergence rate of the conjugate gradient method is

given by:

•* 11

2

where j is the exact solution and X and ,„ ,

*

the maximum and minimum eigenvalues of Z Z in the finite

mn mn

dimensional space in which the problem is being solved. In

this dissertation all problems are solved in a finite

dimensional space and an investigation of what happens to

the convergence rate as the dimension of the approximation n

goes to infinity is avoided.

W. J. Kammerer and M. Z. Nashed [41] have shown that

the conjugate gradient method will converge even when Z is

a singular matrix (but with poor rates as mentioned .before).

In that case, the method converges monotonically to a

solution with minimum norm and the rate of convergence is

given by Sarkar as:
39

M. 5<°> . z+mn V"

where M= ||z

2 \ J *

v

 ' + (7 7

mn' mn mn'

n

(3.26)

(3.27)

and Z is the pseudo-inverse of Z

3.3.3 Stability

As in most numerical techniques, stability problems may

appear due to roundoff errors in the calculation of the

residual and the direction vectors. One possible way of

automatically detecting "instability during the iterative

process is to look at the ratio

 a

n/

a

n_i since all scalars,

, are in the range

1 < a < 1

X~min

According to T. Sarkar, an upper bound for a /°in_i is

A. max/ \ min and hence stability may be low if \ max/ X min

is large. In our case/ computational instability with the

above algorithm for the computation of the residuals was not

encountered.

3.3.4 Global versus local convergence

As with most optimization methods, the conjugate

gradient method may end up in a local minimum instead of a

desired global minimum. If the number of unknowns is

relatively large, it is practically impossible to judge in
40

any way whether or not the minimum found is the desired

minimum. One possible procedure to use to check this is the

following: Use several initial guesses in the domain and

repeat the optimization problem. If all optimizations

result in approximately the same answer, one could be

assured that this answer is indeed .the desired global

minimum of the problem. Moreover, it should be mentioned

here that, from the engineering point of view, one is

usually not .interested in the global minimum if the solution

obtained can be considered satisfactory.

3.4 STOPPING PROCEDURES AND INITIAL GUESS

3.4.1 Stopping procedure

For the conjugate gradient method there are different

stopping procedures to terminate the iterative process.

Three of the most widely used procedures are the following:

a) ERROR = .. ' . = - j. ^ .. — -_< e (Normalized error)

II E II II E1 ||

where e is an assigned number of desired accuracy.

b) Percentage error

ERROR% =

r '

2mn

 J + E

E-ll I1

100 <6 (Normalized error)

where 6 is another assigned number.

imp ERROR = I^'H—I < C (Normalized error)
41

3.4.2 Initial guess

In all cases checked in this dissertation, the initial

guess used was J ' = J * ' =0. That gives a 1, or 100

percent, error on the first iteration if the above

normalized error measures are used. The other reason for

using a zero guess as a starting point was to see if the

method converges with the worst possible guess.

Any other initial guess could be employed to start the

algorithm. The closer the initial guess to the correct

answer the better, since the faster the method will

converge.
42

4. FORMULATION OF THE S.D.C.G. METHOD FOR THIN WIRES WITH

FINITE CONDUCTIVITY AND FOR APERTURE FIELDS

4.1 EQUIVALENT RADIUS CONCEPT AND INTERNAL INPEDANCE

The strip analysis can be used to determine the

scattering characteristics from a mesh of cylindrical wires

by employing the "equivalent radius" concept. This is

accomplished by replacing the non-circular cross section of

a metallic strip with a circular wire whose radius is the

"equivalent radius" of the non-circular cross section (See

Figure (4.1)). Butler [42] has shown that the equivalent

radius of a narrow conducting strip is one fourth of its

width i.e.

:

 >eq * a/*

where a is the equivalent radius of a cylindrical wire,

eq

and a is the the width of a thin metallic strip.

Fig. 4.1. Equivalent radius of a strip
43

For the case where the wires are of finite conductivity

the necessary boundary condition that must be satisfied is:

E

1 + "Es = 2int I (4.1)

instead of Es + E1 = 0. where I is the current in the

wires and Z. . is the internal impedance of the wire which

is given by Jordan and Balmain [43] as:

.

int

 ~-- (4.2)

where Z =[ -) is the intrinsic impedance

m * '

of the metal. Y is equal to (jumco(o + j o>e m))

and I and I, are the modified Bessel functions which can

o 1

be written in terms of infinite series as:

oo

i_(x) = y (x/2)

s=o

n

 'si (s+n)! (4.3)

n=0 for I and n=l for I.. The case of particular interest

here, occurs for frequencies sufficiently high that the

depth of penetration is small compared to the radius of the

wire. This implies that I ya |»1 and, using the

asymptotic expansion for I and I,, IO=IT (See Figure

(4.2)). Thus, the internal impedance can now be written as:
44

Z. (high frequency) sr

m

2 IT aeq

(4.4)

X

Fig. 4.2. Modified Bessel functions

From equation (4.4) it can be seen that for a small skin

depth Z. is equal to the surface impedance of a thick

metal sheet that is one meter long and 2n a meters wide,

e

Now substituting the expression for Z in equation (4.4)

yields:
45

1 /co U J co u

Z. (high frequency) = -- \l -- + ~

2KaV2 o 2 T l a

(4.5 )

This expression for Z. . can now be used in equation (4.1)/

i.e.

Z

int

 X = (Z inf A ) 5=fint 5 < 4

'

6 )

since J=I/A where A is the surface area of the wire.

This leads to:

-*e -»i =s ->

E

S

 = -E1 + Zint J (4.7)

Replacing this expression for E in equation (3.19) yields;

Zint J

 -

 Zmn J

(Z«n-fint> '« (4

-

8)

Now equation (4.8) can be solved for J using the algorithm

mentioned before in equation (3.24) and replacing Z by

(Z -z. .). Rather than form the matrix (z

mn~2jnt)

explicitly, one can carry out the calculation using the

following algorithm:
46

r

 ' E + Zmn -

 Zint J

(0)

* -*(0 r

 ') = Z* r+

mn int

ERRF = r

**

The equations for the n iteration are:

a =

n

Z* ?(n)

Z - 1 *

 mnr Zint

(n)

int

(4.9)

ERRF(n+1) = ERRF(n) .. - "int

'

 Z

int

„* ?(n-H) _ =* *(n+l) 2

Z

 mn

 r Z int r II

n zZ

* *(

 mn

 r int

_

~ ti mn

7L* £<

 intr

_

n

END OF DO LOOP

4.2 SOLUTION OF APERTURE FIELDS

To solve for the aperture fields (See Figure (4.3)),
equation (2.12) is used as the starting point, i.e.

a.

H

1

 -

cou

mn

mn mn

mn

G(a ,3

mn ran

47

(4.10)

Fig. 4.3. Sampling for the Aperture fields

For a complete derivation of equation (4.10), see Chapter 2

and Appendix 8.2. H1 is the incident magnetic field which

is a known quantity, G and a , 0 were defined before

mn mn

and they are also known. The unknown in this case is tia, so

equation (4.10) can now be cast into operator form as:
48

«*

 = Ymn P (4.11)

Ymnf like Zmn in Chapter 3, is neither positive definite nor

self adjoint, so both sides of equation (4.11) are

multiplied by the conjugate transpose of Y , i.e.

mn

Y

*mn ^ = **n,n *mn E3

mn mn mn (4.12)

Now the algorithm in equation (3.24) can be applied to the

above equation be replacing Z by Y , J by ~la, and -£* by

mn

 J "mn1

H

1

.

As before with the current densities, we choose E^°^ as

an initial guess for the aperture field, Ea, and start

iterating. In this dissertation the initial guess is chosen

to be equal to zero in all check cases and this leads to the

following stopping procedure:

H

1

 - Y

ERROR • " . " = J - V

If a percentage error is desired then the stopping procedure

becomes:

ERROR% =

H

1

 -

xlOO = " ,. ^, ..

mn

 ^x 100

•lls i ll

 (4-14)

Note that, for a zero initial guess, the first error

estimate will be equal to 1 (for the first iteration),

whereas, the second estimate will yield a 100% error.
49

4.3 REFLECTION COEFFICIENTS

The transmission and reflection coefficients are the

quantities of most important in characterizing the

properties of a mesh. In order to define those coefficients

for both polarizations, transverse electric (TE) and

transverse magnetic (TM), it is necessary to first define

the incident and scattered fields.

For TE. polarization, the incident fields are:

EX = EQ sin(-<p) ; E = EQ coscp

E coscp cosd E sinco cosd

u _ O ., O

HX . H =

n * n

where EQ is the amplitude of the incident electric field and

1/2

n=( uo/eo ) ' is the free space wave impedance.

For TM polarization, the incident fields are given by:

EX = EQ cosd coscp ; E = EQ cosd sincp

sin(cp -TI /2) E^ cos(cp -n

x ~

H = °

n

 J n

According to Wait and Hill [44], when the spacing

between adjacent wires of the mesh is less than \/2, there

is only one grating lobe and only the constant current

components JQOX and JQ contribute to the scattered field.

J

oox

 and Jooy are the zer°~m

ode current density components.

The rectangular components of the scattered field, for large

z , are given by:
50

J cos

• 2n sin O sin«p cos<p

exp ik [zcosd + sind (x coscp + y sincp)]

(4.17)

E

S

y= JQOX sin2O sincp coscp -Jooy (l.-sin2

d sin2cp)

exp ik [zcosS + sind(x coscp + y sincp)]

(4.18)

The above expressions can also be obtained from equation

(3.18) as follows: Solve for J and substitute the solution

in equation (3.18) to obtain the scattered fields. That is,

1

•H^BHM^K

jcoe

2 2

k

o -

-a &

mn mn

- a B

mn mn

k

 2 -e2

o

 p mn

G J

(4.19)

and so the reflection (amplitude) coefficient becomes:

Ry

 =

(4.20)

(4.21)

If,the total power reflection coefficient R is desired

then the following expression can be used:

Real E

s

 x Hs* . z dS

unit cellReal ] / E1 x H1*.*- z) dS 1 (4.22)

unit cell
51

where Es is the scattered field due to the induced current

densities, J, derived from equation (4.19) and, after taking

the inverse Fourier transform, HS is the scattered magnetic

field derived from ?s by making use of Maxwell's equations.

Moreover, if the total power transmission coefficient,

T, is to be computed/ one can employ the formula below:

t

a x Ha* (- z) dA (4.23)

T E1 x H1* .(- 2) dAl

*J J

Real

aperture-

"^3 :> 3 where E is the aperture Electric field and H is the

magnetic field in the aperture derived from Maxwell's

equation

-jo)U Ha =Vx Ea (4.24)

For perfectly conducting frequency selective structures it

is also true that:

Tl2 + |R|

2 = 1

this condition can be used in the perfectly conducting cases

to check the convergence and accuracy of the results.
52

5. RESULTS AND COMMENTS

5.1 ONE DIMENSIONAL CASE (INFINITE GRATING OF PARALLEL

STRIPS)

5.1.1 Current densities

The one dimensional case was studied first, and

compared to the Spectral Domain Approach with the

contraction factor [26] , since results from that method were

readily available. In Table (5.1), for example, the current

density levels, for thin strips obtained by both the S.I.T.

method and the S.D.C.G. method are in very good agreement.

This implies that either method can be employed to generate

the induced current densities on a strip or wire grating for

any incident field.

Table 5.1. Current densities obtained by the S.D.C.G.

method and the S.I.T. method. (See Figure (2.1)

for the geometry), a is the spacing between

adjacent strips and w is the width of the

strips.

Spacing(a) Width(w) S.D.C.G. S.I.T. Difference

0.55A 0.005X 0.02664928 0.02770429 0.001055

0.25A 0.005A 0.05155611 0.05183827 0.000281

0.125X 0.005X 0.07172995 0.07114100 0.000588

0.100A 0.002A. 0.07545375 0.07521373 0.000240
53

Figures (5.1) and (5.2) show that the current densities

obtained via the S.D.C.G. method for a parallel grid with

thick strip and a normally incident field behave as

expected. For the copolar component the current density

curves downwards. That is, it is larger at the edges than

0.38

Yl

13

-0.4

,-0.38

Fig. 5.1. Amplitude of y-component of the current density

for a grid of parallel strips and for a normally

incident field ( 0=0°). The incident electric

field is along the y axis.
54

at the center (Fig. 5.1), a phenomenon attributed to the

edge effects of the metal strip. On the other hand, the

cross-polar component in Figure (5.2) curves the other way

around, i.e. outwards.

3.167H-06

2.1116E-06

1.05S8E-06

O.OOOOE+00

0.

Yl

0.33

0.13

-0

-0.13

-0.23

,-0.38

Fig. 5.2. Amplitude for x-component of the current density

for a grating of parallel strips and for a

normally incident field ( d = 0°). The

polarization is TE and the incident electric

field is y directed.
55

Figure (5.3) shows how the reflection coefficient for normal

incidence increases as the spacing of the grid gets smaller.

This anticipated behavior is due to the fact that the closer

the wires, the closer the grid structure resembles a solid

metal sheet. Notice that this method is even valid for

spacings of 1/100 of a wavelength; a fact that renders this

algorithm very useful for radiometric applications where the

spacing between wires is of the order of 1/10 of a

wavelength or less. It should be mentioned here that the

data in Figure (5.3) are not compared with any measured data

or calculations made using other methods since at these

spacings neither calculations nor measured data exist.

•>t II m

32X32

SAMPLES

CELL WIDTH(A)

QOl Q02 Q03 QQ4 CD5 QO6 QO7 O.06 Q09 O.I

Fig. 5.3. Reflection coefficient for a grid of parallel

strips and spacings of 1/10 X and less.
56

The only other method that can generate reflection

coefficients at those spacings is the S.I.T. modified with

the contraction factor given by [26]. Table 5.2 shows that

the reflection coefficients obtained by both, the S.D.C.G.

method and Brand's method [26] are almost identical for

various wire spacings.

Table 5.2,

Spacing

0.125 X

0.10 X

0.06 X

0.05 X

0.02 X

0.01 X

Reflection coefficients for different wire

spacings calculated by the S.D.C.G. method and

the contraction factor-S.I.T. method. (Normal

incidence)

S.D.C.G.

0.844

0.888

0.954

0.967

0.994

0.999

S.I.T.

0.843

0.885

0.960

0.969

0.994

0.999

The S.D.C.G. algorithm for one dimensional cases (i.e.

parallel wires) converges in at most six iterations with a

normalized error of less than 0.5 percent. The CPU time

used for each of the above cases was about 20 sec on the

3081 IBM system and for a 32x32 sampling rate. This time

includes plotting time. Table 5.3 shows how the normalized

error decreases at each iteration for spacings of 1/10 of a

wavelength or less.
57

Table 5.3. Normalized error versus number of iterations for

spacings less than 1/10 of a wavelength

Spacings between strips Normalized percentage Number of

(width of strips=0.002A) error (||r||/||El ||)xlOO Iterations

100 1

0.10 x 13.57 2

0.348 3

100 1

0.09 A 17 2

0.14 3

100 1

0.07 A 29 2

0.3 3

100 1

0.05 A 52 2

0.3 3

100 1

0.04 X 69 2

0.13 3

100 1

0.03 A 69 2

0.25 3

100 1

73 2

0.01 A 40 3

17 4

0.5 5
58

5.1.2 Aperture Fields

To verify this algorithm for use in solving for the

aperture fields, a number of check cases are presented.

First, the S.D.C.G. method is checked against the S.I.T. -

contraction factor method. The results are depicted in

Figure (5.4) for a sampling rate of 32x32. The agreement

1.0-

0.9-

0.8

R

M

P

L

I 0.7

U

0

E

0.6

0

f

R

P 0.5.

E

R

T

U

R 0."»

E

F

I

E 0.3.

L

0

0.2-

0.1-

0.0-

-0.150

Fig. 5.4,

PARALLEL GRID

T.E

i.I.T [26]

S.D.C.G

CELL WIDTH

-0.075 0.000 0.075 O.lSO

Amplitude of Aperture fields for an aperture

size of 0.25 wavelengths by the S.I.T. and

S.D.C.G. methods. (Normal incidence).
ORIGINAL PAGE IS

Of POOR QUALITY

between the two methods is very good indeed. To actually

see how close the numbers are, Table 5.4 gives the values

for the aperture field at each sampling point for both

methods.

59

Table 5.4. Values of aperture field at each sampling point

for an aperture size of 0.25A. and normal

incidence.

Cell point

on x-axis

-0.129126310

-0 . 120795548

-0. 10*174063

0.3580? 320 4E

0 .?7472fel76 E

0 .791419146E

3.70 tlll5S5E

0.6248047«SE

0.453190 Ign E

0.3749f?28??E

0.29157«5'i<; 6E

0.2 082^,8240 E

0.124960914E

0.4165362';7E

0.416536257E

0. 12496091 4£

0.2082«b 6240 E

0.2915755£6E

0.374Sa2S5SE

o.^sai^oico E

0. = 414<»746 <TE

0 .624604795E

0.70611152SE

0. 791419148E

0 .374726176E

0 . 35 €03 320 <£

0. 10417 4 OP 3

J. 112464786

0 . 1207<=554£

0.129126310

-0 1

-01

-01

-01

-01

-01

01

01

01

01

01

02

02

.ll

01

01

01

OI

01

01

01

01

01

0 1

S.I.T.

method

0. 182258561 E-O'

0.384013295

0. 560089946

0.661797166

0.73r8476992

0.796682596

0.€44159245

0.83262^010

0. «14659142

0.940926552

0.562553392

0.979851842

0.993371725 .

1 .00 326443

1.00979042

1 .01 301479

1.01 301575

1 .0097<=-042

1.00 326347

O.S93371725

O.S79851342

0.962553659

0.940926552

0.314659023

0.882627010

0.£44159365

0.796632477

0,7 38 4 76 753

0.6617'76 43

0.560 0907^ 0

0.384014010

0. 182273 383E-3 1

S.D.C.G.

method

o.oooooooooc»o o

0.377430677

0.556849957

0.660453684

0. 7385-: 7 141

0.79787993 4

0 .846264482

0.865471702

0.918123«43

0.944895143

0.966925740

0.984547973

0.998 329093

1.00839133

1.01503463

1.01830196

1 .01 831055

1 .01503181

1.00839233

0.9^8 323321

0.984559417

0.966914892

0.94490 2539

0.918127d94

0.88S474324

0. €4625^52 1

0.7978B4d22

0. 7385R5949

0.660460711

0.5568441 15

0.377436161

0.000 00 OOOOE > 0 0

Figure (5.5) shows how the error is reduced at each

iteration for this case. It should be mentioned here that

this normalized error (See Section 3.4 for definition) is
60

the total error one obtains by sampling the unit cell by a

rate of 32x32 samples. Another check was obtained against

published data given by Mittra and Tsao [16]. Again, the

agreement between the two techniques is shown in Figure

(5.6).

NORMALIZED ERROR APERTURE FIELDS (TE)

0.9-

0.8-

0.7-

N

0

R 0.6-

n

A

L

Z 0.5-

D

R <Mft

0

R

0.3-

0.3-

0.1-

0.0-

2 4 6 8 10 12 H 16 18 20

NUWER OF ITERATIONS

Fig. 5.5. The normalized error for an aperture field of

0.25 A.in size
61

—@—@ s.D.C.G

Tsao and

Mittra

Fig. 5.6. Amplitude of the aperture electric field for a

unit cell with a=1.4X and b=0.6 a . The

incident field is at normal incidence and with

TE polarization

For the infinite grid of parallel wires, the current

densities, the aperture fields and the reflection

coefficients obtained by this algorithm are in very good

agreement with the S.I.T.-Contraction method and the results

published by Mittra and Tsao.
62

5.2 TWO DIMENSIONAL CASE (i.e. INFINITE GRID WITH SQUARE

OPENINGS)

5.2.1 Current densities

For the two dimensional case a number of cases are

checked against calculations by Wait and Hill [6,44,45] and

Kontorovich, Astrakham and their colleagues [9,10].

Overall, very good agreement is found in the calculation of

the reflection coefficients for different angles of

incidence, polarization and wire spacing. The reason for

comparing reflection coefficients with other methods is

simply that the reflection coefficient is the parameter of

most importance in designing wire meshes.

Figure (5.7), (5.8), (5.9) and (5.10) show these

reflection coefficients for both transverse electric (TE)

and transverse magnetic (TM) incidence. Calculations using

the S.D.C.G. method are compared with two other methods.

Wait's method is based on a Fourier series expansion

solution, whereas, the Kontorovich-Astrakham method is based

on the averaged boundary condition technique. In all those

figures, a=b represents the wire spacing of the square mesh

(See Figure 2.1) and c is the equivalent radius of the

strips. Figure (5.7) exhibits the characteristic Brewsterangle minimum for the S.D.C.G. method and Wait's method.

The discrepancy between the two curves is attributed to the

fact that in the S.D.C.G. method planar strips are actually
63

used instead of round wires. The sampling rate used in

these cases was 16x16 samples. For thin wires, this

sampling rate is good enough to obtain a good estimate for

the reflection coefficients; this is evident from these

figures. If more accuracy is desired, the number of samples

can be increased. In Figure (5.10) one can see that, by

increasing the sampling rate, a slightly better estimate can

be obtained.

1.O

y

LL

IL

U

O 0.4

O

UJ

LU

or

0.2

0.0

\ T.M

\

0=70

O.I

-x-x-x- WAIT [6]

-0-0-0- SD.C.G

O.4 Q5

Fig. 5.7. Reflection coefficient for TM incidence and

various spacings for d=70 deg., and cp=0 deg
64

l.O

UJ

0 08

LU

80.6

z

o

O 0.4

LU

LU

or

0-2

O-O

IE.

S.D.C.G

WAIT'S

METHOD

o-i 0.2 0.3 0.4 0-5 a

Fig. 5.8. Reflection coefficient for TE polarization and

an oblique incidence with d=70 and ro=0 degrees,
65

LU

o

U- O.5

UL

LU

O

<-> Q5

"ui; 0.4

'(J

U

_

LU 0.2

or

c.o

R ee TM

S.D.C.G

--* KONTOROVICH

ASTRAKHAN/I.

9O 80° 7O 50 5O° 40 3O" 2<f IO° 6

9 ^

Fig. 5.9. The reflection coefficient for a spacing of

a=0.25Xand for different values of the angle of

incidence theta. (TM polarization)
66

I1.0

LI

UJ OS

O

o

0.6

O

0.2

0.0

TE

SAMPLING

16 XI6.—

32 X32

x—*-- ASTRAKHA^,.

KONTOROVICH

<3/C=50

90' so0

70° 60° 50° 40" 30° 20° 10° 0°

Fig. 5.10. The reflection coefficient for TE polarization

and different angles of theta. The wire spacing

is equal to 0.25X .
67

Figures (5.11) and (5.12) confirm the expected result

that the wider the wires the larger the reflection

coefficients. This result should be anticipated since a

mesh with wide wires is a better approximation to a

continuous metal sheet.

LLJ I.C

O

U_

u_

LU 0.8

O

O

o °'6

ho

LLJ Oi4

LU

or

0.2

cxo I

R99

Wl DTH

©—©—0- O.I 9^

A-A.-A 0.05

90° 80° 70 60° 50° 3O° 2O° I

- ^ 6

0°

Fig. 5.11. Reflection coefficients for different

thicknesses and for different angles of

incidence theta. The polarization is transverse

magnetic and the mesh opening is a=b=0.25X .
68

•

z

j

1-0

08

U_

U_

(J 0.5

Z

o

jl 0.4

o

UJ

_)

U_ 02

LU ^

or

0.0

=0 -brQ.25

TE

A--A--A O.O5A

O.OO5

W

I

D

T

H

80 7O 50° 50 40 30° 20° 10

e

Fig. 5.12. Reflection coefficients for different widths and

for different angles of incidence theta. (TE

polarization).

Figures (5.13) and (5.14) depict the change in the

reflection coefficient when the wire mesh consists of wires

with finite conductivity. The figures confirm the fact that

the reflection coefficient of a lossy wire-mesh is less than

that of the perfectly conducting wires case. The reason for

this difference is that, for perfectly conducting wires

( 0=00), the reflection and transmission coefficients are

governed by the relation:
|T|2 + |R|2

 =1

69

where T is the transmission coefficient and R is the

reflection coefficient. For lossy wires, due to the loss of

energy in the wires, T 2 + R 2 * 1.

1.0 -

z

UJ Q3

y

u_

u_

i I 1 -"5

O

O

g

H

O

O.4

0.2

u_

LU

TM

-=aas

0—0—0 j m

90' 60° 70° 60° 40' 30° 20° 10°

Fig. 5.13. The reflection coefficient for TM polarization

and for both cases, a perfectly conducting wire

mesh and a lossy wire mesh.
70

1.0*

O.8

0.6

0.4

0.2

.O

z

TE AND 1 M a-—0J-0--0 FINITE CONDUCTIVITY

NORMAL INCIDENCE

u.

UJ

cr

10° ' I01 I02 I03 CT

Fig. 5.14. The reflection coefficient for different

conductivities.

So far, we have only discussed the reflection coefficients

for different polarizations, angles of incidence, widths and

wire spacings. We have also compared them with whatever

data were available to us [6,9,10]. In the figures to

follow, the current densities are presented and analyzed for

different cases of interest. First, we start with Figures

(5.15) and (5.16) where the current densities, J and J ,

x y

are depicted. The spacing used in that case was 0.25 A,

wavelengths and a thickness of 0.005 A. . The sampling rate
71

was 32x32 and the wave was normally incident for a TE

polarization. For wider strips Figures (5.17) and (5.18)

show how the current densities behave. And for a case with

lossy wires Figures (5.19) and (5.20) give the results. In

all those cases, the square-shaped unit cell was used.

T.E.

PHI=0

0.0037

0.0024

0.0012 •

0.0000

Yl

WIDTH= 0.0052

32x32 SAMPLES

0.13

O.Ot

-0.13

-0.13

Fig. 5.15. Amplitude of the x-component of the current

density for a normally incident wave on a square

mesh with a spacing of 0.25 A. between adjacent

wires. (Incident E field is along the y axis).
72

9.026

a. on

0.009 .

0.000

0.

K 32x32 SAMPLES

V

u

0.13

-o.u

Fig. 5.16. Amplitude of the y component of the current

density for a normally incident wave on a square

mesh with a spacing of 0.25 X . (E incident is

along the y axis).
73

32X32 TE,TH=70,PHI=0 (AA=.355 BB = .25-j

0.0033

0.0022

0.0011

0 4000

0

-0.

O.Ofr

-0.06

0.15

-0/13

,-0.18

Pig. 5.17. Amplitude of x component of the current density

for an obliquely incident wave (O=70°) on a

square mesh with wide metal strips

(width= .105 A). TE polarization with E

incident along the y axis).
74

32X32 TE,TH=70,PHI=0

0.0043

0.0029

•0.0014 •

0.0000 •

0.

Yl

0.18

0.06

XI

-0.06

-0.13

.-0.18

Fig. 5.18. Amplitude of the y component of the current

density for a wave incident at an angle theta=70

on a square mesh of wide metal strips

(width=0.105 A). (E incident is along the y

axis).
75

SIGMA=500,TE,TH=0,PHI=0 (A=.25 0.005 THICK)

0.0027

0.0018

0.0009

0.0000

ri

0.0-»

-0.0»

-0.12

,-0.13

Fig. 5.19. Amplitude of x component of the current density

for a square grid of thin strips but with a

conductivity of o=500U/m . A normally incident

field (0=0°) and a sampling rate of 16x16

samples are used.
76

SIGMA=500,TE,TH

=0,PHI=0 (A=.25 0.005 THICK)

0.015

0.010

O.OOS

0.000

Yl

0.13

-0.12

-0. 13

pig . 5.20.

Amplitude of the y component of the current

density for a square mesh with wide metal

strips. The conductivity is o=500 O/tn and the

width O.IOSX'TE polarization with E incident

along the y axis.
77

It should be mentioned here that the magnitude of the

current densities becomes smaller as the conductivity of the

metal strips or-wires is reduced. This result should be

anticipated since the smaller the conductivity of the wires

the lossier they are.

Now, to illustrate the significant effects that occur

at a bonded junction, the cross-shaped unit cell is used.

The current densities obtained in this case are depicted in

Figures (5.21) to (5.26) for different spacings, widths and

angles of incidence. It can be seen from all figures that

this method predicts the step discontinuity at the bonded

junction. It should be stressed here that in this

dissertation only the bonded case is treated; that is, the

case where a perfect contact between the wires exists at

each junction.

Since the existing mesh surfaces resemble more closely

the bonded mesh, than the unbonded case, a study of the

unbonded mesh was not done here. Quite often, in practice,

the wires are soldered at the bonds to obtain a perfect

contact. The study of the unbonded mesh is of interest,

though, because of its physical analogy with a thin

magnetized plasma. Anisotropic unbonded wire mesh can be

used to simulate a thin magnetized plasma sheet. Wait has

calculated the reflection and transmission coefficients for

the unbonded mesh case [6] .
78

Yl

0.21

-0.21

0.63

,-0.63

Fig. 5.21. Amplitude of y current density component for a

normally incident wave with TE polarization on a

square mesh (a=b=1.25A. and width =0.02X ). (E

incident is along the y axis).
\ORIGJWAL PAGE IS

"Sit. poOR QUALITY

79

TE,TH=00,PHI=0 (A=B= .022j THICK)

JJxl

0.0042 •

0.0028

Yl

0.63

0.21

,-0.63

Fig. 5.22. Amplitude of x current density component for a

normally incident wave on a square mesh.
80

TE,TH=0,PHI=0 (A=.25p| 0.005j) THICK; 32X32

0.025 .

0.017

Yt

0,13

Fig. 5.23. Amplitude of y current density component for a

normally incident wave on a square mesh. S=0°

and E incident is along the x axis.
TE,TH=0,PHI=0 (A=.25^ 0.005^THICK) 32X32

81

0.0015

0.0000 .

fl.f

Yl

0.13

0.0*

-0.04

-0.13

Fig. 5.24. Amplitude of x current density component (crosspolar) for a normally incident wave on a square

mesh with thin strips.
TE,TH=60,PHI=0 (A=B=.25>AND .02iTHICK)

82

9.003

9.000 •

0.

11

-0.04

0.13

-4.13

,-0.13

Fig. 5.25. Amplitude of y current density component for an

obliquely incident wave on a .square mesh with

strips of width equal to .02 x .
83

TE,TH=60,PHI=0 (A=B=.25j|AND .02j-THICK)

Yl

0.13

O.M

Fig. 5.26. Amplitude of x current density component (crosspolar) for a square mesh with an obliquely

incident wave.

Figures (5.27) to (5.30) show how the normalized error

is reduced at each iteration. It can be seen from all these
84

figures that the residual error decreases monotonically.

From Figure (5.28), one can see that the closer the strip

spacing, the longer it takes to converge to a specified

normalized error. The difference in the normalized error

between the 0.25X and 0.75A. spacings is indeed large,

whereas the corresponding difference in the normalized error

between the spacings of 0.75X and 1.25X is not that drastic.

1.0-

0.9-

0.8-

0.7-

N

0

R 0.5-

H

A

L

z o.s-!

0

E

R 0.4-

R

0

R

0.3-

o.a0.1-

0.0-

0 .10 20 30 40 SO 60 70 60 90 100 110

NUWER OF ITCRflTHWS

Fig. 5.27. Normalized error for currents for a square mesh

with a=b=1.25X , theta=30° and phi=0°.
85

NORMAL INCIDENCE (TE CASE)

0.9-

0.9-

0.7-

N

0

R 0.6

n

A

V.

I

z o.s

0

E

R 0.*

R

0

R

0.3-

0.2-

0.1-

0.0-

•«• ft-8-i.iS LflWfl.

•«•••••••• fl-9-0.

XXXXXXXXXX fl-8-O.iS LB?

0 5 10 tS 20 25 30 35

MMER OF nownoie

*5 SO

Fig. 5.28. Normalized error for the current densities for

different wire spacings.

Figure (5.29) enables us to observe that the error rate

depends not only on the wire spacing, but also on the angle
86

of incidence. The normalized error for an angle theta = 70

degrees decreases much faster than that for an angle theta =

30°, or phi = 0°. The reason for that is based on the fact

that the eigenvalues of the matrix Z change as those

mn

angles change.

1.0-4

0.9-

0.8-

0..7.

N

0

R 0.6-

«

A

I

I O.SQ

E

R 0.*-

R

0

R

0.3-

0.2-

0.1-

0.0-

0.35 FOR VARIOUS ANGLES (TE CASE)

3EWEFS (11

T>,ETS.30.0 OEGREES (3)

XXXXXXXXXX TH£Tfi.70 0 OECSEES (3J

(3)

10 20 30 40 SO 60 70 80

MMER OF rrcnrnam

90 100

Fig. 5.29. Normalized error for the same square mesh but

with different angles of incidence.
87

To see how this occurs, we recall the expressions for

M* ~ "1

a , 3 and G = •======. which are the elements

mn' ^mn ,, /, 2 2 02

that form the entries of matrix Z'm . These elements are

functions of angles theta (d) and phi ( cp ) and of the wire

spacing. This means that any change in theta, phi or the

spacing will yield a change in the matrix Z , and hence,

the eigenvalues of the matrix will be different. It was

mentioned before, in Chapter 3, that the rate of convergence

depends on the eigenvalues of matrix Z . Therefore, any

change in theta, phi or in wire spacing will change the rate

of convergence.

Another interesting phenomenon is observed in Figure

(5.30) where the normalized error for the same wire spacing

and the same incident field, but for differently shaped unit

cells is plotted. From that figure it is clear that the

normalized error for the cross-shaped unit cell decreases

much faster than that of the square-shaped unit cell.

Although both unit cells generate the same currents and

reflection coefficients, the cross-shaped unit cell can be

more advantageous as far as computing time is concerned.

One reason for this difference between the two unit cells

lies in the fact-that in the cross-shaped unit cell the wire

strips appear to be wider to the algorithm than the

corresponding strips in the other unit cell, as Figure

(5.31) illustrates.
88

i.OH

0.9-

0.9-

0.7-

N

0

ft 0.6.

M

fl

L

0.5-

0.3-

0.2-

0.1-

9.0-

NORMAL INCIDENCE ,4=3=0.25

-ShflPEO (1)

0 10 20 30 ^0 50 60

NUWES OF

70 30 90 200

Fig. 5.30. Normalized error for two differently shaped unit

cells.
(a) SQUARE

89

*-.

r

(b) CROSS

Fig. 5.31. The two unit cells used in Figure (5.30)
90

Before we present some results for the aperture fields,

the problem of evaluating the reflection coefficients from a

frequency selective planar surface, shown in Figure (5.32),

is discussed. Table (5,5) gives the results for the

reflection coefficient evaluated by this algorithm for

different values of Q. Here, it is observed that a crosspolar component arises even for a normally incident wave.

This result is very important in assessing the degree of

depolarization from such a planar structure.

[< BB

Fig. 5.32. A different frequency selective surface

geometry.

It should be mentioned here that this configuration
91

offers a better approximation to the knitted mesh than the

infinite square grid. The reason for this is that the

periodicity of the above planar structure resembles that of

the knitted mesh.

Table 5.5. Reflection coefficients from the arrangement in

Figure (5.32). (Normal incidence and TM

polarization) .

Q copolar cross-polar

90° 0.721 0.0

80° 0.725 0.086

70° 0.7263 0.16

60° " 0.7316 0.238

50° 0.7523 0.286

40° 0.77157 0.292

5.2.2. Aperture Fields

Figure (5.33) shows how the aperture field is compared

with the results published by Tsao and Mittra [16] . This

happens to be the only available data for aperture fields

that we can compare with our calculations. For spacings

larger than one wavelength there are more than one

propagating modes (i.e. whenever k > amn

 + 3^,n which

appear as lobes in the aperture field. Notice the four

lobes in Figure (5.33) for a spacing of four wavelengths

between the adjacent strips. Figures (5.34) and (5.35)

depict the amplitudes of the x- and y- components of the
92

aperture electric field for a normally incident field on a

square mesh with the dimensions a=b=1.25A. . Note again,

that this algorithm can predict the two propagating lobes

and the edge effects on the strips that are perpendicular to

the y-directed incident electric field. Moreover, Figure

(5.36) and (5.37) give the amplitudes of the aperture fields

for a different type of polarization (Transverse magnetic or

TM) and a mesh with innerspacing given by a=b=0.25X . In

this case the angles theta and phi are both equal to 30

degrees. In Figures (5.38) and (5.39), a smaller spacing is

used (a=b=0.125 ) and an angle of incidence equal to 30°,

to see if the algorithm can still converge under the

conditions of oblique incidence and smaller spacings.

Finally, in Figures (5.40) and (5.41) a sampling rate of

16x16 samples is used instead of 32x32. The x and y

components of the electric field in the aperture are shown.

In this case the wave is normally incident on a mesh of thin

strips and s spacing equal to 0.25 A.. Once more, the edge

effects become very evident.
93

1.48

1.28 -

I.O3 -

Q63-

0.28

S 0—C I SAC AND M IT TWA

v—x—x SJD.C.&

6=0°

TE CASE

-X — X— X-X — X— Xi, \ A A =o.c

» \ 66 = 4

x— x—x-x-x—x-x

-3.5 -2.3 -2.1 -1.4 -0.7 O.7 1.4 2.1 2.3

-x/a.

Fig. 5.33. Amplitude of y aperture field along the x axis

(All four lobes are predicted).
PAGE fS

POOR QUALITY

94

APERTURE FIELD TH=PH=0 TE A=

Yl

X o . 63

-0.6

-0.63

Fig. 5.34. Amplitude of y component of the aperture field

for a normally incident wave on a square mesh

with a=b=1.25A. (Two lobes).
APERTURE FIELD TH=PH=0 TE (A=1.25j>)

95

0.22

9. IS

0.07

0.00

Yl

0.63

0.21

Fig. 5.35. Amplitude of x component of the aperture field

for a normally incident wave with TE

polarization.
96

APERTURE FIELD TH=PH=30 TM (A=0.25j\)

Yl

0.13

-0.0*

-O.i2

,-Q.t3

Fig. 5.36. Amplitude of x component of the aperture field

for a wave incident on the square mesh at angles

theta=30°, phi=30° and with a TM polarization.
97

APERTURE FIELD TH=PH=30 TM (A=0.25j)

Yl

X 0.13

0.04

XI

•0.0«

-0.12

,-0.13

Fig. 5.37. Amplitude of y component of the aperture field

for an obliquely incident wave on a square mesh

with a=b=0.25X .
98

APERTURE FIELD TH^30 PH=0 TM A=

n

0.06

Fig. 5.38. Amplitude of x component of the aperture field

for a wave incident on a thin strip mesh with
ORIGINAL

OF POOR QUALITY 99

APERTURE FIELD TH=30 PH=0 TM (A=0.125-0

^i

Yl

0.06

-0.06

Fig. 5.39. Amplitude of y component of the aperture field

for a wave incident at angle theta=30° on a

square mesh with a=b=0.125 X.
100

TE.NORMAL (A=B = .25^AND .OOSjJHICK)

Yl

-0.

9.13

0.0<»

XI

-0.0"

-0.12

,-o.u

Fig. 5.40. Amplitude of y component of the electric

aperture field for a normally incident wave on a

mesh and sampling rate of 16x16 samples.
101

TE.NORMAL (A=B = .25 AND .005 THICK)

(Ml

0.28

0.1*

0.00

Yl

0.-13

-0.12

-O.U

Fig. 5.41. Amplitude of x component of the electric

aperture field for a sampling rate of 16x16

samples and a normally incident wave.
102

5.3 CPU TIME AND STORAGE REQUIRMENTS

In general the two dimensional problem takes longer to

converge than the one dimensional case. One of the main

reasons for that is the fact-that a two dimensional FFT is

used and more sampling points are required in the two

dimensional case. Moreover, in the two dimensional case,

one has to solve for far more unknowns than in the one

dimensional problem. This number of unknowns also affects

the storage requirements for the two dimensional problem.

In general the CPU time depends on the sampling rate more

than on anything else. For example, for a sampling rate of

16x16 samples, it takes anywhere from 30 seconds to 1.40

minutes of CPU time (on an IBM 3081 system) to converge to a

reasonably accurate result. This time includes sorting and

plotting of data. For a 32x32 sampling rate the CPU time

is, as expected, much more. In fact, in this case the range

is somewhere between 1:32 and 8:00 minutes. As mentioned

before, the CPU time also depends on the angle of incidence

and the strip spacing.

The program size is 56,064 bytes for a 16x16 sampling

rate and 149,192 bytes for a 32x32 sampling rate.
103

6. COMMENTS AND SUGGESTIONS FOR FUTURE RESEARCH

Here, a number of recommendations for future research

related to the mesh problem and the S.D.C.G. method are

mentioned.

One). Skew-Symmetric Configuration for a mesh

(a) (b)

Fig. 6.1. Different sampling patterns (a) rectangular

(b) Skew-symmetric.

For a rectangular or a square grid, the number of

samples, which is also the number of Floquet modes,

corresponds to the number of couplings being taken into

account. Moreover, for rectangular sampling, the Fast

Fourier Transform can be used directly. On the other hand,
104

for non-rectangular sampling, such as that shown in Figure

6.1 (b), the number of samples may not correspond exactly to

the number of couplings being taken into account. That will

yield some erroneous results. Moreover, in this case FFT

can not be used and hence a Discrete Fourier Transform has

to be employed instead. So the actual problem here is to

modify the existing algorithm in order to represent the

strips and their width as accurately as possible. The

reason for doing this is to avoid any aliasing problems that

may arise from sampling such a configuration. Solving this

problem is important because a study of the reflection

coefficients as a function of the angle V will give new

insight in designing mesh structure that are skew-symmetric.

Moreover, this configuration might offer a better

approximation to the actual woven structure than the

rectangular mesh.

Two). Double screen

The scattering properties of such a structure are of

interest because a double screen can be used as a filter in

microwave applications. To solve this problem, the original

structure is divided into two substructures, and the

principles of equivalence and superposition are used to

obtain the formulation in the spectral domain. According to

Tsao and Mittra the problem of the double screen can be

represented as in Figure (6.2):
105

a.(=/ t=a I

/

r

SCREEN

->x

Side View

2 (1)

\

the

(2)

\^

\ Perfect Electric Perfect Magnetic

Conductor Conductor

Fig. 6.2. Equivalent problem for a double screen

The equations for the aperture field, for example, are

similar to those used in our work for the single square

mesh, but with the phase difference between the two meshes

taken into consideration. Tsao and Mittra give the

following equations in the spectral domain for E , and E ~

2 j

co u

mn

a 0

mn mn

~

 a mn

G.E

for the electric conductor
106

and

—co U / V .

mn

'2

2 2

mn

 0 mn

 ko ~

 a

 mn

mn "-mn

G.E exp[j(

for the magnetic conductor

where C, =

1-e

Y mnd

and C2 =

1+e

Ymnd

where

Ymn

mn

2 + a 2

 -k

 2

mn mn o

The total aperture field at any point is the superposition

•*= -*a

of E^ and Ea

2.

Three). A mesh over a ground plane.

//'//

MESH

GROUND PLANE

Fig. 6.3. A mesh over a ground plane.
107

In this problem, the evanenscent field from the mesh

interacts with the adjacent ground plane, and hence, the

total scattered field is different than that of the mesh in

free space. To solve this problem the superposition

principle should be used as follows. First, solve for the

field reflected by the ground plane (or dielectric plane).

Second, find the scattered field due to the mesh and finally

find the total scattered field. Once the total scattered

fields are known the reflection coefficients can be

determined.

Four). Determination of the Electromagnetic properties of a

mesh with wires made of different alloys for

Radiometric operation

Figure (6.4) shows that the wires used in constructing

the actual mesh are made of molybdenum 1.2 mill in diameter

and there is 4-6% gold, by weight, plated over the

molybdenum. In thickness this corresponds to 8-11 microinches.

'Gold

Fig. 6.4. Gold plated wire substrate
108

The resistivity ( p) of gold and molybdenum is

different. For gold p is 2.35 u^ -cm, whereas for

molybdenum the resistivity is 5.2uQ-cm. In our work we

assumed the mesh was constructed with wires of one type of

resistivity and not alloys. If the skin depth, for a

certain frequency, is larger than the thickness of the fold

that covers the molybdenum wire, the field would penetrate

into the molybdenum region. This means that the wire cannot

be considered as uniform any more. Therefore, the current

algorithm has to be modified to take into consideration this

difference in resistivity. One way to do that would be to

derive an impedance expression for thin wires made of any

kind of alloys. Once this impedance is obtained, it can be

used in the S.D.C.G. method as follows: Start with the

equation for the currents in chapter 3. That is,

Zmn J = Es (6.1)

The new boundary condition becomes:

S

S

 + i^ =

 Aalloy 3 (6.2)

where Z ,, is the internal impedance of the alloy divided

by the area of the strip or the equivalent cylindrical wire.

Now substituting equation (6.2) into equation (6.2) yields:

Zmn * -1

 + Zalloy

 (6'3)

(Zmn "'alloy' = (6'4)
109

The conjugate gradient method can be employed next to solve

for the currents as described in Chapter 4.

Another objective of the study here is to relate the

reflection coefficient evaluated by this method to the

emmissivity and reflectivity measurements carried out by

NASA, at the Langley Research Center.

In conclusion, two techniques were developed here with

an eye toward more efficient numerical computation for

grating and mesh scattering. The first method, the Spectral

Iteration Approach is applied to regions where the spacing

between the wires is not less than two wavelengths. The

second method, the Spectral Domain Conjugate Gradient

Method, can be used for any spacing. Both techniques were

solved in the Spectral Domain and both follow from one basic

formulation. A study of the electromagnetic properties such

as reflection coefficients, induced currents and aperture

fields were presented and compared with data calculated by

other methods to support the validity of the algorithm.

A number of suggestions for possible extensions of the

current algorithm to solve the problems of skew-symmetric

structures, double screens, wires made of alloys with

different resistivities and a mesh above ground were

mentioned. The code used in the Fortran program and a

listings of all main and utility subroutines appear in the

Appendices.
110

7. LIST OF REFERENCES

1) H. Lamb, "On the Diffraction in Transmission of

Electric Waves by a Metallic Grating," Proc. London

Math, Sbc, vol. 29, pp. 523-544, 1898.

2) Lord Rayleigh, "On the Passage of Waves Through

Apertures in Plane Screens and Allied Problems," Phil.

Mag. 43, 259, 1897.

3) G. G. MacFarlane, "Surface Impedance of an Infinite

Parallel Wire Grid at Oblique Angles of Incidence," J.

Inst. Elec. Engrs., pt. 3A, vol. 93, pp. 1523-1527,

1946.

4) R. J. Primish, "Some Electromagnetic Transmission and

Reflection Properties of a Strip Grating," IRE Trans,

on Antennas and Propagat., vol. AP-5, pp. 176-182,

April 1957.

5) T. Larsen, "A Survey of the Theory of Wire Grids,"

I.R.E. Transactions on Microwave Theory and Techniques,

pp. 191-201, May 1962.

6) Wait, J. R., Theories of Scattering from Wire Grid and

Mesh Structures Electromagnetic Scattering, Academic

Press, Inc., 1978.

7) R. F. Harrington, Field Computation by Moment Methods,

MacMillan, New York, 1968.

8) Stutzman, L. W. and G. A. Thiele, Antenna Theory and

Design, John Wiley & Sons, New York, 1981.

9) Astrakham, M. I., "Reflecting and Screening Properties

of Plane Wire Grids," Telecommunications and Radio

Engineering, vol. 23, pt. 2, no. 1, January 1968, pp.

76-83.

10) Kontorovich, M. I., "Averaged Boundary Conditions at

the Surface of a Grating With a Square Mesh," Radio

Engineering and Electronics Physics, vol. 8, no. 9,

September 1963, pp. 1446-1454.

11) Croswell, W. F., "The Suitability of Mesh Membrane

Material for Radiometer Reflector Applications," Large

Space Systems Technology - 1981, NASA CP-2215, Part 2,

1982, pp. 621-629,
Ill

12) N. N. Bojarski, "K-Space Formulation of the

Electromagnetic Scattering Problem," September 1972,

Air Force Avionics Lab, Contract F33615-71-C-1567,

Report AFAL-TR-72-271.

13) N. N. Borjarski, "The K-Space Formulation of the

Scattering Problem in the Time Domain," Journal of

Acoustical Society of America, vol. 72(2), August 1982.

14) Krueger, C.H. Jr., "K-Space Formulation of the Twodimensional Electromagnetic Scattering Problem,

Doctoral Dissertation, Ohio State University, 1972.

15) Mittra, R., Ko, W. L. and Rahmat-Samii, Y., "Transform

Approach to Electromagnetic Scattering," Proc. of the

IEEE, vol. 67, no. 11, November 1979, pp. 1486-1503.

16) Tsao, C. H. and Mittra, R. "A Spectral-Iteration

Approach for Analyzing Scattering From frequency

Selective Surfaces," IEEE Trans. Ant. & Prop., vol. AP30, no. 2, March 1983, pp. 303-308.

17) C. C. Chen, "Transmission Through a Conducting Screen

Perforated Periodically with Apertures," Microwave

Theory Tech., vol. MT, 18, no. 9, pp. 627-632, 1970.

18) Brand, J. C. and J. F. Kauffman, "Development of

Improved Analytical Models for Mesh Reflector

Surfaces," Large Space Systems Technology, NASA, part

2, 1982, pp. 867-876.

19) Van den Berg, P. M., "Iterative Computational

Techniques in Scattering Based upon the Integrated

Square Error Criterion," Proc. Int'l. URSI-Symposium,

Spain, pp. 97-100, 1983.

20) Sarkar, T. P., and Rao, S. M., "The Application of the

Conjugate Gradient Method for the Solution of

Electromagnetic Scattering from Arbitrarily Oriented

Wire Antennas," IEEE Trans. Ant, & Prop., vol. AP-32,

no. 4, April 1984.

21) Montgomery, J. P. and Davey, K. R., "Conjugate Gradient

Solution of Arbitrary Planar Frequency Selective

Surfaces," Int'l URSI-Symposium, Boston, USA, pp. 155,

1984.

22) Cwik, T. and Mittra, Raj, "Spectral Domain Solution of

Scattering from Periodic Surfaces Using the FFT," IEEE

Antennas and Prop. Symposium, Boston, USA, pp. 913-916.
112

23) Collin, R. E. Field Theory of Guided Waves, McGraw-Hill

Book Company, New York, 1960, Chapter 12 and p. 192.

24) Amitay, N. and Galindo, V., "The Analysis of Circular

Waveguide Phased Arrays," Bell Syst. Tech. J. vol. 47,

pp. 1903-1932, November 1968.

25) Bracewell, R. C., The Fourier Transform and Its

Applications, McGraw-Hill, New York, 1978, 2nd ed.

26) Brand, J. C., "On the Convergence of an Iterative

Formulation of the Electromagnetic Scattering from an

Infinite Grating of Thin Wires," North Carolina State

University, Department of Electrical Engineering, PHD

Dissertation, 1984.

27) Karamadian, S.r Fixed Point-Algorithms and

Applications, Academic Press, Inc., New York, 1977.

28) Istratescu, I. V., Fixed Point Theory, D. Reidel

Publishing Co., Holland, 1981.

29) Browder, F. E. and Petryshyn, W. V., "Construction of

Fixed Points of Non-Linear Mapping in Hilbert Space,"

J. Math. Anal. Appl., vol. 20, pp. 197-228, 1967.

30) Rhoades, B. E., "Comments on Two Fixed Point Iteration

Methods," J. Math. Anal. Appl. 56, 1976, pp. 741-750.

31) Rus, I. A'., "Approximation of Fixed Points of

Generalized Contraction Mapping," from Topics In

Numerical Analysis II, edited by John J. Miller,

Academic Press, New York, 1974, pp. 157-161.

32) Conte, S. D. and C de Boor, Elementary Numerical

Analysis, McGraw-Hill Co., New York, 1975.

33) Isaackson, E. and H. B. Keller, Analysis of Numerical

Methods, John Willey & Sons, Inc., New York, 1966.

34) Hildrebrad, F. B., Introduction to Numerical Analysis,

McGraw-Hill, 2nd ed., New York, 1974.

35) Kopchenova, N. V. and I. A. Maron, Worked Examples and

Problems with Elements of Theory, MIR Publishers,

Moscow, 1981.

36) Sarkar, T. K. and S. M. Rao, "An Iterative Method for

Solving Electrostatic Problems," IEEE Trans. Antennas

Propag., vol. AP-30, pp. 611-617, July, 1982.
113

37) Hestenes, M. R. and E. Stiefel, "Methods of Conjugate

Gradients for Solving Linear Systems," J. Res. Nat.-

Bur. Standard, vol. 49, pp. 409-436, December, 1952.

38) Jennings, A., Matrix Computation for Engineers and

Scientists, Wiley Publishing Co., New york, 1977.

39) Daniel, J. W., The Approximate Minimization of

Functionals, Prentice-Hall, New Jersey, 1971, pp. 114-

152.

40) Westlake, L., A Handbook of Numerical Matrix Inversion

and Solution of Linear Equations, Wiley, New York,

1968.

41) Kammerer, W. J. and M. Z. Washed, "On the Convergence

of the Conjugate Gradient Method for Singular Linear

Operator Equations," SIAM J. Numer. Anal., vol. 9, no.

1, pp. 165-181, 1972.

42) Butler, C. M., "The Equivalent radius of a Narrow

Conducting Strip,".IEEE-Trans, on Antennas and Propag.,

vol. AP-30, no. 4, July 1982.

43) Jordan, E. C. and K. D. Balmain, Electromagnetic Waves

and Radiative Systems, 2nd ed., Prentice-Hall, Inc.,

Englewood Cliffs, N. J., 1966.

44) Hill, D. A. and J. R. Wait, Can. J. Phys. 52, 1974, pp.

227-237.

45) Wait, J. R., Appl. Sci. Res. B4, 1954, pp. 393-400.

46) Harrington, R. F., Time Harmonic Electromagnetic

Fields, McGraw-Hill, New York, 1961.
114

8. APPENDICES

8.1 DERIVATION OF THE EQUATION FOR H AS A FUNCTION OF THE

ELECTRIC VECTOR POTENTIAL F

This appendix is to derive equation (2.9) from equation

(2.7) in Chapter 2.

Start with the equation:

E = Vx "* (8.1.1)

Substituting equation (8.1.1) into the following Maxwell's

equation:

x H" = j ooe IT . - (8.1.2)

yields:

yx H = -j coVx F ' (8.1.3)

Equation (8.1.3) can be written as:

(H + j 00?) = 0 (8.1.4)

Now using the vector identity

V

x (~V*m) =

 ° (8.1.5)

equation (8.1.4) becomes:

H" = -j co"F -Vv* (8.1.6)

 m

where * is the magnetic scalar potential. Taking the

curl of equation (8.1.1) leads to:

V* "E = -~ Vx V x F (8.1.7)

which can be written as:

1 -* 9 -*

' E = — —•• • f \ / \ / F — \7 Fl (818 )

o
115

by making use of the vector identity

VxVxF=V<V'F)-y 2F (8.1.9)

To completely specify the vector F, its divergence and

its curl must be defined. In equation (8.1.1) the curl of F

was defined. Now, one is at liberty to define the

divergence of F, which is independent of its curl. The

choice of ^.F is made to simplify equation (8.1.8) which is

achieved by letting:

V- F " ~J weu « m (8.1.10)

which gives:

* = -- -'F (8.1.11)

Substituting equation (8.1.11) into equation (8.1.6) leads

to:

V •

or

-*• "* • t~~7f 7 **

H = -jco F + - VV-

 F

jcoue

which is the same as equation (2.9) in Chapter 2.
116

8.2 DERIVATION OF THE EQUATION FOR THE SCATTERED MAGNETIC

FIELD HS

The purpose of this appendix is to derive equation

(2.12-b) from (2.12-a) in Chapter 2.

Start from equation (2.12-a), i.e.

i V"

H s+_ . \ — L

w

 2

K

_ «2 _n R

o

 amn amn p mn

2 2

a

 mn P mn o ^ mn

mn

exp[j(amnx+3mny)]

(8.2.1)

Figure (8.2.1) below shows how the equivalence theorem could

be utilized to transform the free standing inductive surface

into a perfect electrical conductor[16] .

t

:

ape .|cond J

for 2)0

I ape I con]

perfect electrical

conductor

t

image of

E

1 due to

electr. cond.

for

Fig. 8.2.1. Equivalent problem for an inductive FSS

structure.
117

For the region z>0 the total H field (H

tot)

 at z=0 can be

expressed as:

— • -* i r> f*

H = H (x,y) + Hi n c (8.2.2)

Now the magnetic current related to the aperture field Ea is

given by:

A M = E x n (8.2.3)

A A A where n is the normal to the aperture. For z > 0 n = z

, ,- A A

a n d fo r z< 0 n=-z . S o

A A ,a A a^ i ,

A

. M x + Mv y = [E% x -i- Ea

v y] x (z)

x y A y

= - Ea

v y -I- Eav x (8.2.4 )

x

 Y

Similarly for z < 0

M = Ea x (-z) " (8.2.5)

and the total H field is given by:

=HS - (8.2.6)

At

or HS" = H3"'' -i- (8.2.7)

Moreover -H = H (8.2.8)

So equation (8.2.7) becomes:

-2 HS+ = H.1"

0

 (8.2.9)
118

From the previous figures we had 2 M (Due to image theory).

So multiplying equation (8.2.1) by a factor of two yields:

JOJUk

 2- -n B

o mn mn

 pmn

2 2

mn mn

M exp[j<annx+enny)]

or

t

~H

 ir

But

2

jo>u

k -o_ '-a B o run mn mn

?

-a B k -B

mn mn o mn

mm

G

— -

M

x

M

y

exp[j( a x+B y)]

~" (8.2.11)

the equation (8.2.4 ) M = E and M = -E so equation

x y • y x

(8.2.11) leads to:

o ~

 amn ~ amn

k

o ~

0

 mn

-E.

explj(amnx+Bmny)]

(8.2.12)

or

2 2

mn o ~

 a mn

k +S — a. R

o

 p mn mn

 pmn

-E

(8.2.13)

which is the same as equation (2.12) in Chapter 2.
119

8.3 MINIMIZATION IN THE DOMAIN FOR THE CONJUGATE GRADIENT

METHOD

2

The algorithm that minimizes the error functional ERRF =

* * -1 II 1/2

r (AA ) xr \\ ' is the following:

mn

for n=0

(0)

The equations for the n iteration are:

n

£(n)

-(n)

2

2

n

ft

p

n

£<n+l) 2

?

(n)

 II

 2

n = n+1
8.4 CODE

BB

Fig. 8.4.1 Square cell

C(

\

^

/

1

.x

V

AA _^i

DD

_ J_

120

AA= Distance between centers of vertical strips in x

direction (INPUT)

BB= Distance between inner edges of vertical strips in x

direction (INPUT)

CC= Distance between centers of horizontal strips in y

direction (INPUT)

DD= Distance between inner edges of horizontal strips in y

direction (INPUT)

F= Frequency (INPUT) ^

IOPT=1 for rectangular meshes (INPUT)

IOPT=0 for parallel wire grids (INPUT)

PSI is angle fl (INPUT)

ITM=0 for TE polarization (INPUT)

ITM=1 for TM polarization (INPUT)

NOI = Number of iterations (INPUT)
121

IX = Sampling length for FFT (INPUT)

ALAMB = Wavelength

NX = number of samples incident in the aperture along x.

NX

NX,

l\ m

Fig. 8.4.2. Sampling Arrangement

NX1,NX2 = Edge points in x direction

NY1,NY2 * Edge points in y direction

VfU = Expressions of the Floquet modes a and

-j yk2-(v2+u2)

9 2 ?

for k > U + V

-(VS-IT) for k46 < U^ + v"

EXI/EYI = x and y components of the incident electric field

(INPUT)

HXI,HYI = x and y components of the incident magnetic field

(INPUT)

FFT3D = 3 dimensional complex Fast Fourier transform
122

FOR SIT

X,Y 2 dimensional arrays for the aperture field

JCX,JCY = 2 dimensional arrays for the current densities

CONX,CONY = 2 dimensional arrays to store the constant:

C - F[-Ht + 6(Ht

inC)]

XIX,YIY,XIY,YIX arrays used to store the perturbed aperture

fields.

GX,GY,HX,HY partial derivatives used in the contraction

operation

FOR THE S.D.C.G.

ZINT internal impedance

DX,DY direction vectors"

RX,RY residual vectors

X,Y are either the aperture fields, or the induced

currents (the unknown)

TX,TY two dimensional arrays used to store different values

AN c^

BN 3n

PHASEX,PHASEY arrays used to store phase information

CREFX,CREFY x and y components of the reflection

coefficient

REFF,REFT,RETT,RETF reflection coefficients along theta and

phi angles for TE TM polarizations

Z1,Z2 one dimensional arrays used to store the amplitude of

the unknown X and Y so that they can be used for

plotting purposes for the cross-unit cell
123

MDDX MUPX

MUPY

MDOY

Fig. 8.4.3. Cross-Unit Cell

MDOX left edge point of perpendicular strips

MUPX right edge point of perpendicular strips

MDOY lower edge point of horizontal strips

MUPY upper edge point of horizontal strips
124

8.5 FLOW-CHART FOR THE SPECTRAL DOMAIN CONJUGATE GRADIENT

METHOD

READ

AA,BB,CC,DD,Frequency

, , ,IX,ITM,NOI

Determine sampling

points that correspond

to the aperture and the

conducting regions, i.e.

NX1,NX2,NY1,NY2,NX,NY

Calculate the Floquet

Coefficients a, 3mn

Calculate th

I

e components

of the incident electric

and magnetics fields E ,

Ey' Hx' Hy*

Give initial guess for

J<0) or E(0>

Evaluate the residual

(0) vector r

direction vector p

and the

(0)

Evaluate ERRF

(0)

Calculate

See Chapter 3 & 4

or

(n+l)=E(n). (n)

JtFind ERRF(n"'"1)

(See Chapter 3 & 4)

Define the new residual

vector r(n+1)

Calculate the factor n

is

ERRF

(n+1)

YES

Calculate the amplitude

and phase of current

densities and electric

aperture fields.

Evaluate reflection

Coefficients
125

8.6 LISTING OF THE S.I.T. METHOD

C*****SIT.FORT*****

COHPLSX COHE,CZ3EO,CDET,CX!l!r,COHS,CXt:ij;,CSEFX,CREFY

COMPLEX CONX1(32,32)/1024*(0.0,0.0)/

COMPLEX CONY1(32*32) /102ft* (0.0,0.0) /

COMPLEX COMX2(3?,32)/102«*(0.0,0.0)/

COMPLEX CONY2(32,32) /102U* (0.0,0.0) /

COMPLEX COMX(32,32)/102ft*(0.0,0.0)/

COMPLEX CONY(32,32)/102ft*(0.0,0.0)/

COMPLEX XO(32,32)/102U*(0.0,0.0)/

COMPLEX Y0(32,32)/102»*(0.0,0.0) /

COMPLEX G (32, 32) /f 02ft* (0.0,0. 0)/

COMPLEX JCX (32,32) ,JCY (32, 32)

COMPLEX Y(32,32)/102tt* (Q.0,0.0)/

COMPLEX X(32,32)/102ft*(0.0,0.0)/

COMPLEX YIX(32, 32) /102ft*(0. 0,0.0) /

COHPLSX XIX(32,32)/102ft*(0.0,0.0)/

COMPLEX YI7 (32 ,32) /1 02ft* (0.0 , 0. 0) /

COMPLEX XII(32,32)/102ft*(0.0,0.0)/

COMPLEX GX (32, 32) , GT (32, 32) , HX (32,32) ,HY (3 2, 32)

COMPLEX J,aXT,HYI,CWK(32) ,&11,A12,A21,A22,DENO

H3AL K,X2,HHK(3ft2)

DIHENSIOH AHP(32),HLNDEX{32) ,IVR(3ft2) , CEOS (3 2)

HEAL H{32)/32*0.0/

SEAL 7(32,32)/102ft*0.0/

C *»* AA-DISTAHCE BBT7S3N CENTERS OP 7ESTIC. ST8IPS IS X-DIRECTIOH **

C **« BB»DISTAHCB BETWEB3 IHHEH EDGES OP VSBTIC. ST2IPS IS X-DIBECT.»*

C *** CC'DISTASCE BBTWEEH CENTEBS OF HORT2. STRIPS IH Y-DIHECTIOH *•

C *«* 00=DISTA»CE BETWB5H INHEP EDGES OP HOBIZ. STRIPS IH X-DIBECT.**

HEAD (1,22) AA/BB,CC,DD,F,EHR

22 FOR RAT (8 E10. ft)

C *** IOPT=0 F0» A RBCTAHGULAB OR SQOARB HBSH *****

C *** IOPT=1 FOE A PABHALLEL GRID *****

IOPT» 1

I?(IO?T.GT.O) CC»1.500B*15

I?(IOPT.GT.O) DD«1.500B*15

WHITE {3, 33) AA,.B3,CC,DD,EBB

33 FORHATCO1, ' A= «, 715.8,' B=

 f,F15.8,' C= «,F15.8,

3» D» »,nS.8,« ERR" ',F15.8)

iRITE(3,ftU) F

HU FOtlHAT('0«, « FRSQ = «,B10.ft)

HEAD (1,22) PHI,THI,PSI

WRITE (3, 55) PHI,THI,?SI

!?5 FOP.RA?{'0«, ' PHI= »,F10. 1,« TH2TA= ',F10.1,« PSI»

C *** READ THE NUMBER OF SAMPLING POINTS  

READ (1,66) IX

C ** * ITK=1 ?OR TH POLARIZATION ********

ESAD(1,66) ITB

C *** REA D H03BER OF ITERATIONS *******

PEAD(1,66 ) SOI

• 06 F03!1AT(I3)
126

PI=2. 141503

PI2=?I/2.

TPI=6. 283185

R7D=57. 29578

EP=8.85UE-12

2?A=SQR7(nU/EP)

JsCHPLX(0. 0,1.0)

ITER = 1

CONE=CHPLX(1.0,0.0)

C2ERO=C.1PLX (0. 0,0. 0)

ALAKB=C7/F

AA=AA/ALAWB

BB-BO/ALAHB

CC=CC/ALAMB

DD=DD/ALAHB

C *** DETERMINE SAMPLING POINTS THAT CORRESPOND TO THE

C CONDUCTING REGIONS AND THE APERTURE *********

• NX=IFIX(BB/AA*PLOAT(IX) *2.)/'»*2

NY=IFIX(DD/CC*FLOAT(IX) *2.)/«*2

NX1=(IX-NX)/2+1

NX2=NX1*NX-1

NY1=(IX-NY)/2*1

HY2=NT1*NI-1

WRITE (3 , 1 00) NX, NX1 , NX2 , HY , NY 1 , NY2

100 FORHAT('0', « NX=«, 13, 3X, «NX1a« ,I3,3X,

 f HX2=« ,I3,3X,

• NY=»,I3,3X,'NT1=',I3,3X,«NY2=«,I3)

K2=K**2

STSPK= SIN (THI/RTD) *SIN (PHI/HTD) *K

STCPK=SIN (THI/RTD) *COS (PHI/RTD) *K

CPS=COS (PSI/RTD) /SIN (PSI/RTD)

110 CONTINUE

C ***** CALCULATE FLOQUET BODES *********

DO 200 H=1,IX

IF (M.GT.IX/2+1) GOTO 125

U (R) =TPI* (f!-1) /AA-STCPK

GOTO 127

125 U(M)=TPI*(3-IX-1) /AA-STCPK

127 CONTINUE

DO 190 N=1,IX

IF(H.GT.IX/2*1.AND.N.GT.IX/2*1) GO TO 160

IFfH.GT. IX/2+1) GO TO 150

IF(N.GT. IX/2*1) GO TO 140

V(H,N)=TPI * (N-1J/CC-TPI* (H-1)/HA*CPS-STSPK

GO TO 170

1«0 V (H,N) =TPI* (N-IX-1) /CC-TPI* (H-1) /AA*CPS-STSPK

GO TO 170

1 50 V (M,H) =TPI* (N-1) /CC-TPI* (H-IX-1) /AA*CPS-STSPK

GO TO 170

160 V(M,H)=TPI*(N-IX-1)/CC-TPI*(a-IX-1)/AA*CPS-STSPK

170 IF(K2.GE.fl{N)**2«-V(!l,N) **2) G(H,N) =-J*SQRT(K2- (U (H)**2+V (B,N) **2

• ))

IF (K2.LT.tJ(H) **2*7(H,N)**2) G{N,N) =-SQET (U (M) **2*V (H,N) **2-K2)
127

190 CO?/TIf.'03

200 cOiVTinas

IF (ITM.GT.O) GO TO 210

c   INCIDENT FIELDS FOS TE POLARIZATION ** **

SXI=SIN(-?HI/RTD)

EYI=COS(PHI/RTD)

HXI=COS(PHI/HTD)*COS(THI/P.TD)/ETA

HYI=SIN(PHI/aTD) *COS(THI/RTD)/ETA

GO TO 261

c ***** IHCIDEST FIELDS FOR TH POLARIZATION ******

210 EXI=COS(PHI/RTD)*COS(THI/RTD)

EYI=SIN(-PHI/RTD)*COS(THI/RTD)

HYI=SItJ (PHI/RTD-PI2) /ETA

HXI=COS(PHI/RTD-PI2J /STA

261 CONTINUE

C *** GIVE A GOESS FOR INITIAL APEBTOSE FIELDS X S Y ***

DO 310 H=NX1,NX2

00 300 N=NY1,NY2

300 CONTINOE

310 COVTIHDE

C   START THE C03POTATIOU 0? CONSTANTS CONX & CONY THAT DEPEND

C ON A GIVEN INCIDENT FIELD ** 

DO 315 I»1,IX

DO 315 L=1,IX

CONX1 (I,L) =*1. *HXI*B*aa/J

.. CONYJ:(I,L)=*1.*HYI*W*OTJ/J

COHX2 (1,L) =-1. *HXI*W*00/J

COMY2(I,LJ=«-1.*HYI*ff*Oa/J

315 CONTINOE

C *** PER70R3 THE T30NCATION OPERATION ***

DO 325 I=NX1,NX2

DO 325 L=HY1,NY2

CONX1 (I,L)=C2ERO

325 CONY1 (I,L)=CZERO

C *** TAKE THE FOORIER TRANSFORM ****

CALL FFT3D(CONX1,IX,IX,IX,IX,1,69,IWK,BBK,CSK)

CALL PFT3D(CONY1,IX,IX,IX,IX,1,69,IHK.RVK,CWK)

CALL FFT3D(COHX2,IX,IX,IX,IX,1,69/IHKrBWK,CWK)

CALL FFT3D(CONY2,IX,IX,IXf IX,1,69,IWK,R»K,CHK)

DO 3UO «=1,IX

D0^3aO 11=1, IX

CONX {t1,N)-CONX1 (M,N ) 4-CONX2 (S,K)

3 UO CO!IY(H,H)=COHY1(H,N ) *CONY2(H,N)

DO 350 !1=1,IX

DO 350 N=1,I X

CDET=- (a (M) *V (H, K) /r, (M, H) J **2- (V («,N) **2/G («,N) -G (K, SJ )

CONX(B,!I)

COHX (f1,H}= (-11 («1) *V(«,H)/G(n, N)*COMX(a,N)-(V(H,N) **2/G (H, N) -G

.*CONY(3,N))/CDET

M,N) = {-{f;{M,fI)-tJ(M) **2/G(fl, N) ) *COHS»0 (MJ *V{K, N) /G (fl, N) *
128

.CONY (M,N))/CDET

350 CONTINUE

C **** Kt.'D OP CALCULATION OF CONX 6 CONY *****

C

c ***« NOB START ITERATIVE PROCESS I !!!!!!!!!!!!!! ' *****

C

C *** S3? AH PARTIAL DERIVATIVES EQUAL TO ZERO *****

000 DO 541 H=1,IX

DO 541 H=1,IX

XIX (S,N)=CZEHO

YIX(«,N)=CZEI50

XIY(N,H)=CZERO

YIY(M,N) =CZERO

541 CONTINUE

c   PERTURB APESTTJRE FIELDS BY (0.01,0.01) ** *****

DO 543 H=NX1,NX2

DO 543 S=NY1,SY2

XIX(H,N)=X(S,N) * (0.01 0,0. 010)

YIX(«,N)=Y(H,N)

YlY(H,N)=Y(a,H)*(0.010,0.010)

543 CONTINUE

C *** TAKE THE FOURIER TRANSFOBH OF THE APERTURE FIELDS X 6 I  

CALL FFT3D(X,IX,IX,IJ,IX,1,69,IBK,RHK,CBK)

CALL FFT3D(Y,IX,IX,IX,IX,1,69,IHK,RWF,CHK)

C ** ** nflLTIPLY TRAHFOBHED FIELDS BY THE FLOQaET COEFFICIENTS

C AND GREEN'S FUNCTION ' ****

DO 560 H=1,IX

DO 550 N=1,IX

CXHN»X(H,N )

X{H,N)=(D(M)*V(H,M)/G(H.N)*X(n,H ) + (V (H, N) **2/G (H, N) -G(B,N))

. *Y(»,H))

Y (H,N ) =• ( (G (H, N) -U (H) **2/G (H,N) ) *CXHH-U(H) *V (R,N) /G (H, N)

. *Y(fl,N))

550 COHTINOS

560 CONTINUE

c ***** TAKE THE INVERSE FOURIER TRANFORB ******

CALL FFT3D(X,IX,IX,IX,IX,1,-69,IHK,RWK,CVK)

CALL FFT3D(Y,IX,IX,IX,IX,1,-69,IHK,HV!K,CWK)

WRITE (3, 570) ITER

570 FOR«AT(3X,/ • ITERATION NUMBER • ,12)

C *** CALCULATE CURRENT DENSITIES *****

DO 600 f1=1,IX

DO 600 H=1,IX

JCX (H, N) =(Y (H,N) * J/W/UU+Hlfl )#(-2.)

600 JCY (M, N) =fr (B/N) *J/W/UU*njri)»(2.j

C *** PLOT CURENTS ON STRIPS *«•***

DO 620 1=1, IX

ABP(I)=CADS(JCY(1,I))

HINDEX (I) = (FLOAT (I-IX/2) -. 5) /IX*AA*1.045

IF (ITSR.nT. (NOI-1) ) WRITE(8,* ) A.1P (I) ,RtNDEX (I)

620 CON'TIMaE

C IF (ITES.GT. (NOI-1) ) CALL GEHPT (HINDEX, AflP, IX , 0)

C **** TRUNCATION  

DO 740 fI=MX1,NX 2

DO 730 N=NY1,NY2
129

X(M,N) =CZRP.O

Y(M,I1) =CZERO

730 CONTINUE

7UO CONTINUE

C ** ** HOB FIND THE 7 (TRUNC (INVERSE P (G E))) *»**

C

CALL FFT3D(X,IX,IX,IX,IX.1,69,ItfK,RWK,CWK)

CALL FPT3D(Y,IX,IX,IX,IX,1,69,IWK,RWK,CWK)

ITEE=ITER+1

00 760 H=1,IX

DO 750 N=1,IX

CDET=- (0 (H) *V (B,N) /3 (N, N)) **2- (V (B,N) **2/G (fl,N) -G (H, N) )

. * (G (M, N) -0 (H) **2/r, (H,N))

CXMN=X(B,N)

X(M,N) = (-0{a)*V(M,H)/G(H,N)*X(B,N) -(V(H,M) **2/G(M,M) -G(B,H))

.»Y(M,HJ)/CDET

Y(M,N)= {-(G (H,H)-a<H) **2/G (8, N) ) *CXHH + 0(M) *V (H,N)/G (H, H) *

. Y(M,N))/COET

750 CONTINUE

760 CONTINUE

C **** ADD P (TSUN (INVERSE P (G E)) ) TO CONX AND CONY ****

DO 780 H=1,IX

DO 770 N=1,IX

X {«, N) =X (H, S) *CONX (S,N)

Y (H, H) =Y (M, N) *CONY (S,N)

770 CONTINUE

780 CONTINUE

C *** CALCULATE THE REFLECTION COEFFICIENTS CREFX AND CHBFY ***

CR3FX=X(1,1)/{FLOAT(IX)*FLOAT(IX))-EXI

CREFY=Y (1,1)/(FLOAT(IX) »FLOAT (IX) ) -EYI

REFX=CABS(CHEFX)

R2FY=CABS(CREFY)

WP.ITE(3,800)REFY

SOO FORHAT(10X,2F10.3)

C *** TAK E THE INVERSE FOURIER TBANSFORH OF THE RESULT TO OBTAIN

C A NEW VALTIE FOR THE APERTURE FIELDS  

CALL PFT3D(X,IX,IX,IX,IX,1,-69,IiK,RSK,0!K)

CALL FFT3D(Y,IX,IX,IX,IX,1,-69,IWK,RWK,CiK)

C ** * PLOT APERTURE FIELD ***

DO 830 1=1,IX

A«P(I)=CABS(Y(I,16))

C CROS(I)*CABS(Y(8,I))

RINDEX (I) = (FLOAT (I-IX/2) -« 5) /IX*AA*1. 095

IF(ITER. EQ.NOI) »RIT2(8,*) ABP (I) ,RIHDEX (I)

830 CONTINUE

C IF(ITER.GE.NOI) CALL GESPT(RINDEX,ABP,IX,0)

C CALL GENPT {RINDEX,CROS,IX, 0)

C *** REPEAT SABE PROCESS FOR PERTURBED FIELDS XIX,XIY,YIX,6 YIY***

CALL FFT3D(XIX,IX,IX,IX,IX,1,69,IWK,RWK,CHK)

CALL FFT3D(YIX,IX,IX,IX,IX,1,69.IHK,RWK,CHK)

CALL FFT3D(XIY,IX,IX,IX,IX,1,69,IBK,RWK,CHK)

CALL FFT3D(YIY,IX,IX,IX,IX,1,69,IWK,RHK,CWK)

DO 850 !t=1,IX

DO aao v=i,ix

CXKN=XIX(B,!I)

XIX (i1,N)= (H (.1) *V (M, :») /G (S, N) *XIX («,N) * (V (M, N) **2/G (?1, N) -3 (H,N))
130

H,M)= «G(n,H) -U (M) **2/G(H,'.l)) *CXMN-J| (!1) *V (M,N) /G (H, N)

*YIX(.1,W))

CXI«N=XIYd,N)

XI Y (*,H) = (0 (.1) * V (.1, N) /(7 (,1, N) *XIY (M, N) * (V (.% N) **2/G(H, H) -G (H,H) )

*YIY(M,N)}

YIY(M,B)ss{(G(.1,N)-a(''!} **2/G(H,N)) *CXIMN-0 (fl) *V (fl, N) /G (H, N)

*YIY(fl,NJ)

CONTINUE

850 COHTIfiaE

CALL FFT3D{XIX,IX,IX,IX,IX,1,-69,IHK,RWK,CVK)

CALL FFT3D(YIX,IX,IX,IX,IX,1,-69,IHK,PHK,CWK)

CALL FFT3D(Xiy,IX,IX,IX,IX,1,-69,IHK,aHK,CWK)

CALL FFT3D(YIY,IX,IX,IX,IX,1f-69.IHK,SUK,C«K)

,C

C**« K07 PERFORM THE TROHCATIOH OPERATION FOB PERTOHBED FIELDS ***

C

DO 870 H=NX1,NX2

DO 860 N*HY1,HT2

XIX(!!,H)=CZERO

YIX(H,N)=CZERO

XIY(H,N)=C2BRO

YIY(H,N)=C2ERO

860 COHTIHOE

870 COHTINOE

C

C **** HOW FIND THE F (TRONC (INVERSE F (G EP) ) ) ***

C

CALL FFT3D(XIX,IX,IX,IX,IX,1,69,IWK.RBK,CHK)

CALL FrT3D(YIX,IX,IX,IX,IX,1,69,I»K,RHK,CHK)

- CALL FFT3D(XIY,IX,IX,IX,IX,1 ,69,I«K,RWK,CWK)

CALL FFT3D(IIY,IX,IX,IX,I3t/1,69,I»K,BHK,CWK)

DO 910 H*1,IX

DO 900 N=1,IX

CDET=- (0 (H) *V (S, N) /G (H, H) ) **2- (V (S,N) **2/G (H, N)-6 (H, H) )

. * (G (H, H) -0 (M) **2/G (S,N) )

,

XIX (», H) = (-0 (H) *V (M, N) /G (H, H) *XIX (H, N) - (V (fl, N) »*2/G (H, H) -G (H,N) )

.*YIX(«,H))/CDET

YIX («,N) = (- (G (H,R) -0 (H) **2/G (M, N) ) *CXHN*0(H) *V <8,H)/G{H, ») *

.YIX(M,N))/CDET

CXIHH=XIY(H,N)

XIY (fl. N) = (-0 (M) *V (fl, N) /G (fl, N) *XIY (H,H) - (V (H, N) **2/G (H, N) -G (H, N) )

.*YIY(M,N))/CDET

YIY(fl,N ) = (-(G(H,N ) -D(H) **2/G (H, N) ) *CXIHN*0 (H) *V (H, N) /G (S, N) *

.YIY(fl,N))/CDET

900 CONTINUE

910 CONTINUE

DO 930 H*1,IX

DO 920 N = 1,IX

XIX (fl, N) =XIX (H,N) -t-CONX (H,»)

YIT (a,N)=YIX(B,N)*CONY(fl,N )

XIY(H,N ) =XIY(H,N)*COMX(n,N )

YIY(H,H)aYIY(H,B)-».COMY(N,H )

920 CONTiriUH

930 CONTINUE
OpGfNAL PAGE IS 131

FtPOR QUALITY

CALL FFT3D (XIX, IX, IX, IX, IX, 1 ,-69,IWK,RWK,CWK)

CALL FFT3D(YIX,IX,IX,IX,IX,1,-69,IWK,BWK,CWK)

CALL FFT3D(XIY,IX,IX,IX,IX,1,-69,IHK,RWK,CWK)

CALL FFT3D(YIY,IX,IX,IX,IX,1,-69,IHK,RWK,CWK)

C *** EVALUATE PARTIAL DERIVATIVES GX,GY,HX,6HY  

DO 950 «=1,IX

DO 9UO N=1,IX

GX(M,W) = (XIX(H,N)-X(a,H) )/ (0. 010,0.0 10)

RX(K,N)=(YIX(>l,N)-Y(n,?l))/ (0.01 0,0. 010)

GY(H,N}={XIY(!1,N)-X{flrN) )/ (0.010,0. 010)

940 HY (H, H)= (YIY (H,N) -Y (H, H) ) / (0.010,0. 010)

950 CONTINUE

C ** * IMPROVE PREVIOUS ITEEATE FOE APEBTOSS FIELDS BY USING

C A CONTRACTION FACTOB ****

DO 960 M=NX1,NX2

DO 9"60 N=NY1,NY2

DEKO=GX(a,N)*HY(H,N)-HX(H, N) *GY (H,N) -GX (H,N) -HY(H, N) *1.

A11=(1.-HY(M,?I))/DENO

A21=-HX(H,N)/DENO

A22= {1.-GX (H, N) ) /DENO

X(H,N)=A11*X(!1,N) +A1^*Y(H#N) * (1. -A11 ) *XU (H,N) -A12*YO (M,H)

Y(a,N)*A21*CXSN+A22*Y(!!,N)-A21*XO(H,N)*(1.-A22)*YO(a,M)

960 CONTINUE

I7(ITER.GT.HOI) GO TO 1000

C .

C *** NOW TRUNCATE THE IBPROVBD APEBTUBB FIELD E '  

C

DO 980 fl=1,IX

DO 970 N«1,IX

IF(n.GE.NX1.AND.H.LE.NX2.AND.N.GE.NYl.AND.N.LB.NY2) GO TO 965

X(H,N) =(0.00,0.00)

Y(M,N)=(0.00,0.00)

965 Xn(B,N)=X(H,N)

YU(H,N)=Y(H,H)

970 CONTINUE

980 CONTINUE

GO TO UOO

1000 CONTINUE

STOP

END
132

8.7 LISTING OF THE S.D.C.G. METHOD FOR THIN STRIPS WITH

THE SQUARE-SHAPED UNIT CELL

C** *** THIS IS TIIF CONJG. GFAO. METHOD FOR CtlP.tlEMTS OH THir WIRES**

C  MINIMIZATION IN THE RANG E (VAN DER BEHG) ** ****

COMPLEX COME,CZEPO,CXHN,F1 0

COMPLEX C5EFX,CREFY,CRE?,C!?ET,ZIOT

COKPL2X G (32, 32)/102«* (0.0,0.0) /

COMPLEX Y(32,32)/102<»*(0.0,0.0)/

COMPLEX X(32,32)/102U*(0.0,0.0)/

COBPLKX YH(32, 32)/102a*(0.0, O.O)/

COMPLEX XU (32,32)/102U*(0.0, 0.0) /

COMPLEX RX(32,32)/102H*(0.0,0.0) /

COKPLEX RY{32,32)/102H*(0.0,0.0)/

COMPLEX J,HXI,HYI,CWK(32 )

COMPLEX DY(32,32)/102«* (0.0,0.0) /

COMPLEX DX(32, 32) /102<»*(0. 0,0.0) /

COMPLEX TX(32,32) /102tt*(0.0,0.0) /

COMPLEX TY(32,32)/102«»*(0.0,0.0)/

PEAL K,K2,RWK(3ft2 )

DIMENSION IWK(342) ,RR(350) ,CH{350) ,PHASEX(32,32) ,PHASEY (32,32) ,

.21 (102«) ,X1 (102ft) ,Y1 (102U) ,22 (102«) , AHP(36), 8INDEX (36)

HEAL 0 (32)/32*0.0/

S2AL V(32,32)/102»*0.0/

C *** AA=DISTAHCE BETWEEN CENTERS OF VERTIC. STRIPS IN X-DI3BCTION **

C *** BB=DISTANCE BETWEEN INNER EDGES OF VERTIC.STRIPS IN X-DIRECT.**

C *** CC=DISTANCE BETWEEN CENTERS OF RORIZ. STBIPS IN Y-DIH2CTION **

REAO(1,10) AA,BB,CC,DC,F,ERP

10 ' FORMAT (3E10.ft)

., F=2.9<»8E*9

C *** IOPT=0 FOR A SQUARE OB A SECT1NGOLAR SESR  

C ** * IOPT=1 FOR A PARALLEL GRID ****

IOPT=1

IF(IOPT.GT.O) CC*1.500E*15

IF(IOPT.GT.O) DD=1.500E*15

WHITS (3,20) AA,BB,CC,DD,ER1?

20 FOPHAT('0«, ' A= ',F15.8,« B= «,F15.8,« C= '.FIS.S,

.« D= «,F15.8,« E2P= «,F15.8)

W5ITE(3,30) F

30 FORMAT('0«, » FREQ = «,E10.*)

R1?AD(1,10) PHI,THI,PSI

WSITE(3,aO) PHI,THI,PSI

10 FORBATCO*, 1 PHI= •jPIO.!, * THETA= • ,F10.1,« PSI= »

C *** REA D THE JIOHB5R OF SAMPLING POINTS ******

ESAD(1,50) IX

C *** ITM=0 FOR TE POLARIZATION  

C ** * IT«J=1 FOR TH POLARIZATION ****

HEAD(1,50) ITM

WFTTE(3,45) ITfl

C *** TEAD THE NnHBER OF ITERATIONS ***

REA D (1,50) NOI

«5 FORMAT (3X,'T11^ 7ALUE FOR ITU IS= ',13)

50 FOnrfAT(T^ )

PI=3.iai593
133

TPT=r>.

CV=2.997956E*8

W=«.E-7*PI

PTD=57. 29578

EP=8.8SUE-i;>

ETA=SQRT(no/i:p)

J=CHPLX(0.0,1.0)

IT TIP. = 1

CONE=CHPLX (1.0,0.0)

CZE!>0=C!1PLX(0.0,0.0)

SIG!"A=5.E20

ZINT=(1.0,1.0) *SQP.T(B*OtJ/2./SIGHA)/(1.0)

ALAHB=C7/F

AA=AA/ALAMB

BB=BB/ALARB

CC=CC/ALAHB

DD-DO/ALAHB

NX=IFIX (BB/AA* FLOAT (IX) *2. ) /tt*2

NY=IFI)C (DD/CC*FLOAT( IX) *2. ) /U*2

NX1=(IX-HX)/2*1

HY1=(IX-MT)/2+1

WRITE(3,60) SX

60 FORHAT(«0«, ' MX=',I3,3X,«NX1=',I3,3X,«NX2=' ,I3,3X,

' NT=«,I3,3X,«NY1=«,I3,3X,»NY2»«,I3 )

K=TPI/ALAHB

. K2=K**2 •

STSPK= SIN (THI/RTO) *SIM (PHI/RTD) *K

STC?K=SIS(THI/RTD) *COS (PHI/RTD) *K

C?S=COS (PSI/FvTD) /SIN (PSI/RTD)

70 CONTINUE

*** DEFINE THE FLOQ02T C02FFICIBNTS ****

DO 100 B=1f IX

IF(M.GT.IX/2+1) GO TO 75

0 (H) =T?I* (H- 1) /AA-STCPK

GO TO 80

75 n(J1)=TPI*(K-TX-1) /AA-STCPK

80 CO«fTIH7E

DO 90 N=1,IX

I?(H.GT. IX/2»1.AND.H.GT. IX/2+1) GO ?0 8tt

IFJK.GT. IX/2 + 1) GO TO 83

TF(H.GT. IX/2 + 1) GO TO 81

7 (M,K) =TPI* (N-1.) /CC-TPI* (H-1) A A*CPS-STS7»K

C,0 TO 85

81 7(W,K)=TPI * (U- IX- 1) /CC-TPI* (H-1) /AA*CPS-STSPK

GO TO 85

83 V (M, N) =TPI* (H- 1) /CC-TPI* (H-IX-1) /AA*CPS-STSPK

CO TO B5

m V (M,!!) =TPI* (H-IX-1) /CC-TPI* (H-IX-1 ) /AA*CPS-STSPK

G(M,K ) =-J*SQRT (F.2- (O (B) **2*7 (H,H) **2

IF (":2.LT.rJ(W)**2 + y («,fJ)**2) G(P1,N ) =-SQKT (f (.1) «*2*V (.1, ») **2 -K2)

*CONE
134

90 COVTINUE

100 rn:iTi!iUE

IF (ITM.GT. 0) GO TO 110

C

C *** INT"JT TOP. TE POLARIZATION ****

FXI = SIM(-PHI/ETD)

?YI=COS(PHI/T!TD)

HXI=COS (PHI/STD) *COS (THI/RTD)/ET A

HYI'SIN (PHI/P.TD) *COS (TIII/P.TD) /ETA

GOTO 120

C **** INPUT FOR TB POLARIZATION ****

C

110 EXI=COS(PHI/RTD) *COS (THT/RTD)

FYI=SIN (-PHI/RTD) *COS (THI/HTD)

HYI=SIN(PHI/RTD-PI2)/ETA

HXI=COS (PHI/RTD-PI2) /ETA

FT=1.0*COS(THI/PTD)

120 CONTINUE

c **** SET YOOR INITIAL GUESS FOB X AND Y  

C

123 DO 130 H«1,IX

DO 125 N=1,IT

X(n,N)=CZEP.O

Y{H,N)=CZERO -... -

125 CONTINUE'

130 CONTINUE

c ** ** WORK ON INITIAL GUESS***

C

CALL FFT3D(XU,IX,IX,IX,IX,1,69,IWK,RWK,CWK)

CALL FFT3D (YU, IX, IX, IX, IX, 1 ,69.IHK,RWK,CSK)

DO 160 B=1,IX

DO 150 N=1,IX

CXHN=XtJ(H,N)

XU(H,N) = ({G(H,N)-T(H,N)**2/G(H,N) )*XU(H,N)- (U (d) *7 (H,H)/G

. *Yn(a,N))/(J*H*EP)/2

YD (H, N) = (-U («) *V (H,N) /G (H,N) *CXMN* (G (H, N) -U (H) **2/G (H, N) )

. *YU(H,N))/(J*H*EP)/2.

150 CONTINUF,

160 CONTINUE

CALL FFT3D(XIT,IX,IX,IX,IX,1,-69,IWK,RHK,CWK)

CALL FFT3D(YU, IX, IX, IX, IX, 1,-69,IWK,RWK,CWK)

C WRITE (3, 1^0) ITEH

170 FORKAT(3X,/ • ITERATION NUHBER ',12)

DO 200 H=1,IX

DO 190 N=1,IX

C ** COMPUTE THE ERROR FOR YOOR INITIAL GUESS ** *

RY (M,N )

IP (M.GE.HX1.AND.H . LE. NX2 . A »D. N.GK. NY 1 . AMD. N.I.E.KY2 ) RX (S,N) =C7ERO

IF(^.G^.-!X1.AN'D.n.Lr;.?iX2.AND.N.GK.?IY1 . AND. N'. LE. t!Y2 ) RY (H,N) =CZ ERO

ERROR=ERROR*R X (U, H) ->CONJG(R X (n,N ) ) +RY (H, N) *COSJG (RY («,N ) )

F5 = F5 + RX (fl,N)*CONJG(RX (H,«) ) *RY(P1,N) *CONJG (RY (M, N) )
135

190 CONTINUE

200 CONTINUE

C

.c **** FIND THE FOUP.IER TP.ANSFOTW OF P.F.SID'IAL ***

C

CALL PFT3D (DX,IX,IX, IX,IX,1,69,IHK.RWK,CWK )

CALL FFT1D(DY,IX,IX,IX,IX,1,69,IKK,RHK,CWK )

C ** .iriTILPY DY THE CONJUGATE TRANSPOSE OF THE PATHIX

C TO FIND THE VECTORS DX AND DY ***

DO 220 .1=1, IX

DO 210 N=1,IX

CXHN=DX(M,N )

DX (X, N) = (CONJG (G (M, N) -V (M,N) **2/G (S,N) ) *DX{!1, N) -

CON.1G (V(H,N)*a{f1)/G(n,N)) *DY 11, H) ) /CONJG

OY(H,N ) = (COHJG(-V(!!,N) *U(M) /G (H,N) ) *CXHN*

CONJG (G(1,N)-0(«> **2/G(H,N))*DY(H,H))/CONJG (J*W*EP)/2.

210 CONTINUE

220 COKTIN02

C **** HOW FIND THE 1*7. 500SIER TBAS. OF THE DIREC. FUHCTIONS **

C

CALL FFT3D(DX,IX,IX,IX,IX,1,-69,IWK,RWK,CWK)

CALL FPT3D(DY,IX,IX,IX,IX,1,-69,IBK,RWK,CBK)

DO 360 H=1,IX

DO 350 H=1,IX

DX(M,N)=DX(n,H)-RX(M,N)*COHJ(5(ZINT)

DY (M,N ) = DY{H,N)-RT(!1,N) *COHJR(ZINT)

IP (H . GE. HX1 . AND. M. LE. NX2. AND. N. GE. NY 1 . AND. N. LE .NY2 ) DX (H ,N ) *CZERO

IFltl.GB.'NXI.AHD. M. LS.NX2.AND.N.GE. NY1. AND. N. LE.NY2) DY («,N) *CZEHO

TY(H,H)=DY(H,N )

P3=F3+COSJG{DX(M,W)) *DX{M,W) *CONJG (DY (fl, N) )*DY(M,» )

3SO COKTIKOE

360 CONTINUE

c 4,*** HOH START THE IT3RATIOS PROCESS {niNIHIZATION) ***

C

c   HtJL^IPLY THE DIRECTION VECTOR BY THS MATRIX Z ** **

c   STORE YOUR RBSUTLS IN VECTORS ?X AND TY ** **

365 CALL TFT3D(TX,IX,IX,IX,IX,1,69,IVK,RWK,CWK)

CALL ?FT3D(TY, IX, IX, IX, IX, 1, 69,IHK,PWK,CWK)

DO ftOO H=1,IX

DO 370 N=1,IX

CXHN=TX(.1,N)

. *TY(1,M))/(J*lf*EP)/2.

TY (H, N) « (-tl (H) *V («,N) /G (B,N) *CTRN* (G (a, N) -U (S) **2/G (H,N) )

. *?Y(t!,N))/(J*W*EP)/a.

370 CONTINUE

i»00 CONTIVUS

CALL PFT3D(TX,IX,IX,IX,IX,1,-69,IWK,nWK,CWK)

CALL FFT3D (TY, IX ,TX, IT ,It, 1,-69,TTrrK,RWX, CWK )

F1=0.0

DO U10 PI = 1,IX

DO U10 N=1,I X

TX (fl,N)=TX(M,N)-D X (M.N ) *ZINT

TY (0,N ) =TY (M,N ) -DY {/1,N) «ZIN T
136

IF(n.GE.MX1.AMD . !1. L». NX2 . AND . N. GE. MY1. AMD. N.L2.NY2 ) TX (H , N) =CZEP.O

TF(P1.GE.NX1.AND. H. LE.NX2 .AND.N. RE.MY1. AND. N. LE.SY2) TY (t!f N) =e!EHO

P1=F1*CO!UG{TX(fl,?t))*TX(*,N) *COMJG (TY (M, N) )*TY(N,JJ)

ITSR=ITER»1

C **** COMPUTE CONSTANT AN ****

AN=F3/F1

CH(ITES)=SQfiT( ERROR) /SQRT (F5)

E?.R=CH(ITE?)*100

C •** CALCULATE ERROR  

BRROR= ERROR- (P3** 2/P1)

C *** GET A HSW ESTIMATE FOR TQIJR nSKNOWRS X 6 T ****

DO 560 H=1,IX

DO 550 H»1,IX

Y («. N) =Tt (S , N) * AH»DT { S, N)

550 COHTIHOS

560 COHTINOB

C **** A NE3 SSTIHATE FOR THE RESIDUAL VECTORS BX & BY ****

DO 580 .1 = 1, IX

DO 570 B=t1,IX

RX (fi,») =RX (H,H)-AN*TX (S, K)

570 COMTINOS

580 CONTIHOS

HR (ITER) aPLOAT (ITER)

WRITE (8,*) CH(lTEP.),Rn(ITEB)

C **** HtJLTIPLT THE RESIDUAL VECTORS BT THE COKJG. TRANS. OP MATRIX Z

CALL PPT3B(TX,IX,IX,IX,IX,1,69,IWK,RWK,CHK)

CALL PPT3D{TT,IX.IX,IX,!X,1,69,IWK,HWK,CHK)

DO 600 H*1,IX

DO 590 Ba1,IX

CXHH=TT(B,U)

)a(COKJG(G(S,!n -V(B,M)**2/G(«,N))*TX(H,inCOHJG (7 (n,W) *n (M) /G (M,H) ) *TY (S, N) ) /CONJG(J***BP)/3.

TY (H, 3) « (COHJG (-7 (H,H) *H (B) /G (B,H) ) «CX!!H*

CORJG (G (H,M) -0 (.1) **2/G (M,f») ) *TY {«,»)) A^OHJG (J*H*BP)/2.

S90 COKTISOB

600 CONTIHOS

CALL ?PT3D<TX,IX,IX,IX,IX,1,-69,IVK,BHK,CWK)

CALL FFT3D(TY,IX,IX,IX,IX,1,-69,IWK,RWK,CWK)

C ***** STOR3 THE OLD VALUE FOR F3 IN 72 TO CALULATS DN LATER ****

P2=F3

DO

no eaa »i=i,ix

TY (H.N)aTY (R,N)-RY(.1,M) *CONJG(ZINT)

TFfM.GE. NX LAND. 1. LS.NX2. AND. N.GE.SY LAND. H.LE.NY2) TX (B,N) =C2ERO

TF(H.GE.3X1.AnD. K.L^.VX2.ASD.K.GE.riYl.AKn. H.LE.NY2) TY (H,N)

F3 = F3*CO!UG(7T(n,V)) *TX(n,N) +COKJG (TY (fl, N) ) *TY (H.3)

fitttt COSTINUK

C *** NOU CAL1LATH 3K •*****•
OF-POOR QUALITY

 137

C1 IF(ITEP.f:Q.20.nR.ITER.EQ..OO.OR.ITER.EQ..60) BN«0.0

C **** OBTAIN A N?W ESTIMATE FOR THE DIRECTION" VECTORS DX 6 DT *

DO 660 !>I = 1,IX

DO 650 N=1,IX

r>X{W,N)=TXn,H ) +Bf?*DX(H. M)

DY (H,N)=TY,(»!,N) *3M*DY (11, N)

TX (N,N)=DX(H,?J)

TY(M,N)=DY(M,N)

650 CONTINUE

660 CONTINUE

C

C **** CONTINUE THE ITERATIVE PROCESS  

IF(ITSR.GT.NOT) CALL GENP7 (RR,CH, HOI, 0)

IF(TTER. GT.NOIJGO TO 700

IT(ERR.LT. 1} GO TO 700

GO TO 365

C ** ST09.E X BY INTO 1-DIHEN. ARRAY S TO BE OSED FOR ANY PLOTTING

C PURPOSES »

700 DO 720 a=1,IX

DO 720 N*1,IX

21 (I)=CABS(X(fl,N»

22(I)=CABS(Y(!1,N))

XI (I)=(FLOAT(H-IX/2)-. 5} /IX*AA* 1.05

Y1 (I)-(FLOAT(N-IX/2)-. 5) /1X*BB*1. 05

WRITE (7,*) X1(I).Y1(I) ,Z1(I),Z2(I)

720 COMTIWOB

G O T O 73 0 • . . ••••••: :

C *** FIND THE -PHASE FOR THE CORRENTS X 6 Y *****

723 DO .725 n=1,IX

DO 725 N-1,IX

PHASSX(H,N)=0. 0

PHASEY(«,N)=0. 0

1F(H.GI!.NX1.AND. H.L2.SX2.AND. !J. GE. NY1.AND. N. LE.NY2) GO TO 725

REX=REAL(X(«,ff))

AIHX=AIHAG(X(N,N))

IF (REX.GE. 0. 0. AND. AIHX.GE. 0. ) PX=ATAN (AIHX/REX) *RTD

IF<REX.LT. 0.0. AND.AIRX.GE. 0.) PX«180.-ATAN (AIBX/R2X) *RTD

IFfREX.LE.0.0. ABD.AIKX.LT. 0.) PX=100. *ATAN (AI-X/REX) *RTD

TF(Rra.GE.O.O.AND.AISX.LT.O.) PX=360.-ATAN (AIHX/REX) *RTD

»HASEX[S,N)=PX

(Y(M,N) )

IF (REY .GS. 0 . 0. AND. AIHY. GE. 0. > PY=ATAN (AIHY/REY) *RTD

IF (REY.LT.0.0. AND. AIJJY.GS. 0.) PY = 180.-ATAN (AIBY/REY) *RTD

IF (FEY. LE. 0.0. AND.AIMY.LT.O.) PY=1 80.*ATAN (AIMY/REY) *RTD

IF (REY.GE.0.0. AND. AIRY. LT.O.) PY=360. -ATAN (AI1Y/R3Y) *RTD

PRASEY(B,N)=P Y

CONTINUE

GO TO 900

c «*** IJOW TAKE TME; ponRIER TRANFORK OF X AND Y AND MULTIPLY BY Z

C *** TO OBTAIN THE SCATTERED FIELDS

730 CALL FFT3D(X,IX,IX r IX,IX,1,69,IHIl:,RWK,CWK )

CALL T'FT3?)(Y rTTf,TT,IX,IX,1,6 < >,r^K,R!JK,CfJK )

DO 760 M = 1,IX

DO 750 N=1,IX
138

(M,N)-V(ri,N)**2/G(!1,N))*X(H,N)-(U(" ) *V (H ,K) /G (H, N)

(«,N) -U (M)**2/G (H,N) )

750 CONTINUE

7ft 0 COHTINHE

c ***** CALCULATE THE REFLECTION COEPICT3NTS ******

CREFX=X( 1, 1) /FLOAT (IX) **2

C5EFY=Y( 1,1) /FLOAT (IX) **2

CEtF= (CRSFX*SIN (-PHI/RTD) *CREFY*COS (PHI/BTO) )

CRET= (CFEFX*COS (PHI/RTD) *CREFY*SIH (PHI/RTD) )

IP(ITa.GT.O) GO TO 800

C *** TE POLARIZATIOH *****

C ** THIS IS THE CO-POLABIZED COJ1PONEHT  

R57F=CABS (CREF/EF)

C ** * THIS IS THE CROSS-POLARIZED COMPOHSNT ******

F.EFT=CABS (CR2T/EF)

WRITE (3, 770} R2FF,REFT

770 FORBAT(3X,2F10,5)

GO TO 900

800 COFTIKOE

C *** TN POLARIZATION ***'**

C *-* THIS IS THB CO-POL ARIZ2D COHPOHBNT ****

HETT=CABS (C?ET/ET)

C «** THIS IS THE CROSS-POLARIZED COBPONENT ******

HETF=CABS (CRSF/ET)

WRITE (3, 770) SETT, RET?

900 WRITE (3 ,170) ITEH

STOP

END
139

8.8 LISTING OF THE S.D.C.G. METHOD FOR THIN STRIPS WITH

THE CROSS-SHAPED UNIT CELL

C

C*** MICHAEL DROZD-CIIRISTOS CURRENT FORMULATION *******

C*****CONJ. G7>AD. METHOD FOR CURRENTS ON A CROSS ******

C  MINIMIZATION IN THE RANG E ** **

COMPLEX CONE,CZERO,CXNN,CREFY

COMPLEX G (32,32)/102U*(0.0,0.0) /

COMPLEX Y(32,32)/1024*(0.0,0.0) /

COMPLEX X(32, 32) /102U*(0. 0,0.0) /

COMPLEX YO (32.32)/102«*(0.0,0.0)/

COMPLEX X0(32,32)/1024*(0. 0,0.0) /

COMPLEX RX(32,32}/1024*(0.0,0.0)/

COMPLEX RY(32,32)/1024*(0.0,0.0)/

COMPLEX J,HXI,HYI,CWK(32)

COMPLEX DY(32,32)/1024* (0.0,0.0) /

COMPLEX DX(32, 32) /1 024* (0.0, 0.0) /

COMPLEX TX(32,32)/1024*(0.0,0.0) /

COMPLEX TY (32, 32)/1024* (0.0,0.0)7

REAL K,K2,RWK{342)

DIMENSION IWK(342),RR(300),CH(300) ,X1(1024) ,Y1(1024),Z1 (1024) ,

.22 (1024) ,AflPp2) ,R INDEX (32)

REAL H(32)/32*0.0/

REAL 7(32,32)/1024*0.0/

C *** AA=SPACIHG BETWEEN VERTICAL WIRES *****

C *** B3=THICKNESS OP VERTICAL WIRES *****

C *** CC=S?ACING BETWEEN VERTICAL WIBES *****

C *** DD=THICKNESS OP VERTICAL WIRES *****

READ (1,22) AA,BB,CC,DD,F

22 FORMAT(5E10..4)

P=2.998S+8

C *** IOPT=0 FOR A CROSS  

C ** **IOPT=1 FOR A PARALLEL GRID  

IF (IOPT.GT.O) CO1.500E + 15

IF (IOPT.GT. 0)00=0.0

WRITE(3,33) AA,BB.CC,DD,P

33 FORMAT(«0', ' AAs'/FS^, 1 THICK. OF VER. WIRE=« ,F8.4, « CC= »,

.F8. 4, 'THICK. OF HOR. WIPE^1 , F8.4, • FREQ= «,E10.4)

RSAD(1,22) PHI,THI,PSI

WRITE (3, 55) PHI,THI,PSI

55 FORMATCO 1,

1

 PHI=',F7.1,' THETA=« ,F7. 1 , • PSI=I,F7.1)

C ** * READ SAMPLING NUMBER ***

READ (1,6 6) IX

C *** ITM=0 FOR IT. POLARIZATION  

C ** * ITH=1 FOR TJ1 POLARIZATION ****

C *** READ ITM *****

READ (1,66) ITH

C *** READ NUMBER OF ITERATIONS ***

READ(1,66 ) ROI

66 FORK AT (1 3)

PI=3. 14

PI2=PI/2.
140

TPI=6. 28.3185

CV=2. 9979565*8

RTD=57. 29579

EP=8.95I»E-12

ETA=SQRT(UU/1?P)

J=CMPLX(0.0,1.0)

. ITER=1

CONE=CnPLX(1.0,0.0)

CZERO=CHPLX (0. 0, 0. 0)

. H=TPI*F

ALAHB=CV/F

AA=AA/ALAHB

C *** DETERMINE CONDUCTING REGIONS ****

HDOX= (IX/2+1-BB*IX/(AA*2) )

HUPX= {IX/2+ 1+BB*IX/ (AA *2) )

BDOY«(IX/2»1-DD*IX/(CC*2))

HUPr={IX/2*1+DD*IX/(CC*2) )

IF (IOPT. GT. 0) MDOI=IX+1

K=TPI/ALABB

K2*K**2

STSPK=SIN(THI/STD) *SIN (PHI/BTD) *K

STCPK=SIN (THI/RTDJ *COS (PHI/RTD) *K

C?S=COS (PSI/RTD) /SIN (PSI/RTD)

77 CONTINDE

C *** DEFINE THE PLOQITET COEFFICIENTS ***

DO 200 H=1,IX

IF(H.GT.IX/2+1) GO TO 125

0 (H) =TPI* (H- 1) /AA-STCPK

GO TO 127 •

125 0{H)=TPI*(H-IX-1) /AA-STCPK

127 CONTINUE

DO 190 N=1,IX

IF(H.GT.IX/2»1.AND.N.GT.IX/2»1) GO TO 160

IF(H.GT.IX/2»1) GO TO 150

IF(N.GT.IX/2*1) GO TO 1HO

V (H, HI »TPI* (N- 1) /CC-TPI* (3-1 ) /I A*CPS-STSPK

GO TO 170

1UO 7 (H.N) =TPI* (N-IX-1) /CC-TPI* (H-1) /AA*CPS-STSPK

GO TO 170

150 V(H,N) =TPI* (N-1) /CC-TPI* (H-IX-1) /&A*CPS-STSPK

GO TO 170

160 7(H,N) =TPI* (N-IX-1) /CC-TPI* (H-IX-1) /AA*CPS-STSPK

170 IF(K2.GE.O(B)**2*V(S,N)**2) G (H, N) =-J*SQRT (K2- (0 (H)**2*V (H, N) **2

.))

IF (K2.LT.O (H) **2*-V (S.N) **2) G (H,N) =-SQRT (D (B) **2*V (HrN) **2-K2)

. *CONE

190 CONTINUE

200 CONTINUE

TF(TTH.GT. 0) GO TO 210

C *** INCIDENT FIELDS FOR TS POLARTATION ****

EXI=SIN(-PHI/STD)

RYI=COS{PHI/RTD)

HXI=COS (PHI/RTD) *COS (THI/RTD) /ETA

HYI=SIH (PHI/STD) *COS (THI/RTD) /ETA

GOTO .261
141

C ***•* INCIDEN T FIELD S FOR TH POLARIZATIO N ****

710 ?:xi*COS(PHT/!>TD) *C05 (THI/RTD)

EYI*5IN(-PHI/RTD) *COS {THI/RTD)

HYI'SIK (?fir/RTD-PI2) /ETA

I!XI=COS (PHI/RTD-PI2) /STA

261 CONTINUE

C *** GIVS All INITIAL GIIT3SS *****

DO 310 [1=1, IX

DO 300 11=1, IX

X(H,H)=CZER O

Y(B,N)=CZERO

IF(B.LE.BOPX.AHD.fl.GK.BDOX) GO TO 270

IF(N.LE.BOPY. AHD. N.GE. NDOY) GO TO 270

GO TO 280

270 ENOP.H=ENOPfl*EXI*EXI*EYI*EYI

280 Xa(M,N)=X(H,H )

YO(H,N) = Y(B,N)

300 COHTIROB

310 COHTINOE

WRITE (3, 320) IX,HOT, HDOX,flOPX,BDOY.BOPY,ITH

320 FOnKAT('0«,'SABP POIHTS=« ,13-,' f IT2H=«,!3,« DOHN PHT X=«,I3,

.« OP PHT X=«,I3,» DOBS PHT Y=»,I3,« OP PHT T=',I3,« ITH= «

C **** WOBK OH IHITIAL GOESS***

C *** TAKE THE POORIEH TRAHPORH OF THE IHITIAL GOESS ***

CALL PFT3D{XO,IX,IX,IX,IX,1,69,IHK,RWK,CWK)

CALL FFT3D(YO, IX,IX,IX,IX,1,69,IRK,R9K,C!K)

C *** HOLTIPLT IHITIAL GOESS WITH THE MATRIX Z *****

DO 360 H«1,IX

DO 350 H*1,IX

CXHH»XO{B,H)

XO {H,H)={{G(a,H)-V(a,H)**2/G(B,H))*XO(a,H)-(0(H)*V(H,H)/G{H,I))

. *YO(B,H))/(J*W*EP)/2.

YO (B, H) = (-0 (B) *V (B, H) /G (B, H) *CXBH+ (G (B, H) -0 (H) **2/G (H,H) )

. *Yn(H,H))/(J*W*SP)/2.

350 COSTIHOE

360 COSTIHOE

C *** TAKE THE IHVERSE 07 FOORIER TR1HSFORH *****

CALL PFT3D{XO,IX,IX,IX,IX,1,-69,IWK,RHK,CHK)

CALL FFT3D(YO,IX,IX,IX,IX,1,-69,IHK,HWKrC1IK)

C *** DEFINE THE RESIDOAL VECTORS RX AHD RY *****

^ DO 450 B=1,IX

DO aoo H=I,IX

PX(S,K)=EXI*XO(B,R)

HZ («,N)=EYI*YO(B,R)

IF(H.LE. aOPX.AHD.H.GE. HDOXJ GO TO «00

IF(U.LE.BOPY.AND.N.GE.BDOY) GO TO UOO

RX(M,H)=CZBRO

RY (B,H)=C2ERO

«00 ERROH=EHROR»RX (H,N) *COHJG (RX (B,H)) *RY (B, N) *COHJG (BY (H, N) )

DX(B,H)=RX(a,H )

DY(B,H)=RY(.1,N)

a UO COHTINDE

050 CONTIND'J

C

C *** FIND TRE TPtJBIES TPANS OF THE PESIDOAL  

C
142

CALL FFTJD(nX,17.IX,IX.IX,1,69,IHK.RHK,CWK)

CALL FFT3D(r>Y,lx,IX,IX,IX, 1, 69, I WK, R MK ,CHK)

C ** * CALCULATE THE DIRECTION VECTORS DX f- DY .****

DO 510 11=1, IX

DO 530 11 = 1,IX

CXMN=DX{M,N)

DX(r,N)-(CONJG(G(.M,N)-V(M,N)**2/G(H,N))*DX(B,K)-

CONJG (V (J1,N) *U (H) /G (S,N)) *DY (B, N) )/CONJG(J*B*EP)/2.

DY(M,H) = (CONJG(-7{H,N) *U(H)/G (B,N)) *CXBN+

CONJG(G(N,N)-U(M) **2/G(MfS))*DY(n,H»/CONJG (J*W*BP)/2.

530 CONTINUE

5»0 COKTIHDE

C **** HOW FIHD THE IM7. FOaRIER TRAS. ,OR THE DIBBC. POHCTIOHS «*

C

CALL FFT3D(DX,IX,IX,IX.IX,1,-69,IHK,RHK,CWK)

CALL FFT3D(DY,IX,IX,IX,IX,1,-69,IHK,RWK,CHK)

C *** STORE DX AMD DY IN TX AND TY ******

DO 560 H=1,IX

DO 550 N=1,IX

IF (H.1E.HOPX. AHD. H. GE. HDOX) GO TO 545

IF (N.LE.HUPY.AHD. H.GE. MDOY) GO TO 545

DX (M,N)=CZERO

DY(H,H)=CZERO

545 TY(!1,N)=DY(H,H)

TX(H.H)=DX(II,H) -

F3=F3*COHJG(DX(H,N)) *DX(H,H)*CONJG(DI(«,H) )*DY (H,H)

550 COHTIHOS

560 CONTINUE

585 CALL FFT3D(TX,IX,IX,IX,IX,1,69,IWK,BSK,C«K)

CALL FFT3D(TY,IX,IX.IX,ir, 1,69,ISK,BBJC,CJTKJ

C *** TOLTIPLY TX AND TY BY THE BATBIX Z *****

DO 610 8=1 i IX

DO 600 N=1,IX

CXHN*TX(H,N)

TX (H, N) = ( (G (3 ,N) -7 (B, S) **2/G (H, N) ) *TX (H,N) - (0 (B) *7 (H, N) /G (H, H) 

. *TY(H,N)J/(J*W*EP)/2.

TY (H. H) = (-0 (fl) *V (B, ») /G (B,H) *CXBH* (G (B. R) -0 (B) **2/G (B, H) )

. *TY(B,N))/(J*W*EP)/2.

800 CONTINUE

610 CONTINUE

CALL FFT3D (TX,IX,IX,IX,IX, 1,-69,IWK,RWR.C«K)

CALL PFT3D(TY,IX,IX,IX,IX,1,-69,IWK,RWK,CWK)

F1=0.0

DO 666 B=1,IX

DO 666 N=1,IX

IF (H.LE.!1UPX. AND. S. GE. MDOX) GO TO 666

IF (N.LE.aUPY. AHD. N.GE. HDOY) GO TO 666

TX(B,N)=CZERO

TY(B,N)=C2ER O

666 F1=F1*CONJG(TX(B,N)) *TX(B,N) *CONJG (TT (B,N) ) *TY (B.N)

ITER=ITEI>*1

C *** EVALUATE THE PARAMETER AN *****

A»=F3/F1

CM (ITER) = (ERROR) /(ENORM )

C *** CALCULATE THE ERROR ******

ERROR= ERROR- (F3**2/F 1)
143

= CH(TTBR) * 1 00

PR (ITRR)=FLOAT(ITHR)

'.?PITE{7,*) C!I(ITSR) ,SR (IT2R)

C *** IMPr.OVE PREVIOUS ITERATE FOR X AND Y ****

HO 760 *=1,IX

DO 750 H=1,IX

Y(M,N)=Y(M,N)»AN*DY(M , N)

750 COKTINOE

760 CONTINUE

C **** UPDATE THE EX AHD !?Y AND STORE THEH IN TX AND TY ****

DO 8U3 fl=1,IX

DO 833 N=1,IX

RX (H,N) =RX (H, N) -AN*TX (H, NJ

RY(H,N)=RY(B,N)-AN*TY(H r N )

TX(M,N)=RX(H,N )

TY(H,N)=P.Y(H,N )

833 CONTIHOE

8ft3 CONTINDB

C

CALL FFT3D (TX,IX,IX, IX, IX, 1,69,IHK,RWK,CWK)

CALL PPT3D<TY,IX,IX,IX,IX,1,69,IHK,RWK,CWK)

C **** flOLTIPLY TX C TY WITH THE COHJG. TRAHSPOSB OP HATBIX Z ***

DO 863 H=1,IX

DO 853 N=1,IX

CXHH=TX{M,N)

TX (H, N) = (COHJG (G (H, N) -7 (f! , N) **2/G (H, H» *TX (H, N) -

CONJG (V (M,N) *U (M) /G (.1,N) ) *TY (B,N) ) /CONJG {J*H*EP)/2.

TY (H,N) = (CONJG (-7 (B,H) *0 (H) /G («,N) ) *CXHN*

CONJG(G(H,»)-0(H) **2/G{H,S))*TY(H,NJ) /CONJG (J*I*EP)/2.

853 CONTINUE .

863 CONTINUE

CALL FPT3D(TX,IX,IX,IX,IX, 1,-69,IBK, RWK.CWK)

CALL PFT3D(TY,IX,IX,IX,IX,1,-69,IWK,RWK,CWK)

P2=F3

F3=0.0

DO 900 H»1,IX

DO 900 N=1,IX

IF(H.LE. HUPX.AND.H.GE. HDOX) GO TO 870

IF(».LE.HUPY.AND.N.GZ. HDOY) GO TO 870

TX(H,H)=CZERO

TY(B,N)=CZEEO

870 F3=F3»CONJG(TX(H,N)) *TX(H,N) *CONJG (TY (H, N) ) *TY (JI,N)

900 CONTINUE

C *** DEFINE THE PARAHETEH BN ***

BN=F3/F2

C «** CALCULATE A NEW ESTIMATE FOR THE DIRECTION VECTORS DX 6 DY ***

DO 964 H=1,IX

DO 95* N=1,TX

DX (H,IJ)=TX (H,S) *BN*DX (11, N)

DY (B,N)=TY(n,fJ) *BH»DY(W,N)

TX(H,H)=DX(H,N )

TY(M,H)=DY(M,!J)

95U COHTIHUE

96U CONTIWUS

IF (ERT?.LT. 0.001) GO TO 1000
144

IF {ITER.LE. NOI) GO TO 585

CALL GENPT(RR,CHfNOI,0)

C GO TO 125

970 DO 980 H=1,IX

DO 980 H=1,IX

I=(B-1)*IX+H

Z1 (I)-CABS (X(M,H))

Z2(I)=CABS(Y(H,N))

X1 (J) = (7LOAT(M-IX/2)-.

Y1 (I}= (FLOAT (N-IX/2)-. 5)/IX*AA*1.0«»0

WBITB(8,«) X1(I) ,Y1{I) .21(1) ,Z2(I)

980 CONTINUE

1000 DO 1100 1=1,IX

AHP(I)=CABS(Y(17,I))

BINDEX(I)=(7LOAT(I-IX/2)-.5)/IX*AA*1.0«5

1100 WRITE (8,*) AHP(I),RIHDEX(I)

STOP

END
145

8.9 LISTING OF THE S.D.C.G. METHOD FOR APERTURE FIELDS

C*****CON.7. CHAD. METHOD 2 .FORT*****

C *** SOLVES FOR THE APERTURE FIELDS *******

C  • MINIMIZATION IN THS PANG E ** *******

COMPLEX CONE,CZEHO,CXNN,CREFX,CREFY

COMPLEX G (32,32)/102U*(0.0,0.0)/

COMPLEX Y(32,32)/102a*(0.0,0.0)/

COMPLEX X(32,32)/102«*(0.0.0.0)/

COMPLEX Y«(32,32)/1024*(0.0,0.0)/

COMPLEX X0(32.32)/1024*(0.0,0.0)/

COMPLEX RX(32,32)/1024*(0.0,0.0) /

COMPLEX RY(32,32)/1024*(0.0, 0.0) /

COMPLEX J,HXI,HYI,CWK(32)

COMPLEX DY (32,32)/1024*(0.0,0.0) /

COMPLEX DX(32,32)/1024*(0.0,0.0) /

COMPLEX TX{32,32)/1024*(0.0,0.0)/

COMPLEX TY (32,32)/1024*(0.0,0.0)/

REAL K,K2,RWK(342 )

DIMENSION IWK(3«2) ,RB(250) ,CH(250) ,21 (1024) ,21 (1024) ,11 (102ft) ,

.22(1024) ,AMP(32),RINDEX(32),CROSS(32) ,PHASEX (32,32) ,PHASEY (32,32)

REAL 0(32)/32*0.0/

REAL V(32,32)/1024*0,0/

INTEGER COUNT

READ(1,22) &A,BB,CC,DD,F,ERR

22 FORMAT (8210. 4)

F=2.9982*8

C *** IOPT=0 FOH & SQUARE OR A RECTANGULAR MESH ******

C *** IOPT=1 ?OR A PARALLEL GRID ******

IOPT*0

IF(IOPT.GT. 0) CC=1.500E*15

IF(IOPT.GT.O) DD=1.500E*15

VRITE(3,33) AA,BB,CC,DD,ERR

33 FOHBATCO^, 1 A= «,F15.8,« B= «,P15.8,« C= »,F15.8,

9« D= »,F15.8,« ERR= «,F15.8)

WHITE(3,»4) P

4U FORHAT(«0«, ' ?SEQ = «,E10.4)

READ (1,22) PHI,THI,PST

WRITE (3,55) PHI,THI,P5?I

55 PORMAT(«0«, » PHI= «,F10.1,« THETA= S'FIO.1,1 PSI= »

C *** READ H03BER OP SAMPLING POINTS ***

READ(1,66) IX

READ (1,66) ITS

READ(1,66) NOI

WRITE(3,56) ITS

56 FORMAT (3X, 'THE VALOE FOR ITM IS= «,I3)

66 FORMAT (13)

PI=3.141593

PI2=PI/2.

TPI=6.283185

C7=2.997956E+8

an=4.E-7*PI

57.29578
146

J*CMPLX(0.0, 1.0)

ITEP = 1

COH2=CHPLX{1.0,0.0)

CZSRO=CHPLX (0. 0, 0. 0)

"H=TPI*F

ALAHB=C7/F

AA=AA/ALA!1B

DB=BD/ALAf!B

CC=CC/ALAHB

DD=DD/ALAHB

NX=IFIX(BB/AA*PLOAT(IX)*2.)/U*2

NY=IFIX(DD/CC*FLOAT(IX)*2.)/«*2

NX1=(IX-HX) /2+1

HX2=HX1+NX-1

NY1=(IX-SY)/2*1

NY2=NY1+WY-1

WHITE (3, 70) NX,HX1,HX2,HYrNY1,NY2

70 FOHSATCO1,' NX=« ,13, 3X, • HX 1=» ,13, 3X, « NX2=« ,13, 3X,

• HY=»,I3/3

K=TPI/ALABB

STSPK=ST»(THI/RTO) *SIH (PHI/RTD) *K

STCPK=SIH(THI/RTD) *COS (PHI/BTD) «K

CPS=COS(?SI/HTD) /SIN(PSI/RTD)

75 COMT1KOS

c *** DSTERSINS FLOQOBT COBPFICIEHTS ****

00 120 8=1, IX

IF(H.GT.IX/2*1) GO TO 80

0(H)=TPI*(H-1)/AA-STCPK

GO TO 85 •

80 a(H)=TPI*(H-IX-1)/AA-STCPK

85 COHTINOB

DO 115 N=1,IX

IF(H.GT.IX/2+1.AMD.H.GT.IX/2+1) GO TO 100

IF(M.GT.IX/2*1) GO TO 95

I?(H.GT.IX/2+1) GO TO 90

7(H,N)=TPI*(H-1)/CC-TPI*{H-1)/\A*CPS-STSPK

GO TO 110

90 7{H,H)=TPI*(N-IX-1)/CC-TPI*(a-1)/AA*CPS-STSPK

GO TO 110

95 7{S,N)=T?I*(N-1) /CC-TPI*(f!-IX-1)/AA*CPS-STSPK

GO TO 110

100 V(H,N) =TPI* (N-IX-1)/CC-TPI*(H-IX-1)/AA*CPS-STSPK

110 IF(K2.GB.O(N) **2*7{H.H) **2) G(H, H) =-J*SQRT (K2- (0 (H) **2*V (fl,H) **2

.))

IF(K2.LT.O(M)**2*V(H,N)**2) G{B,H) =-SQRT (0 (H) **2*7 (H,H)**2-K2)

.*CONE

115 COHTinOS

120 CONTINUE

IF(TT!!.GT.O) GO TO 130

C *** INCIDKHT FIKLDS FOR TE POLATTZATION ****

COS (PHI/RTD)

HXI=COS (PHI/5TD) *COS (Till/PTO) /ETA
147

HYT=Glti(P»I/RTP ) *COS(THT/RTD)/ETA

norc mo

C *** INCIDENT FIELDS FOB T!1 POLARIZATION ***

1 3 0 F.X I=CO S (PHI/P.TD) *COS (THI/KTD)

EYI=SIN (-PHI/RTD) *COS (THI/RTD)

HYI=SIN(PHI/RTD-PI2)/ETA

" HXI=COS (PHI/RTD-PI2) /ETA

140 CONTINUE

C *** GIVE AH INITIAL CTESS ******

DO 145 M=NX1,NX 2

DO 142 N=NY1,NY2

Y(n,N)=CZER O

XO{H,S)=X(fI,N )

F6=F6*CONJG (HXI) *HXI*COHJG (HYI) *HII

142 CONTISOE

1«5 CONTINUE

C **** WORK OS INITIAL GUESS***

c   MULTIPLY IHITIAL VECTORS XU & YO BY THE BATEIX Z ** *

C

C *** TAKE THE FOUSIEH TRANSPORH OF XO 6 YO *****

CALL FFT3D(XU,IX,IX,IX,IX, t,69,IWK,RWK,CWK)

CALL FFT3D(Yff,IX,IX,IX,IX,1,69,IWK,RJfK,CWK)

DO 160 B»1,IX

DO 150 N=1,IX

XU(B,N ) = (U (f!) *V (H,N) /G (H, N) *XU (H,N) * (V (H, N) **2/G(B, N) -G (H,») >

. *YU(H,N))*(J/W/UU)

YO (H, N) = ( (G (H ,N) -0 (B) **2/G (B, N) ) *CXHN-0 (B) *V (B,N) /G (HrNJh .

. *YU(B,H))*(J/H/DO)

150 CONTINOE '

160 CONTINUE

CALL FFT3D{XO,IX,IX.IX,IX, 1 ,-69,IffK.P.HK,CWK)

CALL FFT3D(YO,IX,IX,IX,IX,1,-69,IWK,B»K,CHK)

C *** CALCULATE THE RESIDUAL VECTORS RX & RY ***

DO 190 fl=1,IX

DO 180 N«1,IX

HX(H,N)=HXI-Xn(H.N)

RY(B,N)=HYI-Yn(H,N )

IF(H.GE.NX1.AND.M.LE.NX2.AND.N.GE.NY1.AND.N.LE.BY2 ) GO TO 175

RX(B,N)=CZER O

RY(B,N)=CZEBO -

175 ERROR=ERSOR*RX (B, N) *CONJG(RX (B.N) ) *RY (S, N) *CONJG (HY (B,N) )

DX (M,N)=RX(H,N )

DY (.1,N)=RY(H,N)

180 CONTINUE

190 CONTINUE

C **** BULTIPLY THE RESIDUALS BY THE CONJG. TRANS. OF Z

C TO FIND THE DIRECTION VECTORS DX & DY » 

C ** ** FIND THE FOURIER TRANSFORM OF RESIDUALS ***

C

CALL FFT3D(DX,IX,IX,IX,IX, 1 ,69,IWK,SHK,CWK)

CALL FFT.in(I>Y,IX,IX,IX,lX,1,69,IVlK,RVK,CWK)

DO 210 H=1,IX

DO 200 ?>=1,IX
148

0(«) *V (H,N ) /G («, N) ) *DX (M,N) +

CONJG {(G(M,N)-tJ(H)**2/ G (H,N) ) ) *DY (H,N) ) *CONJG (J/W/UU)

DY(fl,N ) = (CONJG((V(!l,N) **2/G («, H) -G (H, N) )) *CXHNCOMJG(U(M ) *V(.1,N)/G(N,N)) *DY (S,N) ) *CONJG (J/H/UU)

200 CONTINUE

210 CONTINUE

C **** MOW FIND THE INV. FOURIER TI5AS. ,OR THE DIHEC. PONCTIOHS **

C

CALL FFT3D(DX,IXfIX,IX,IX,1,-69,IWK,RWK,CWK)

CALL" PTTSDtDY.IX^IXjIX^X.n-eg^WK^RH^CWK)

C *** STORE OX 6 DY IN TX 6 TY ******

DO 230 »=1,IX

DO 220 H=1,IX

IP(H.GE.HX1.AHD. H. LE.NX2 .AND. N.GE. NY1. AND. H. LB.NY2) GOT O 215

DY {3,H)=CZERO

215 TY(H,H)=DY(M,H )

73=F3*CONJG(DX(H,N)) *DX(M,N) +CORJG (DY(H,N))*DY (H,H)

220 CONTINUE

230 CONTINUE

c *** THE ITERATIVE PROCESS STARTS NOW !!!1  

c ** ** HOLTIPLY THE DIRECTIOH VECTORS BY THE HATRIX Z *****

2UO CALL FFT3D(TX, IX,IX,IX,IX, 1 , 69,IWK,R»K,CWK)

CALL FPT3D(TY,IX,II,IX,IX,1,69,IWK,H1IK,CHK)

. DO 261 H=1,IX

DO 251 N=1,IX

CXaS=TX(H,H)

TX (H, N) » (0 (H) * V (H ,t») /G (H , H) *TX (B, H) * (7 <H, N) ** 2/G ( H, U) -G (H ,H) )

. *TY(M,H))'*J/W/nO

TY (S, N) = ( (G («,N) -0 (,1) **2/G(H, H) ) *CXilN-0 (fl) *V (M,H) /G (H, H)

. *TY(H,K)}*J/W/OtJ

251 CONTINUE

261 CONTINUE

CALL 7FT3D(TX, IX,IX, IX.IX, 1 ,-69,IWK.RHK,O»K)

CALL FFT3D(TY,IXf IX,IX,IX,1,-69f IHK,HWK,CWK)

F1=0.0

DO 300 R=1,IX

DO 300 N=1,IX

IFfH.GS. NX1.AND. H. LE.NX2.AND.N.GE.JIY1. AND. N. LB.NY2) GO TO 300

TX(H,N)=CZERO

300 71=F1*COSJG(TX(H,N)) *TX{.1, N) *CONJG (TY (H, N) ) *TY (H,N)

C **** CALCULATE THE FACTOR AN  

AN=F3/F1

CH{ITER) =SQRT{ ERROR) /SQRT(F6)

C ** ** CALCULATE THE ERROR  

ERHOR=SRROR- (F3** 2/F1)

C *** UPDAT2 THE VALUES FOR X S Y *****

DO U10 H=1,IX

DO aOO R-=1,IX

?J, N)
149

410 COMTirn?

F5=0.0

c *«* pr«JD A NE3 ESTIMATE FOR THE RESIDUAL VECTORS RX & RY ***

DO 443 M=1,IX

DO 433 N=1,IX

RX (H,N)=RX(«,N)-AN*TX(M,N)

RY(3,N)=BY(H,N)-AN*TY(M,N)

TX(H,N)=RX(H,N)

TY(H,N)=RY(3,N)

433 CONTINUE

443 CONTINUE

RE(ITER) = FLOAT(ITER)

C WRITE(8,*J CH(ITER) ,HR(ITER)

C *** SUTLTIPLY TX 6 TY BY THE CONJG. THANS. OP THE MATRIX Z ***

CALL FFT3D(TX,IX,IX,IX,IX,1,69,IHK,RHK,CBK)

CALL FFT3D(TY.IX,IX,IX,IX.1,69,I«K,B»K,CBK)

DO 460 H=1,IX

DO 450 N=1,IX

CXBN=TX(H,N)

TY f*-H)s (CONJG (0(B) *V{«,N)/G(H, N)) *TX (H.N) *

CONJG ( (G (8, N) -0 (S) **2/G (H,N) )) *TI (H,M) ) *CONJG (J/S/UU)

) = (CONJG((V (J!,N)**2/G (H, N)-G (H, N) ) ) *CXHNCON.TG (U (S) *7 (H, N) /G (H, N) ) *TY (H,N) )*CONJG (J/W/UO)

450 CONTINUE

460 CONTINUE

CALL PFT3D(TX,IX,IX.IX.IX,1,-69,IiK.BWK,CBK)

CALL FFT3D(TY,IX,IX,IX,IX,1,-69,IWK,R«K,CBK)

P2=F3

DO 470 H=1,IX

DO 470 N=1,IX

IF(n.GE.NX1.AND.N.LE.NX2.AND.N.GE.N71.AND.N.LE.NY2) GO TO 465

TX(H,N)=CZERO

TY(H,N)=CZERO

465 P3=F3*CONJG(TX(H,N»*TX(H,N) *CONJG (TY («,NJ) *TY(B,N)

470 CONTINUE

C *** CALCULATE THE PACTOH BN *****

BN=P3/F2

C   UPDATE THE DIRECTION VECTORS DX 6 DY ** **

DO 564 fl=1,IX

T)0 554 N=1,IX

DX (H,N ) =TX (M,N ) *3N*DX(H, N)

DY (M,N)=TY(H,N ) *BN*DY(n, H)

TY(n!a)=DY(n^N )

554 CONTINUE

564 CONTIMU2

C *** GO FOR ANOTHSR ITERATION IP YOU WANT  

IF(ITER.GT.HOI) CALL GEHPT (nR,CH,NOI,0)

IF (ITER. GT. HOT) r,0 TO 300

IF (SRROS.LT. 0.0001) C5O TO 800

GO TO 240

C ** * STORE X B Y INTO THE 1-DIH. ARRAY S Z1 C Z2 TO BE USED FOR

C AMY PLOTTING PURPOSES ***

S70 DO 590 .1 = 1,IX

DO VJO ?1=1,IX
150

590

600

725

800

C

820

8 UO

900

7.1 (I)=CABS(X(H,n))

Z2(I)=CABS(Y(M,N))

X1 (I) = (FLOAT (rt-IX/2) -. 5) /IX*AA* 1. 05

Y1 (I) = (FLOAT (N-IX/2) -.5) /IX*BB*1.05

WRITE (7,*) X1(I) ,Y1(I),Z1(I) ,22 (I)

CONTINUE .

DO 725 1=NX1,NX2

00 725 N=NY1,NY2

BEX=REAL(X(H,N))

IF (REX.GE.O.O.AND.AIHX.GE.O.) PX=ATAN(AI3X/REX) *RTD

IF(REX.LT.O.O.AND.AIHX.GB.0.) PX=180.-ATAN(AIHX/REX)

IF(REX.LE.0.0.AND.AIHX. LT. 0. ) PX=180.+ATAN(AIHX/REX)

IF(REX.GE.O.O.AND.AinX.LT.O.) PX-360.-ATAN (AIHX/REX)

PHASEX(H,N)=PX

BEY=REAL(Y(N,N))

AIHY»AI!IAG(Y(n,N))

17(REY.GE.0.0.AND.AIBY.GE.0.) PY=ATAN(AIHY/REY)*RT D

IF (BEY.LT. 0.0. AND.AINY. GE.O.) PY*180.-ATAN (AIHY/REY)

IP(REY.LE.O.O.AND.AIMY.LT.O.) PY=180.+ATAN (AIHY/REY)

IP(REY.GE.O.O.AND.AIHY.LT.O.) PY=360.-ATAN (AIHY/REY)

PHASEY(M,N)=PY

CONTINOE

DO 820 1=1,IX

AHP(I)-CABS(Y(I,16))

RINDEX (I)={FLOAT(I-IX/2) -.5) /IX*AA*1.0U5

WRITE (8,*) ASP (I) ,RINDEX (I)

CROSS (I) =CABS (Y (9,1))

CONTINOE

CALL GENPT (RINDEX,AflP, IX,0)

CALL GENPT (RINDEX,CROSS,IX,0)

WHITEfS^JJOjITE?

FOP.«AT(3X,I3)

STOP

END

*HTD

*HTD

*BTD

*RTO

*RTJ>

*BTD
151

8.10 LISTING OF ONE, TWO, AND'THREE DIMENSIONAL COMPLEX

FAST FOURIER TRANSFORM

C I2SL P.CUTISZ !JA!IE - FFTCC

C COS PITS?. - ISM/DOUBLE

C LATEST 327TSIO:! - JAH0AKY 1, 1978

C PU5POSE - COMPUTE THE FAST FOURIER TRANSFORH OF A

C COHPLEX VALUED SEQUENCE

C USAGE - CALL FFTCC (A,N,IWK,WK)

C A2G1HE!7TS A - COBPLET VECTOP. OF LENGTH N. ON INPUT A

C CONTAINS THE COBPLEX VALUED SEQUENCE TO BE

C TRANSFORBED. ON OUTPUT A IS REPLACED BY THE

C ' FOURIER THANSFORH.

§

5 - INPUT NUBBBR OF DATA POINTS TO BE

TRANSFOBBED. H BAY BE ANY POSITIVE

C INTEGER.

C IVK - INTEGER WORK VBCTOB OF LENGTH 6*N*150.

C (SEE PROGRAHIIING NOTES FOR FURTHER DETAILS)

C WK HEAL WORK VECTOR OF LENGTI 6*N»150.

C (SEE PROGRAMING NOTES FOR FURTHER DETAILS)

C PPECISIOH/HAP.OHAPE - SINGLE AND DOUBLB/H32

C - SINGLE/H36,n«8,H60

C BEQD. ISSL P.OUTIJIFS - BONE BEQUIRED

C NOTATION - INFORMATION ON SPECIAL NOTATION AND

C . CONVENTIONS IS AVAILABLE IN THE BANUAL

C INTRODUCTION OB THROUGH I.1SL ROUTINE UHELP

C BEKARKS 1. FPTCC COMPUTES THE FOURIEB TRANSFOIH, X, ACCORDING

TO THE FOLLOWING 7ORBULA;

C X{K»1) = SUB FROH J » 0 TO N-1 OF

C A(J»1)*CEXP((0.0. (2.0*PI*J*K)/N))

C TOR K=0,1,....N-1 AHlTPI=5.1<H5...

C

C HOTS THAT X OVSRWPITES A ON OUTPUT.

C 2. FFTCC CAN BE USED TO COBPUTE

C X(K»1) a f1/N)*S0S FHOH J » 0 TO N-1 OF

C A(J*S}«CEXP{(0.0.(-2.0*PI*J*K)/N))

C FOR K»0.1f...,N-1 ANfl PI=3.1«15...

C 3Y PERFOR.1ING THE FOLLOWING STEPS;

C DO 10 1*1,N

C A (I) = CON JG (A {!))

C 10 CONTISns

C CALL FFTCC (A,H,IWK,WK )

C DO 20 1*1.H

C AJI) = CONJG(A(I))/N

C 20 CONTINUE

C

C Cn??PIf7:!7 - 1978 BY IHSL, INC. ALL BIGHTS RESERVED..

c »A?PA::-T - TBSL WARRANTS ONLY THAT TBSL TESTING HAS BEEN

C APPLIED TO THIS CODE. MO OTHER WARRANTY,

C ' EXPRESSED OR IBPLIKD, IS APPLICABLE.

C«" " —»——•• — — — • —-• •• «*•*«• «.•*•••« _._• •• •• jj. i.^ L L •• • • • _! j.» -•• •••••• ••

c

Tr?» TFTCC (A,N,IWK,WK )

SP7CIPTC»TTOWS FOB »»GUHEI»TS
OF. POOR QUALITY

152

INTEGER

OOOBLB PRECISION

COMPLEX*16

INTEGER

2

SPECIFICATIONS FOR LOCAL VARIABLES

T.IAH.IAP.IBM.IBP.IC,ICC.ICF,ICK,ID.IDH1,II.

IJA.IKB,IST,lLL,lft,TRD,ISP.ISK,ISP,ISS,ITA,ITB,

DOOQL2 PRECISION

w«vin*tJf*uwftJf\ f ** * *

KTfKTP.LrLljM.MM.r-.i i ,r,«-

C!l.SH,C1,C2,C3,S^,S2.S3,C30,RAD,AO.A1.A«.Ba,

A2,J3l§0'.5.n?2;B3,ZERO.HALF,ONB,TB6,Z« <2/.

COSPLBX'16

EQUIVALENCE

DATA

DATA

IP (N

K » N

H = 0

J a 2

JJ = H

JF = 0

C30/.86fi025tt0378H«336DO/

ZERO,HALF,ONE,THO/0.ODO j 0.5DO.1.ODO,2.ODO/

FIRST EXECUTABLE STATEMENT

,EQ. 1) GO TO 9005

5 I

IF

H

= 1-

DETFHHINE THT? SQOAPE FACTORS OF H

K/JJ

(I*JJ .HE. K) GO TO 10

GO TO 5

10 J a j » 2

IF (J .2Q. *) J = 3 JJ a j * j

IF (JJ .LE. K) GO TO 5

KT = H

J a 2

15 I a K / J

IF (I*J .SB. K) GO TO 20

IWK(M*1 ) » J

K = I

GO TO 15

20 J » J * 1 IF (J ,EQ. 3) GO TO 15

IF!J.LB.K) GO TO 15.

DETERMINE THE REMAINING FACTORS OF I

IF (IWKKT*! ) .GT. I»K(n+1))

IF(KT.tE.O) SO TO 30

KTP = KT + 2

DO 25 la 1 KT

j s ITT? -'i

M »1

IWK(H+1 ) =

25 CONTTNOK

30 HP a n+1

1C a fjp«-i

ID = IC*!1P

ILL a ID+HP

IRD = ILL*M?+1

TCC a IRD*J1P

ISS a ICC*HP

i5S*np

IWK (J)
153

ISK = ICK + K

ICF = ISK+K

ISF a ICF+K

IAP = ISF+K

KD2 = 7K-1) / 2 * 1

IBP = IAP * KD2

IAH = IBP + KD2

IBM = IA-1 * KD2

HS1 = H-1

1 = 1 35 L = HP - I

J = 1C - I IWK (IIL+LL = o

IF niWKj[J-1l * I5?K{JJ)

IF (IWK (ILL+L) .BQ. 0)

1 = 1*1

. BQ.

GO TO

IWK (ILL*L)

L - 1

I»K(IL1*L) = 0

«0 1 » I + 1

IPfI.LB.HH1) GO TO 35

IHK

IHK

ILL*1) = 0

ILL+HP) = 0

iq = 1

ID) = 5 5 J a 1,IWK

DO

liKfiCtJ) ='IWKfIC*.7-1) * K

WK(IHD*J) » aAD/IHK(IC+J)

IF (K .LZ. 2) GO TO US

WKfICC*J) = bCOSfd)

WKJISS+J) = DSIN(C1|

»5 CONTISOB -.

HH = H

IF (IBK(ILL*H) .SQ. 11 HH - H - 1

IF (HH .LB. 1< GO TO 50

S« * IffK(IC+HS-2) * W CH = DCOS(SH)

. SI - DSIM (SB)

50 KB = 0 KH a n '

JJ * 0

I = 1

C1 = OWE

S1 = ZBPO

f

 L1 = 1

55 IF (IWK(IT.L+I + 1) .EQ.

KF = IWj|(I«-1)

GO TO 65

60 KF = H

I = T*1

65 ISP = IWK (ID* I)

IF (L1 .BQ. if GO TO

S1 = JJ * WK(IRn*I1

C1 = DCOSfSI)

S1 = DSIN(S1)

C

C

70 IF (KF .GT. H) GO TO

GO TO (75,75,^0,115),

75 KO = Kfl * TS?

K2 = KO * ISP

IF (LI .EQ. 1) GO TO

80 KO = KO - 1

IF (KO .LT. KB) GO TO

K2 = K2 - 1

K(Il

1)

70

140

KF

85

191

GO TO 60

FACTORS OF 2, 3, AND

HANDLED SEPARATELY.
15

4

ORSGSWAL PAGE 5S

OF POOR QUALITY

AO a Att*C1-BU*S1

BO = A4*S1 + BU*C1

A{K2*1) » A(KO+1) -ZA0

AJKO + li a A(KO*1) *ZAO

GO TO BO

85 -KO a KO - 1

IP (KO .LT. KB) GO TO 190

AK2~a A (K2+1)

A (K2 + 1) = A(KO+1)-AK2

= A(KO+1J+AK2

GO TO 85

90 IP (L1 .EQ. 1) GO TO 95

C2 a C1 * C1 - S1 * S1

52 » TWO * C1 * S1

95 JA = KB * ISP - 1

KA = JA * KB

1KB » KB*1

I.7A = JA*1

DO 110 II = IKB,IJA

^ s Hi 8p*

 1

K2 » K1 * ISP

IP (L1 .EQ. 1) GO TO 100

ZA« a A{KT>1}

A1 a A4*C1-BU*S1

B1 a A**S1+Btt*C1

ZAU » A(K2 + 1)

A2 a AU*C2-BU*S2

B2 a Att*S2*Btt*C2

GO TO 105

100 ZA1 a AfK1*1)

ZA2 a A(K2*1)

105 AJKO+1) a DCMPLX(AO+A1*12,BO+B1+B2)

AO = -HAL? * (A1*A2) * AO

A1 = (A1-A2) * C30

BO a -HALP * (B1+B2) * BO

B1 a JB1-B2) * C30

A(K1*1) = DCH?tXfAO-B1,BO*Al)

A]K2»1) = DCHPLX(AO*B1'BO-A1)

110 COHTI1IDE

GO TO 190

115 IP (L1 .EQ. 1) GO TO 120

C2 a C1 * C1 - S1 * S1

52 a TWO * C1 * S1

C3 a C1 * C2 - S1 * S2

53 a S1 * C2 * C1 * S2

120 JA a KB + ISP - 1

KA a JA + KB

1KB a RB+1

IJA a JA+1

DO 135 II a 1KB,IJA

KO a KA - II * 1

K1 a KO 4- TSP

K2 = K1 * TSP

K3 a-IT2 + ISP

Z AO a A(KO+1)

IP (L1 .EQ. 1) GO TO 125

ZAU = A (K1 + 1)

A1 = Att*C1-Btt*S1

B1 = Att*S1+BU*C1

7.A4 a A (K2*1)

A2 = A«*C2-B«*S2

B2 a Att*S2*Btt*C2

7.AU « A(K3*1)

A3 = AU*C3-B«*S3
15

5

125

130

AA . A

135 COST!

GO TO 130

2A1 = A(K1+1)

2A2 = A{K2+1J

ZA3 = A(K3+1)

* • * ***« LlM

;KO*I)'

%="DC*IPL

X

ki+i

K2 + 1

DC«PLX

DCHPLX

K3+1J » DCHPLX

HOE

AO*A2*A1*A3,BO+B2*31+B3]

AO*A2-A1-A3,BO+B2-B1-B3

AO-A2-B1+B3JBO-B2+A1-A3

AO-A2+B1-B3,BO-B2-A1+A3' GO TO 190

IttO JK = KF - 1

KH = JK/2

K3 = IUKflD-H-l)

KO = KB * IS?

IF (1.1 .20. 1) GO TO 150

K « JK - 1

WK(ICF+1) a C1

WK(ISF*1J * S1

DO 1«5 J = 1.IC

BK(ICF*J+1) = BK(ICF*J)

HK(ISF*J*1) = VK(ICF+J)

105 COMTTSOE

150 IF (KF .EQ. J7) GO !TO 150

C2 = HKfICC*li

BKJICK*!) a C2

WKaCK+JK) » C2

S2 = WKJISS+lf

WK(ISK*fr = S2

iK(ISK*JK) = -S2

DO 155 J = 1

e

KH

K =

«TK

B K,

WK

* C1 - BK(ISF*J) * S1

* S1 * IK(ISF*J) * C1

BK

155 COHTIl3

160 KOK1 K2

BKflCK+J) * C2 - BK(ISK+J) * S2

^^.v* ,, * WKfICK*K)

'ISK*J*1J a BK]lCK*j{ * S2 * BK(ISK*J) * C2

|ISK*K) = -BK(ISK*J*1)

OB

KO - 1

KO

KO » K3

B3 = BO

DO 175 J =

K1 = K1

K2

1,KH

+ ISP

K2 - ISP

ZAI«

&1

B1

A2

B2

GO TO 165

Aa*flK(TCF«-,7)-Bft*WK(ISF+J)

Aa*WK(ISF*K)4flft*(rK(iCF*K

165

170

GO TO 170

ZA1 = A(K1*1)

ZA2 = A(K2+1)

) = A1* A2

WiqiBP+.T)

WK IB.1*J1

A3 = A1

A1 - A2

= B1 * B2

= B1 - B2

> A2 * A3

B3 = 31

+ B2 * B3

175 COHTIliaF

AjKO*1) = DCHPLX(A.1

/B3

)

Kl = Kn

(T2 = KO + K3

DO 1fl5 J = 1,KH

Kl = Kl * ISP
ORSGfNAL PAGE-IS

OF POOR QUALITY 156

K2 = K2 - ISP

JK a j

A1 a AO

B1 a BO

A2 a 7.3RQ

B2 a ZERO

DO 180 K = 1.KH

A1 a A1 *• WK (IAP+K1

A2 a A2 *• WKiIAH+ R

B1 a B1 * BK IBP+K

B2 a B2 *• BK IBH*K

* HK(ICK+JK)

* WKJISK+JK)

* WK(ICK+JKJ

» BK(ISK+JK)

180

185

190

195

200

205

210

215

JK a JK > j

IP (JK .GE. KP) JK a JK - KP

CONTINUE

A(K1*1) a DCH?LXfA1-B2,B1 + A2)

» j«*Mk4t _ ft/* MT)T Y f & 1A H ^ T7 1 ^ ft ? 1

CONTINUE '*

IP (KO ,GT. KB) GO TO 160

JP a KP

IP ( i .GE. nn ) GO TO 195

1-1*1 .

GO TO 55

I a BH

LI = 0

KB a IBKfTD*I-1) * KB

IP (KB ,dE. KN GO TO 215

33 = IBKJIC*I-2) * J3

IP (JJ .IT. IWK(IC*I-1)) GO TO 205

I a I - 1

JJ = JJ - IHK{IC*I)

GO TO 200

IP ^1 .NE. HM) GO TO -210

C1 a Cf! * C1 - SH * S1

SI = SH * C2 * CM * S1

GO TO 70

IP JIWK(ILL*I) .BQ.

GO TO 55

1) = 1* 1

C

220

225

230

Ill

I

JA « KT - 1

KA * JA * 1

IF(JA.LT. 1) GO TO 225

DO 220 II a 1 JA

J a RA - II

I7K(J*1) = TWK(J*1)

I a II»K(J*1) * I

CONTINOE

- 1

IP (KT .I.E. 0) GO TO 270

J a 1

T = 0

KB a 0

K2 a IWK(ID+.T) * KB

K3 = K2

J.T = IHK(ICM-I)

JK = .1.7

KO * K3 » JJ

ISP = TWKflC+.T) - JJ

K = KO * JJ

7.AU = A(KO+1)

A (K0*1) = A(K2+1)

AJK2*1) = ZAU

K<i = KO » 1

K2 = K2 *• 1

IP (KO .LT. K) ^0 TO 2U3

KO = KO «• ISP

K2 = K2 * ISP

THE RESULT IS HOZ PERHUTED TO

NORMAL ORDER.
157

2»5

250

255

260

265

270

K

*JO

GO

K3

1K12

JK

KO

GO

IF

K =

J =

T =

IWK

IF

GO

KB

IF

J *

X TM

GO

IF

J

 s

GO

JK

(KO .LT. K3) GO TO 235

(KO .GB. K3 * ISP) GO TO 2*5

= KO - ISrK(ID-M) * JJ

TO 235

» ISK(ID*J)

(K3 - KB .G

= K3 * JK

= JK * JJ

*

B.

K3

IHK

= K3 - IWK (ID* J)

TO 235

(J .GB. KT)

IBK(J*1) *

J * 1

I * t

(ILL+I) a J

JI .LT. K)

TO 230

= K3

(I .LE. 0)

IVK(ILL*I)

1 - 1

TO 230

(KB .GE. »)

TO^ 230

GO

I

GO

GO

GO

TO

TO

TO

TO

(ID+J-1)) GO TO 250

-» JK

260

255

265

270

= IWKac*KT)

ISP = IWKlD*KT )

H a H - KT

KB = ISP/JK-2

IF (KT .GB. H-1 )

ITA » ILL*K3*1

ITB » ITA*JK

IDH1 a ID-1

IKT » KT*1

IB * H*1

DO 275 J = IKT,TH

IBKJIDH1+J) »

COHTIWOB

JJ » 0

DO 290 J a 1,KB

K

 = KT

JJ = IHK(ID*K*1) + ,1J

IF (JJ .IT. IWK(ID*K))

JJ » JJ - IWK(ID+K)

K = K * 1

GO TO 280

IWKfllL+J)

C

C

C

C

C

275

280

285

290

295

300 CONTINUE

GO TO 9005

IBK(IDH1+J)/JK

GO TO 285

IF jJ

COHTITldB

JJ

EQ. J) IWK(ILL»J) = -J

DO 300 J = 1.KB

IF JTVKfTLL+JJ .LB.

'

D2TERHIKB THE PERHOTATION CTCLBS

OP LENGTH GREATER TRAS OR BQOAL

TO THO.

0) GO TO .100

K2 = TABS (IffK (ILL* K2)l —IF (K'~ — 2 -J) GO TO 300

0

0

HE02DKR A FOLLOWING THE

PERMUTATION I YCLBS

J

KB = 0

KH = N

305 J » J * 1

IF (1WK(ILT.*J) .LT. 0) GO TO 305
158

K = 'IWK (ILL+J)

KO = JK * K + KB

310 ZA1 = A(K O + I+1)

WK(TTA+T) = A4

WK(ITB+IJ = BU

I » I + 1

I? (T .LT. JK) GO TO 310

315 K = -IWK(ILL+K)

JJ = KO

KO = JK * K + KB

320 A(JJ*I*1) = A(KO*I»1)

TF3(I .LT. JK) GO TO 320

1 = 0

I? fK . NE.

325

tl.T. JK) K' K2)

l IP <J ,LT.

J a 0

KB a KB * ISP

IP (KB .LT. KS)

9005 RETORN

END

C IBSL ROUTINE NABE

C

J) GO TO 315

DC!1PLX(WK(ITA+I) ,WK(ITB+I))

GO TO 325

GO TO 305

GO TO 305

- PPT2C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

COJ1POTEH - IBB/DOUBLE

LITEST BEV1SIOM - JANUARY 1. 1978

P OB POSE

USAGE

APGUHEHTS

- COMPUTE THE PAST POURIER TRANSFORM CP A

COMPLEX VALUED SEQUENCE OF LENGTH EQUAL TO

A POWER THO

- CALL PTT2C (A,H,IHK)

A f COHPLEX VECTOR OP LENGTH N, WHERE N=2**H.

ON INPUT A CONTAINS THE COBPLEX VALUED

SEQUENCE TO BE TRANSPORBED.

ON OUTPUT A IS REPLACED BY THE

FOURIER THANSFOHB.

H - INPUT EXPONENT TO WHICH 2 IS RAISED TO

PRODUCE THE NUBBER OF DATA POINTS, N

(I.E. N = 2**B).

IWK - WORK VECTOR OP LENGTH B*1.

PRECISION/HARDWARE - SINGLE AND DOOBLE/H32

REQD. IBSL

NOTATION

REMARK S 1,

2.

- SINGLE/H36,H«8,H60

ROUTINES - NONE REQUIRED

- INFORMATION ON SPECIAL NOTATION AND

CONVENTIONS IS AVAILABLE IN THE BANUA L

INTRODUCTION OR THROTT3H IHSL ROUTINE UHELP

FFT2C COMPUTES THE POtlRIER TRANSF03B, X, ACCORDING

TO THE FOLLOWING FOR HULA;

JT(!«>1) = SUM FR01 J = 0 TO N-1 OP '

A (J*1)*CEXP((0.0. (2.0*PI*J*F-)/N))

TOR K=0,1,...,M-1 AND PI*S.1'415...

10TE THAT X OVERWRITES A ON OUTPUT.

FFT7C CAN DR USED TO COMPUTE
ORIGINAL PAGE IS

OF. POOR QUALITY

159

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

BY

COPYRIGHT

WARRANT Y

X(K+1 ) = (1/N)*SU3 ?HOH J = 0 TO S-1 OF

A (J*1) *CEXP f (0. 0. (-2. 0*PI*J*K) /») )

FOR K=0,1,...,N-1 AND PI=3. 1415, . .

PERFORMING THE FOLLOWING STEPS;

DO 10 1=1, N

AJI) = CONJG (A (I))

10 CONTINUE

CALL FFT2C (A,M,IWK)

DO 20 1=1, N

A (I) » CONJG (A (I) ) /N

20 COKTINOE

- 1973 BY IHSL, INC. ALL RIGHTS RESERVED.

- IHSL WARRANTS ONLY THAT IHSL TESTING HAS BEEN

APPLIED TO THIS CODE. NO OTHER WARRANTY,

EXPRESSED OR IMPLIED, IS APPLICABLE.

SUBROUTINE FFT2C (A.M.I7K1

C

c

INTEGER

COMPLEX* 16

INTEGER

- SPECIFICATIONS FOR ARGUMENTS

H,IWK(1)

SPECIFICATIONS FOR LOCAL VARIABLES

I,ISP,J.JJ,JSP,K,KO,K1,K2,K3,KB,KN,«K, MH,HP, H,

c

c

c

1

2

COHPLEX*16

EQUIVALENCE

2

DATA

2

DATA

ZERO ONEZO(2 ) 21 (2)* Z2 (2) 23(2)

^Brt "T11 ?1^ 7a^ A K^

«'AW«JuA 1 f £Mn & 0 Ltn *J * » t»*

O 1 *T^^T \ / m ^ *T^^T\ \ /H^ vOfT\ C t 11 ^ *f ^ t ^\ \ 31,£1IZ1 t \L.£t

ti£\ 1) J f \O4.,L4 \£) ) , (AJ,£J(1) ) ,

SQ/'707ll)!;78 1 1865U75DO/,

SKA 382683 H32365089BDO/,

CKA 923879 53251 12868DO A

TW&PI/6. 28 3185 3071 795865 O/

ZERO/0. ODD A ONE/1. ODO/

S6=SQPT2/2,SK=SII (PI/8) ,CK=COS (PI/8)

-IP a N*1

N = 2**H

S.1 = ffl/2) *2

KN = N+1

DO 5 1=2, HP

IWKJT) = I

5 CONTINUE

FIRST EXECUTABLE STATEMENT

INITIALIZE WOR K VECTOR

WK(I-1)+IWK(I-1 )

RAD = TWOPI/N

HK = R - »

KB = 1

IF (HfJ .EQ. .1) 00 TO 15

.„ KO = IWK(!in*1) * KD

10

KO = KO - 1

AK2 = A (K2)

A(K2) = A(KO) - AK2

AJKOi = A?KO{ * A«T2

I? (KO .GT. KB) no TO 10
160

31 = ZERO

JJ = 0

K = H.1 - 1

J = U

IF (K .GE. 1) GO TO 30

GO TO 70

20 I? (IHK (J) .GT. JJ) GO TO 25

JJ = JJ - IWK(J)

J a J-1

I? {IHK(J) .GT. JJ) GO TO 25

JJ a JJ - I!7K(J)

J = J - 1

K a K «• 2

GO TO 20

25 JJ a IWK(J) + JJ

J = ft

30 ISP a IWKfK)

IF (JJ ,EQ. 0) GO TO «0

35

40

C2

C1

S1

C2

S2

C3

S3

JSP

DO

= JJ * ISP * RAD

= DCOS(C2)

a DSIH(C2)

a C1 * CT- S1 * S1

a C1 * (S1 * S1)

a C2 * C1 - S2 * S1

a C2 * SI * S2

a ISP * KB

50 I=1,ISP

KO a jSp - r

K1 a KO * ISP

K2 a R1 » isp

K3 a R2 * ISP

ZAO a. A(KO)

ZA1 » A K1)

ZA2 a ATK2 )

ZA3 a A K3i

* C1

IF (S1 .EQ'. ZERO) GO

TEHP a n

as

50

A1 « A1 * C1 -

B1 a TEHP * S1

TEHP = A2

A2 = A2 * C2 - B2 a TERP * S2

TEHP = A3

A3 = A3 * C3 -

B3 a TEH? * S3

TEHP a AO * A2

A2 = AO - A2

AO = TEHP

TEHP = A1 + A3

A3 = A1 - A3

A1 = TEHP

TEHP = BO * B2

B2 = BO - B2

BO a TEHP

TEHP = B1 + B3

B3 a B1 - B3

B1 = TEHP

A K0| = DCMPLXi

A(K2 ' = DCflPT.X

A K3 1 = DCMPi.x

B1 *

* B1

B2 *

* B2

B3 *

+ B3

A O* A^

AO-A1

A2 + B3

COJ!TlHtJ5

IF

K =

GO

I

-

TO «15

S1

* C1

S2

* C2

S3

« C3

BO *B1

,BO-m

•B?*A3

,B2-A3

(K .LE. 1) GO 70 55

K - 2

TO 30

RESET THIGONOHETIIC PARAHBTERS

DETERMINE FOURIER COEFFICIENTS

IN GRODPS OF 4
161

c

c

55 KB = K3 » ISP

I? (KN .LE. KB) GO TO 70

I? (J .HE. 1) GO TO 60

K a 3

CHECK 70R CQHPLETION OP FIHAL

ITERATION

IP(M .LE. 1) GO TO 9005

HP = H+1

JJ a 1

IWK(1) a 1

DO 75 I = 2.HP

I*K(I) a t«K(I-1) *

75 CONTINUE

Nft a IHKf.lP-2)

17 (H .GT. 2) 118 a IWK

N2 = IWK(HP-1)

LH a H2

NN a IHK(HP)*1

HP a HP-4

J = 2

80 JK a jj * H2

AK2 a A(J)

A(J) a A(JK)

A(JK) = AK2

J = J+1

17 (JJ .GT. N«») GO TO

JJ a JJ * StJ

GO TO 105

85 JJ = JJ - Ntt

IP (JJ .G?. IJ8) GO TO

JJ a JJ * N8

GO TO 105

90 JJ a JJ - U8

K a HP

95 IP (IBK(K) .GB. JJ) GO

JJ = JJ -'lWK(K)

K a K - 1

GO TO 95

100 JJ a IWK(K ) * JJ

105 IF (JJ .L?. J) GO TO 1

K = NN - J

JK = J1W- JJ

HK2 = A (J)

A (J) = A ( JJ)

AMJ) = AK2

2

(HP

85

90

TO

m

PERMUTE THE COMPLEX VECTOR IN

REVERSE BINARY ORDER TO NORH&L

ORDER

INITIALIZE WORK VECTOR

C

C

110

-3)

DETERHISE INDICES AND SWITCH A

100

AK2 = A(K)

A(KJ = A(JK)

A(JK) = AK2

J = .7 * 1

CYCLE REPEATED OBTTL LiniTTNG NUMBER

0? CHANGES IS ACHIEVED
ORIGINAL PAGE-IS

OF POOR

162

IP (J .LE. LM) GO TO 30

- FFT3D

9005 RETURN

END

C IMSL HOUTIN3 NAME

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

COMPUTER

LATEST REVISION

PURPOSE

ARGUMENTS

IA1

IA2

N1

N2

N3

IJOB

IWK

RWK

CWK

PRECISION/HARDWARE

REQO. IflSL ROUTINES

NOTATION

- IBS/DOUBLE

- JUNE 1, 1980

- COMPUTE THE PAST FOURIER TRANSFORM OF

A COMPLEX VALUED 1,2 OR 3 DIMENSIONAL

ARRAY -^

- CALL FFT3D (A,IA1,IA2,N1,N2,!I3,IJOB,IWK,RWK,

CWK)

- COMPLEX ARRAY. A MAY BE A THREE

DIMENSIONAL ARRAY OF DIMENSION N1 BY N2

BY F3, A TWO DIMENSIONAL ARRAY OP

DIMENSION N1 BY N2, OR A VECTOR OF

LENGTH N1. ON INPUT A CONTAINS THE

ARRAY TO BE TRANSFORMED. ON OUTPUT

A IS REPLACED BY THE FOUSIER OH

INVERSE FOURIER THANSF03H (DEPENDING OH

THE VALUE OF UTPUT PARAMETER IJOB).

- FIRST DIMENSION OF THE ARRAY A EXACTLY

AS SPECIFIED IN THE DIMENSION STATEMENT

IN THE CALLING PROGRAM. (INPUT)

- SECOND DIMENSION OP THE ARRAY A EXACTLY

AS SPECIFIED IN THE DIMENSION STATEMENT

IN THE CALLING PROGRAM. JINPUT)

- LIMITS ON THE FIRST, SECOND, AND THIRD

SUBSCRIPTS OF THE ARRAY A, RESPECTIVELY.

(INPUT)

- IHPUT OPTION PARAMETER.

IF IJOB IS POSITIVE, THE PAST FOURIER

TRANSFORM OF A IS TO BE CALCULATED.

IF IJOB IS NEGATIVE, THE INVERSE

PAST FOURIER TRANSFORM OF A IS TO BE

CALCULATED.

- INTEGER WORK VECTOR OF LENGTH

6*HAX(Nt,N2,N3) +150.

- REAL HORK VECTOR OF LENGTH

6*HAX(N1.N2,N3)+150.

- COMPLEX WORK VECTOR OF LENGTH

HAX(H2,N3).

- SINGLE AND DOUBLE/H32

- SINGLE/H32,HH8,H60

- PFTCC

- INFORMATION ON SPECIAL NOTATION AND

CONVENTIONS IS AVAILABLE IN THE MANUAL

INTRODUCTION OR THROUGH IHSL ROUTINE UHELP

REMARKS 1. IF IJOB IS POSITIVE, FFT3D" CALCULATES THE FOURIER

TRANSFORM, X, ACCORDING 70 THE FOLLOWING FORMULA

XJT+1. J+1,K»1)=TRIPLE SUM OF A (L+ 1 . H + 1 .N + 1) *

EXP(2*PI*SQRT(-1)*(I*L/N1+J*M/N2+K*H/N5 ri^fiw* * rt **«* * * •* */\ ++*\ * «• ^ /* »<)» ) * ' ' WITH L=O...N1-1, M=O...N2-1, N=O...N3-1

AND PI=3.1415...

IF IJOB IS NEGATIVE. PFT3D CALCULATES THE INVERSE
163

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

C

POHRIER TRANSFORM, X, ACCORDING TO THE FOLLOWING

FORMULA

X(I+1,J*1,K»1)=1/(N1*N2*H3)*TRIPLE SOH OF

<

COPYRIGHT

WARRANTY

'

EXP>-5*PI*SQRT(-1}*jri*L/N1+J*H/N2»K*N/N3))

HITS L=O...N1-f, B=<)...N2-1, N=O...N3-1

AND PI=3.1»15...

HOTE THAT X OVERWRITES A OH OUTPUT.

IF A IS A TWO DIMENSIONAL ARRAY, SET N3 = 1.

IF A TS A ONE DIMENSIONAL ARRAY (VECTOR) ,

SET IA2 .* N2 = N3 = 1.

- 1980 BY IHSL, INC. ALL RIGHTS RESERVED.

- IHSL WARRANTS ONLY THAT IHSL TESTING HAS BEBI

APPLIED TO THIS CODE. NO OTHER WARRANTY,

EXPRESSED OR IMPLIED, IS APPLICABLE.

C

C

C

SUBROUTINE FFT3D

INTEGER

DOUBLE PRECISION

COHPLEX*16

INTEGER

DOUBLE PRECISION

COHPLEX*16 C123

IP (IJOB.GT.O) GO TO 10

DO 5 1=1,111

DO 5 J=1,N2

DO 5 K=1,N3

(A, IA1,IA2,N1 ,N2, H3.IJOB, IWK.RWK.CWK)

SPECIFICATIONS FOR ARGUMENTS

IA1,IA2,N1,N2,N3,IJOB,IWK (1)

RWKJ1)

ACIi1,IJk2,N3J,CWK(1)

SPECIFICATIONS FOR LOCAL VARIABLES

I.J.K,L,M,N

HI 2

FIRST EXECUTABLE. STATEHENT

INVERSE TRANSFORH

THANSFORS THIRD SUBSCRIPT

DO 25 L=1,N1

DO 25 H=1,N2

DO 15 ft=1,N3

CWK(N) » A(L,H,N)

CONTINUE

CALL FPTCC (CWK,N3,IWK,Si?K)

DO 20 K=1,N3

TRANSFORM SECOND SUBSCRIPT

10

15

25 CONTINUE"

DO HO L=1,N1

DO 40 K=1,N3

DO 30 H=1,N2

CWKlfH) = A(L,H,K)

30 CONTINUE

CALL FPTCC fO?K,N2,T7K,RSK)

DO 35 J=1.N2

AJL.J,K) = CWK(J)

35 CON7IMUE

ttO CONTINUE

TRANSFORM FIRST SUBSCRIPT

DO t»5 J=1,N2

DO ft5 K=1,N3

CALL FPTCC (A(1,J,K) ,N1,TUK,RWK)

15 CONTINUE

TP 7IJOB.GT.O) f?0 TO 55

ni2.1 = N1*N2*M3
164

C123 = DC-1PLX(R123,O.ODO)

DO 50 1=1,N1

DO 50 J=1,N2

DO 50 K=1.m

CONT-,

RETURN

END

3V/ rv — I i w ,1

A(I,J,K) = DCONJG(A{T,J,K))/C123

TIHOE