I NBS

PUBLICATIONS

NBSIR 81-1643 AlllOt. 171223

EFFICIENT COMPUTATION OF THE FAR FIELD

RADIATED BY AN ARBITRARY RECTANGULARAPERTURE DISTRIBUTION

COMPUTER PROGRAM DOCUMENTATION

Richard Lewis

Electromagnetic Fields Division

National Engineering Laboratory

National Bureau of.Standards

Boulder, Colorado 80303

. U56

81-1643

1981

c.2

March 1981

NBSIR 81-1643

EFFICIENT COMPUTATION OF THE FAR FIELD

RADIATED BY AN ARBITRARY RECTANGULARAPERTURE DISTRIBUTION

COMPUTER PROGRAM DOCUMENTATION

RATIONAL flUKSAU

O?1 standards

LIBRA RT JtiN i Q 195|

no-i aCC -G/C

Gluoo

.056)

C.9-

Richard Lewis

Electromagnetic Fields Division

National Engineering Laboratory

National Bureau of Standards

Boulder, Colorado 80303

March 1981

Q V

v 'W0* *

U.S. DEPARTMENT OF COMMERCE, Malcolm Baldrige, Secretary

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director

TABLE OF CONTENTS

Page

List of Figures..... iii

I. INTRODUCTION 1 II. PROBLEM DESCRIPTION .1

III. USER INPUT TO THE PROGRAM 3

IV. PROGRAM OUTPUT ..6

V. PLOTTING SUBROUTINES 8

VI. DIMENSION STATEMENT SPECIFICATIONS 9

VII. PROGRAM VERIFICATION 9

APPENDIX A. Computation of the Far Field due to a Tangential Electric-Field Distribution over an Aperture Plane using Fast Fourier Transform

Processing; Theoretical Analysis........ 11

Appendix A.l Far-Zone Electric-Field Components

Expressed in Terms of Plane-Wave

Spectrum Components 11

Appendix A. 2 Finite Fourier Transform Representation

of the Plane-Wave Spectrum Integral 14

REFERENCES 18

APPENDIX B..... Computer Program Listing for Computing Azimuth

and Elevation Vector Components Along User-Specified

X- and Y-Plane Cuts in the Far Field Arising from a Prescribed Rectangular-Aperture Field Distribution 26

APPENDIX C Computer Program Modifications to obtain Theta and

Phi Vector-Component Output Instead of Azimuth and

Elevation Vector-Component Output 36

APPENDIX D Computer Program for Computing

Principal-Plane-Cut Far Fields 37

LIST OF FIGURES

Figure 1. Antenna coordinate system using A and E spherical angles

with y as the polar axis 19

Figure 2. Antenna coordinate system using a and e spherical angles

with x as the polar axis 20

- i i i -
Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Problem geometry, showing the rectangular aperture within an aperture plane and the far-zone position coordinates. ......... .21

The principal plane cut <p = PHI -component amplitude

plot for a 5-wavelength-wide slit aperture

(varied data-point-spacing increments) 22

The principal plane cut 4> =

-J-. PHI -component amplitude

plot for a 5-wavelength-wide slit aperture

(varied Fourier integration ranges).. 23

Antenna pattern magnitude for a uniformly excited circular

aperture. Y-plane cut with constant elevation angle E = 10°...24

Antenna pattern magnitude for a uniformly excited circular

aperture. Y-plane cut with constant elevation angle E = 30°...25

-iv-
Efficient Computation of the Far Field

Radiated by an Arbitrary Rectangul ar-Aperture Distribution

Computer Program Documentation

Richard Lewis

This report contains the computer documentation for calculating

the far-zone electric field due to a user-prescribed electricfield distribution within a rectangular aperture. The far-field

output is computed along two arbitrarily selected, perpendicular,

spatial -frequency plane cuts. Program execution time is minimized

by the use of fast Fourier transform (FFT) processing. The

program was designed so that the required far-field output is

obtained by processing only two, vector, one-dimensional FFTs.

The far-field results are obtained in the form of elevation and

azimuth vector components and electric-field-vector magnitude. A

complete analytical discussion of the problem is presented, along

with sample graphical output to illustrate how aliasing and output

resolution limitations effect the graphical results.

Key words: Algorithm; antenna; aperture; computer program;

electromagnetic; far-field; FFT; Fourier transform.

I. INTRODUCTION

This report constitutes the computer-program documentation for

calculating the far-zone electric field, due to a prescribed electric-field

distribution within a rectangular aperture, along two perpendicular spatial

frequency plane cuts. Program execution time is minimized by use of fast

Fourier transform processing. The analytical basis for the computations is

presented in Appendix A. The following text describes the computer program,

listed in Appendix B, that has been developed for implementing these

computations. The same nomenclature is used by the text as in the program

listing for the variable parameters; consequently, reader familiarity with the

input parameter list, as specified by the program comment cards, should prove

beneficial

.

II. PROBLEM DESCRIPTION

We divide space into two regions by an infinite plane surface (the

aperture plane) containing the rectangul ar-aperture surface. The region in

which the far-field computation is desired is assumed to be homogeneous and
source free. The theoretical basis for calculating the far field due to a tangential electric-field distribution over an aperture plane is well known

[1,2,3]. For physical antenna apertures that are large compared with a wavelength, the electromagnetic field amplitude will be small everywhere on the aperture plane except on the aperture surface. As a first approximation

the far-zone radiation field may be found from the field on the aperture

surface alone [1]. The aperture field may be determined either through direct

near-field measurements or by suitable approximate analytical methods, such as ray tracing. It is presumed that the computer-program user can specify the

tangential, aperture, electric-field vector (both amplitude and phase) at

equidistant points on a rectangular grid over the aperture surface, while the

field elsewhere on the aperture plane is assumed to be zero. The far-field

output is obtained as a function of the spatial -frequency components k

x and

k . Here, we compute the far-field output along user-chosen k

x

- and k

y

-plane

cuts using fast Fourier transform (FFT) processing. We will refer to the k x =Q

and k =0 plane cuts as the principal -plane cuts. With just one k -plane cut

and one ky-plane cut far-field computation, the computer program herein

described was designed so that the input data set would only be read in once while the far-field output would be obtained by just carrying out two vector

one-dimensional FFTs. The result is an extremely efficient program, in terms

of computation time and computer storage, for computing the far fields arising

from a rectangul ar-grid aperture distribution along two perpendicular plane

cuts

The x-dimension of the rectangular aperture over which the user specifies

the tangential electric-field vector is a; the y-dimension is b. The aperture

plane is given by z=0, while the assumed aperture radiation travels in the

positive z-direction. The number of data points within the aperture along the

x-direction is NX, while the number of data points along the y-direction is

NY. Let us presume that a matrix of x-component data and a matrix of

y-component data is available for specifying the aperture field, and that each

component value at a given data point is specified as a complex number,

| E | e

1 ^

= | E | cosip + i

j E { s i nip

j E j being the amplitude and ij; being the phase of a particular datum. The

dimensions of these two vector-component matrices would each be given by

2
NY x NX, so that each matrix row corresponds to a constant value of the

y-coordinate, while each matrix column corresponds to a constant value of the

x-coordi nate. Thus, the row index r of the data matrix is proportional to the

y-coordinate, while the column index c of the data matrix is proportional to

the x-coordi nate. For each row or column, a data point at the center of the

respective x-dimension or y-dimension of the aperture shall be specified,

along with an equal number of data points on each side of the midpoint.

Consequently, NX and NY are both odd numbers. The spacing between adjacent

data points in the x-direction is 6 = a/NX, while the spacing in the

X y-di recti on is 6^

= b/NY. It is presumed, by these equalities, that the edges

of the rectangular aperture are located half a data-point spacing beyond the

farthest data points within the aperture.

III. USER INPUT TO THE PROGRAM

The program user provides a subroutine, GETARAY, that obtains or calculates one row of x-coordinate aperture-field data, correspondi ng to a particular value of the y-coordinate and to a particular vector component,

each time it is called. The subroutine arguments are NX and STOR, where NX is

the number of data points in the row and STOR is a complex array of dimension

NX for storing the aperture-field data. The value of NX may not be changed by

the user's subroutine, GETARAY. This subroutine is called NY times for the

x-component data, followed by another NY times for the y-component data. Each

successive call represents a lower row (greater row-index number) in the data

matrix of each vector component, or a larger value of the y-coordinate.

Within the row of data supplied by GETARAY, each successive datum represents a larger value of the x-coordinate. In those regions of the rectangular

aperture that are beyond the boundary of the user's physical aperture, the user is expected to provide zero-fill via subroutine GETARAY.

It may be observed that spatial integration of a physical electric field

over the aperture plane will give rise to a function of spatial frequency k x and ky that is virtually bandlimited [2]. According to the sampling theorem

[4], a band-limited two-dimensional function with bandlimits K

x and Ky may be

completely recovered using sampling intervals 6 and 6 that satisfy the

x y

relations.

3
K = f- and K - f- . x 6 X y 6 y

A virtual bandlimit only slightly greater than the spatial propagation

constant k = 2tt/a ordinarily prevails, where A is the wavelength in the

medium, so that an aperture data-point spacing of 6 - 6 < A/2 generally is

x y

~ adequate. However, in order to minimize aliasing errors in the far-si delobe

region, closer data-point spacings may be required (cf. Appendix A, eq. (16)

and the preceding discussion, for analytical details.) Consequently,

computer-program flexibility is maintained with the specification

6 x

£ and 6y

= A-.

Here, S * SX and S = SY, where SX and SY are dimensionless numbers

x y

specified by the program user. There is only a small difference between S

x and SX and between S and SY, which the computer program can automatical ly

compute. Thus, the user need only be concerned with the order of magnitude of

SX and SY. With a uniformly-excited-aperture test case, highly accurate

comparisons with theoretical results required values of SX, SY > 9 due to the

presence of evanescent modes. If evanescent modes are not present, values of

SX, SY = 2 should suffice for most applications. Values of SX, SY < 2 should

be used cautiously. In addition to SX and SY, the user must specify the

variables ALPHA = a and BETA = 3, where a and 3 serve to specify the aperture

dimensions normalized to the wavelength. Thus,

a . b a e j

and 3 = “•

As a consequence of these definitions, the values of the aperture-spacing

normal ization constants, Sv and S.., are given by

a y

S

x NX

and c NY

y

' 3*

In practice, the computer program makes the initial computations, NX = a*SX

and NY = 3*SY, and then checks to see if the values of NX and NY so computed

are odd. If not, then unity is added to the result so that the final values

of NX and NY will be odd. Then, the exact values of S

x and Sv are computed.

Although only the field within the aperture region will contribute to the

4
integrand, with FFT processing the resolution of the far-field output will

depend on the size of the area of integration. As an example, if a uniform

aperture excitation is specified, and if the area of integration just equals

the aperture area, then the program would produce far-field output at the

points k = 0, + 2ir/a, + 4ir/a, + bn/a, etc., corresponding to the exact

X (elevation-plane) far-field pattern function.

E - sin (k

In other words, except for one point, one would just obtain output values at

the nulls of the far-field pattern function, which would provide extremely

limited information. A more detailed picture would be obtained by increasing

the output resolution so as to produce output values at intermediate points

between the pattern-function nulls. By extending the area of integration and

adding zero-fill to the FFT data input, one can increase the output

resolution.

The extent of the area of integration is increased by increasing the

number of data points supplied to the FFT processing routine. The data

amplitude values at each additional data point are automatically set equal to zero by the computer program. The number of data points supplied to the FFT

processor is designated NNX in the case of x-direction FFT processing and NNY

in the case of y-direction FFT processing, where NNX > NX and NNY > NY. In the case of the FFT processor used by the computer program, the number of

input data points must equal a power of two. Thus, NNX = 2**MMX and

NNY = 2**MMY. The program automatical ly chooses NNX to be twice the largest

power of two that is just greater than NX, and similarly for NNY and NY.

However, these automatically generated values may not suffice for all

applications; consequently, the automatically generated values of NNX and NNY are multiplied by 2**NZX and 2**NZY, respectively, to produce updated values

of NNX and NNY. The two constants NZX and NZY are specified by the user,

where NZX, NZY >_ -1. The recommended specification is NZX=NZY=1 or greater,

although specifying NZX or NZY equal to zero may be adequate if the aperture

field tapers to near zero at the aperture boundaries. With a uniformlyexcited-aperture test case, the values NZX=NZY=2 gave smooth plotting

resolution with a microfilm plotter subroutine. It is not recommended that

5
the -1 specification for NZX or NZY be used unless NX or NY is just slightly

larger than some power of 2, as otherwise there will not be adequate zero-fill

added to the input data. It should be recognized that the larger that NNX and

NNY become, the longer the FFT computation time becomes and the larger the

computer-storage requirement becomes.

In addition to the input variables already described, the user must

specify the azimuth angle A

0 and the elevation angle E

c corresponding to the

particular plane cuts desired. These input quantities are best discussed

within the context of the next section.

IV. PROGRAM OUTPUT

The program computes the azimuth (E^) and elevation (E^) vectorcomponents of the far field, along with the antenna pattern or electric-vector

magnitude (the program listing for azimuth and elevation vector-component

output is as given in appendix B, while a listing of substitute cards for use when theta and phi vector-component output is preferred may be found in

appendix C). The program computes both output vector components along a constant ^-coordinate plane cut and along a constant ky-coordinate plane cut

as respective functions of elevation (e) and azimuth (A). That is, the

spatial -frequency component that is not held constant becomes the sole outputfunction variable. Note that both the elevation (e) and azimuth (A) angles

become equal to the polar angle theta in the event that principal -plane cuts

are selected (users who only require principal -pi ane-cut computations may wish

to utilize the simpler version of the program given in appendix D). The

relationships between the various spherical -coordinate-system angles that have

been mentioned are shown in figures 1 and 2, with supplementary discussion

found in [5]^.

At this point it would be appropriate to note that the final computed

values of the far-zone electric fields are renormal ized, such that the factor

'The notation used here to designate the alternative azimuth and elevation

angles is "a" and "e", whereas in [5] these alternative angles are designated

a and p. The notation A and E normally is used with az/el antenna mounts,

while the a and e designations apply to el/az antenna mounts. Our A and E angle designations agree with the notation used in [5] for these quantities.

/

6
— e^r is di vided out.

r Corresponding to each start-up of the program, a complete set of input

x-component and y-component aperture data is requested from GETARAY by the

program, and one k -plane cut and one k -plane cut far-field computation is

carried out. The particular k x -plane cut generated is selected by specifying

the azimuthal angle A

Q

, at an elevation angle of zero degrees, that

corresponds to the particular value of k x desired. Similarly, the particular

ky-plane cut generated is selected by specifying the elevation angle EQ

corresponding to the particular value of ky desired (the "o" subscripts are used here to distinguish the user-specified azimuth and elevation angles that

define the plane cuts, whereas the omission of this subscript implies the

general angular-coordinate variables). The relationships between the farfield spatial -frequency variables, k x and k^, and the spherical -coordinate

angles, depicted in figures 1 and 2, are given by,

k = k Cos E Sin A = k Sin 0 Cos <j> = k Sin a X

ky = k Sin E = k Sin 6 Sin <j> = k Cos a Sin e

From the preceding, we see that the desired k - and k^-pl ane cuts are

X Jr given by k /k = Sin A and k /k = Sin E , respectively, where the angles A

X 0 J/ 0

and EQ are specified by the user. Upon completing the far-field computations, the azimuth and elevation

vector-components and the antenna pattern are plotted as a function of

elevation (e) and azimuth (A), respectively, for the k

y -pi ane cut and the

ky-plane cut. The plotting is limited to angles in the visible range (i.e.,

real angles within the spherical spatial -frequency region of radius k, such

that k = /k2 - k

2 - k2 is real). This visible range not only excludes

z x y

complex evanescent angles, but also limits the range on k^ and k x to

k^ <_ k Cos A

q

and k x

_< k Cos EQ

, respectively, in the case of the k x

- and

ky-plane cuts. Since the k x -plane cut is plotted versus elevation (e), while

the ky-plane cut is plotted versus azimuth (A), one obtains a full ± 90°

plotting range. However, it must be recognized that rapid variations, as respective functions of k or k x

, will appear greatly expanded when A

Q or E 0

do not equal zero.

7
V. PLOTTING SUBROUTINES

The subroutine PAMPIOT that is supplied with the program for setting up

the plotting arrays permits the user to restrict the plotting range through

the specification of maximum and minimum abscissa values. Both relative

amplitude in dB and phase in degrees are plotted for each vector-component,

with just relative amplitude in dB plotted for the antenna pattern. This

routine labels each plot with the peak amplitude level in dBL^ and with a self-generated caption that automatically indicates the vector-component

plotted; phase, amplitude, or magnitude plot; and the value of either the

azimuthal angle (A

Q

) or the elevation angle (E

Q

) corresponding to the

particular k x

- or k^-plane cut that has been generated. Also, a user- specified 40-character graph legend is placed on each plot. In addition to

plotting far-field output, subroutine PAMPLOT allows the user to obtain a printout of every jth data point, j=l,2,..., etc. (Note that the reciprocal

of j times 100 gives the percentage of output points to be printed.) If j=0,

printout is suppressed. This printout gives the abscissa point, absolute

amplitude, and phase in degrees.

The actual plotting routine supplied with the program is limited to point

plotting on the user's printout. The user can elect either to point-plot the

amplitude and phase of each vector component separately, or else the amplitude

and phase curves can both be point-plotted on the same graph. The latter

plotting mode has the advantage of placing two closely related curves in

proximity to each other. The disadvantage is that ordinate-axis labeling of

the phase plot will be suppressed by ordinate-axis labeling of the amplitude

plot. A microfilm plotting routine is highly desirable. Subroutine PAMPLOT

contains an illustrative calling statement to the dummy routine CRTPLOT, which

may be replaced by the user's own pi otting-routine call. The dummy subroutine

that is called here is SUBROUTINE CRTPLOT (XARRAY, YARRAY, XMAX, XMIN, YMAX,

YMIN, N, LABEL, NO, NOMAX), where XARRAY contains the N abscissa points,

YARRAY contains the N ordinate points, XMAX is the greatest abscissa value

plotted, XMIN is the least abscissa value, YMAX is the greatest ordinate

‘Here, dBL is taken to mean 20 Login(E), where E designates the quantity whose

amplitude level is desired, (cf. IEEE Standard Dictionary of Electrical & Electronics Terms under "level".)

8
value, YMIN is the least ordinate value, N equals the number of data points

per curve, LABEL is a nine-word graph legend, NO is the number of the curve

being plotted, and NOMAX is the number of curves per graph (1 j< NO _< NOMAX <2 )

.

VI. DIMENSION STATEMENT SPECIFICATIONS

The program comment cards supplied with each routine should enable the user to specify the correct input to run the program. Particular attention is

called to the array-dimension specifications in the main program.

Insufficient dimensions (as specified by the user's data statement) will cause

the program to abort. Note that the dimensions of arrays THETA and DATA are based on the final computed value of NNX, while the dimension of array El

equals twice the sum of the final computed values of NNX and NNY. The

dimensions of arrays EXPX and STOR are equal to NX, while the dimension of

array EXPY is equal to NY.

VII. PROGRAM VERIFICATION

In order to verify that the computer program was working properly, it was tested to see how accurately it could compute the far fields of uniformly

excited apertures, for which exact formulas are available [7, pp. 332-336].

Good agreement was obtained between the theoretical results and the computer

program output, particularly as the data-sampl ing interval was decreased. As

mentioned earlier, very short data-point spacings were required due to the

presence of evanescent modes in these hypothetical aperture distributions.

Tests were carried out both for a rectangular aperture and for a circular

aperture with zero fill to the boundary of the enclosing square aperture. In

order to facilitate numerical comparisons, computations with the exact

formulas were carried out at the same far-field coordinates that the computer

program used to calculate the far-zone fields.

In figure 4 we show some pri ncipal -pi ane-cut results for a uniformlyilluminated slit aperture five wavelengths wide, in which we compare results

from computations with a data-point spacing of A/5 and x/3 to the

correspond!' ng exact theoretical curve. The solid-line curve in figure 4

corresponds to the A/5 data-point spacing curve, while the dots and crosses,

respectively, correspond to the theoretical curve and to the A/3 data-point

9
spacing curve. Just those portions of the latter two curves which do not

overlap the solid line curve are shown in the figure. Thus, the only

differences shown are in the far side-lobe region where the overlapping ends

of the replicated far-field function contribute the most to the aliasing error

(refer to Appendix A, eq. (15) for analytical details). As may be anticipated

from this figure, we have verified that the computed curve does get closer and

closer to the theoretical curve as the data-point spacing gets smaller and

smaller.

To illustrate the effect of adding zero fill to the FFT to increase the

output resolution, in figure 5 we show a comparison between two curves having

the same data-point spacing but different specified values of NZY. Here, the

solid line curve corresponds to NZY = 0 while the dotted line curve corresponds to NZY = 3. Adjacent points on each curve are connected by a straight line with no curve fitting between points. Both figures 4 and 5 correspond to the same plane cut and the same aperture excitation, but the

solid line curve of figure 4 was generated with NZY = 4; consequently, the

nulls shown in figure 4 are uniformly deeper than the nulls shown in the

dotted line curve of figure 5. However, the computed values of figure 5 are the more accurate, as they were obtained using a data point spacing of A/7.

Finally, in figures 6 and 7 we show some antenna pattern magnitudes for a uniformly excited circular aperture three wavelengths in diameter. Here, we were interested in comparing the computer program output against the

theoretical model for some arbitrarly selected plane cuts. Consequently, we show two ky-plane-cut antenna-pattern curves plotted versus azimuth, in

figures 6 and 7, corresponding respectively to a 10° elevation angle and a 30°

elevation angle. The peak amplitude in figure 7, incidently, is about 15 dB

below the peak amplitude in figure 6. In order to achieve highly accurate

program output, a data point spacing of was used, which resulted in the

generation of numerical dB values for the furthest side lobes that agreed with

the theoretical results within eight-tenths of one percent. This compares to

just five or ten percent agreement between the furthest side lobes of the

theoretical model and the far side lobes of the computed curves plotted in

figure 4. A value of NZX = NZY = 2 was selected for our circular aperture

computations, resulting in just moderately smooth plotting resolution. It might be remarked that the aperture excitation was assumed to be diagonal to

the cartesian coordinates of the aperture, resulting in the asymetrica! k,,-

plane-cut patterns shown.

10
APPENDIX A

Computation of the Far Field due to a Tangential Electric-Field

Distribution over an Aperture Plane using Fast Fourier Transform Processing;

Theoretical Analysis

A.l Far-Zone Electric-Field Components Expressed in Terms of Plane-Wave

Spectrum Components

An expression for the far-zone electric-field vector transverse to the

z-direction is readily obtained by integrating the near-zone tangential field over an aperture plane [1], [2]. Thus, we have the far-zone electric field

transverse to the z-direction expressed as

£t (£) ~

~ lk

fr

° S6

£> ( 1 )

where the aperture-plane integral is defined as

k • » -A R-P'

B(pR)=/ / £t (P') e

r dx'dy'.

—CO —OO

( 2 )

Here, P' = x'a , + y'a. and R = xa + ya , where a Y = a Y. and a u — = a yt are —x —y — —x —y —* —* —y —y

unit vectors along the x and y coordinate-system axes, and where x and y

denote far-zone position coordinates while x' and y' denote aperture-pl ane position coordinates. The far-zone radial coordinate is r = /x2 +y2 +z2 . These geometrical relationships are depicted in figure 3. The propagation

constant is k = 2tt/A, where X is the wavelength. An e -iwt time dependence is

implicitly assumed by (1); thus k = — where c is the velocity of light. The

spatial -frequency components k Y and k v are related to the far-zone position

coordinates by the relations, k = k— and k = k—. E. (P 1 ) is the tangential

x r y r l electric-field vector in the aperture plane. The z-component of the far-zone

electric field is also required, but it can be obtained from the expression

B (— R) =

7X— R-B(— R), (3)

z V —y r Cose v r —

which results from the fact that each plane-wave spectrum vector,

k k B(— R) + a B (— R), must be orthogonal to the corresponding propagation

11
direction k_ = ^

(_R + za_j

Consequently, it readily follows from (1) and (3) that the theta and phi

components of the far field can be expressed as [1]

E eW Ui

2tt

ikr

[BL (-

x'r

R) Cos<j> + B (- R) Si n<j>] (4a)

Ik e

2tt

ikr

r [-B

x (i R) Sin* + B

y

(iR) Cos*] Cose (4b)

In place of the theta-phi component description of the far field, the

alternative elevation and azimuth component description is frequently chosen

(cf. figure 1). These far-field components are given by

E E (r)

-ik

2n

B

y

(^R) CosA (5a)

ikr

E A(r) ~ ^ [B

x (^ R) CosE + By

(| R) SinE SinA] (5b)

The radial component of the far field, which is perpendicular to each of the

components in equations (4) and (5), is equal to zero.

The elevation (E) and azimuth (A) angles, shown in figure 1, are obtained

from

SinE = Sine Si n<j>, TanA = Tane Cos*}), (6)

whereas the usual spherical-coordinate angles theta (e) and phi (<f>) may be

expressed in terms of the far-field spatial -frequency components k y and k u a y

usi ng

k /k2 + T2 Tan4

> = Sine = —~

--

. (7)

K K x

The alternative azimuth (a) and elevation (e) angles, depicted in figure 2, are obtained from the expressions

Sin a = Sine Cos<f> = CosE SinA (8a)

12
and

Tan e = Tane Sin<j> = TanE SecA. ( 8 b)

We can express the magnitude of the far-zone electric-field vector,

|EJ 2 = E^E*, where * implies the complex conjugate, as

lEl 2 = iE

e

l 2 + |E |2 = |e e

|2 + |E a

|2.

Substituting from either (4) or (5) into the above expression results in

I

E

| 2 ~ -7^—5

- { |B | 2 Cos 2 E + |B | 2 Cos2 a + (B B* + B*B ) Sin E Sin a}, (9) — 4it 2 r

2 x y x y x y

k k k

where the argument, — R, of B (— R_) and B (— R) is understood.

1 A I Jr I Finally, we give some explicit relations between the vector-components

which hold on the principal -pi ane cuts. Thus, for the principal plane <j> = 0 we have:

E = e = 0, A = a = e —R = r sine a—x

and

E A(r) = E,(r)

Ik P

fkr

k S’ f- B x<7

E F (r) = E (r)

t — <j> —

ik e 1krk sVv^1 Cose - ( 10 )

On the principal plane <j\> = j

the relationships become:

A = a = 0, E=e=e

JR = r sine a^

and

E

A

(r) - E (r)

<P “

ik e

"2tT “

ikr — B

x (^ R) Cose

13
(11)

.. ikr , E

E (I) - Ee (r)~ VBy4 R>

It may be noted that a jump change In sign of the 0 and <j> components, due to

the jump change In 4

> by ir radians upon passing through the coordinate-system

origin, has been suppressed in writing expressions (10) and (11).

A. 2 Finite Fourier Transform Representation of the Plane-Wave Spectrum

Integral

We next require computational expressions for evaluating the aperturek plane integral defined by (2). If we presume that JB(— R) is a bandlimited

function of the spatial -frequency components k Y and k u , having bandlimits

a y

-K <_ k < K and -K <. xxx k < K , then the aperture-pl ane integral (2) can be

y y y

expressed as [2]

BCp R) I I £t <tJ 6 6 e XV L J m=-°° n=-<»

•4

(12)

tore, Pmn

- n*5 x

a x

+ n6

y

ay

, where « x- and S y

- f- • x y

Furthermore, if the aperture-pl ane tangential electric-field vector is

restricted to a rectangular aperture of dimensions a = NX6 and b = NY6 , x y

then the aperture electric field will be non-trivial only over the

1-MY NX-1 1 NY NY-1

range —

y

_< m _< - y , y _< n <_ 311—± , so that the doubly infi nite series

in (12) may be replaced by two finite sums. Moreover, we can sample the

k

transverse plane-wave spectrum vector J3(— R) at the equally spaced

points k x

= J 7

and ky

= £ £

A

y

, where ~ a

x and ~ A

y are the spatial -

frequency sampling intervals, and where j and z are integers. The total

number of x-coordinate spatial -frequency sampling points is N

x

, and the total

number of y-coordinate spatial -frequency sampling points is Ny

, where N

x

> NX

and Ny

> NY. The range on the integers j and ji is

N N N N

- < j < ! and --£ <*<-£- 1 , (13)

so that the bandlimits are

14
K ( 14 ) x 2r

A

x

N

x and K

y

- ^ ay

N

y

The numbers N„ and N„ respectively correspond to the quantities NNX and

NNY described in the computer-program documentation. For the sake of computer

program compatabil ity with existing software, the number of input data points

should equal the number of output data points. This may be accomplished by

zero-filling the input data array correspondi ng to data points outside the

rectangul ar-aperture boundaries.

Although (12) is an exact expression if the aperture-field integration

results in a bandlimited function, our restriction of the input data to a rectangular aperture is not compatable with this condition. Consequently,

there will result an aliasing error with the f i nite-Fourier-transform

relation, so that the actual finite-Fourier-transform relation becomes [6]

-j* l l B(

r

R

. N N ) p=_oo q=_oo ° r x ^ y

6 6 x y

2

1

m= -

h _1

2

1

l N n= -

-H7in e 4VAnn

(15)

where the range on the integers j, i is given by (13), and where we have

written r n for E+ (Pmn ). Also, R„ . = sa a + tA.a , where s,t take on the —him —

l

—»iir i —s,i x—x y—y

k values indicated in (15). From (15), we see that R. ) is replicated as

^ J ^ B o+nw ) for P >(1 = etc.; thus J3. is formed, within the

r ^ ^P^x ^y

range on j and i given in (13), by summing contributions from the primitive

function B (— R. ) to those non-negl igible contributions, from the overlapping

• vJ X/ ends of these replicated functions, that arise from spatial -frequency values

outside the assumed bandlimits K

x and Ky. Here, of course, Kx and Ky simply

correspond to the quantities tt/ 6 and n/6 , rather than to actual bandlimits.

x y

should be noted that an analytical transition from the right-hand side of (12)

to the right-hand side of (15) also involves a replication of _Ej. (_Pmn ) , so that

both the input and output functions in the finite Fourier transform relation

(15) are doubly periodic with periods N

x and Ny However, there is no overlap

of the replicated input functions within the integration interval.

It

15
Although the fini te-Fourier-transform relation (15) does not give a perfect evaluation of the aperture integral (2), it can be recognized that

aperture-pl ane integration of a physical electric field will produce a virtually bandlimited function. This implies that the approximation

B(— R . ) * B . - r -jr -j£ < i < (16)

will be valid provided K

x and are large enough. This in turn implies that

and 6 must be sufficiently small,

x y

J

One difficulty that arises with utilizing (15) is that fast Fourier

transform (FFT) subroutines require non-negative integer values for the

indices j, £, m, and n. Fortunately, this difficulty is readily resolved

using the periodic character of B. and L . Thus we have

B.

6 6 —JA

x y

i E-mn

-I 2tt ( jm/N

x

e

+ £ n/N ) y

N

x

-1

I m=0

V1

l £

n=0 mN

x T’

iWj+Jt) -i

tt ( jm/N +tn/N ) e e

^

n- y"7

Finally, we obtain

V1 = e iiT(j+£)

j x

“y

j+ —^-,£+ —^

m=0

y

I e n=0

-iir(m+£ ) E m-

-i2ir( jm/N +tn/N ) ^ e x y

"T

(17)

This result enables us to integrate over the aperture surface,

- T

X £ X ' i V

“ T

y £ y' £ T

y

> where T

x

= \ Vx and T

y

= 7

Hy6y

» US1

‘ n 9

standard FFT-program format. The output of the FFT routine is a function of j

and £, which range over the values

16
0 < j < Ny - 1, 0 < £ < N - 1. x y

Consequently, since B. is periodic in j and £ with periods Ny and Nv , respectively, by using (17) the FFT output can be plotted directly as though

the range on j and £ were as given in (16).

Finally, we write down expressions that are compatible with the problem

of computing the far field on the plane cuts. Thus, for the plane ky = const,

we have

i \ . £ + -p N -

"c— SinE , y $

y

(18)

where S„ = 7— and E assumes a fixed value. We can now collapse the summation

V 0 y

over y (i.e., reformulate the problem so as to just carry out an FFT on data

formed by summing the input-matrix columns) to obtain

1 B

6 6 - N N xy

j + 4*+ -^

. . v1

. V1

= e

l7rJ' y e“

l7rm

y m=0 n=0

(Em- ^n3

2

* 2

-i (^ Sin E)(n-

_"j

2

tt jm/N

e y }e

(19)

Similarly, for the plane k y = const, we have

1 N

x . J +2 N

x = “Sin a , ( 20 )

where S = — and a assumes a fixed value. We now collapse the summation over x

x 6 x (i.e., reformulate the problem so as to just carry out an FFT on data formed

by summing the input-matrix rows) to obtain

1 B = e V1

I

e- n=0

V1

l m=0

. /2tt N

" N

x m- — ,n-n —

-1

( $

Sin a ) (m2

) _i 2-n^n/N

e

x

} e y . (21)

17
REFERENCES

[1] Collin, R. E., Zucker, F. J., Antenna Theory Part I, McGraw Hill (1969).

[2] (a) Kerns, D. M. Correction of near-field antenna measurements made with

an arbitrary but known measuring antenna. Electronics Letters 6, pp.

346-347 (1970).

(b) Kerns, D. M. Plane wave scattering matrix theory of antennas and

antenna-antenna interactions, NBSIR 78-890 (June 1978) (forthcoming

as NBS Monograph 162)

[3] Yaghjian, A. D. Upper-bound errors in far-field parameters determined

from planar near-field measurement. Part I: Analysis, NBS Tech Note 667

(Oct. 1975).

[4] Landau, H. J. Sampling, data transmission, and the nyquist rate, Proc.

IEEE 55, #10, 1701-1706 (Oct. 1967).

[5] Newell, A. C. Planar near-field measurement techniques, (lecture notes

for the NBS short course. Antenna Parameter Measurement by Near-Field

Techniques, 7-11 July 1975).

[6] Cooley, J. W. , Lewis, P. A. W. , Welch, P. D. The fast Fourier transform

algorithm and its applications, IBM Research Paper RC-1743, (Feb. 9,

1967).

[7] Johnson, C.C., Field and Wave Electrodynamics, McGraw Hill, 1965.

18
Y

Figure 1. Antenna coordinate system using A and E spherical angles

with y as the polar axis.

19 .
Y

Figure 2. Antenna coordinate system using a and e spherical angles with x as the polar axis. H

20 .
ro|

fu

Figure 3. Problem geometry, showing the rectangular aperture within an aperture plane and the far-zone position coordinates.

21 .
Relative

Amplitude

in

dB

0.0

" 6.0

12.0

18.0

- 24.0

- 50.0

- 36.0

* 42.0

- 48.0

-54.0

1 LEGEND

-

® ® • XXX

LMLl

A/3 spacing

A/5 spacing

A

A IT !

' 1

vyr

/

1

*

1

\

1

/

1\\

j

1

1

- 60.0 Li - 90.0 - 72.0 - 54 .0 -36.0 -18.0 0,0

Theta in degrees

.0 36o0 54,0 72.0 SC

Figure 4. The principal plane cut <j> = PHI -component amplitude

plot for a 5-wavelength-wide slit aperture

(varied data-point-spacing increments).

22 .
Relative

Amplitude

in

dB

0.0

-6.0

12.0

-18.0

-24.0

-30.0

•36.0

-42.0

-48.0

54.0

-60.0

-90.0 - 72.0 -54.0 -36.0 -18.0 0.0 18.0 36.0 54.0 72,0 90.0

Theta in degrees

Figure 5. The principal plane cut <\> = J. PHI-component amplitude

plot for a 5-wavelength-wide slit aperture

(varied Fourier integration ranges).

23 .
Relative

Amplitude

in

dB

Figure 6. Antenna pattern magnitude for a uniformly excited circular

Q

aperture. Y-plane cut with constant elevation angle E = 10 .

24 .
0.0

- 6.0

- 12.0

18.0

- 24.0

- 30.0

- 36.0

-

\l

.

-42.0

- 48.0

-64.0

-60.0

- 90.0 - 72.0 -54.0 -36.0 -18.0 0.0 18.0 36.0 54.0 72.0 90.0

Azimuth in degrees

Figure 7. Antenna pattern magnitude for a uniformly excited circular

aperture. Y-plane cut with constant elevation angle E = 30°.

25 .
APPENDIX B

Computer Program Listing for Computing Azimuth and Elevation Vector

Components Along User-Specified X- and Y-Plane Cuts in the Far Field

Arising from a Prescribed Rectangular-Aperture Field Distribution

PROGRAM F AREA ( INPUT ?0UTPUT,T4 p £i=TNPUT)

DIMENSION S T 0P(1C24» ?EXPX C 1 C 2

)

? EXPY (1024) * El (4396)

DIMENS ION THETA(lC24) t DATA(2343>

COMPLEX STOR,EXPX?EXPY? EX?ZZ

EQUIVALENCE ( T MET A, STOP ) , (EXPX ? DATA)

DATA <NMAX=1024) ? ( NOT M FN = 40 95

)

94 FORMA T(4F7. j, 1X21 3, 2FE. j, 5XA10 ) 95 FORMAT (1H1, 2X52 HA LP HA =X -COORDINATE APERTURE 0 TME NS I ON IN WAV ELENG T

0HS?9X‘- 1HBETA=Y-C00R3INATE APERTURE DIMENSION IN WAVELENGTHS/ 3X57 HS

PX=RATI0 OF WAVELENGTH TO X-COO RDTNAT E DAT A-POTNT SPACING, 4X5ZH3 Y = R

QATIO OF WAVELENGTH TO Y-COORDINATE DATA-pQTNT SPACING/3X49HNZX-X-C

ROCROINATt ZERO-FILL- REGION DOUBLING F ACT OR, i ?X 49HN7Y= Y - C OORO I N A T

t

SZERO-FILL-REGION DO JB LING FAC r OR/3X64HA=AZTMUTH ANGLE (AT ZERO D EG

TREES ELEVATION) OF A7IMUTH-PLANE CUT ? 1 4X 40 HE= EL I V A T ION ANGLE O

r EL

UEVATION-PLANE CUT/3X 3 9H3ANNER=10- CHARACTER PLOT IDENTIFICATION? 4X

V45HNNX, NNY=NUMBER OF TERMS IN X OR IN Y FF T S UM S , 6X 25HNNX=2 *

A

X W, NNY=2**MMY////BH ALPHA = , F8 . 3, 4X6HBET A = , F 8. 3 ? 4 X4H SX =,R7.3?4X

X4HSY =,F7.3*4X5HNZX =,I3,4X5HN7Y =?I3?4X3HA = ? F 5 . 1 , 4X 3H£ =,F5.1»4X

Y7HBANNER. » A10? 1H. //2LX5HNNX = , 15 ? 10 XEHMMX = 13? 30X5 HN NY =,I5?10X

75HMMY = ? I 3///

)

96 FORMAT (1H;., 25X7 tH ARRAY DIMENSIONS TOO SMALL OR ILLEGAL OR ELSE E 0

1R A EQUALS 9 2 DEGREES// IX 5HNNX = ,15? 5X5HMMX = ? 14 ? 10 X5HNNY =,I5?5XS

2HHMY =,I4? 10X23HCOMPUTEO El DIMENSION =,I6,1!X^HE = »F C . 1 ? 5X 3HA = t 3F5.1)

C

C COMPUTATION OF THE FAR FIELD ALONG SPECIFIED PERPENDICULAR

0 AZIMUTH- AND £LEVA r TON -PLANE CUTS? WHERE T H E FAR-FIEL0

Z EXCITATION CORRESPONDS TO A GIVEN NEAR-FTELQ DISTRIBUTION

C IN A SPECIFIFO APERTURE. IN ADDITION , THE FAR-FXELQ

Z ANTENNA PATTERN ALONG THESE SPECIFIED CUTS IS ALSO C0M a lJTC[).

C

c

c

c

cc

-9cc

fc

c

O’c

THE O A

~ A SPACING BETWEEN X- OR Y-COOR DINA TE ELEMENTS IS ASSUMED TO BE FIX

NM A X = D I ME NS ION OF ARRAYS S

T OR AND ^ H E T A 8 DIMENSION 0- ARRAY

DATA SHOULD EQUAL 2*NMAX. NCIMEN IS THE DIMENSION OF ARRAY El.

NDIMEN SHQU„ n EQUAL T WICE THE SUM 0~ "HE X - AND Y»-Ff DIMENSIONS

T H E MAXIMUM USEFUL VALUE OF NDIMEN IS 4 *NMAX.

INPUT P ARA M ET E p S • ® ALPHA = A PEP TUBE WIDTH IN WAVELENGTHS f X-COOROE NATE

)

BE T A = APEPTJRE LENGTH IN WA V EL r NGTHS f Y-COORD T NAT E)

5 X = p A T I O OF WAVELENGTH to X-COOPDINATF QATA-pDINT SP6CT

SY= RATIO OF WAVELENGTH TO Y-CODRDINATE OATA~POXN T SRACI

NZ X = NUMBER OF ExTP A DOUBLINGS OF THE ZERQ-FIL. REGION

FOR TH" X-FAP-FIELO-COOPDINATF FFT. (N0®M4LLY=i»

N Z Y = N UM 3ER, OF EXTRA DOUBLINGS OF t H E ZERO-FIL. REGION

FOR the Y-FAR-FIELD-COOPOINATE FFT. < N 0 ? MAL L Y = 1

)

A= AZIMUTH ANGLE (AT ZERO DEGREES ELEVATION),

IN DEGREES? OF AZIMUTH-PLANE CUT.

E= r LEVA T TON ANGLE? IN DEGREES? CF ELE VATION-PLANE CUT

BANNER = 10-CHARACTER GRAPH IDENTIFICATION

26 .
THE INPUT DA T fi IS ASSUMED to PE SUrn L

T ED AS ROWS O p X-? OCR 01 NATS

DATA, EACH ROW CCRGESPQNDTNG t 0 A FIXEP VALUE OF THE Y - COOROT MA

r . ALL OP THE E-SUB-X PA~A T p TO BE SUPPLTE n FIRST, ANO Tjr N this

INPUT DATA I r TO BE FOLLOWED P Y ALL O r T HF E-SUB-Y DATA.

THE PROGRAM OBTAINS HE IN p UT DATA THROUGH REPEATED CAlLS TO

THE USER-SUPPLIED SUBROUTINE G

E

" A R A Y

EACH CALL T 0 SUBRCLTJNE GFTA PA Y RFSULTS IN A ROW O r NX n OMPL EX

DATA POINTS, T A CH SUCCESSIVE CALL CORRESPONDING TO A ROW

SUCCESSIVELY FURTHER BELOW THE TOP n F

T HE A P ER T UPE. HE'E, NX

IS E DUAL T 0

T HE FFODUC T SY * ALPHA. ~wr NUMBER OF ROWS. 0°

THE NUMBER OF CALLS TO GETARAY, IS GIVEN BY NY = SY * B

r

A

,

NUMBER OF DATA CAPOS USED BY PRO GRAM* SUBROUTINES = 3.

FOLLOWING t HF COMPLETION OF A GIVEN PAR-FIELD COMPUTATION ALONG

SPECIFIED PE D PENC ICULA c A7IMUTH AND ELEVATION PLANES, ~HE PROGRAM

RETURNS TO THE 8EGINING TO READ A M"W DATA CA = 0 AND START a LL

0 V E p AGAIN ON A NEW FAR-FIELC COMPUT ATIDN. THIS time, HOWEVER,

THE PLOTTING ROUTIN'" PARAMETERS AND L

rG E N 0 THA T WERE Rr AD IN, ON

THE FIRST TIME ’HROUGH ''HAT SUBp CU T IN r , A p S SIMPLY RECALLED AND

USED 0 VEF AGAIN. THE USER MUST MAM- SUBROUTINE G

T TAPAY CAN

SUPPLY A NEW c E

T C c N

r A p - F T c

L 0 DATA UPON RESTARTING T H~ p^OG^AM.

PARAMETER p E AD IN

1

READ 94, AL p HA,8FT A,SX,SY,N7X,N7Y, A,£,FANNE r IF ( EOF ( 1 ) ,NE. G ) GO to 36

CALCULATION OF NX AND NY

NX = 0 . F*SX* ALPHA $ NX=2*NX4l 3 NY = C . D *SY* P ETA $ NY=7*NY*-1

SX=NX/ALPHA 3 S Y = N Y / 3 E T A CALCULATION of NNX and NHY, ALONG with MMX AN~ MMY

MMX=MMY=NNX=NNY=i

2 MMX = MMX 4 1 $ NNX=2*NNX 3 IF ( NNX . GT ..NX ) 3,2

3 NNX=2*NNX $ IP(NZX) -4,6,4

44 NNX = NNX/ 2 3 MMX=MMX-1 3 GO TO 6

4 CD F M = 1 , N 7 X * MMX = MMXf!

5 NNX = 2 * NNX

6 M M Y = M ?' Y + 1 3 NNY=2*NNY * I F ( NMY . GT . NY ) T ,6

7 NNY=2*NNY 3 IP(MZY) -5,lC, a 45 NNY = NNY/ 2 $ MMY=MMY-i * GO 1G

8 DO 9 M = i , N 7 Y 3 M M V = M M Y 4i

9 NNY=2*NNY

10 M M = M A X C ( M M X , MM Y ) * N!#M AX 0 ( NNX , MNY 1 * ND = NMX 4 MN Y S 1=?*N0

IF ( J .

T

. NDI MEN. CP . NN , GT . NM AX .0 R . *M . GT . 14 ) 11,12

11 P- I NT 96 , NNX, MMX, NNY. MMY, J,F, A 3 CALL EXIT

INPUT PA p AMPTE p PRINT-OUT

12 P p I NT 95 , ALPHA, BETA, SX , SY ,N2X , N7Y , A , E,BA NN-R, NNX , mmx,NNY ,M m Y COMPUTATION OF A P EB T UR E -0 A T

A PCS TTT ON WITHIN INTEGRATION IN T F-7 VAL

F ACT0R=1. I/ (SY*SY) S IG X= (NNX- NX41) /2 3 IOY = ( NNY -NY 41 ) / 2 «-NN X STNE=SIN(C. "17 4 5329 2519 9- 7* c ) 3 <? IN A =S IN ( 0 . 0 1 5

7 29 251 a 943* A ) COMPUTATION OF APER T UR.F 0 AT A-PCTN T SPACTNGS

OELX = SX/NNX 3 LLX = SQR T ( 1. G-SIM r **2)/D c LX 3 IF(LLX.ED.:> GO TO 11

DEL Y = G Y/ NNY 3 LL Y = SO ~ T ( 1 .

] -<~TM A * * 2 > / n ELY 3 TF (LL Y . EQ . !? ) GO T 0 11

27 .
C RE -ZEROING OF INTEGRAL FIELD

DC 13 M = 1 , J

13 E I ( M ) = ( 0 • , G • )

C COMP IT A TT ON OF CONSTANT EXPONENTIAL MULTIPLICATION r AC T OP - Y-PLAME CUT

ZZ=CEXP (CMPLX (0

, 6 . 23 31 85 3C 71 7 96*SINE/SY )

J=(NY-l)/2 * K = J + l 3 EXPY (K) = (1.C,3.

DO 1C M=1,J 3 L= K-M 3 LL=K+M

EXPY (L) =ZZ*FXPY (L+i J 1C EXPY (LL)-CONJG (EXPY (L) )

C COMPUTATION OF CONSTANT EXPONENTIAL MULT I PL ICAT I ON RACTQF - X- D LANE CUT

ZZ=CEXP( CMPLX (0. ,6. 23 31 8 5 30 71 7 96*SIN A/SX )

J=(NX-l)/2 3 K = J l $ EXPX(tO=(l.C*a.)

DO 15 *=i,J 3 L= K-M 3 LL=Y+M

EX PX (L) =ZZ*EXPX (L +1)

15 EXPX (LL )=COMJG (EX p X (L )

C DO LOOP T Q COMPUTE TWO PLANE-CUT INTEGRALS OV r R APEP^U RE-FIELD VECTOR

DO 21 1 = 1,2 3 MX=NO*(I-l> 3 IX=MX+I0X $ I Y = MX +T H Y

C SET-UP OF INTEGRAND - COLLAPSE ON V- AND Y - P L A N

7 CU T S

DO 16 K=1,NY 3 J=K+IY

CALL GETARAY (STOP ,NX)

DO 16 L= 1 , NX 3 M = L I

EI(M)=EI(M)+EXPY(K)^TTOP(L)

16 ET( J)=EI( JJ +EXPX(L) *STOR(L)

C INITIALIZE X - C OMPONE

T PARAMETERS FOP FF T I NT r GR A I ON

N V Y = N X 3 M M = M M X 3 NN= NNX 3 J=

17 1 2= X X + NX Y 3 I1=IX+1 3 I

p (M00(NXY/2, 2) .EO.O) 11 = 11 + 1

C PRE-FFT MULTIPLICATION Q Y (-1>**N

DO 16 M=I1,I2,2

18 EI(M)=-EI(M)

C FFT COMPUTATION

CALL CCFFT2 (El CMX+1) ,MM ,-FACTQP,-i)

1 1 = M X + 2 3 I 2= MX +NN

C POST -FFT MULTIPLICATION BY

DO 19 M=X1,I2,2

19 El C M » = - E I ( M

IF( J.EO.O) 20,21

C INITIALIZE Y-COMPONENT PARAMETER 7 FOR FFT I NTGRA7

X ON

20 NXY = NY 3 MM=MMY 3 NN= NNY 3 mx=MX+NNX 3 IX= T Y 3 J = 1 3 GO T 0 17

21 CONTINUE

C INITIALIZE Y-PLANE-CUT PARAMETERS FOR VECTO=-7 OM°OM r NT AND PATTERN COMP*

NN= N NX 3 LL =LL X 3 DEL = DEL X 3 MX=J=3 3 AN GL c = r V - S I N E 3 V V = V

* * 2 3 C 7 p E = S C PT ( 1 . 0 - VV »

C RESTRICT COMPUTATIONS TO p OL AR SPATIAL-FREQUENCY ANGLES IN VTSIBL" D A N G

22 NXY=NN/2+l 3 IF ( NXY , LT . LL +2 ) ?3,?C

23 Il-MX + 1 3 1 2=MX + NN 3 GO TO 25

2C I1 = NXY-LL + MX 3 I2=Xi+2*LL

25 K=I1-1 3 NXY=NXY+ M X

28 .
C VEC'CR -COMPONENT COMPUTATION

DO 33 M = 1 1 , 12 $ L=M -K 3 LL=M+NP 3 U= (H-M XY) * CEL f COST=U**2+VV

IF(J.EQ.O) 26,27

26 CC?E=CO°E $ COST=SQPT (l ,r -COS" ) /COSE 3 STNT =U*V/OOSE 3 GO T 0 3

27 IF(U.EQ.l.G) 28,29

28 CCST = i.: 3 SINT = COSF=C • t IF(V.EQ,E,) 30,11

29 COSE = SQPT (1. C-U**2) 3 COST=SQ° I (1,0-COST I/O O rr 3 S T NT=U *V/CO ^E

30 TFE T

A <L) =57,2957795130 8 2* A c IN ( U/CCP r ) El (M)= (C • ,1. 0) *<EI(M) *COSE+EI ( LL)*SINT)

33 El (LL) =CMPLX (C . ,COST) +E1 ILL) <3 LL=Ii+ND

C VECTOR-COMPONENT PLOTTING

CALL PA MPLO T (L » El ( II ) , THETA ,D A TA , 2* NN , BA NN r R , A NGLI)

CALL PA'IPLOT (L , El (LL) , T EE T A , OA T A , 2* NN , PA NNF° , 2, A NGL r >

C ANTENNA -PA T T-RN COMPUTATION

DC 34 1=11,12 3 L = m-K 3 L L= M+M C 34 DATA (L) = SQ°T (CEAL (El ( M) ^CCNJG ( =T (M) ) fCT ( ll> *CONJG (EKLL ) ) ) )

c ant enna-pattepn plotting

CALL PAM PLOT ( L » 0 • ,THF T A,DATA,L , BANNER ,1, ANGL r ) IF(J.EQ.O) 35,1

C INITIALIZE X-PLANE-CUT PARAMETERS FOR VEC T OP-CO M PON ir NT AND PATTERN COMP,

35 NN = NNY $ LL =LL Y 3 DEL= n ELY 3 MX=NNX $ J=1 I A NGL E= A V=S I N A 3 VV=V**2 $ CO P E= S 0= T ( 1 . C - VV ) 3 GO ’D 22

36 C0N T INUE

END

SUBROUTINE PA MPLOT ( M, £=>, XVALUE , 0 A T A , N2, BANNER ,MnnE, ANGLE)

DIMENSION E

c (N) ,XVALUE(N) ,C£TA (N?)

COMPLEX E P

DIMENSION HE A D(9) ,YMAX( 2) ,YMIN (2) , A M P ( 4 ) , C A M { 2)

DIMENSION BP (4) ,ELA7( 2) , CCMPff

)

DATA (COMP(M) ,M=1,6')/1GH AZIMUTH ,5HCOMR . , l'"H c L r VATION ,5HCOMP.,10

AH AN r ENNA PA ,5HT7ERN/ , ( AMP (M) , M =l,fc| /1CH AM P L TTUD E, 6H PHASE, 1 0 H A

M

D R/FHASE, II H MAGNITUDE/, (EL A7(M) ,M = 1,2) /IS HEIEV

A

T 1 0N=,1DH AZIMU T H = /

DAT A (CAM ( i)

i:h R r A L C A )

,

(CA M ( 2 ) = 1 CHCO MPL r X DA)

DATA ( K = G ) , ( I A = 1 ) , (I Q =i)

99

FORMA" ( 2X11 , 3F7 .1 ,I T ) 100 FCPMA t (4HAMP=,G8,2,A1E,F3.0,A:'!,A5)

101 FORMAT (-Aid

)

102 FORMAT ( 1HG , 1GX1H* ,9 A IP , 1H*)

103 FORMA' (1H1, 15X25HPLO" POU T INE JOB COMPLETE, 1 G Y3? HNUM9E P O p NON-TRI

ZVJAL DATA P

T

S = , I G//2 6X 31 HNUMB

C

p DATA points ppINTED =,I c //>

104 FCRMATdHl, 10X2GHFRROR EXIT - PAMp|_OT//>

105 FCPMA T (1HQ, 3X7HNO=,I?,^X5HXKIM = ,F7» ?,7X"HX''AX = ,Ft.S, 3X«HYMTN<1) = 1,F7. 3, 3X8HYMAX ( 1 ) = , F 7 . 3, 3X8HYM IN ( 2) = , F7

.

7 ,~X8 HY M AX < 2) = , F7 . 3 , 3X, i*H

2 PRINT SPACING=,I3//)

1 C 6 FORMAT ( 1H*//25X44HALL INPUT p ATA VALUES tq pA

x*PL CT AR p TRIVIAL'/)

107 FORMAT (5 (F10«4,G10*4, F7 .3)

)

108 FCF.MA r ( lH'J, 25X24HMAXIMUM CA'A A MPLI T UPE =,G1 7 .£//)

109 FCRMA T (50X, A1G ,16HT A PLOTTING MODE//)

11C FCRMA T (8 (F8.2,G9. 3)

)

29 .
C ROUTINE TO SET UP ARRAYS FO° FLO T TTNG AMPLITUDE AND PHASE

OF THE COMPLEX AFRAY E R . DIMENSION 0F ARRAY E c EQUALS N,

C DATA PLOTTEO VERSUS ABSCISSA ARRAY XVALLJE.

C

c

c

c

ccc

cc

c

c

c

cc

cc

c

c

c

c

ccc

c

oc

cc

1/c

cc

cc

c

c

c

cc

cc

c

c

c

ccc

c

c

c

cc

COMPLEX INPUT DATA ASSUMED IF MODE- 2 • IF MQOE=i, HOWEVER,

IT IS ASSUMED THAT ONLY REAL-AMPL I

r UPE DATA IS TO BE PLOT'FD.

T HI S F E AL -A MPLT UCE D A " A WILL BE A

c SUM !r D TO BE IN ARFAY DATA*

ARRAY DATA IS ALWAYS CHANGED BY the PROGRAM, WHILE APRAY c p

IS LEFT AS I T WAS ON ENTRY TO THE PROGRAM, IF MOOE=i, AR° A

ER IS IGNORED.

ROUTINE ALSO PRINTS OUT THE AMPLITUDE, PHASE, AND ABSCISSA

VALUE AT A SELECTED NUMBER OE DATA POINTS. USEP SPECIFIES

THE POINT SPACING BETWEEN THF DA 'A POINTS TO BE PRINTED AND

THEN DATA POINTS WILL BE SELECTED FOP PRINTING UNIFOPMLY

DISTRIBUTED OVER THE ABSCISSA RANGE X MIN TO X M AX CONSISTENT

WITH THE DAT£-PCTNt SPACING SELEC t ED. NOTE THAT THE

PECTPPOCAL OF THE PA t A-POTNT SPACING pAPAMETE 0 EQUALS T H

PERCENTAGE OF DATA FOIN T S TO PE PL OT T r D

AMPLITUDE OF ARRAY EP (DATA) IS PLOTTER TN D8 BELOW MAXIMUM ARRAY

VALUE. THE DB RANGE IS SPECIFIED BY YM TN (i ) (MUST BE NEGATIVE),

PHASE OF ARRAY E

r IS PLOTTED TN DFG'EER r FT W E r N YM AX (?) AN'*' YM I N ( 2

BOTH PHASE AND AMP PLOTTED ON SANE PLOT IF NO=2.

PHASE AND AMP PLOTTED c EPEP ATELY T r WO=l.

A U SE C -SUPPLI EO 33-CHAPACTER GRAPH LA ocr L IS PCIN T EO ON EACH PLO".

T HIS CAPTION IS OBTAINED BY READING ONE CAPTION CARO THE ^IRST

TIM- THE SUBROUTINE IS PALLED, WHILE T HE SAME 30-CHAPACTFR LEGEND

WILL SIM C LY BE USED OVER AGAIN ON SUBSEQUENT CALLS TO T H E ROUTINE

ALSO, A TEN-CHARACTER US CC MESSAGE -PANNES- IS PRINTED ON C A CH

G D A P H , ALONG WI’H THE MAXIMUM DB -LEVEL O r T HE n ATA AND A 4Q«

CHARACTE 0 DESCRIPTION OP WHAT TYPE OF GRAPH IS BEING PLOTTED.

THIS LATTER DISCRIPTION IS OBTAINED FROM THE PRESUMED FALLING

SEQUENCE OF T HE MAIN PROGRAM, FA PEA • ^H r VALUE OF THE CONSTAN T

AZIMUTH 0- EL E V

T ION ANGLE -ANGLE- FQP t HF GRAPH IS ALSO PRINTED.

A TO T AL LEGEND OF cr SHAPACTEPS IS °P T N TE D ON EACH GFARH,

USER SUPPLIES DUMMY A ° R A Y DATA FQP INTERMEDIATE STOP AG r WHEN

MODE = 2. T N

r HI S CASE, T HE DIMENSION OF ARPAY RATA EQUALS N2=2*N»

WHEN MOCE = i, THE INPU" IS ASSUMED '0 B

r T M ARRAY QA T A AND ITS

DIMENSION IS N2=N.

DATA CARDS READ ON FIRST CALL '0 SUBROUTINE ONLY

FIRST INPUT CARP.. NOO=NUMBER OR CU°VES CP R GRAPH (i Qp 25

XMTN, XM AX = SMALLEST AND GREA T EST ABSCISSA POINT

YMIN(i)=AMo L XTUDE FLO t LOWEST DB VALUE

J=POIN? SPACING BETWEEN ADJACENT PRINTED DATA VALUES*

(IF J= 0 NONE OF THE D A T A VALUES ARE PRINTED)

SUBROUTINE ALSO READS A 30-CHARAC T EP L r GEND CARO.

IF(K 8 EQ,3) 13, 14

13 READ 99,N00,X M IN,XMAX ,YMIN(1) ,

14 YMIN(2)=S. S YMAX ( 2) = 36 u • YMTN ( 2) =PHAS E PLOT SMALLEST NU MPE° OF DEGREES (C.)

YMAX (2) =PHAS E PLOT LARGEST NUMBER OF DEGREES (362.)

IF (N2.LT.2»N.AN0.M0DP,EQ.2.0R. N2.LT.N. AND.MQDE.EO.i) 1,2

1 PRINT IS 4 t CALL EXIT

30 .
2 If(XMftX.LE.XMIN.O r '.YMINdJ.GE. " . .OR. YMAX

(

2) . L " . YMIN ( 2 > ) 1,3

3 HEAD (9) CANNES $ Y M 4 X (II = 0 * $ 1 = 0 t KL=KI = 1 * ft NORM - 0 , $ NP=C

P c I NT 1C 5 ,NOO,XMIN,XMAX, YMTNtl ) , YMAX ( 1 } , V M

T N ( ?) , Y M ft X ( 2 ) , J

PRINT 1C9,CAM(M0DE)

IP ( MOD E. c C» 1 ) GC T 0 22 t I F ( MO OS . NE

2) GC T 0 1 T NO = NOO f LL = 0 DO 6 M=i,N

IF(XVALU£ (M> •LE.XMIN) KI = M $ I F ( XVAL U£ (M > . LE . XM A X ) KL = M

OATA (M) =C'ABS C£R(M)

)

IF<DAPA(M) . PQ. C.) 4,5

4 DATA (M+N) =YMIM (2) $ GO TO 6

5 1=1+1

IF (DATA (M I . GT. ANOPM) AN CRM = DAT A ( M)

9= A TAN 2 ( AIM AG ( EF ( M) ) , PEAL

(

rr (M ) )

)

IFtB.LT.j .) B = e + 6.28318T3il 7 lTa<:

DA T A (M+N) =57,295779513182*8

6 continue

IF(I.EQ.O) 7,6

•7 PRINT 10 6 T GO “0 21

6 P = INT 1 C 0 , A NO r M

IF(J,LE.C) 29,25

25 IF(J.IO.l) 26,27

26 I1 = KI $ 1 2= KL S NP = KL - K I + 1 $ GO to 28

21 M=KL-KI+1 « L = ( M - 1 ) /

NP=L/J % L = J * N P 3 I2=M/2+KI $ Ti=I2-L $ I2=T?+L * NP=?+NP+i

26 IF (MOTE. EG. 1) 3C,31

30 PFINT 110, ( XVALUE (M) , DATA (M),M=I1,I2, J) f GO X D 29

31 P C I NT 1Q7, (XVALUE(M) , DATA (M), DATA (M+N) ,M = I1,X2,J)

29 A N O 9 M = 2 3 . 0* ALOG1C (ANC D M)

ENCODE (u 3 ,1 CO , SB > ANOPM, ELAZ(T°) , ANGLE,COM D ( T A) , COMP ( I ft +D

DECODE iu ; , 1C1,BB) (

r AO ( M ) , M=t , A

)

DO 11 M=KI,KL

IF(DA T A(M).EQ.C. ) 9,10

9 OAT A (M) =YMIN (1 ) % GO TO 11

10 DATA (M) = 2C , C *ALOGlC (DATA («)) -ANOPM

11 CONTINUE

IF(K.EQ.O) 15,16

15 K = 1 T READ 101, (HEAD(M) , V = 6 , 0

)

16 HEAD (5)=AMP (3) % L=1 S 9=1H, T TF(NC-2) 10,19,1

17 P^INT 102, (HfAC(M) , M = 1 , 9

)

18 HEAD (5 )=AMP (L+LL) 3 B=1H+

19 JC = (L-l) *N + KI $ N = L t tp (NO, EG, 2) M=MOD(L,2)+l t T2= (M-i) *N+KI

CALL CPTPLOT (XVALUE (KI ) ,C

T A( J?) ,XMAX, XMIN,

! MX, YMIN,KL-KT+i ,HE AD,

X L , NO)

CALL PLT120R ( XVALUE <*1 ) , CATAf I 2) ,XM«X,XMIN, Y«AX (M) , YMIN (M) , KL-KI+l

Y , 8 , L , NO ) ? IF (L.EQ.2.0F. NODE. EQ.l) 21,2

C

20 L = 2 % B= 1 H+ T IF(NO-i) 1,17,19

21 PRINT 1C2, (HEAD (M) ,M= 1,9)

PF'INT 13 3, I. NP

GC TO 32

22 DC 24 M= 1 ,

IF(XVALUE(M) .LE.XMIN) KI = M $ I F ( XVALUE (M ). LE . XMA X ) «L= M

IF ( D AT A ( M ) . EG . C . ) 24,27

23 1=1+1

IF ( DATA (M) , CT. ANOPM) ANCRM=OA^A(M)

24 CONTINUE

NC=1 S LL =3 $ IF(T.EO.I) 7,8

32 IF (MODE. EQ, 1) 33, *4

33 IP = MOD(IB ,2) +1 $ I A = 1 * GO x O

34 I A = I A + 2

35 RETURN

END

31 .
SUBROUTINE PL' 1 2C E ( X , V, XNAX, XMIN, YMAX, Y M T N » LAST, I SYMBOL , NOPL“,'I 1, MOST) PL t I

C MOO I FT ED li/4/68 P L T i

DIMENSION Y(l), Y ( 1 ) , 7X(13), GRAPH(12i, 51) PLTi

INTEGER GRAPH, COLUMNS, BLANK, BORDER PL t 1

DATA (LINES = 51), (COLUMNS = ill) PL T 1 YLAR = YM AX PLTi

Y e MA = Y M IN PLTi

YSCALI = ( YL A c - YSMA) / (LINE P - 1.) PL T 1

IF (NO ,N£. 1 ) GO TO 190 PLTi

KM A X - COLUMNS / 10 + 1 PLTi

X L A R = XMAX P L T i

XSMA = XMIN PL T 1

EOF 0 E

p = i H I PLTI

BLANK = 1H PLTI

M

A

T R I X = COLUMNS 3,1 LINES PLTi

IF ( MA T p

I

X . L

~ . 1) GH TQ IRC PL^l

dc i:: T = i, mat ix plti

liO G 5 A P H ( I ) = BLANK PLTI

120 CONTINUE PL r l IF (LINES .LT. l) GO TO 140 Pl T l

DO 13.. * = 1, L I N - S PLTi

130 G

t A P H ( 1 , I) = G P A PH (COLUMNS, I) = BORDER PLTI

1-

0 CONTINUE PI ^

i

IF (COLUMNS .LT. 1 ) GO T 0 i€G PL T i

CC 15 1 1=1, COLUMNS PLTI

l': C GRA PH ( 1 , 26) = 1H. PLTi

160 CONTINUE PL T i

XSCALE = ( XL AS - XSMA) / (COLUMNS - 1.) PL T 1 IF ( KM AX . L ” « 1) GO 'O 18? PLTi

DO 1 ?.: K = 1, KM A X PLTi

170 ZXCKI = 1 C. * FLOAT (K - 1 ) * XSCALE X CM A PLTI

180 CONTINUE PLTi

190 IF ( L AS T .LT. i) SO ”0 257 PL T 1

DO 24 E I = 1 , LAS' P L T

i

IF (X(I) .GT. XL A - .CP. X(I) . L” • XS^O) GO ’S 24" °IT

i

IF ( Y ( I ) .G T

. Y L A , 0 P . Y(T) ,L T . YS M A ) GO ”0 ?4 3 PL~i

IX = ( X ( I ) - XSMA) / XSCALE + 1.5 PLTi

IY = ( Y ( I ) - YSMA > / Y SCALE .5 PLTi

IY = LINES - T Y PL T

i GRAPHdX, IY) = I Y M p 0 L PLT

i

2-

0 CONTINUE PLTI

250 CON T I N UE PLTi

IF (NO . NE . MOST) s p TURN PL T i

PRINT 1G0C PL T i

YFS = YLAR + YSCALB PLTi

IF (LINES .LT. i) GO ”0 27 : »L T i

00 26. I = 1, LINES PLTi

YES = YES - YSCALI PL^l

PRINT i E i D , YES, ( G R A P H ( J , I), J = 1, COLUMN*) PLTi

260 CONTINUE PL'5'! 2 <’ 0 CONTINUE P L T

i

PRINT 1520 PL T 1 P c I NT l

c 32, 7 X PL r l

RETURN D LT 1

151 £! FCPMA T (lHi,9X,2-(5HT....)iHT» °LT i

1510 FORMAT (1H ,

E

8 . 2 , 1 X , 1 2 i A 1 ) D L T 1

1570 FORMA” (1H , PX , 2 4 ( 5 H T . . . .

)

1 H* 5 PL T i

15 3 Q F CP M

A

T (1H ,?X,13(1X , F 9 • 3)

)

e L7i

END °LTi

32 .

•tv

n-

-x.

n

•»

n

.n,

^

?

,• .

rf

4-' v

rE°A

f

-A

a

A

a

a

a

um

aA

i\>
SUBROUTINE CCFFT 2 CC , M , S C , NX ) 1 DISCRETE COMPLEX P AST FOURI c R TRANSFORM. 2

L. DAVID LEWIS, NOAA-SEL, T 3!)31E 3 CALL CCFFT 2 ( C , M , SC , N X ) 4

COMPLEX INPUT C(J) IN NORMAL SEQUENCE. ' COMPLEX OUTPUT C ( K ) IN NORM AL SEQUENCF. f> SEQUENCE LENGTH IS N = 2**M 7

SC IS PEAL SCALING MULTIPLIES. « NX IS THE SIGN OF THE EXPONENT IN T H r TRANSFORM orriNTTION. 0

USES CFFTPCO AND PEVBIND, 0. V. 11

OTMENSION C(2) 11

CALL REV3IND(C(i),CC?),M,2) 12

CALL CFPTPCC (C (1) ,C (?) ,M,SC,NX ,2) IT

RETURN 14

END

SUBROUTINE REV RIND f A , R , NM , NO I L

CALL REVBIND(A,B,M,ND)

REVERSIBLE PERMUTATION OF ARRAYS A AND D FROM NORMAL SEQUENCE TO PRVFRSE BINARY SEQUENCE,

OR VICE V E~SA.

ND IS SU P SCR IRT INCREMEN T FO^ A, P.

SEQUENCE L

r NGTH IS N = 2**M

WRT t teN BY L. DAVID LFMS AND MARIE WEST, r SSA. . MODIEIED F D OM , OR INSPIRED BY THE ALGOL P c OCcOURE

REVERSEBINAP Y , BY R. C. cr NGL r TQN, C RI.

SUBROUTINE CFFT RCD ( A , 0 , MM , SCAL E , N’EX P , NOEL)

DISCRETE COMPLEX FAS T FCURIER TRANSFORM.

CALL CFFTRCD(A,B,M,SC,NX,ND ) OUTPUT Am 4 I*B(K) IN NORMAL S

r QUENCS.

ND IS SUBSCRIP’ INCREMENT ^0* A, P,

SEQUENCE LENGTH IS N = ?**M

SC IS REAL SCALING MULTIPLES.

NX IS THE SIGN OF t H£ EXPONENT IN T HE TRANSFORM DEFINI t ION.

INNER LOOP SINES AND COSIN'? CCMPU T ED

RECURSIVELY BY SINGLETON*? 2ND-D 1RFF P “NC F ALGORITHM,

INITIALIZED FROM A BATA TABLE.

DIMENSION A (163E4) ,B(16T6R) 1

COMMON /RFTCC/ M, JD (1R) , ST( IF) 12

DIMENSION JO ( l 5

, 7 i CALL c OLLCALL UeHc/4/71 ) 14

M^MM % CALL FRTC T IF(*.L£.i) "TUEN 1 7 x y no=no-:l

00 13 : L C = 1 , 1

I’

JC(LC) = J Q (LC) * N 0 1

* N - J C ( M 4 1

) - N 0 + 1 S NP = N4 1 1 B K-l S I = N 0 4 1 $ JsN'-NO ? 1 LC=M 21

K = K+ JC ( L C ) T JC t L C) =” JC (LC) 2 ° IFCJC(LC) • L T • C ) GO TO 4 2 3

IFCLC.EQ.2) RETURN 24

L C = L C - 1 3 GO t o Z IF(K.LE.I.OF, J.LT.K) GO T C E SB

T = A ( I ) 3 A ( I) =A <K ) 3 A m=T 7 7

T = B(I) 3 B ( T ) = B ( K ) $ B(K)=T S 3 IFCJ.EQ.K) GO TO E 29

kk=np-k 3 1 T=A

KK ) * A ( KK ) = A ( J ) 3 A( J) =7 7 1 T-B(KK) 3 B(KK)=8(J> 3 B(J)=t 3 2 I ~

1 4 N 0 3 J=J-ND 3 GC T 0 2 3T

END 7 4

1

->7

B

7

s9

11

12

33
WRITTEN BY L. DAVID LEWIS A NR MARIE WEST, ^SS A

1 3 modified FROM, OR INSPIRED q Y THF ALGOL P r OCEPU Q E 1 4 REVERSE r OURIERC

,

BY P. C. SINGLETON, SRI. 13

DIMENSION A (16384) ,B ( 16384) 16

COMMON / FFTCC/ M, JD ( 1 5 ) , S ( IS) l 7

CALL ROLLCALL (48H6/4/71 ) 19

M = MM % CALL FFTC 19

N C=NOEL 20

N = JO (M + l

21

K = N/4 7 O

NO = K 9 1 JSP AN=NQ 24

ND2=2*ND ?E

NLIM=i4(N-i )*ND $ NN= NLI M -1 ££

SC = SCALE 2

» IF ( ARS( P C-1 . ) ,LT. 1. E-l 1 ) GOTO 7 2 9 DO 5 JC=l,MLIM,NO 29

A(JC) = SC * A ( JC

T 0 B ( JC) = SC * B ( J C

)

31

CONTINUE -? 2

IFCM.EO.C) RETURN 33

CO 10 KK = i, NLIM, ND2 ? 7 : K c = KK + NO ; RE = A ( K K ) « A ( K S

To

A ( KK ) = A ( KK } A (KS) 7 7

A(KS) - PE 3 3 FIM = B ( KK ) - q ( KS ) T9

B(KK) = 9 ( K K ) f B(K3) '* 3 B ( KS ) = FIM 41

CONTINUE 4 2 IF(M.EQ.i) R E

r U F N U T EX D S = ISIGN(1 ,N'XP) 7 4

CC 90 J

q = 2,M 46

sr = -sud -i> /. 7 CD - 2. * S ( JB) * q ( J q ) !i9

F = ” 2 e * CD 4 9

CM = 1. F 1 c m = •3 i

SN = "

F 2

JJ - 2 3 3

KK = 1 E 4 i

12

20

SM - 4- E XPS

J S ° A NH = JS ° AN

JSPAN = JSPAN 4- ISP A

1 KS = KK 4 JSPAN

RE = ON * A ( KS

- SN

FIM = SM <• A ( KS ) 4 C' A ( KS) -

A ( K K) - q 7

A ( KK) = A ( KK ) q R B (KS) - B(KK) - FIN

B( KK) = BCKK) 4 FIM

KK = KK JSPANH

KS - KS 4 J ; PA N

F T M = S M * A ( KS ) 4 Cl RE = CM * A(K p ) - SN

A ( KS) = A ( KK) - RE

A ( K K ) -

A ( KK) 4 9 E B (KS) B ( K

47

) - FIM

B (KK) = B ( KK) 4 FIN

* 3 ( KS

)

! * q ( KS

q(t<C)

9 (KS)

7- >

3 7

3 9 ~ 3

s 7 w -i

5j 7

5 7

3 A

5 3 5 6

^

•*»

C, A 3 9

7

3 7

1 7 7

34 .
3C

8C

9C

KK = KS + JSPANH 7 7

IF( KK.LT. NLIM) GO T 0 20 ' *« KK = KK - NN 75

JJ = JJ K

T 5

IF(JJ.GE.NQ) GO TP a: 7*’ CD = P * CN + CD 78

CN = CD + CN ”9

SN= CN* EXPS 8 0 SD = R * CM + SD

CM = SD + CM 3 7 SN=-CM*EXPS 8 T GO TO 20 94

K = K/2 97

CCNTIMUE 97

RETU C N 83

EMC 39

SUBROUTINE FFTC 1? COMMON SUBROUTINE FOR FFT SUBROUTINES. 7 JC IS POWEFS-OF-TWO ARRAY. . JC(M)=2*MM -1

)

4

ST IS SINE ARRAY. . $ T ( M ) =S I N ( PI/ ( 2**M )

)

5

'1 IS T ES TEH FOR PROPER INPtJT RANGE, ! ,L~ .M.LE . 1A. 6

7 COMMON /FFTCC/ M, JC (15 ) ,S 7 (13) S DATA ( J C= 1, 2,4,8,16,32,54,128, 256,512,1924, *; 48, A" 96, 8192,1638-) 9

data (ST = 1.390000 3 90GCE+ 000, 7 . -: T

i oe73ii8?E-aai, 1 0 3, 82683A37765E-U 01, 1, n ~ 1 5

r 7 27 r 16E-9-31, 11

9 , S CITING 3 296E- G 3 2, 4. ° ’ 6‘7 6’T 432 T 4E- 3 D 2, 12

2. 4541228- 229E-G C 2, 1. t??1E7828E7 c -1j2, 1 7

6 . i 3 5 531* 0 - B1EE“ ’0 3, 3 . ’ 6 9 5 6 7

6

2 9 7

£

3 j 7 , 14

i . 5 339 fl r 18 629”-

I

r 3, 7. 66 95 031 374 3£- 3 04, 15

3

. *34351*7 571E-C 04, i. 91747597311^-3.84, 16

9. 5873799* B 6 f E - C 05) 17

IF(M,LT,0 .OR.M.GT.l^) CALL c Y7 T EW(52,iGHM T L L EG A L .

)

1 8 IF(M.LT. J.OR.m.GT.14) CALL 080 ERPCR ( 0,10 HM TLL“G AL. ) 18

RETURN 19

ENO 2

1

35 .
APPENDIX C

Computer Program Modifications to obtain Theta and Phi Vector-Component

Output Instead of Azimuth and Elevation Vector-Component Output

I. Required changes to VECTOR-COMPONENT-COMPUTATION section of program FAREA

replace DO loop, DO 33 M-Il, 12 thru statement 33, with the following:

VECTCF-CCMP ON::n t ccm d UT A T I 0 N

DC 33 M= Hi,12 $L=M -K ILL - M + NO tu= ( M

• N XY) *0 r L 7TF(VV*tO.O.

46 CC2T=tJ** 2 T S I NT H = U % IF ( SINT H . eo. , ) 26

26 IF ( J. EQ. G) 27, 2P

27 ccs p=i .

$ SINF=C

* GO T 0 32

28 SIN°= 1 . c $ C 0 c P = 0 . $ GO T C 3?

47 C0ST=tJ** 2 *VV * SIT. TH = SQP T

('' 0 S

?

) 29 IF ( J .='0. : > 30,31

30 CCS P = U/S INI' H 7 SIN D = V/f I K *H S GO T 0 3 2

31 S IN P=U/S TNTH t ro~ P = V/SI N

32 THETA (LI = 3’ .295773 El 3 0 3 2

4? AS IN f U/OCPP ) Z7= (G . » 1 . : i

I

(EX (Ml *C OS P + r KLU * S I N p ) El

LL) =C MPLXC. ,SQ PT (1.0 - COST) > * (-EI ( M) * STNP+TI f i_

> *COSP)

33 F. T < N ) = Z 7

* LL= 1 1 + ND

II. Subroutine PAMPLOT DATA-statement change; replace the first DATA

statement card with the following:

OCT A (C0 M P (M) ,M = 1, - ) A 0 OMP , F H ON ON H ’HI 00MO,gH0NFhr, i"
APPENDIX D

Computer Program for Computing Principal-Plane-Cut Far Fields

(This program ignores the jump change in sign of the electric-field vector when

passing through the coordinate-system origin. Refer to equations (10) and (11)

of Appendix A for details).

PROGRAM FARE UNPUT, OUTPUT # TAPE1=I NPUT)

DIMENSION ST0R( 1024) ,EU096) , THETA (1024) * DATA (2043)

EQUIVALENCE (THETA, STOR)

COMPLEX ST OR ,

£

DATA (NMAX=1024) , ( NOIMEN=409fc>

94 FORMAT (4F7 «Q,2XA10, 1X212)

95 FORM AT ( 1H1 * 2X5 2H ALP HA=X -COORDINATE APERTURE DIMENSION IN WA VELENGT

T HS , 9 X51HBE TA=Y-C00 ROI NAT E APERTURE DIMENSION IN W A V ELENGTHS/ 3X57HS

UX=R AT 10 OF WAVELENGTH TO X-COORDINATE DATA-POINT SPACING, 4X5 7HSY=R

V AT 1 0 OF WAVELENGTH TO Y-COORDINATE DATA-POINT SPACING/3X39HBANNER*

W10-CHARACT ER PLOT IDENTIFICATION, AX 39HNZ X, NZ Y=ZERO-FILL AREA DOUBL

XING FACTORS, 4X44H NNX , NNY=NUM6ER OF TERMS IN X- AND Y-FFT SUMS///8H

Y ALFhA =, F6.3,4X6H8ETA = , F8 ,3 , 4X4 HSX = ,F7 .3, 4X4HSY = ,F7. 3 ,4X7 HQ ANN

ZER.t A 10, 1H . ,*X4HNZX=,I2,3X4HNZY=, I2,7X5HNNX =,I5,4X5HNNY =,I5//)

96 FORMATC1HO .25X37HARRAY DIMENSIONS TOO SMALL OR ILLEGAL//1X5HNNX *,

AI5, 3X20HNNX=2**MMX , MMX = , I 3 , 1 5X 5HNNY =, I5,3X20HNNY*2**MMY , QMMY = , 13 , 1 5X22HC0MP UT ED E-DIMENSICN =,I6)

97 FORMAT (25X5 HMMX = , 13, 25X5HMMY =,I3//)

COMPUTATION OF THE FAR-F I ELD PRINCIPAL-PLANE CUTS OF THE

THETA- AND PHI-COMPONENTS OF THE ELECTRIC -FI ELD VECTOR OUE

TO A GIVEN NEAR-ZONE ELECTR IC-FIELD DISTRIBUTION IN A

SPECIFIED APERTURE, ALONG WITH COMPUTATION OF THE FAR-FIELD

PATTERNS FOR THESE PRINCIPAL-PLANE CUTS.

THE DATA SPACING BETWEEN X- OR Y-COOROINATE ELEMENTS IS ASSUMED TO BE FIXED

NM AX=C IMENS I ON OF ARRAYS STOR ANO THETA. DIMENSION OF ARRAY

OATA SHCULD EQUAL 2*NMAX. NOIMEN IS THE DIMENSION OF ARRAY E.

NDIMEN SHOULD EQUAL TWICE THE SUM OF THE X- ANO Y-FFT DIMENSIONS.

THE MAXIMUM USEFUL VALUE OF NOIMEN IS 4*NMAX.

INPUT PARAMETERS.. ALPHAS APERTURE WI0TH IN WAVELENGTHS (X-COOROINATE)

B £T A= APERTURE LENGTH IN WAVELENGTHS (Y-COORDINATE)

SX = RAT 1 0 OF WAVELENGTH TO X-COOROINATE DATA-POINT SPACING

SY=RATIC OF WAVELENGTH TO Y-COURCINATE DATA-POINT SPACING

BANNER = 10-CHA RACTcR G l-APH IDENTIFICATION

N ZX=NUMBER OF EXTRA DOUBLINGS OF THE ZERO-FILL REGION

FOR THE X-FA R-F IELO-C CORDI NA TE FFT • (NORMALLY*!)

NZY=NUHB£R OF EXTRA DOUBLINGS OF THE ZERO-FILL REGION

FOR THE Y- FA R-F IELO-COOROINATE FFT. ( NORMA LL Y = l)

THE INPUT DATA IS ASSUMED TO BE SUPPLIED AS ROWS OF X-COORDINATE

DATA, EACH =0W CORRESPONDING TO A FIXED VALUE OF THE Y-COORDINATE.

ALL OF THE E-SUB-X DATA IS TO BE SUPPLIED FIRST, AND THEN THIS

INPUT OATA IS TO BE FOLLOWED 8Y ALL OF THE E-SUB-Y DATA.

THE PROGRAM OBTAINS THE INPUT OATA THROUGH REPEATED CALLS TO

THE USER-SUPPLIEO SUBROUTINE GETAFAY.

37 .
C EACH CALL T 0 SUBROUTINE G^ARAY RESULTS T N a ROW OF NX -DUPLEX

C DATA POINTS, r CC H SUCCESSIVE CALL CORPUS ENDING TO A ROW

C SUCCESSIVELY FURTHER BELOW THE T OP or T H

r APERTURE. HERE* NX

C IS EOUAL TO THE FFODUC" SX * ALPHA. 0 HE NUMBER O

c ROWS. 0°

C the NUMBER OF CALLS TO GFTARAY, I c GIVEN BY NY = SY * BETA.

c

C NUMBER OF DATA CAROS USEO BY PROGRAM* SUBROUTINES = 7

,

C

C UPON COMPLETING A GIVEN S E

T OF F AR-FT rL 0 D RINCIPAL-PLA M E

C COMPUTATIONS, T HE COMPUTER PROGRAM WILL t hEN peTU°N t C

t HE

C BEGINING OF THE PROGRAM TO Ps fAO A W pata CARO AND START ALL

C OV E ° AGAIN ON A NEW F A F-FTELO COM°UT A "'TON. thts t IME, HOWEVER,

C THE PLOTTING ROUTINE PARAMETERS A NO L EG r NO TH A T WERE P r AD IN, ON

C THE FIRST tine t HROUGH r H

A

t SUBPOU~TN r , ARE SIMPLY RECALLED AND

C USED OVER AGAIN. THE US r R MUS T MA*S SURE SUBROUTINE G

r T AR A Y CAN

C SUPPLY A NEW SE T O

c NEAR-FIELC! DATA UPON RESTARTING T H

r ppOG c AM.

C

c PARAMETER °E AQ-IM

1 ROAD Hi. , AL°HA, BETA, SY,SY, BANNER, N7X.N7Y

IF ( EOF ( 1 ) .ME. C ) GO to ?g

C CALCULAT T ON CP NX A NO NY

NX = I

. "»0X*ALPHA 3 NX-2*NX*1 * NY= : . 5 * SY* B

E

’ A * NY=2*NY+1

sx= nx/ alpha $ SY = ^Y/ n EiA

C CALCULATION OF NMX AND NNY, ALONG WITH Mvy AND mmY

MMX=MMY=NNX=NNY=i

2 MMX = M M X + 1 3 NNX=2*NNX 3 I F ( NNX .

G

T . NX ) 3.2

3 NNX = 2*NNX 3 IP(N7X) - u , « , a 4P NNX=NNX/2 3 M M X = M N X - 1 3 GO T 0 F 4 00 5 M=i,N7X t MMX = MMX4-I

5 N NX = 2 * NN X

e MMY = MMY + i 3 NNY=2*NNY 3 IF ( NNY . C-^.NY ) 7 ,6

7 NNY = 2* NNY 3 ir(N7Y) - - , 1 u , o U5 NNY = NNY / 2 3 MMY = MiiY-i 3 GO T o 1”

g nr 9

v = 1 , N 7 Y 3 mMY=MMY *1

9 NNY =2* NNY

1C MM = M AX C ( MM X ,MMY) I NN = MAXC (NNX. NNY) 3 MP = N M Y + NN Y ? J=2*N0

If ( J . GT . NO I MEN . OR . NN. GT . NM AX. 0 p , nm

.

0- , i A ) 11,12

11 P C IMT 9E,fJNX,MMX,NNY, MM Y, „ * CALL 'XI T

C INPUT RAG A M ETEP FRTN'-OU T 12 Pf- I NT 9P , ALPHA, BETA ,SX,SY ,BANN-P,NZX ,N7Y , KNX, NNY

P

r I NT 9 ^ ,M m X,MMY

C COMPUTATION OF A P E F T UR E DA T A-°OIN'p SPACING:

D t L X = S X /N NX 3 LLX = 1, /O ELY 3 0 EL Y = S Y/NNY 3 L L Y= l . o / 0 E L

Y

C COMPUTATION OF A P£P T UPE -0 A T A PO c I

r TON WI t HTN ~H r INTEGRATION IN^povAL

FAC T OR= 1. C/ «SX*SY ) 3 I 3X= (NNX-NX + 1) /2 3 J 3 Y= ( NNY - N Y *1 ) / 2 *MNX

C PR E “ZERO I NG CF INTEGRAL FIELD

CC 13 M = 1 ,

J

13 E

( M ) = (*:., C . )

C DO LOOP “0 COMPUTE T W0 PLANE - r UT T N T EG C ALS CV’jr- AP£

t URE-F TELO VECTOR

o° 2 f i=i ,2 * mx=nd*(T-d t ix=mx+i:x 3 TY= m X+t:y

C SET -UP OF INTEGRAND - COLLAPSE ON X- AND Y- PL IN'

- OlJ~ c DC 14 K - 1 , N Y * J=R*IV

38 .
CALL GETARAYCFTCR.NXI

DO 14 1=1, NX f f‘ = L + TX

E(M)=E(M) + c T 0 p ( L ) 14 E(J ) = -. < J) +STOR (L)

C INITIALIZE Y-PLANE-CU? PA D A NRT-oc

NXY = NX t MM=MMX 3 NN=NNX f n EL=OELX ? L L = L l X FACTO~=-FACTOF $ J=C

15 NMIQ=NN/2+l T I2=IXfNXY 3 T1=IX+1 $ IF (MOD (NX Y/ 2 , 2) . EG . 0 ) 11=11+1

C PRE-FFT MUL T IPLICATION D Y (- 1 )**M

DC 16 M= 1 1 , T 2 ,

2

16 E (M) = - E

(

M )

C FFT COMPUTATION

CALL CCFFT 2 (E(MX+1) ,MM,FACTn P, -D

IL=MX+2 T I 2 =MX+NN % JJ=KCD ( J+ T -1 ,2)

C POST-FFT multiplication p Y ( - 1 ) * * J DC 17 M=I1,T2,2

17 E(M)=-E(M)

C REST-ICT COMPUTATIONS r O REAL POL A© SPATIAL -F D E n U

c N CY ANCLES

IF (LL. G” . NMID-2) 18,19

18 II = M X + 1 T J 2= MX + NN S KK=C T on T 0 2

19 KK=NMID-LL 3 I 1 = KK+MX S I2=Ii+?*Ll * KK= KK1

20 NMTD = NMI r

5 + HX * K=KK+UX

C VECT CP -COMPONENT COMPUTATION

00 23 M=I1,I2 $ L=M-K 3 LL = M-^t<

F IN T

H

T A = (M-NMTO) -PEL % -Hc t A ( L )c 7. 2957705^ 3

«2* A SIN (STNTHT A

)

IF(JJ.EQ.O) 21,22

21 E (LL > = (.., 1 . T

» *E ( ) * GO ~0 2 7 22 E (LL) =C^PLX ( C , ,S0- T

( 1 . : -SIN7H' 5**2) ) * E (“

)

23 CONTINUE

C VEC T 0P-C0MP0NENT PLOWING

CALL PA PLO T (l,E (MY+ 1 ) , T 4 FT A ,0 A T A , 2*N N , BA NN 7 =• , 2 ) IF(J.EQ.C) 24,25

C INITIALIZE X-PLANE-Cir P A R A 6 E ~ E

3 C

24 IF(I.EQ.l) FACTOR = -FAC'r OF * MY-^y+NNX S J =1 3 LX=L 3 KX = I 1 -N^I ", -1

NXY=NY 3 MM=M M Y 3 NN=NNY T $ LL=LLY * TX=IY * GO TO 1

E

25 LY=L

Il=NNX+l 3 1 2 = N N X + L Y

C ANT£NNA-PA'"r E

p N COMPUTATION

DC 26 M= 1 1 , I 2 $ J = M*\'Q

26 DA T A(M)=SQP'T

( c E A L (E ( M ) * C C N J G ( ~ f M > ) + E ( J

) * COM J G

(

T ( J )) )

)

C ANTENNA-PATTERN PLG T ~ING (X-PLANE CU T

) CALL PA PL OT ( L Y , G . , TH pt a , p A

T A ( 1 1 > ,LY , B A MNE R . 1

>

DC 2 7 M=i,LX 3 J=U+N9

SINTH T A= ( M + KX) *CELX 1 T HE T M8 ) = r T,2CC?7qr 1 ^o?«,/i r IN (SJMTHT A)

27 DATA(MI=SQRT(FEAL(E(M)*CCNJGC-:(M))+E ( J) » COM JO (r ( I) ) n

C A NT ENM A -P A T

r £F N PLOTTING (Y-FLAN- CUT)

CALL PAPL0T(LX,C*, T HETA,0*TA,LX,PANNEP,1)

GC TO 1 28 CCNTjnuENP

39 .
Appendix D: Modifications to subroutine PAMPLOT for

compatibility with principal-plane-cut program.

I. SUBROUTINE name change and SUBROUTINE argument-list change:

SUBROUTINE P A PLOT ( N , £ & , X V AL UE , D

T A , N 2 , B A NNE~ , NO ° E

II. Subroutine DATA-statement change:

DATA COMP ( M ) ,M=l,o) /lCHTHrlTA C0 M P , 5 H CN EN T , 1 1 H

~ 3 i>i a

, 5 P0NEM T, 1 ] A H A N T £ N N A P A , 5 H T T £ F N / , ( A M P ( m ) , M = 1 , A ) / 1 0 H A * P L T ' U 0 E . £. H F H A S 5 , 1 G H A M P /

13 P h A ' I E , 1 : ^MAGNITUDE / , ( EL A7 ( M) ,* = 1, 2 ) / P^T - £

, ^ = PJ/ 2/

III. FORMAT-statement change and change in ENCODE statement (the latter

statement is located just following statement 29):

1 l £ FORMAT ( H A MC = » G 9 . ^,2XA3,lXAi£ , A<r , 1 X )

ENCODE UI , 1C 3 , BB) A NORM , E LA 7 ( I B ) , COM P ( I A ) , COMP ( I A + 1

IV. Replace statements 32 thru 35 (located at end of subroutine) with the

following six statements:

3? IF (IA . rn— ' 9 1) 33,36

3 3 IF (13. EO. 1.) 34,3 6

3^ I A = 3 * I 3 = 2 T GO

3 5 I A = 5 $ GO TO 7 7 36 I b = M 0 0 (IB ,2) +1 $

37 TURN

E Q . 2 ) I Ari j jr (n'ODE. EO. 1 ) 1^1

40 .
HBS-114A (REV. 0-78)

U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET

1. PUBLICATION OR REPORT NO. NBSIR 81-1643

2. Gov't m No,

-A

.

'

, , -Ms--

'

Z, AccossIor No.

4. TITLE AND SUBTITLE

Efficient Computation of the Far Field Radiated by an Arbitrary

Rectangular-Aperture Distribution

5. Publication Date March 1981

6. Performing Organization Code

""ST

7. AUTHOR(S)

Richard Lewis

8. Performing Organ. Report No.

9. PERFORMING ORGANIZATION NAME AND ADDRESS

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, DC 20234

10.

11. Contract/Grant No.

12. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street, city, state, ZIP)

Electromagnetic Compatabil ity Analysis Center

Annapolis, Maryland 21402

13. Type of Report & Period Covered

*4. Sponsoring Agency Code

.!» 15. SUPPLEMENTARY NOTES

| | Document describes a computer program; SF-185, FIPS Software Summary, is attached.

16. ABSTRACT (A 200-word or less factual summary of moat significant information. If document includes a significant bibliography or

literature survey, mention it here.) This report contains the computer documentation for calculating the far-zone electric

field due to a user-prescribed electric-field distribution within a rectangular aperture

The far-field output is computed along two arbitrarily selected, perpendicular, spatial

frequency plane cuts. Program execution time is minimized by the use of fast Fourier

transform (FFT) processing. The program was designed so that the required far-field

output is obtained by processing only two, vector, one-dimensional FFTs. The far-field

results are obtained in the form of elevation and azimuth vector components and electric

field-vector magnitude. A complete analytical discussion of the problem is presented,

along with sample graphical output to illustrate how aliasing and output resolution

limitations effect the graphical results.

17. KEY WORDS (six to twelve entries ; alphabetical order; capi talize only the first letter of the first key word unless a proper name

;

separated by semicolons)

Algorithm; antenna; aperture; computer program; electromagnetic; far-field; FFT;

Fourier transform.

18. AVAILABILITY [~X1 Unlimited

| | For Official Distribution. Do Not Release to NTIS

| | Order From Sup. of Doc., U.S. Government Printing Office, Washington, DC

20402, UXOrder From National Technical Information Service (NTIS), Springfield, VA. 22161

19. SECURITY CLASS

(THIS REPORT)

UNCLASSIFIED

20. SECURITY CLASS

(THIS PAGE)

UNCLASSIFIED

21. NO. OF

PRINTED PAGES

45

22. Price

$6.50

USCOMM-DC
-