FINAL REPORT 

TO 

GODDARD SPACE FLIGHT CENTER 

ON 

FREQUENCY SELECTIVE REFLECTION AND TRANSMISSION 

AT A 

LAYER COMPOSED OF A PERIODIC DIELECTRIC 

b 

PREPARED UNDER 

GRANT NAG 5-714 

(NASA-CB- 180274) E3ECUENCY SELECIIVE N67- 19554 

REFLECTl.CN AtiL! !IBANSflISSlCE A'1 A LAYER 

CCIYPOSED OF A PERIODIC DIELEClfIC Final 

&€Fort, 1 Octo 14E5 - 1 Cct, 196k Onclas 

{Eolytechnic Inst, of New Ycrk, Erocklyn.) G3/32 43544 
Polytechnic University 

Department of Electrical Engineering and Computer Science 

Brooklyn, NY 11201 

Final Report 

to 

Goddard Space Flight Center 

on 

Frequency Selective Reflection and Transmission 

at a 

Layer Composed of a Periodic Dielectric 

Prepared Under 

Grant NAG 5-714 

Henry L. Bertoni 

Li-hsiang S. Cheo 

Theodor Tamir 

Oct. 1, 1985 - OCt. ' 1, 1986 
li 

Abstract 

This report examines the feasibility of using a periodic 

dielectric layer, composed of alternating bars having dielectric 

constants el and e2, as a frequency selective sub-reflector in 

order to permit feed separation in large aperture reflecting 

antenna systems. For oblique incidence, it is found that total 

transmission and total reflection can be obtained at different 

frequencies for proper choice of el, e2 and the geometric 

parameters. The frequencies of’ total reflection and transmission 

can be estimated form wave phenomena occurring in a layer of 

uniform dielectric constant equal to the average for the periodic 

layers. About some of the frequencies of total transmission, the 

bandwidth for 90% transmission is found to be 40%. However, the 

bandwidth for 90% reflection is always found to be much narrower: 

the greatest value found being 2.5%. 

I 
I. Introduction 

Separation of the feed structures for different frequency bands in 

large reflecting antennas has been achieved using sub-reflectors whose 

transmission and reflection coefficients are frequency dependent. For 

frequencies in one band, the sub-reflector acts as a perfect reflector, ’ 

while for frequencies in another distinct band the sub-reflector is 

transparent to the radiation, thus permitting direct illumination of the 

main reflector by the feed. To date, periodic arrays of conducting 

plates, or apertures in a conductive screen, have been used as the frequency selective surface [l-41. Typically, the conductors are placed on 

a dielectric layer that provides mechanical support. 

Use of a dielectric layer with periodically varying dielectric 

constant has been suggested as an alternative way to obtain a frequency 

selective surface. As considered here, the layer is composed of alternating strips of two materials having different dielectric constants, as 

shown in Figure 1. At mm frequencies, such dielectric layers offer the 

advantage of low absorption loss as compared to metallic screens. Since 

the layer thickness is on the order of a wavelength, the amount of material required would not be excessive at these high frequencies. 

This report describes a theoretical study of frequency selective 

reflection and transmission at dielectric layers of the type shown in 

Figure 1. Because this effort was intended as a limited feasibility 

study, we consider the case when the plane is incident perpendicular to 

the strips, and assume the electric field to be polarized along the 

strips, as in Figure 1. The layer is found to exhibit the desired frequency selective properties. It is also found that the approximate 

1 
cu 

w 

3r 

cd 

'k 

k 

d 

0 

.d 

a 

0 

.d 

k 

Q, 

a 

d 

4-r 

0 

a 

a, 

o 

0 

a 

0 

Q, 

0 

cd 

w 

k 

8 

1 

(Il. 

rn 

2 
locations of reflection and transmission bands can be predicted from wave 

properties of a layer having uniform dielectric constant equal to the 

average of that in the periodic layer. 

The significant wave properties are discussed qualitatively in 

Section 11. In Section 111, the necessary mathematical analysis is 

carried out to permit numerical evaluation of the reflection and transmission coefficients. Numerical results are presented in Section IV. 

3 
11. Wave Mechanisms for Freuuency Selective Behavior 

In this section, we describe the wave phenomena that can provide 

frequency selective reflection and transmission at a periodic dielectric 

layer. The description given here is intended to clarify the nature of 

the subsequent analysis, and to provide a context for discussing the 

numerical results that have been obtained. 

Consider first a dielectric that is periodic along x but infinite 

along z, as shown in Figure 2. For two dimensional propagation in the 

(x,z) plane, the dielectric will support an infinite set of modes with 

different wavenumbers, rrn along z, but each having the same Bloch wavenumber % along x [5,6]. At low frequencies, only the lowest n = 0 mode 

will have a real wavenumber xO, while all other modes will be cut off 

imaginary or complex). At somewhat higher frequencies, the n = -1 mode 

will also propagate (K-~ real), while higher modes remain cut off. 

Further increase in frequency will result in more propagating modes. 

(K~ 

Consider now a semi-infinite, periodic dielectric illuminated by a 

plan wave incident from vacuum, as shown in Figure 3. The incident wave 

will excite all of the modes of the periodic structure. At a low enough 

frequency fl, only the n = 0 mode will propagate along z, as suggested in 

Figure 3a. Higher modes will decay exponentially away from the surface z 

= 0. 

propagate. If el and e2 are not close to the dielectric constant of free 

space, two modes can propagate in the dielectric even for the periodicity 

However, at a higher frequency f2, the n = 0 and n = -1 modes can 

.d small enough compared to the free space wavelength Ro so that no grating lobes are present in the field reflected into the region z < 0. In 

this case, the reflected field propagates only at the specular angle, as 

4 
a, 

Ll 

a 

.d 

cv 

a, 

k 

3 

M 

*d 

Erc 

5 
d 

I 

0 

c 

H 

P 

W 

C, 

eo Q, 

Q, II c 

a EC, 

d c 

c, oc d 

-d E 0 

0 cd 00 

Xdd II 

w am c 

6 
indicated in Figure 3b. 

ness h, as in Figure 4, the propagating modes will excite a plane wave in 

the vacuum region z > h below the dielectric. 

If the periodic dielectric is of finite thickAt low frequencies fl, only one mode propagates along z with real 

wavenumber no, so that the layer acts approximately as if it had a uniform dielectric constant equal to the average of that for the periodic 

layer. Thus the transmission properties will be similar to those of a 

uniform layer. In particular, the reflection coefficient will vanish at 

about the frequency for which Koh = 7~. In this case, the dielectric 

acts as a half-wave window, and there will be total transmission of the 

incident plane wave, as suggested in Figure 4a. 

At a higher frequency f2, both the n = 0 and n = -1 modes will 

propagate along 2. 'These modes are excited at the top surface of the 

layer by the incident wave. When each mode reaches the bottom surface, 

it excites both modes traveling back to the top, as well as a transmitted 

plane wave in the air. Because of the phase-matching conditions at the 

top and bottom surfaces, the transmitted plane wave in the air propagates 

in the same direction as the incident plane wave. 

The modes in the layer that are traveling back towards the top 

surface excite a reflected plane wave in the air above the layer, as well 

as being scattered back into the layer, as suggested in Figure 4b. 

Repetition of this scattering process establishes the total field in the 

layer, and the total reflected and transmitted plane waves in the air. 

For some frequency f2, the phases of the two modes in the layer will be 

such as to add destructively in producing the transmitted plane wave, 

while constructively adding for the reflected plane wave, thereby produc7 
m 

n 

P 

v 

cd 

k 

b 

d 

cd 

c, / 0 

B 

n 

cd 

v 

8 
ing the desired frequency selective property. However, at other frequencies in the range where the n = 0 and n = -1 modes propagate along z, the 

phase's of these two modes may be such as to add for the transmitted plane 

wave and cancel for the reflected wave. Thus it is possible to have 

multiple frequencies for which total reflection and total transmission 

take place. 

The foregoing behavior is known in other related geometries to be 

associated with the excitation of waves guided along the layer [7,8]. To 

understand the connection with the guided waves, consider a layer of 

uniform dielectric constant equal to an average permittivity of the 

periodic layer defined by 

This layer will support guided waves whose fields vary sinusoidally in 

the layer, and decay away from the layer in the air [9]. For TE guided 

wave modes, the normalized wavenumber Pgh with g=0,1,2 is plotted along 

the horizontal axis in Figure 5 versus the normalized free space wavenumber koh = oh/c, which is plotted along the vertical axis for ea = 2. 

Because ,9, is greater than co, there are no angles of incidence 8 at 

which a plane wave can satisfy the phase-match condition kosin8 = pg for 

direct excitation of the waveguide modes. However, excitation is possible if the dielectric constant of the layer is a periodic function of 

x. In this case, the fields of each waveguide mode will consist of a 

series of space harmonics, one of which has fields that are very similar 

to those of a mode in a uniform layer having the same average dielectric 

constant ea. This space harmonic is designated m = 0, and has wavenumber 

9 
0 

rl 

10 
fig, along x that is close to B for the uniform layer. Other space harg 

monics have wavenumbers along x given by 

While Bgo = fig > ko, it is possible to choose the periodicity d such 

that IOg,-1 I < ko for m = -1. In this case, a plane wave incident at an 

angle 0 = sin'( I fig,-l I /ko) will couple to the space harmonic, and 

through it excite the waveguide mode. Once excited, this mode will reradiate plane waves into the air regions above and below the layer 

through the same space harmonic. The process of excitation and re-radiation is depicted in Figure 6a for the case when 

This condition is sufficient to guar'antee that only one space harmonic 

will give rise to a plane wave propagating away from the layer, and also 

implies that the plane wave propagates backward with respect to the 

direction of the waveguide mode. Guided waves that radiate some of their 

energy as they propagate are known as leaky waves[lO]. 

leo + fig > 2n/d > 8,. 

The same physical processes hold if the incident wave is from the 

left, as shown in Figure 6b, except that the direction of propagation 

along x is reversed for the waveguide mode. The re-radiated plane wave 

above the layer adds to the reflected plane wave generated directly at 

the top surface of the layer to give the total reflected field. When the 

two components are in phase, strong reflections take place. However, 

when they are out of phase the reflected field is small and strong 

transmission occurs. Because the phases are frequency dependent, the 

overall reflection can have the desired frequency selective behavior. In 

11 
Re-Radiated 

Reflected Wav Plane Wave 

Directly Reflected 

X 

Re-Radiat ed 

Transmitted Wave 

Incident \ 

\ 

Guided Wave 

Total Reflected 

Wave 

7 

Plane Wave 

Guided ' Wave \,\- 

Transmitted 

Wave 

(b) 

Figure 6. 

X 

Plane wave excitation and re-radiation of the wave 

guided by a dielectric layer with periodically 

modulated dielectric constant for the case: a) plane 

'wave incident from the right; and b) plane wave 

incident from the left. 

12 
Section IV, it is shown that total reflection can be achieved, and that 

- the frequencies of total reflection can be predicted from the properties 

of the leaky wave modes. 

13 
111. Formulation of the Reflection Problem 

In this section we develop the mathematical formalism for the reflection and transmission coefficients at a layer of a periodic dielectric in a way that can be implemented on a computer for numerical evaluation. 

= 0) but it can be readily extended to the TM polarization. Expressions 

The analysis is restricted to the case of TE polarization (Ex = E, 

for the fields in the infinite periodic medium of Figure 2 are first 

derived. Subsequently the scattering matrix for a single surface normal 

to z, as shown in Figure 3, is found. Finally, network concepts are used 

to join th'e scattering matrices for the two surfaces normal to z that are 

shown in Figure 4, so as to obtain the reflection and transmission coefficients for the layer. 

A. Fields in an Infinite Perlodic Medium 

In studying the fields in the infinite periodic medium we use the 

approach of Collin [ll] and Lewis and Hessel 1121 which expresses th'e 

fields as a superposition of modes having different wavenumbers along z. 

Assuming a time dependence exp(-iot), the fields of the n-th mode have 

dependence along z given by exp(+ixnz). Thus the only non-zero component 

of electric field Ey and the z-component of magnetic field H, can be 

written as the sum of modal fields in the form 

Here en(x) and hn(x) are the x-dependent mode functions while An and B, 

are the modal amplitudes for waves propagating in the +z and -z direc14 
tions, respectively. 

Because each slab of dielectric has a homogeneous E, the mode func- 

. tions can be expressed in trigonometric form. We first define the wavenumbers along z in the two dielectrics as 

where k12 = ko2el and kZ2 = ko2e2. 

admittances for the dielectrics by 

We further define the impedances and 

Referring to the coordinate system in Figure 2, the mode functions 

for -dl, < z < 0 are given by 

e,(x) = V, cos(sx) + i ZlnIn sin(u,x), 

hn(x) = I, cos(unx) + i YlnVn sin(ii,x). 

In the range 0 < z < d2, the mode functions are 

en(x) = V, COS(V,X) + i Z2,1n sin(v,x), 

hn(X) = I, COS(V,X) + i YZnVn sin(v,x). 

The constants V, and I, in (6) and (7) are defined using the Floquet 

condition discussed below. 

The Floquet condition requires that the fields at one end of a 

period (z = -dl) differ from those at the other end (z = d2) by at most a 

phase factor. Writing the phase factor as exp(iSod), the Floquet condi15 
tion is 

en(d2) = en(-dl) exP(iSod), 

(8) 

hn(d2) = hn(-dl) exp(iSod). 

Substituting from (5), (6) and (7) into (8) gives two homogeneous equations in the two unknowns Vn and In. . In order to have a non-trivial 

solution of these equations, it is necessary for un and vn to satisfy the 

secular equation 

COS Sod = COS(Undl) cos(vnd2) 

Since and vn are functions of K, (9) can be viewed as a relation 

between K~ ar,d So. 

angle of incidence in Figure 1. Thus (9) serves as an equation whose 

roots are the allowed values of 'c~, and represents the dispersion equation of the periodic medium. For large Inl, the roots are well approximated by 

As will be seen later, So = ko sine where 8 is the 

Kn 5 [ko 2 ea - (S0+n2n/d) 2% I (10) 

which holds even for relatively small values of Inl. From expression 

(lo), it is seen that only for small values of In1 will K~ be real, 

whereas higher-order solutions will be below cutoff along z. 

When the dispersion equation (9) is satisfied, the ratio In/Vn can 

be determined form (8). After some manipulation, it is found that 

whereas In can be found from (11) if Vn is known. The value of Vn is 

16 
itself arbitrary and is usually obtained by normalizing the mode functions. In this analysis, we use the normalization 

d2 

-d1 

f le(x)I2 dx = d. 

This normalization is carried out numerically during the computations, as 

described subsequently. 

Because of the Floquet condition, the mode functions en(x) and hn(x) 

are periodic functions of x multiplied by the phase factor exp(iS,x). 

Thus we may write en(x) as the Fourier sum 

where 

Alternatively, the expansion coefficients anq can be found from the 

integral 

1 d2 f en(x) exp(-iSqx) dx. 

=nq = ;i 

-d1 

Substituting (6) and (7) into (15), we obtain after m 

tion that 

(1) (2) 

a nq = Vn [Jnq - Jnq It 

where 

ch manipi 

(15) 

lawith 

17 
iunvn [cos(vndz) - exp(iSgd)cos(undl)] 

unsin(vndp) + vnexp(iSod)sin(undl) 

Yn = 

(2) (1) 

The quantity Jnq in (16) is of the same form as that of Jnq with un 

replaced by vn, and dl replaced by -d2 in (17). 

The normalization (12) is equivalent to requiring 

00 

lanqi2 = 1. 

q=-OD 

From (16) it is seen that the normalization condition (19) implies 

Substituting expression (13) into (3) for Ey(x,z) and changing the 

order of summation gives 

00 OD 

When applying boundary conditions at the surface z=O in Figure 2, it 

is necessary to consider the x component of magnetic intensity Hx(x,z). 

For the TE polarization, and using the Maxwell curl equations, it is 

readily seen that Hx can be found from the derivative with respect to z 

of Ey. The resulting expression is 

18 
Expressions (21) and (22) are used in the next section to find the scattering matrix for a single interface z = 0. 

B. Scattering at a Sinale Interface 

A single interface at z-0 is depicted in Figure 3. The field in the 

air region z < 0 consists of an incident plane wave propagating at an 

angle 8 with respect to the z axis, and reflected plane waves corresponding to the specular and higher space harmonics. We define 

So = ko sin e, 

C, = ko cos 8, 

and assume the incident electric field to be polarized along y with 

amplitude Eo. The electric field of the incident wave is then given by 

E, exp (isox) exp(iCoz). (24) 

The electric field due to the reflected wave is the sum of the space 

harmonics and takes the form 

where Sq is given by (IC) and 

Cq = (ko2 - Sq2)JI. 

For the conditions of interest here Sq2 2 kO2 for all q # 0 so that only 

the fields of the specular (q = 0) space harmonic propagate away from the 

interface, while the fields of all other space harmonics decay. 

amplitude coefficients Rq have yet to be determined. 

The 

From (24) and (25) the total electric field in the air region is 

19 
seen to be 

where dqo is the Kronecker delta. The x component of magnetic intensity 

in the air can be found from the derivative with respect to z of (27), 

which yields 

The boundary conditions at z=O require that Ey and Hx be continuous 

there. Equating the pair (21), (27) and the pair (22), (28), and making 

use of the orthogonality of the functions exp (iSqx) over a period, one 

obtains 

Qo 

In (29) and (30), n ranges over all positive and negative integers, so 

that these equations represent two infinite sets of equations. 

We wish to solve (29), (30) for the amplitudes Rq and An of the 

waves traveling away from the surface in terms of the amplitudes Eo and 

Bn of the incident waves. In this way, we obtain the scattering 'matrix 

of the surface. The waves incident from the periodic medium arise from 

reflection at the second surface, as suggested in Figure 4. For the 

20 
conditions of interest, only one or two waves are propagating in the 

periodic medium, while all other waves are cutoff along z. It can be 

shown that the two propagating waves correspond to the indices n = -1, 0. 

Since the two surfaces are separated by at least one half-wavelength, the 

fields of the cutoff waves excited at one surface are exponentially small 

at the other surface. 

0 for n # -1, 0. 

Hence, to a good approximation we may assume Bn = 

While all the higher modes are excited in the air and in the periodic medium, the interaction at the two surfaces and radiation into the 

air are described by the amplitudes Ro, A,, A,1 of the propagating waves. 

Thus we ultimately need only the 3x3 portion of the full scattering 

matrix relating Ro, Ao, A,1 to E,, Bo, B,l. This scattering relation 

takes the form 

To solve for the elements Sa ,9 in the scattering relation (31), we 

first multiply (29) by Cq and then add (29) and (30). Since Bn = 0 for 

n # -l,O, the resulting equation may be written in the form 

Expression (32) represents an infinite set of equations with index q in 

an infinite number of unknowns A,. 

The terms (xn + Cq) aqn can be viewed as elements of a matrix. With 

21 
this view, it can be shown that the elements decrease as one moves away 

from the main diagonal. Thus it is reasonable to solve (32) by using a 

finite truncation of the summation and a corresponding limitation on the 

. number of values of q considered. For the parameters chosen in this 

study, it was found to be sufficient to allow q, n to range over the 

integers -3, -2, -1, 0, + 1, + 2. With this truncation, (32) is solved 

for A, with n- -3 ... +2. Returning to (29), we then compute Ro from 

2 

Collecting the resulting values of Ao, A,1 and Ro due separately to Eo, 

Bo and B,l gives the values of the scattering matrix in (31). 

C. Scatterinu From A Periodic Laver 

Having found the scattering matrix for a single surface (z = 0), 

network concepts can be used to treat the interation with a second surface at z = h. As discussed previously, the interations between the two 

surfaces are essentially due to the propagating modes in the periodic 

medium. For our case, these are the n = -1, 0 modes, which are shown in 

the transmission line model for the interation shown in Figure 7. At the 

surface z = 0 the incident waves Eo, Bo and B,l couple to the scattered 

waves ROD A. and A,1. 

At the surface z = h, the incident waves in the layer are 

A, exp(irroh) and A,1 exp(iz,lh), while no wave is incident from the air 

side. 

B,l exp(-iK,lh). The relation between scattered and incident waves is 

again given by (31), which for the foregoing conditions takes the form 

In this case, the scattered waves are To, Bo exp(-ix0h) and 

22 
?7 N N 

2 ,+ 

43 4 

0 

N 

3 

QO w 

NO 

Y 

.rl 

- No 

Y 

Nd .rl 

Qb 

% 

CplL 

I 

Y 

d 4 

I 

k 

: 

I 

A 

v 

0 

a 

5 

23 
1. 

2. 

3. 

4. 

5. 

(34) 

[i:exp(-ixoh) ] = [ :::: s1,o So,o ::I:] 

B-lexp(-iK,lh) S-1,l s-1,o S-L-1 

For the problem of scattering by a layer, the field Eo is known and 

one wishes to solve for the reflected and transmitted wave amplitudes Ro 

and To, respectively. To this end, (31) and (34) can be viewed as six 

inhomogeneous equations in six unknowns Rot Ao, A-1, Bo, B-l, To. Assuming Eo=l, the solution of these equations for Ro and To give the reflection and transmission coefficients of the layer. This approach has been 

used as the final stage in our computer program, as discussed below. 

D. Computer Program for Ro and To 

A program has been written in the PL-1 language to compute Ro and To 

by the methods derived above. The listing of the program is given in 

Appendix A. An outline of the program is given below. 

Given input frequency o, angle of incidence 8, geometric parameters 

of the layer dl, dZ, h and the electrical parameters el and e2. 

For integers n between -3 and 2, compute K~ from (9) using Newton's 

method with starting value given by (10). 

For integers n, q between -3 and 2 compute A using (16)-(18) and 

(20). 

Using (32) and (33), solve for An (-3 n 52) and Ro for Eo = 1, Bo 

= B-l = 0. This gives the elements Sl,l, So,1, S-l,l in (31). 

Repeat for Bo = 1 with Eo = B-l = 0, and then for B-l = 1 with Eo = 

Bo = 0 to get the remaining scattering coefficients in (31). 

Using (31) and (34) with Eo = 1, solve for Ro and To. 

nq 

24 
Several checks were carried out to ensure that the program was 

working properly. Choosing c1 and e2 very close to each other, we computed Ro and To. As expected, they were very close to the values for a 

homogeneous dielectric layer of value ea. For values of el and e2 used 

subsequently, it was found that the scattering matrix in (31) conserves 

power. Finally it was observed that lROl 

as required by power conservation. 

+I TOI * is very close to unity, 
IV. Numerical Studies of Frequency Selective Reflection 

Using the computer program described previously, we have carried out 

numerical studies for several examples. The purpose of these studies is 

to gain insight into the frequency selective behavior that can be expected for Ro and To, and to relate this behavior to wave processes in 

the periodic layer. We have therefore arbitrarily chosen the angle of 

incidence 8 = 45O and set dl = dp = d/2. 

For the initial studies, we have assumed el = 2.56 and e2 = 1.44, 

which are realistic values for low-loss plastics. With these choices, 

the average dielectric constant is ea = 2. The perodicity d must now be 

chosen such that, at the high frequency of interest, two modes (with n = 

-1, 0) propagate-in the dielectric layers. Furthermore, only the specular (q = 0) space harmonic must propagate in the air. 

A. Choice of d, h and frequency 

The restriction on d needed to insure that only the q = 0 space harmonic propagates in the air can easily be interpreted with the help of 

Figure 8a. The incident wave has wavenumber So = ko sin 0 < ko along x. 

Wavenumbers Sq of other space harmonics lie at a distance q(2x/d) away 

from So, as shown for q = -2, -1 and + 1 in Figure 8a. Provided that d 

is small enough so that 

ko sin 0 - 2n/d < - ko, (35) 

the q = -1 space harmonic will lie outside the visible circle defined by 

S2 + C2 = ko2. 

harmonics will also lie outside the circle, so that Cq defined by (26) is 

imaginary for q # 0, and the space harmonics decay away from the layer. 

Provided that the modulation (e1 - e2)/ea is not too large, the concept 

It is further seen that, if (35) holds, all other space 

26 
S =k sin9 

00 

S =k sine 

00 

Figure 8. Vissible circle for determining the propagating space harmonics 

in: a) air; and b) a medium with dielectric constant 

'.- 

27 
of the visible circle can also be used to estimate the number of propagating waves in the layer. When the modulation is small, (10) can be 

used as an estimate for Kn, in which case the visible circle is given by 

Sn2 + irn2 = ko2ea, as shown in Figure 8b. 

that, for the two modes n= -1, 0 to propagate in the layer, d must satisfy 

It is seen from this figure 

ko sin 8 - 4z/d < -ko& < ko sin 8 - 2z/d, 

koFa < ko sin 8 + 27c/d. 

Conditions (35) and (36) can be rearranged into the following inequalities: 

d/A < 1/(1 + sin e), 

l/(G + sin 8) < d/A < 2/(& + sin e), (37) 

d/A < 1/(& - sin 8). 

Assuming ea = 2 and 8 = 45O, these inequalities are d/A < 0.586, 

0.471 < d/A < 0.943 and d/A < 1.414, respectively. To satisfy these inequalities, we have chosen d/A = 0.54. 

Initially a value of layer thickness h at which total reflection 

will occur was obtained by computing Ro for various values of h/R. 

Subsequently, it was found that values of h and d for total reflection 

could be related via the conditions for guidance of a wave by a layer of 

uniform dielectric constant “a. Whereas our initial approach gave us the 

value of h/R = 0.925 for sample calculations, it is the subsequent interpretation that is discussed below. 
B. Variation of Ro With Frequency 

Computations of the frequency dependence of Ro have been made assuming that h = 0.925 and d = 0.54 for el = 2.56 and e2 = 1.44. This choice 

produces total reflection for a frequency f such that R = c/f is about 

unity. Note that, if h and d are scaled by A, then total reflection can 

be obtained at any desired frequency. The results of the calculation for 

lROl are depicted in Figure 9, where we have used the normalized frequency variable k,h = 2xfh/c, and have plotted up to the value koh = 6.30 at 

which the q = -1 space harmonic in air switches from cutoff to propagating along 2. 

For koh < 5.12, only the n = 0 mode in the periodic dielectric is 

propagating along z. In the frequency range 0 < koh < 5.12 for single 

mode propagation, Ro vanishes at the two frequencies koh = 2.56 and 5.03, 

at which KOh = 3.145 and 6.355. Thus, frequencies of total transmission 

occur when the layer thickness is close to a multiple of one half the 

effective wavelength along 2, as predicted in Section 11. The difference 

between the values of Koh for total transmission and n, 2n are due to the 

non-zero phase of the transmission and reflection coefficients at the 

individual surfaces z = 0 and h, which result from excitation of higher 

cutoff space harmonics. 

For k,h < 5.12, where the n = -1 mode also propagates along z in the 

layer, total reflection takes place at two frequencies given by koh = 

5.32 and 5.83. In the vicinity of these frequencies for total reflection, the variation of lR~l is that associated with resonances wherein a 

frequency dependent function has a real axis zero and a nearby pole at a 

complex location. 

29 
0 

4 

00 

0 

to 

d 

dc 

0 

cv 0 

0 

to 

u3 

cil 

W 

30 
To examine the bandwidth of the total reflection resonance, we have 

plotted I RO I on an expanded scale in Figure 10. Since IRo I + I To I = 

1, the curve can also be used to determine IT0I2 by using the vertical 

scale to the right of the plot. 

centered at koh = 5.83. For this peak, the fractional bandwidth between 

The wider of the two peaks of IR0I2 is 

the frequencies at which IR0l2 = 0.9 is 0.86%. 

The narrower of the two peaks in Figure 10 is shown further expanded 

in the insert. For this peak, the fractional bandwidth between the frequencies at which IR0I2 = 0.9 is less than 0.04%. 

wider bandwidths are found for total transmission when the n = -1 mode in 

the layer is cut off. 

Figure 9, a 40% bandwidth is found between the frequencies for which 

IR0l2 =I 0.1. 

C. Prediction by Means of Gulded Waves 

By comparison, much 

For example, in a region about kod = 2.56 in 

The location of the frequencies of total reflection can be predicted 

from the properties of the waves guided by the periodic layer. Consider 

first the case of a wave guided along a uniform layer having dielectric 

constant ea. The normalized propagation constant Bgh of this guided wave 

is plotted horizontally in Figure 5, versus the normalized frequency koh, 

which is plotted vertically. The lowest g = 0 guided wave mode starts at 

the origin and becomes asymptotic to the wavenumber koTa of the layer. 

Higher guided-wave modes start at points koh = gn(ea-l) along the 45O 

line, where g = 1, 2, 3 .... 

In Figure 11, we have repeated the plot of Figure 5, and have added 

the dispersion curves for waves propagating in the negative x direction 

(6 < 0). We have also plotted as a broken line the transverse wavenumber 

31 
0 

- d d 

. 

In 

co 

d 

00 0 

0 d 

4 

m 

In 

m 

m 

In 

d 

d 

m 

0 00 

0 

d 

00 a 

cv d - 0 

- Eo 

d( 

0 

32 

cv 

0 

0 

0 

cd 

k xu 8% 

i!k 

c,d 

'5 dg 0 

3c 
'\.\ \' a 

1. n 

N 

I I I I I I I 1 8 I I I 1 I I I 

I I 

1 0 

0 

d 

33 
Soh = koh sin 6 of the incident plane wave. 

broken line is never close to the dispersion curves for the guided waves. 

Except at koh = 0, the 

As a consequence, an incident plane wave cannot couple to the guided 

waves on a uniform layer. 

If the dielectric layer is made periodic along x, the field of each 

guided wave mode becomes a sum of space harmonics. For the guided waves 

traveling in the -x direction, the wavenumber of the q = -1 space harmonic is (-fig + 2n/d). 

this space harmonic has the same form as -Bgh versus koh, except for a 

shift 2zh/d to the right. While finite modulation affects the value of 

When normalized by h, the dispersion curve of 

19, for the guided wave, for small modulation of the dielectric constant 

8, is close to that for a uniform layer. 

Dispersion curves for the q = -1 space harmonics of the guided waves 

in the small modulation limit are shown in Figure 12 for h/d = 

0.925/0.54. 

koh. 

strong coupling between an incident plane wave and guided waves through 

the q = -1 space harmonic. Note that, above the dashed line having an 

angle of -45O, the q = -1 space harmonic in the air propagates along z. 

Thus for reflection and transmission of a single space harmonic, the 

operating point along the Soh line must be kept below the dashed line. 

For the parameters used in drawing Figure 12, this condition implies that 

koh F 6.30 for one propagating space harmonic in air. 

We have also drawn a broken line representing Soh versus 

Intersection of the Soh line with the dispersion curvers indicates 

In Figure 12, the line Soh intersects the dispersion curve for the g 

= 0 guided wave at koh = 5.27, and for the g = 1 guided wave at koh = 

5.72. These values are close to the values koh = 5.32 and 5.83 for total 

34 
0 

rl 

I I I I-- - 

c 

cn 

k 

0 

c 

M 

0 

ho 

00 (D 

c 

24 

0 

d' 

35 

cv 

0 
reflection (IROI=l) obtained from Figure 11. The deviation between the 

values of koh obtained from Figures 11 and 12 is thought to result from 

the fact that the finite modulation of the dielectric constant of the 

layer alters Bg from the value obtained for a uniform slab. Thus, decreasing modulation should bring the values closer together, while increasing modulation should result in greater deviation. This latter 

condition is shown subsequently. 

To further demonstrate the relation between the frequencies of total 

reflection and the guided waves of the layer, we have considered a layer 

of increased thickness h = 1.1, but the same periodicity d = 0.54. The 

dispersion curves of the q = -1 space harmonics of the first three guided-wave modes are shown in Figure 13 for the limiting case of small 

modulation. The broken line giving Soh = koh sin 8 is seen to intersect 

the three dispersion curves at koh = 6.22, 6.67 and 7.30. Our model predicts that total reflection should take place at normalized frequencies 

close to these values. 

A plot of lROl versus koh for h = 1.1 and d = 0.54 is shown in 

Figure 14. From this plot, total reflection is seen to occur at the 

three frequencies koh = 6.25, 6.78 and 7.42, which are close to those 

predicted by the small modulation theory. The bandwidth over which IR0I2 

> 0.9 about each frequency of total reflection is seen to increase as koh 

approaches the value 7.50 where the line Soh crosses the dashed line, 

above which the n = -1 space harmonic propagates in air. The bandwidth 

about the lowest of these frequencies is only 0.01%, while that of the 

highest is 0.7%. 

36 
I I I I I I I I I 

D 

c 

M 

k 

0 

c 

ln 

0 

M 

ldc 

d 

00 

cu 

d 

cil 

d 

0 

d 

00 

i9 

dc 

cu 

0 

d dc cil 0 

37 

c oal 

a 

mc, 

k 

ale 

-4 m 

aal 

3 

2 

xg 

c, 

a 

aa, 

ea 
I I I I 

I I I I I I I 

/ 

I I 

a3 (0 d cv 

d d d d d 

38 

0 
C. Influence of Modulation 

To explore the influence of modulation, we have computed the reflection and transmission coefficients for a layer with el = 3 and e2 = 1. 

This layer has average dielectric constant ea = 2, as before. 

assume that d = 0.54 and h =0.925, as in the case of the results presented in Figures 10-12. A plot of IRoI versus normalized frequency koh 

is shown in Figure 15. The variation of lROl is seen to be qualitatively 

the same as that of Figure 10. The increased modulation is seen to shift 

the first frequency of total reflection (IROI = 1) to koh = 5.45 and the 

second to koh = 6.12, which are farther from the respective values 5.27 

and 5.72 predicted by small modulation theory. 

We further 

Besides shifting the frequency of total reflection, the modulation 

influences the bandwidth. At the first total-reflection frequency, the 

bandwidth for I RoI 0.9 is 0.001%. However, at the higher total reflection frequency, the bandwidth is 2.5%. The modulation is also seen 

to have a small effect on the frequencies of total transmission (IROI=O). 

39 
0 00 

rl 0 

CD b( 

0 0 

cv 0 

0 

CD 

c 

X 

0 

In 

dc 

i?Y 

cv 

II 

w 

cv 

k 

0 

w 

c 

X 

0 

WK 

cdc 

40 
V. Conclusion 

It has been shown that frequency selective reflection and transmission takes place at a periodically modulated dielectric layer. Frequencies of total transmission and total reflection were found, and they can 

be related to'various wave phenomena. In the limit of small modulation, 

these frequencies can be estimated from the appropriate wave phenomenon 

in a uniform layer having dielectric constant equal to the average of 

that in the periodic layer. 

In the range of low frequencies where a single space harmonic propagates along z in the periodic dielectric, total transmission occurs when 

the layer thickness h is one half the effective wavelength along z, i.e., 

when h = n/rrg. For small modulation, K~ 

average dielectric constant. Total reflection can be achieved at those 

higher frequencies for which two space harmonics propagate along z in the 

periodic dielectric. These frequencies of total transmission are associated with the excitation of leaky waves guided by the dielectric layer. 

In the limit of small modulation, the frequency of total reflection can 

be approximated from the dispersion characteristics of waves guided by a 

uniform dielectric layer. 

( eako2-So2)k, where ca is the 

In the examples treated, the bandwidth over which IRo I 2 0.9 about 

the frequency of total reflection was found to be small. The largest 

bandwidth obtained was 2.5%. While angle sensitivity was not computed, 

the narrow frequency bandwidth suggests that, at the frequency of total 

reflection, IRo12 will be sensitive to the angle of incidence 8. 

Whereas the study was carried out only for the TE polarization, the 

form of the results have implications for the TM polarization. We expect 

41 
that the frequencies of total reflection for the TM polarization are also 

associated with the excitation of the leaky waves guided by the periodic 

layer. However, the dispersion characteristics of the leaky TM waves 

will differ from those of the TE polarization. As a result, it is expected that incident plane waves of the TE and TM polarizations will, in 

general, experience total reflection at different frequencies. Hence, 

the periodic dielectric layer is expected to be polarization sensitive. 

42 
References 

1. V.D. Agrawal and W.A. Imbriale, "Design of a Dichroic Cassegrain 

2. B.A. Munk, R.G. Kouyoumjian and L. Petris, "Reflection Properties of 

1 Subreflector," IEEE Trans., AP-27 (1979)8 pp. 466-473. 

Periodic Surfaces of Loaded Dipoles," IEEE Trans., AP-19 (i971), pp. 

612-616. 

3. E.L. Pelton and B.A. Munk, "Scattering from Periodic Arrays of 

Crossed Dipoles," IEEE Trans., AP-27 (1979), pp. 323-330. 

4. W.V.T. Rusch, "The Current State of the Reflector Antenna Art," IEEE 

Trans., AP-32 (1984), pp. 313-329. 

5. R.S. Chu and T. Tamir, "Guided-Wave Theory of Light Diffraction by 

Acoustic Microwaves,tt IEEE Trans., MTT-18 (1970), pp. 486-504. 

6. S.T. Peng, T. Tamir and H.L. Bertoni, "Theory of Periodic Dielectric 

7. A. Hessel and A.A. Oliner, "A new theory of Wood's Anomalies on 

Waveguides," IOEE Trans., MTT-23 (1975), pp. 123-133. 

Optical Gratings," Appl. Optics, 3 (1965), pp. 1275-1297. 

8. T. Tamir, "Nonspecular Phenomena in Beam Fields Reflected by Multilayered Media," J. Opt. SOC. Amer. A, 3 (1986), pp. 558-565. 

9. R.E. Collin, Fjeld Theory of Guided Wave (McGraw-Hill, 1960), Ch. 

11, p. 453. 

10. T. Tamir and F.Y. Kou, Varieties of Leaky Waves and Their Excitation," IEEE J. Quant. Electron, QE-22 (19861, pp. 544-551. 

11. R.E. Collin, "Reflection and Transmission at a Slotted Dielectric 

Interface," Can. J. Phys., 34 (1956), pp. 398-411. 

12. L.R. Lewis and'A. Hessel, "Propagation Characteristics of Periodic 

Arrays of Dielectric Slabs," IEEE Trans., MTT-19 (1971), pp. 276- 

286. 

43 
._. ... 

Appendix A: Listing of Computer Program 

1. 

2. 

4. 

4.1 

4.2 

5. 

5.2 

5. ,7 

5.4 

5.5 

5.6 

5.7 

6. 

7. 

0. . 

9. 

1 0. 

11. 

12. 

13. 

15. 

16. 

16. (31 

16. 02 

16.03 

16.1 

16.2 

16.3 

16.4 

17. 

18. 

19. 

20. 

21. 

22. 

23. 

24. 

25. 

26. 

27. 

29. 

30. 

31. 

32. 

33. 

34. 

35. 

3s. (31 

35.02 

35. 0.3 

35. 04 

35. (:I4 1 

.--a . (35 

35.06 

35.07 

35.08 

35.09 

- 

.:I . 

5. a 

28. 

-c 

//TEMNNF'OP JOE (F1(39602,CHEO! , 'L. S. CHEO' 

// EXEC PLIXCLG 

//PLlL.SYSIN DD * 

TEMODE: PROCEDURE OPTIONS(MA1N): 

/* FINAL FROGRAM FOR TE MODE WITH ONE AND TWO * / /* PROPAGATION MODES ic/ 

/* GAMA(-1) CAN BE EITHER REAL OR IMAGINARY */ /~*+**************+*~*****************~**~*****+*~~*~*

/* INSTRUCTION FOR DATA ENTRY : ie / 

/* THE FIRST DADA CARD ENTRY ORDER IS : */ 

/+ D,DL,D2,El,E2,LAMTA,DEGREE,H *i 

/* M=-6,-5,-4,-3,-2,-1,0,1,2,3,4,5 * i 

/*+****~*+*****+*******************~****~*~***ic ~** / 

/* THE SECOND DATA CARD IS FOR VALUES OF M AS : it/ 

/* CONSIDER V (-6: 5) ,U (-6: 5) , SN (-6: 5) AS ARRAYS */ /* XI = U/V, U IS REAL , V CAN BE REAL OR IMGGINGRY */ 

/* */ 

/* F(U) = 0 BY NEWTON'S METHOD , TO SIX DIGIT *i 

/* GCCURACY. +/ 

/* */ 

/* FIND A ROOT NEAR UO-SO OF THE FUNCTION ic /. 

/* THIS PROGRAM COMBINES TEMNNPO, RTLCOMP,AND RTLGGUS */ 

/'* AA IS THE COEFFICIENT MkTRIX FOR GAUSS EL-IMINATION */ /* AND X IS THE SOLUTION VECTOR OBTAINED */ /* FROM THE GAUSS ELIMINATION METHOD PROCEDURE */ 

DCL (AA(I:E,I: 11) ,x(i:s,i:a)) FLOAT ; 

/* */ /* MATRIX MN(M,N), EXCEPT PRINTING PART FOR ALPHA*/ 

/* SUMJREAL,SUMJIMAG, SUMJlJZ ARE ELIMINATED. '**/ 

/*+* ***/ 

DCL (S,SS,U, (Ul,U2) (-6:s) ,EPSILON) FLOAT DEC; 

DCL (V,SM) (-6:51 FLOAT , 

DCL (UDl(-6:5),RTRD2,D2,Dl,D,SQHTR) FLOAT ; 

DCL (EI,E2,KO,SOD,UR,RU,SUMl,SUM2,~(-6:~)) FLOAT j 

DCL (A,E,TERM,PI,SMD(-6:5),GAMAIN,GAMA(-6:5)) FLOAT ; 

DCL (Y,YR,YI,YI,YZ,Y3,Y4,YMAG) (-6:5,-6:5) FLOAT ; 

DCL (J1R1,J1R2,J1H3,JlR4,JlR) (-6:5,-6:5) FLOAT; 

DCL (J111,J112,J113,J114,J11) (-6:5,-6:5) FLOAT; 

DCL (Jl,J2) (-6:5,-6:5) FLOAT; 

DCL (J2Rl,J2R2,J2R3,J2R4,J2R) (-6:5,-6:9) FLOAT; 

DCL (J211,J212,J213,J214,J21) (-6:5,-6:5) FLOAT; 

DCL (SUMJlJ2,SUWJR,SUMJI) (-6:5,-6:5) FLOAT; 

DCL (~MN,AMNR,AMNI,MN,MNR,MNI) (-6:5,-6:5) FLOAT: 

DCL SUML (-6: 5) FLOAT ; 

DCL (ALPHA,CN,COSVD2,SINVD2) (-6:s) FLOAT; 

DCL (UMDl,SND,SNDl,SNDZ,VMDZ) (-6:s) FLOAT; 

DCL (COSHVD2,SINHVD2) (-6:5) FLOAT; 

DCL. (TR,TI) (3,-2: 1) FLOAT(6); 

DCL. (SCATH,SCATI) (-2: 1,-2: 1) FLOAT(6): 

DCL !ROR,ROI) (1:.3) FLOAT (6) ; 

DCL < LRMTA, DEGREE, THETA, KO2, RAT IO ) FLORT : 

DCLiN, 1 ,J,F) FIXED(5,O) : 

DCL (AB(12) ,CD(12)) FLOAT j 

DCL GAMAO-1H FLOAT ; /* GAMA (0) *GAMA (-1 ) +H */ 

DCL GAMAQ2H FLOAT; /+ 2+GAMA(O)*H */ 

DCL GAMA 12H FLOAT: /* 2*GAMA(-l)*H */ 

(K:,KMCiX) FIXED DEC; 

DCL M FIXED (4,Q)INITIAL(-5); 

I')CL (RTLR(4,4) ,RTLI (4,4)) FLOAT: 

44 
35.1 

3J. 11 

.36. 

36.1 

36.2 

37. 

37.1 

37.2 

37. .3 

37.31 

37.32 

37.4 

37.5 

37. b 

37.7 

37.8 

37.9 

38. 

39. 

.39.1 

40. 

41. 

42. 

43. 

44. 

45. 

46. 

47. 

47.1 

48. 

49. 

JCL. 

51. 

52. 

53. 

54. 

55. 

56. 

57. 

58. 

59. 

59.1 

60. 

61. 

62. 

63. 

64. 

65. 

66. 

67. 

60. 

69. 

70. 

71. 

72. 

73. 

74. 

75. 

76. 

77. 

'7 13 . 

79. 

80 . 81. 

82. 

-c 

c 

Q3. 

DCL GAMAZH FLOAT :"I* GRMA(!:))*H 

DCL GAMA-114 FLOAT: /* GAMfi(-l)+H */ 

DCL (LL) FIXED(S,(:)): 

DCL FLAG2 FIXED(2,O): 

DCL~RREAL,KIMAG,RFKIME~~~~RItl~~~ FLflAT: 

DCL (CDISC,GAMAD) (-6:s) FLOAT: 

DCL (!32R,B2I,B3R,P31) (-6:s) FLOAT: 

DCL H FLORT ; 

DCL(MM,NN,PIVOT,RS) FIXED !5,0) : 

*i 

DCL ( TOR , TO I , TM I NUS 1 R , TM I NUS 1 I ! 

DCL(RABS,RPHASE,RFRAtiS,RP"ASE) FLOAT: 

FLOAT : 

/* MM=NO. OF COLUMNS b "=NO. OF ROWS IN MN (N, M) .* / /* MATRIX IN PROGRAM TEMNNPO +/ 

/* IF PIVOT=l THEN THE PIVOT IS SUBSCRIBED, +/ 

/+ RS=NO. OF RHS OF THE AUGMENTED COEFFICIENT MATRIX+/ /* IN THE PROCEDURE GAUSSiELIMINATION METHOD * / 

/* */ /* DEFINE THE FUNCTION F(Xj 

F : PROCEDURE(U,M) j 

DCL U FLOAT ; 

DCL M FIXED(4,O) j 

CALL CHECER(U,M) : 

S=COS(UDl (MI ) *A-0. S*SUMl*SIN (UDL (Mj ) *E+COS (SOD) 

END F; 

/* */ /* DEFINE THE 1ST DERIVATIVE OF F(X) */ 

FPRIME: PROCEDURE(U,M): 

; 

RETURN (S) : 

DCL U FLUAT ; 

DCL M FIXED (4,O): ' 

CALL CHECKR (U , M) j 

DCL (TERMl,TERM2,TERM3 ) FLOAT; 

TERM1 =-Dl+SIN(UDl(M))*A-D2*COS(UDl(M))*B+IJH : 

TERM2=-O.S*SUM2*SIN(UDl(M))*~ ; 

TERM3~-O.~*SUMl*~D1*COS~UD1o)+E+TERM); 

SS-TERM1 + TERM2 + TERM3 ; 

END FPRIME; 

/* PROCEDURE TO CHECK THE SIGN OF R=U**2-K0++2+ (EI-EZ! */ */ 

RETURN (SS ) j 

/* 

CHECKR : PROC(U,M); 

DCL U FLOAT ; 

DCL M FIXED(4,O); 

UDl(M) = U * D1 ; 

R(M) =U+*2-K!3*+2*(El-E2) ; 

SQRTR = SQRT(ABS(R(M))) ; 

HTRD2 = SQRTR*D2 ; 

RU = SQRTR/U ; 

UR = U/SG!RTR ; 

SOD = SO * D ; 

IF R(M) ). 0 THEN 

DO ; 

A = COS (RTRD2) j /* 

Ec = SIN(RTRD2) j 

SUMl = RU + UH : 

SUM2 = 2/SQRTH-RIJ/U-lJR**2/~~~T~ ;

TERM = DZ+SIN(UDl (M) )*A : 

END ; /* El'JD OF R:::,i) CASE */ 

ELSE 

DO : / .+ Wt-IEN ii'.::: */ 

A = COSH(RTRD2) : 

E = SINH(RTRD2) : 

SUMl = UR - RU ; 

SUM2 = 2/SQRTR+RU/U + UR**2/SG!RTR ; 

TERM = -D2*SIN(UDl (M) )*A : 

END : I.* END OF R4:O CASE */ 

*/ 

45 
84. 

85. 

86. 

87. 

88. 

89. 

90. 

90.1 

90.2 

90.21 

90.3 

90.301 

90.31 

90.31 1 

90.32 

90.33 

90 I 34 

90.35 

90.5 

90.6 

90.7 

90.71 

90.9 

90.91 

90.91 1 

90.92 

90. 93 

91. 

92. 

93. 

94. 

95. 

96. 

96.1 

97. 

98. 

99. 

100. 

101. 

102. 

103. 

104. 

105. 

106. 

107. 

108. 

109. 

110. 

111. 

112. 

113. 

114. 

115. 

116. 

117. 

118. 

119. 

133. 

121. 

122. 

123. 

124. 

125. 

126. 

127. 

)I: 

END CHECKR j /** END OF CHECKR PROCEDURE **./ 

/+ u/ 

FLAG=l; 

ON ENDF I LE (SYS I N) FLAG=(); 

/* MAIN PROCEDURE + i 

PI = 3.14156 ; EPSILON=l.OE-06; EMAX=%) : 

GET LIST(D,Dl,D2,El,E2,LAMTA,DEGREE,H): 

K.0 = 2+PI/LAMTA : 

PUT SKIP: 

PUT SKIP EDITI'TE MODE WITH TWO PROPAGATION MODES') 

IF LAMTA=l. C) THEN 

(X (5) ,A) ; 

ELSE 

DO : 

PUT SKIP ; 

PUT SKIP EDIT('TE MODE WITH ONE PROPAGATION MODE') 

END : 

(X (5) ,A) : 

PUT SKIP EDIT (REPEAT( '*',40> ) (X (5) ,A) ; 

PUT SKIP EDIT('D=',D,'Dl=',Dl,'D2=',D2,'THETA=',DEGREEj 

(X(5),4 (A,F(9,.3),X(4))): 

PUT SKIP EDIT('El=',El, 'E2=',E2, 'LAMTA=',LAMTA, 'b::Q=' 7t:::Q

(X(5) ,4 (A,F(8,5) ,XC3))) j 

PUT SKIP; 

PUT SKIP EDIT('H=',H) (X(5),A,F(8,4)): 

PUT SKIP EDIT(REPEAT('+',SO)) (X(5) ,A): 

PUT SKIP j 

CALL INPUT(M) ; 

LOOP: DO WHILE (FLAG=l) : 

CALL CALCULATE (M) ; 

CALL PRINT (M) : 

FLAG2= 1 ; 

CHECKl: DO WHILE (ABS ( (U2 (M)-Ul (M) ) /U1 (M) ) >.=EPSILON 8< 

FLAG2= 1 ) ; 

IF K<=KMAX THEN DO: 

U1 (M) =U2 (M) ; 

CALL CALCULATE (M) : 

CALL PRINT (M) ; 

END; 

ELSE DO; 

PUT SKIP(2) EDITt'FAILS TO CONVERGE')(A): 

FLAG2=0 ; 

END; 

END CHECK1; 

PUT SKIP(2); 

PUT SKIPj2) EDIT('R=',R(M)) (X(5) ,A,E(14,6) ); 

IF R(M) (1 THEN 

PUT SKIP(2) EDIT('V IS IMAGINARY') (X(5) ,A); 

V(M) = SQRTR ; 

PUT SKIP(2) EDIT('V(M)=',V(M)) (X(5) ,A,E(14,6) ); 

GAMAIN = K0**2+Ei-U2(M)++2 ; 

PUT SKIP(2) EDIT('KO+*2+El-U2+*2=',GAMAIN) (X(S),A,E(14,6 

IF GAMAIN *< 0 THEN 

PUT SEIP(2) EDIT('GAMA IS IMAGINARY') (X(5) ,A): 

GAMA (M) = SQRT (ABS (GAMAIN) ) ; 

PUT SKIP(2) EDIT('GAMA=',GAMA(M)) (X(5) ,A,E(l4,b)): 

IF FLAG2=1 THEN CALL OUTPUT(M) j 

CALL INFUT (M) : 

END LOOF: 

/+ */ /* TO COMPUTE UM.*D 1 , VM+D2 I SN*D 1 , SN*D2, SN+D **/ /* TO COMPUTE COS(VMD2),COSH(VMDZ),SIN(VMD2),SINH(VMD~) ++/ /******+*)c+***+*****+**+**~++*+*~*******+*+*+*~~+~.~**~*~**/ 

LOOPM : DO Mx-6 TO 5 : 

UMDl(M) = U2(M) * D1: 

46 
8 

128. 

129. 

1 30. 

131. 

132. 

133. 

134. 

135. 

131. 

137. 

138. 

139. 

140. 

141. 

142. 

143. 

144. 

145. 

146. 

147. 

148. 

149. 

150. 

151. 

152. 

153. 

154. 

155. 

156. 

157. 

158. 

159. 

160. 

161. 

162. 

163. 

164. 

164. 1 

165. 

166. 

167. 

168. 

169. 

170. 

171. 

172. 

173. 

174. 

175. 

176. 

177. 

178. 

178.1 

179. 

180. 

101. 

182. 

182.1 

183. 

181. 

185. 

186. 

187. 

188. 

189. 

190. 

VMD21Mj = V(M) * D2 : 

SNDl(M) = SM<M) * D1 : 

SND2(M) = SM(M) * 02 ; 

SND(M) = SM(M) * D j 

SINVD2 (M! = SIN (VMDZ (M) ) ; 

SINHVD2(M) = SINH (VMD2 (M) ) ; 

COSVDE(M) = COS(VMD2(M) > : 

COSHVD2(M) = COSH (WID2 (M) ) j 

END LOOPM : 

/ * * ** * * * *?% ** ** * ** * * *S *.E*** * ** ** * * 9 * .Y * * * * * *. * *.* * * * .iC * * ?t+ *.* * * ic * / 

/* */ /* TO COMPUTE Y (M, N) ++ic+.~*+***sC**********.Y~*~ic*ic**~**.~S~~; 

LOOPYM : DO M = -4 TO 5 ; 

LOOPYN : DO N = -6 TO 5 ; 

Y1 (M,N) = COS(UMD1 (M) )+SIN(SND(N)): 

Y4 (M,Nj = SIN (UMD1 (M) ) + SIN (SND (N) ) /U2 (M) : 

IF RCM) >. (3 THEN 

DO j 

Y3 (M,N) =SINVD2 (M) /V (M) +SIN(UMDl (M) 1 *COS (SND (N) ) /U2 (Mj : 

Y2(M7N) zCOSVD2 (M) -COS (UMDl (M) ) *COS (SND iN) ) j 

END; 

ELSE 

DO ; 

Y3(M,N)=SINHVD2(M) /V(M)+SIN(UMDl (M) )*COS(SND(N) ) /U2(Mj ; 

YZ(M,N)=COSHVD2(M)-COS(UMDl (M) )*COS(SND(N) j : 

END ; 

/+******.E*******+*******+******+ +** **.E**********S*~*****~/ /* TO COMPUTE ABS(Y(M,N)**Z) ** / /* TO COMPUTE REAL AND IMAGINARY PART OF Y: YR, YI **/ 

YMAG (M , N) =Y3 (M , N) **2+Y4 (M , N) +*2 ; 

YR (M, N) = (Y 1 (M, N) *Y3 (M, N) +Y4 (M, N) *YZ (M , N j j lYMfiG (M , N) : 

Y I (M, N) = (Y2 (M, N) *Y3 (M, N) .-Y 1 (M ,N) *Y4 (M, N) ) /'<MAG (M, N j : 

END LOOPYN : 

END LOOPYM : 

/+**+Y**+Y*+*+*++**+******+*+**+*++****************~~**********/ /* TO COMPUTE Jl(M,N) */ /* REAL AND IMAGINARY PART OF JI : JlR, J11 */ 

PUT SKIP ; 

LOOPJlM : DO M = -4 TO 5 j 

DO N = -6 TO 5 ; 

JlRl (M,N)=YI (M,N)+(l-COS(UMDl (M) )*COS(SNDl (N) )! 

JlR2(M,N)=(SM(N)+Yh(M,N) )*COS(UMDl (M) )*SIN(SNDl (N) ) ; 

JlR3(M,N)=(U2(M)+SM(N)*YR(M,N) /UZ(M) )+COS(SNDl (N) ) 

JlR4(M,N)=SIN(SNDl (N) )*SM(N)*YI (M,N) /U2(M)j 

J 1R (M, N) = (J 1R1 (M, N) -J lh2 (M, N) +SIN (UMD 1 (M) ) * 

JlIl (M,N)=(SM(N)+YR(M,N) )*(COS(UMDl (M) )*COS(ShlDl (N) )-i) j 

; 

(JlR3(M,N)-JlR4(M,N)) )/(U2(M)*+2-SM(N)++2); 

Jl IZ(M,N) =YI (M,N) *COS (UMD1 (M) ) +SIN (SND1 (N) ) : 

JlI3(M,N)=SM(N)*YI(M,N)+COS(SNDl (N) )/U2(M); 

JlI4(M,N)=SIN(SNDl (N) )*(U2(M)+SM(N)*YR(M,N) /U2(M) ) j 

J1 I (M,N)f(Jl I1 (M,N)-Jl I2 (M,N) +SIN(UMDl (M) ) * 

(J 1 I3 (M , N) +J 1 14 (M , N) ) ) / (U2 (M) +*2-SM (N) ++2) ; 

END ; 

END LOOPJlM; 

/Y~*+*+**~+Y~~*Y~~****+~***++ +** ******+~*+****~~~***~*~.E*****/ 

/* TO COMPUTE J2(M,N) */ 

/* REAL AND IMAGINARY PART OF 52 : JZR, J2I */ 

PUT SKIP; 

LOOPJ2M : DO M = -5 TO 5 : 

DO N = -5 TO 5 : 

J2H4 (M , N ) -SM (N) +Y I (M , N ) *SIN ( SND2 ( N ) ) /V 11 ; 

J2IJ(M,N)==SM(N)*YI (M,N)*COS(SND2iNj j iVi<M) : 

IF R(M) 30 THEN 

DO; 

J2Hl (M ,N) SY I (M, N) * (COSVDZ (M) *COS (SND2 (N) ) -1 1 ; 

,12R2(M.N)=(SM(N)+YR(M.N) )*COSVD2(M)*.SIN(SND2(N) ): 

47 
191. 

192. 

193. 

194. 

195. 

196. 

197. 

198. 

199. 

2(j(:). 

20 1 . 202. 

203. 

204 I 

2(:)5. 

206. 

207. 

2(:)8. 

209. 

210. 

211. 

212. 

212.1 

253. 

214. 

215. 

216. 

217. 

218. 

219. 

220. 

221. 

222. 

223. 

224. 

225. 

226. 

227. 

228. 

339, 

230. 

231. 

232. 

233. 

234. 

235. 

235.1 

236. 

261. 

262. 

262.1 

262.2 

262. ,7 

263. 

264. 

265. 

265.2 

266. 

267. 

268. 

268.1 

268.2 

268.3 

)j 

265.3 

M) ) 

J2H3 (M, Nj = (V (M) +SM !N) *.YR (M, 1\11 /V (Mj ) *COS (SND2 (N) j : ' 

J2R (M , N! (J2R1 (M, N) -J2RZ (M, N) +SINVD2 (tl) * 

( J ZR3 ( M , N ) + J 2R4 ( M , N ) ! ) / ( V ( M ) -%+t2-SM C N ) * *2 j ; 

JZI 1 (M , N) = (SM (N) +YH (M , N) ) * (COSVD2 (M) *COS iSND2 (N) ) -1.) ; 

J2IZ(M,N)=COSVD2(M)+SIN(SND2(" )*YI (M,N): 

J214 (M, N) =SIN (SND2 (N) ) * (V (M) +SM (N) *YR (M , N) /V (M) ) : 

J2I (M,N)=(-J211 (M,N)-J212(M,N)+SINVD20+ 

(J213(M,N)-J214(M,N) 1 )/ (V(M)**2-SM(M)+*2) j 

END : /*END OF R(M) > 0 FOR. 52 COMPUTATIONS ***I 

ELSE 

DO ; 

J2R1 (M, N) =Y I (M , N) * (COSHVD2 (M) *COS (SND2 (N) ) -1 ) j 

J2R2 (M, N) = (SM (N) +YR (M, N) ) *COSHVD2 (M) *SIN (SND2 (N) ) : 

J2R3(M,N)=(-V(M)+SM(N)*YR(M,N)/V(M) )*COS(SND2(N) ) j 

JZR (M , N) =- (J2R 1 (M , N) -J2R2 (M , N) +S INHVD2 (M) * 

J2I 1 (M, N) = (SM (N) +YR (M, N) ) * (COSHVD2 (M) *COS (SND2 (N) ) - 1 ) : 

5212 (M, N) =COSHVD2 (M) *SIN (SND2 (N) ) *Y I iM , N) : 

J214 (M, N) =SIN (SND2 (N) 1 * (-V (MI +SM (N) *YR (MI N) /V (M) ) j 

J2I (M,N)=-(-J21 1 (M,N) -JZI2 (M,N) +SINHVD2 (M) *(J213 (M,Nj 

(J2R3 (M , N) +J2R4 (M, N) ) ) / (V (M) **2+SM (N) **2) j 

-J214(M,N)) )/(V(M)**2+SM(N)**2); 

END j /* END OF R(M) e:= 0 FOR 52 COMPUTATION ***/ 

END ; 

END LOOPJZM j 

................................................................ 

/* */ /* TO COMPUTE SUMJlJ2 */ /* REAL AND IMAGINARY PART OF SUMJlJ2 : SUMJR , SUMJI */ 

SUMJM : DO M = -5 TO 5 ; 

SUMJN : DO N -5 TO 5 j 

SUMJR (M, N) =J 1R (M, N) +J2H (M, N) : 

SUMJI (M,N)=JlI (M,N)+JZI (M,N) : 

SUMJ 1 J2 (M, N) =SUMJR (M, N) **2+SUMJ I (M , N 1 **2 : 

END SUMJN; 

END SUMJM ; 

/****** TO COMPUTE ALPHA(M), M=-6 TO 5 ***************/ 

/* TRUNCATE INFINITE SUM TO SUM SUMJlJ2(M,N) */ /* */ 

ALFA : DO M -5 TO 5 ; 

SUMLL : DO LL = -5 TO 5 ; 

SUML(M)=O : PUT SKIP : 

SUML(M) = SUML(M) + SUMJlJ21M,LL) j 

END SUMLL; 

I***** TO COMPUTE ALPHA(M) ********/ 

/* */ 

ALPHFI (M) = SQRT (SUML (M) ) : 

PUT SKIP EDIT('ALPHA(',M,')=',ALrHA(M)) (X(5) ,A,F(2,0) ,A,E(12,5) 

END ALFA; /** TO COMPUTE CN(N),GAMA(M),AND MN(N,M) M,N = -6 TO 5 **/ 

/ ** 

LOOPCN: 

GRMADM: 

**/ 

PUT SKIP; 

PUT SKIP EDIT(REPEAT('*',40)) (X(3) ,A); 

PUT SKIP; 

DO N = -5 TO 5; 

CDISC (N) = K0**2-SM (N) **2 j 

CN (N) = SQRT (ABS (CDISC (N) ) ) : 

PUT SKIP EDIT ( 'CN ( ' ,N, ' ) = ' , CN (N) , 'CD ISC ( ' , N, ' ) = ' , CDISC (N! ) 

END LOOPCN : 

DO M = -5 TO 5 : 

G~MA(M)=S~RT(AbS(GAMAD(M))): 

PUT SKIP EDIT('GAMA(',M,')=',GAM~(M),'G~MAD(',M,'!=',GAMAD( 

(X(3),2 (A,F(3,(3),A,E(12,5) ,X(2))): 

GAMRD (M) = t:::C)**2itE1-U2 (M) **2 : 

(X (3) ,2 (A.F(3.0) .A.E(12.5) .X (2) ) ) : 

48 
. 

269. 

270. 

271. 

272. 

273. 

274. 

275. 

276. 

277. 

278. 

279. 

280. 

281. 

282. 

283. 

284. 

285. 

286. 

287. 

288. 

289. 

290. 

291. 

292. 

293. 

294. 

295. 

295.1 

295.2 

296. 

297. 

298. 

299. 

300. 

301. 

302. 

303. 

304. 

305. 

305.1 

306. 

307. 

308. 

309. 

310. 

311. 

312. 

313. 

314. 

314.01 

314.02 

314.03 

314.04 

314.05 

314.06 

314. 07 

3 14. 08 

314.09 

314. 1 

314. 101 

314. 102 

314. 1(:13 

314.104 

314.105 

314. 106 

9 

. .. . i. :<'. I, . 

pub& ,@jAnm 

LDCIPMNA : DO N = -5 TO 5 ; .. 

. ,, .. END GANADM ; ' 

LOOPAMN: DO M = -5 TO 5; 

AMNR(M,N)=SUMJR(M,-N) /ALfiHA(M) ; 

ANN1 (M,N)=SUMJI (M,N)/ALPHA(M): 

IF CDISC (N) > 0 8c GAMAD iM) ? (3 THEN 

CASE1 : DO ; 

MNR (N , M) = (CN (N) +GAMA (M) ) *AMNR i M, N) j 

MNI (N,M) = (CN(N)+GAMA(M) )*AMNI (M,N) : 

END CASE1 ; 

ELSE 

DO ; 

CASE2 : IF CDISC(N) < 0 b GAMAD(M) > 0 THEN 

MNR (N, M) = GAMA (M) *AMNR (M, N) -CN (N) *AMNI (M, N) : 

MNI (N,M) = CN(N)*AMNR(M,N)+GAMA(M)*AMNI (M,N) : . 

END ; /*** CASE2 ****/ 

ELSE 

CASE3 : IF CDISC(N) ? 0 $4 GAMAD(M) .:E (1) THEN 

DO ; 

MNR (N, M) =CN (N) *AMNR (M, N) -GAMA (M) *ANN I (M , N) 

MNI (N, Pi) =GAMA (MI *AMNR (M , N) +CN (N) *ANN1 (M, N) ; 

END ; /*** CASE 3 ***/ 

ELSE 

CASE4 : DO ; 

: 

MNH(N,M)=-(CN(N)+GAMA(M) )*AMNI (M,N) ; 

MNI (N,M)=(CN(N)+GAMA(M) )*AMNR(M,N) : 

END CASE4 ; 

END LOOPMNA ; 

END LOOPAMN ; /** TO PRINT REAL AND IMAGINARY PART OF MATRIX MN(M',N) **/ 

/** : MNR(M,N) AND MNI(M,N) : M , N = -6 TO 5 ** / /*** *+*/ 

PUT SKIP; 

PUT SKIP EDIT(*PRINT VALUES OF MNREAL AND MNIMAG ') 

PUT SKIP ; 

PUT SKIP EDIT(REPEAT('*',SS)) (X(3) ,A): 

(X (6) ,A) ; 

PUT SKIP EDIT(REPEAT('*',55)) (X(3) ,A); 

MNPRINT: DO N = -5 TO 5 ; 

PUT s#Ir; 

PUT SKIP EDIT( 'N=' ,N) (X (7) ,A,F(3,0) ) ; 

PUT SKIP EDIT (REPEAT('*',lO)) (X(5),A); 

PUT SKIP ; 

PUT SKIP EDIT('M','MNREAL','MNIMAG') (XG) ,~,2 (x(~),A(~o))) 

PUT SKIP i 

NMPRINT: DO M = -5 TO 5; 

PUT SKIP EDIT(M,MNR(N,M) ,MNI (N,M)) (X(2) ,F(2,0) ,2 E(19,5)); 

END NMPRINT; 

END MNPRINT ; /** TO COMPUTE RIGHT HAND SIDES VECTORS FOR T(II) AND *******/ 

/+* T(I11) ,B2 AND 83 ,REAL AND IMAGINARY PARTS *****/ /* B2(N)=(GAMA(O)-CN(N) )+AMN(O,N), N=l,Cl,-l,-2 * * ** *'/ 

/*+**** *****/ 

/* E3 (N) = (GAMA (-1 ) -CN (N) 1 *ANN (-1, N) , N=l ,(I, -1, -2 *****/ 

I-OOPB2 : DO N = 1 TO -2 BY -1 j 

IF CDISC(N) :* 0 THEN 

DO : 

B2R (N) ziGAMA ((5) -CN (N) ) +AMNR ((),N) 

E21 (N) 22 (GAMA (0) --CN (N) ) +AMNI (O,N) 

: 

; 

END : 

DO ; 

B2R (N) =GAMA (0) *AMNR (0, N) +CN (N) *AMNI ( 0, N) 

B2I (N) =GAMA (0) *AMNI (0, N) -CN (N) +AMNR ((1, N) 

END : END LOOPBZ i 

ELSE 

; 

; 

49 
a 

.z 14. 107 

314. 108 

314. 109 

314.11 

314.12 

314.13 

314. 14 

314.15 

314.16 

314.18 

314.181 

314.19 

314.191 

314.192 

314.193 

314. 194 

.314.195 

314.196 

314.197 

314.198 

314.199 

314.201 

314.21 

,314.22 

314.23 

314.24 

314.25 

314.26 

314.27 

314.28 

314.29 

314.3 

314.31 

314.32 

314.33 

314.34 

314.35 

314. 36 

314. 361 

314.37 

314.38 

314.39 

314.4 

314.41 

314.42 

314.43 

314.44 

314.45 

514.46 

314.47 

314.48 

314.49 

314.5 

314.51 

314.52 

314.5.3 

314.54 

314.55 

314.56 

314.57 

315. 

316. 

314.2 

3i7. 

318. 

J19 . Y 

.. 

LOOPH3: DO N=l TO -2 BY -1 ; /* I32 (tdj FOR REAL 4. IMOGINARY GGMA i--S j */ 

IF CDISC (N) :> 0 Pc GAMfiD (-1 j 1.:. Ci 

'THEN DO ; 

B.~R(N)=!GAMn(-l)-CN(N) )*RMNFr(-l ,N) : 

B3I (N) =(GAMA (-1 j -CN (N) ) *RMNL (-1 ,N) ; 

. END ; /* C(N) GAMA(-l). ARE BOTH REAL *./ 

ELSE 

IF CDISC(N) <: 1:) t GAMAD(-l) 1:. 0 

THEN DO j 

P3H(N)=GAMA(-l!+AMNH(-1 ,N)+CN(N)*fiMNI (-1 ,Nj; 

PSI (N) =GAMA (-1 ) WMNI (-1, N) -CN (N) *AMNR (- 1 , N) 

END ; /* END OF C(N) IMAGINARY % GAMA(-I) REAL */ ; 

ELSE 

IF CDISC(N) >, 0 & GAMAD(-l) .:; 0 

THEN DO ; 

B3H(N)=-(GAMA(-l)+AMNI (-l,N)+CN(N)*AMNR(-l ,N)) j 

B31(N)=GAMA(-l)*AMNH(-l,N)-~N(N~~A~NI (-1.N): 

END ; /* C(N) REAL $4 GAMA(-l) IMAGINARY */ 

ELSE 

DO j 

b.~R(N)=(CN(N)-GAMA(-l) )*AMNI (-1 ,N) ; 

B3I (N)=(GAMA(-l)-CN(N) )*AMNR(-l,Njj 

END ; /* C(N) t GAMF\(-l) ARE BOTH IMAGINARY */ 

END LOOPBS ; /* END COMPUTING B2 AND B3 */ /** TO PRINT EZR(N) ,b2I (N) ,BSR(N) ,RND B.31 (N) , N=l,O,-1,-2 +*/ 

/** ** / 

PUT SKIP ; 

PUT SKIP ; 

PUT SKIP EDIT(REPEAT( '*' ,55) ) (X (3) ,A) 

PUT SKIP EDIT('PR1NT REAL AND IMAGINARY PART OF B2 &B3') 

(X (6) ,A) ; 

PUT SKIP EDIT(REPEAT ('*',SEI)) (X(3),A) j 

PUT SKIP ; 

PUT SKIP EDIT('N','E~REAL','B~IMAG','EGREAL','B.~IMAI~') 

PUT SKIP ; 

(X(S),A,4 (X(3) ,A(lC?))) ; 

PRINTE: DO N = -2 TO 1 ; 

PUT SKIP EDIT(N,E2R(N) ,821 (N) ,B3R(N) ,E31 (N)) 

(X(2) ,F(2,:3> ,1 E(15,5)! : 

END PRINTB; /** END PRINTING E3 AND b3 **/ /** PRINT AMNR(M,N),AMNI(M,N): M,N = -6 TO 5 *****+*i~*/ 

/ ** ** / 

PUT SKIP; 

PUT SKIP ; 

PUT SKIP EDIT(REPEAT( '*' ,SO) 1 (X (.I) ,A) ; 

PUT SKIP EDIT('PR1NT VALUES OF AMNREAL AND AMNIMAG') 

PUT SKIP EDIT (HEPEAT( '*' ,55)) (X (3) ,A) j 

PUT SKIP ; 

PUT SKIP EDIT('M=',M) (X.(7) ,R,F(3,0)) : 

PUT SKIP EDIT (REPEAT('*',lO)) (X(S),A) ; 

PUT SKIP; 

PUT SKIP EDIT( 'N' , 'AMNREAL', 'AMNIMAG' ) 

PUT SKIP; 

PRINTAN: DO N = -5 TO 5 ; 

(X (6) ,A) ; 

PRINTAM: DO M = -5 TO 5 ; 

(X(3),A,2 (X(EI),A(lO))); 

PUT SKIP EDIT(N,AMNR(M,N) ,AMNI (M,N) 1 

END PRINTAN : 

(X(2) ,F(2,0) ,2 E(15.5)) : 

END PRINTAM ; 

/ *. *.*.** * * 4t ** Y * * *++ * * * * * *+** * ** ** * Y * * ******.***-*I. * * ., * .* * * * * * * * * * * / 

/* SUBROUTINE TO INPUT DATA +/ 

INPUT: PHOC(M) ; 

GET LIST (MI: 

DCL M FIXED (4,O); 

50 
a 

--I - .-. L 'J . 321 I 

322. 

323. 

324. 

325. 

326. 

.4L7. 

,328. 

329. 

330. 

331. 

332. 

.-.a4 I 

rJ35. 

.336. 

337. 

-7.7 

77-7 .-#34 . 7- 

7- .>. a. 1 

339. 

340. 

341. 

342. 

343. 

344. 

,345. 

346. 

347. 

349. 

350. 

351. 

352. 

353. 

354. 

355. 

446. 

357. 

358. 

359. 

360. 

361. 

362. 

363. 

364. 

365. 

.366. 

367. 

368. 

369. 

369.01 

369.02 

369.03 

369.04 

369.04 1 

367.042 

3.59. 043 

369.044 

367.045 

367.046 

369. 05 

369.06 

369.11 

369.12 

369.13 

348. 

7C 

/* /* 

ORIGINAL PAGE IS 

OF POOR QUALITY 

*/ 

*/ 

PUT SKIP (5) EDIT ( 'U1 ' , 'U2', 'F (U1) ' , 'F' ' (U1) ' , :U2-U1 : ' , 'COUNT') (X (3) ,n,X (12) ,A, X ( 10) ,A, 

X(7) ,A,X(8) ,A,X12) ,A); 

END j 

RETURN: 

END INPUT; 

/* /* SUBROUTINE TO PERFORM CALCULATION 

CALCULATE: PROC (M) : 

DCL M FIXED (4,O) ; 

U2 (M) =U1 (M) -F (U1 (M) ,M) /FPRIME (U1 (Mj , M) : 

K=K+1: 

RETURN; 

END CALCULATE ; 

/* */ /* SUBROUTINE TO PRINT TABLE */ 

PRINT: PROC (M) ; 

PUT SKIP EDIT(U1 (M) ,U2(M) ,F(U1 (M) ,M) ,FPRIME(Ul (M) ,M) 

DCL M FIXED(4,O); 

ABS(U2(M)-Ul (M) 1 ,K) 

(5(E(12,5) ,X(l)) ,F(2,0)); 

RETURN ; 

END PRINT; /* */ 

/+ SUBROUTINE TO PRINT FINAL RESULTS */ 

OUTPUT: PROC(M); 

DCL M FIXED(410) ; 

PUT SKIPCS) EDIT('APPHOX1MATE ROOT U2= ',U2(M), 

RETURN; 

'F(U2)= ',F(U2(M) ,M)) (A,E(14,7) ,X(5) ,A,E(14,7) ); 

END OUTPUT: 

TO COMPUTE SCATTRING MATRIX SCAT */ 

CONSTRUCT COEFFICIENT MATRIX AA FROM MATRIX MN(N,M) */ 

MM= 8; 

NN= a 

AAIl: DO I = 1 TO NN/2 j 

AAJl: DO J = 1 TO MM/2; 

AA ( I, J ) =MNK (2-I,2-J) : 

AA(I,J+MM/2)=-MNI (2-1,2-J) j 

END AAJl : 

END RAIl ; 

PIVOT= 1.C) : 

AAI2: DO I = NN/2+1 TO NN ; 

AAJ2: DO J = 1 TO MM/2 ; 

RS- 3 ; 

AA ( I. J 1 =MNI (NN/2+2-I. 2-J ) : 

51 
369.14 

309.15 

369.16 

,369.2 1 

369.22 

,369.23 

369.24 

369.25 

369.29 

369.31 . 369 32 

369.33 

369.34 

369.35 

369.36 

369.37 

369.372 

369.38 

369.39 

369.4 

400. 

404. 

407. 

4O8. 

409. 

4 1 0 . 456. 

457. 

459. 

460. 

461. 

462. 

463. 

464. 

465. 

466. 

467. 

468. 

469. 

470. 

477. 

478. 

479. 

490. 

491. 

492. 

483. 

484. 

48s. 

486. 

487. 

488. 

489. 

490. 

491. 

492. 

493. 

494. 

495. 

496. 

497. 

498. 

499. 

5 0 (3 . 501. 

502. 

44 ( I, J+MM/2) :;:MNR (NN/2+2- I , '2-J) : 

END A&J2 ; 

END AAI2 : 

RA (I ,MM+1) =(I ; 

END RHSl : 

RA (2,MM+lI =2*CN (0) : 

RtlS23R: DO I = 1 TO NN/2 ; 

RHS1: DO I = 1,3 TO NN : 

i+ END OF RHSl COLUMN ' */ 

AA(1 ,MM+2)=B2R(2-I) : 

AA ( I , MM+3) =B.3R (2- I ) ; 

END RHS23R ; 

RHS23I: DO I = NN/2+1 TO NN : 

/* EN,D OF RHS2 COLUMN OF AUGMENTED MATRIX IN GAUSS PROC.*/ 

AA(I,MM+2)=E2I(MM/2+2-I) ; 

AA(I,MM+3)=B3I (MM/2+2-I) j 

END RHS23I ; 

/* END OF RHS3 COLUMN OF AUGMENTED MATRIX IN GAUSS FROC. .*/ 

/* INVOKE GAUSS ELIMINATION TO COMPUTE SCATTERING MATRIX +/ 

CALL GAUSS(AA,MM,NN,PIVOT,~S?X) j 

GAUSS: PROC(AA,M,N,PIVOT,RS,X) ; * i 

*/ 

*/ 

/+ 

/* GAUSSIAN ELIMINATION WITH OR WITHOIJT PIVOTING. /* ANSWERS ARE THEN SUBSTITUTED BACK INTO THE 

/* ORIGINbL EQS. WITH MULTIPLE RHS VECTORS. */ 

/* */ 

DCL (AA(*,*) ,X (*,it)) FLOAT ; 

DCL (M, N, PIVOT, RS) F I XED (5,O) ; 

START: BEGIN; 

DCL (AC(M,N+RS) ,BB(M,N+RS) ,HOLD(N+RS) ,SUM! FLOGT(6), 

. XX(RS,N) FLOAT(6) INIT( (F:S*N)(:I): 

/+ INPUT AUGMENTED MATRIX * i 

CALL INPUTl; /* CONVERT TO UPPER TRIANGULAR MATRiX 

CALL UPTRI; 

/+ BACK SURSTITUTE 

CALL BACKSUB; 

CALL OUTPUT1; 

/+ PUT ANSWERS BACK IN ORIGINAL EQUATIONS 

CALL TESTi; 

/+ SUBROUTINE TO INPUT AUGMENTED MATRIX 

INPUT1 : PROC; 

*/ 

*/ 

*/ 

/* */ */ 

PUT PAGE EDIT('GAUSS1AN ELIMINATION') (X(28),A); 

IF PIVOT=l THEN PUT SKIP EDIT('W1TH PIVOTING') 

ELSE PUT SKIP EDIT('WITH0UT F'SVOTING') (X(30),A); 

PUT SKIP EDIT('F0R ',M,' BY ',N,' MATRIX') 

(X(31) ,a); 

(X(29) ,A,F(2,0) ,A,F(2,0) ,A): 

PUT SKIP EDIT('WITH',RS,'RIGHT HAND SIDES') 

(X (29) ,A,F(3,0) ,X (2) ,A) : 

PUT SKIP(5) ; 

DO 1=1 TO M; 

DO J=1 TO N+RS; 

AC(I,J)=AA(I,J) ; 

PUT EDIT(AA(1,J)) (X(1) ,F<El,T):l: 

BB( I, J) =AC (I, J) : 

END; 

PUT SKIP; 

EhlD: 

RETURN: 

EPJD T NF'UT 1 : 

/* 

/* SUBROUTINE TO PRINT MATRIX 

PRINT: PHOC; 

PUT SKIP(5) : 

52 
fC13. 

504. 

505. 

506. 

507. 

508. 

509. 

510. 

511. 

512. 

513. 

514. 

515. 

516. 

517. 

518. 

519. 

520. 

521. 

522. 

524. 

525. 

526. 

527. 

528. 

529. 

530. 

531. 

532. 

533. 

J44. 

535. 

536. 

537. 

538. 

539. 

540. 

541. 

542. 

543. 

544. 

545. 

546. 

547. 

548. 

549. 

550. 

551. 

552. 

553. 

554. 

555. 

dab. 

557. 

558. 

559. 

560. 

561. 

561.1 

562. 

569. 

564. 

565 

566. 

567. 

e-7 JL3. 

ccc 

DO I=l TO M; 

DO J=1 TO N+KS: 

END: 

PUT SKIF; 

FUT EDIT(AC(I,J!) (X(.l)TF(8,.Z))j 

END: 

RETURN: 

END PRINT; /* */ 

/+ SUBROUT. CONVERTS MATRIX TO UFF'ER TRIANGULAR *./ 

UPTR I : PKOC jDO K-1 TO M-1; 

IF PIVOT-1 THEN CALL PIVOT1; 

DO I=K+l TO M; 

RATIO = AC(I,K)/AC(K,K) 

DO J-K TO N+RS: 

AC ( I , J ) SAC ( I , J ) -RAT I O* AC ( K J ! : 

END : 

END ; 

END; 

RETURN; 

END UPTKI: /* SUBROUTINE TO USE PIVOTING */ 

PIVOT1 : FROC; 

P=K: 

DO I=K+1 TO M; 

END : 

IF P."=K THEN DO; 

IF ABS(AC(P,K)) e: ABS(AC(1,K)) THEN F' = I: 

DO J=1 TO N+RS: 

HOLD (J 1 =AC (K, J) ; 

AC (K , J ) =AC (P , J ) j 

FIC (P. J) =HOLD (J 1 ; 

END; 

END: 

RETURN; 

END PIVOT1; 

/* 

/+ SUBROUTINE 

DO K: = 1 TO 

DO I=N TO 1 

SUM=O; 

BACKSUB: PROC: 

TO BACK SUBSTITUTE 

RS ; 

BY (-1) ; 

DO J=I TO M; 

SUM=SUM+XX (K , J) *AC ( I, J ) j 

END; 

XX (K, I ) = (FIC ( I, N+K) -SUM) /AC ( I, I ) ; 

END; 

END ; 

RETURN; 

END BACESUB; /* /* SUBROUTINE TO PRINT ANSWERS 

OUTPUT 1 : PROC ; 

PUT SC::IP(5) EDIT('ANSWERS') (X(34) ,A); 

DO J = 1 TO.RS ; 

PUT SKIP EDIT ('SET',J) (X(20),A,F(3,0)); 

PUT SKIF; 

DO 1=1 TO N; 

END; 

END ; 

RETURN; 

END OUTPUTI: 

*/ */ 

*/ */ 

53 
568. 

569. 

570. 

571. 

572. 

573. 

'574. 

575. 

576. 

577. 

578. 

579. 

580. 

581. 

582. 

583. 

584. 

585. 

'586. 

587. 

588. 

589. 

589.1 

590. 

591. 

591.1 

591.2 

592. 

592.1 

592.2 

592.3 

592.4 

592.5 

592.6 

593. 

594. 

595. 

596. 

597. 

598. 

599. 

600. 

601. 

602. 

603. 

604. 

605. 

606. 

607. 

608. 

609. 

610. 

611. 

612. 

613. 

614. 

615. 

616. 

617. 

618. 

619. 

620. 

621. 

622. 

623. 

624. 

/* */ 

/* SUBROUTINE TO PUT ANSWERS EtACC; IN ORIGINAL- 

/* EQUAT IONS ic; 

TEST1: PROCj 

PUT SKIP(S) EDIT 

('ANSWERS PUT IN ORIGINGL EBUATIONS') (X(21),H); 

DO K = 1 TO RS ; 

FlJT SKIP EDIT ( 'SET ' ?k:) (X(Y?) ,A,Fi2,(3j ) ; 

PUT SKIP; 

DO 1=1 TO M: 

PUT SKIP; 

SUM=(:) j 

PUT SKIP EDIT(' ') (X(1) !A); 

DO J=l TO N; 

SUM=SUM+BB( I, J) *XX (t, J) : 

PUT EDIT(HE(I,J),'X(',J,') ') 

IF J':N THEN PUT EDIT( '+ ' j (A) : 

ELSE IF J=N THEN PUT EDIT( '= ') (4) : 

(F(9,3) ,A,F(l,O) ,Ai: 

END; 

PUT EDIT(SUM) (F(9,3)): 

END ; 

RETURN: 

END TEST1 : 

END j 

END START ; 

END GAUSS j 

/* */ 

CALL LOOPT; /* PROCEDURE TO COMPUTE T VECTORS */ 

/* TO PRINT OUT T VECTORS PROCEDRE CRLL */ 

CALL PRINTT; 

CALL INNER : /* FROCEDRUE TO COMPUTER DOT PRODUCT */ 

CALL SCATTER; /* PROCEDURE TO CONFUTE SCATTERING MATRIX*/ 

CALL RTL ; /* PROCEDURE TO COMPUTE RTL MATRIX */ /* PROCEDURE TO COMPUTE T VECTORS */ 

/+ TR AND TI ARE REAL AND IMAGINARY PART OF T VECTORS +/ 

/* */ 

LOOPT: moc; 

DO K = 1 TO S ;- 

DO N = 1 TO -2 BY -1 ; 

TR (E, N) =X (E, 2-N) j 

TI (K,N)=X(K,6-N).: 

END: 

END: 

END LOOPT : 

RETURN; 

/* PRINT T VECTORS, TREAL AND TIMAGINARY */ /* */ 

PRINTT: PROC; 

PUT SKIP; 

PUT SKIP EDIT(REPEAT('*',55)) (X(S),A); 

PUT SC:IP EDIT('PR1NT T VECTORS FOR N = 1,0,-1,-2 ') 

PUT SKIP EDIT(REPEAT('*',S5)) (X(3),A): 

PUT SKIP : 

DO K = 1 TO 3 : 

PUT SKIP: 

PUT SKIP EDIT('K=',K) (X(10) ,A,F(.T,O)): 

PUT SKIP EDIT (REPEAT( '*', 10)) (X (5) ,A) : 

PUT SKIP: 

PUT SKIP EDIT('N','TREAL'.'TIMAG') (U(7) ,n.2 (X(4j .Ail(:))!); 

F'IJT SKTP: 

DO N -2 TO 1 : 

FWT SKIP EDIT(N,TR(K,N) ,TI (K,N)) (X(5) ,F(3,0) ,2 E(15,5) ); 

PUT SKIP: 

END : 

(X (7) ,A); 

54 
625. 

626. 

627. 

628. 

629. 

631. 

632. 

633. 

6.34. 

635. 

636. 

637. 

638. 

639. 

640. 

641. 

642. 

643. 

644. 

645. 

646. 

647. 

648. 

649. 

65O. 

651. 

652. 

653. 

654. 

655. 

656. 

657. 

658. 

659. 

660. 

661. 

662. 

663. 

664. 

664.1 

665. 

666. 

667. 

668. 

669. . 

670. 

671. 

672. 

673. 

674. 

675. 

676. 

677. 

678. 

679. 

680. 

601. 

682. 

683" 

c,o/l. 

oE35. 

686. 

607. 

688. 

689. 

6.3 (11 . 

END: 

END PRINTT; 

/+ PROCEDURE TO COMPUTE INNER PRODUCT OF AMN AND T MkTK'ICES-~./ 

/* +/ 

INNER : PROC; 

DO K= 1 TO 3: 

ROR(K) = 0; 

ROI(K) = (:) ; 

DO N = -2 TO 1 : 

ROR (K) =ROR (t:) +AMNR (N, (:)) +TR (K:, N) -AMNI i~, (11 *-r I (K , N) : 

ROI (K) =Ha1 (K) +AMNI (N, (:I) +TR (K , N) +AMNR (N, o) *TI i~? N) : 

ROR (1) =ROR (1) -1 j 

END ; 

END; 

ROR (2) =ROH (2) +AMNR (0,O) j 

KO1 (2) =ROI (2) +AMNI ((),(I)) : 

ROR(3)=ROR(3)+ANNH(-l ,O) : 

ROI (3) =ROI (3) +AMNI (-1 ,(:I) ; 

/* PRINT INNER PRODUCT OF AMN(N,O)AND T(b::,N) */ /* */ 

PUT SKIP; 

PUT SKIP EDIT (REPEAT( '*' ,55) ) (X (7) ,A) : 

PUT SKIP EDIT(REPEAT('*',55)) (Xi7),A); 

PUT SKIP EDIT('R0,INNER PRODUCT OF A h T') (X(5) ,A): 

PRINTRO: DO I = 1 TO 3 : 

PUT SKIP EDIT ( 'ROREAL ( ' , I , ' ) = ' , ROR ( I ) , 'ROIMAG ( ' , I , ' ) = ' , ROI (I) 1 (X(10) ,2 (X(2) ,A,F(2,0) ,A,E(12,5)) ): 

END PHINTHO ; 

END INNER j /* END ON INNER PRODUCT PROCEDURE */ /* PROCEDURE TO COMPUTE SCATTERING MATRIX */ 

SCATTER: PROC; 

LOOPSCl : DO J = -1 TO 1 j 

SCATR ( 1, J)=ROR (2-J) 

SCATI ( 1, J) =RO I (2-5 ) : 

END LOOPSCl; 

; 

LOOPSC2: DO I = 0,-1 : 

DO J = -1 TO 1 : 

SCATR ( I ,J> =TR (2-5, I ) : 

SCATI ( I ,J)=TI (2-5, I) : 

END ; 

END LOOPSCZ; /* PRINT SCATTERING MATRX SCATREAL, SCATIMAG PARTS */ 

/* */ 

PUT SKIP; 

PUT SKIP EDIT(REPEAT('*',55)) ()s(3),A); 

PUT SKIP EDIT('PR1NT SCATTERING MATRIX SCAT') 

PUT SKIP EDIT(REPEAT('*',55)) (X(3) ,A); 

PUT SKIP; 

DO K = -1 TO 1 ; 

PUT SKIP; 

PUT SKIP EDIT('K=',E) (X(10) ,A,F(Z,O)); 

PUT SKIP EDIT(REPEAT( '*' ,lO) ) (X (5) ,A) ; 

PUT SKIP; 

DO N = -1 TO 1 : 

PUT SKIP EDIT ( 'SREAL ( ' ,N, ' , ' , t:::., ' ) =' , SCATR (N, t:::) , 

(X(5) ,2 

(X (10) ,A) ; 

'SIMAG ( ' ,N, ' , ' , t:: , ' ) = ' ,SCAT1 (N, K:) ) 

(X(2) ,A,F(2,0) ,A,F(2,0) ,A,E(12,5) 1 ): 

END : 

END: 

RETURN: 

END SCATTER; 

/+ COMPUTE COEFFICEINT MATRIX RTL, AND THE RHS VECTOR+*/ /* OF THE SYSTEM OF EQUATIONS TO SOLVE R,R',T(O), AND */ 

/* T(-l) */ 
a 

690. 

691. 

692. 

693. 

694. 

695. 

696. 

697. 

698. 

699. 

700. 

702. 

703. 

704. 

705. 

706. 

707. 

708. 

709. 

710. 

711. 

712. 

713. 

714. 

715. 

716. 

717. 

728. 

719. 

731). 

721. 

722. 

72s. 

724. 

725. 

726. 

727. 

728. 

729. 

730. 

731. 

732. 

733. 

734. 

735. 

736. 

737. 

738. 

739. 

740. 

741. 

742. 

743. 

744. 

745. 

746. 

747. 

748. 

749. 

750. 

751. 

752. 

753. 

754. 

755. 

701. 

END : 

ELSE' 

DO : 

RTLR ( 1 ,3) =CD ( 1 1 *SI N (GAMAOZH) -Ab ( 1 1 *COS (GAMA02H ) 

+EXP (-GAMA (-1 ) ) * (CD (2) *SIN (GAMAOH) -Ab (2) it 

COS (GAMAOH ) ) : 

RTLI ( 1,s) =-(AB ( 1 ) *SIN (GAMA02H) +CD ( 1 ) *COS (GAMA02H) 

+EXP(-GAMA(-l) )*(CD(2)+SIN(GAMAOH)+AB(2)*CO~(GAMA(:)H) ) ) 5 

END; 

/* TO COMPUTE RTL(l,4) */ /* RECIL s( IMAGINARY PARTS OF S(l,O)*S(O,-l) & S(l,-l)*S(-l,-l)*/ 

AB~3-I~~SCATR~1,I~*SCATR~I,-l~-SCATI~l,I~*SCATI (1,-1); 

CD (3-1 ) =SCATR ( I, I )*SCAT1 ( I ,-I ) +SCAT1 ( 1, I )*SCATH ( I, -1 ) ; 

DO I a 0 TO -1 BY -1 ; 

END ; 

/* RTL(l,4) */ 

IF GAMAD (-1) :> 0 THEN 

DO /* RTL(1,4) WHEN GAMA(-l) IS REAL */ 

RTLR ( 1 ,4) =- (AB (3) *COS (GAMA(>- 1H) -CD (3) *SIN (GAMA(>- 1H) 

RTLI ( 1 ,4) =-(AB (3) *SIN (GAMA(:)- 1H) +CD (3) *COS (GAMA(:) 1H) 

+AB (4) *COS (GAMA-lZH) -CD (41 *SIN (GAMA-12H) ) : 

END : 

56 
I 

756. 

757. 

758. 

759. 

76O. 

761. 

762. 

763. 

764. 

765. 

766. 

767. 

769. 

770. 

771. 

772. 

773. 

774. 

775. 

776. 

777. 

778. 

779. 

780. 

781. 

782. 

783. 

784. 

785. 

787. 

788 

789. 

790. 

791. 

792. 

793. 

794. 

79s. 

796. 

797. 

798. 

799. 

800. 

801. 

802. 

803. 

804. 

805. 

806. 

807. 

808. 

809. 

810. 

811. 

812. 

813. 

014. 

015. 

816. 

817. 

818. 

819. 

3 

748. 

) 

786. 

IF GAMAD(-l) :. il THEN 

RTLR (2,rl) =- (SCATR ( 1, -1 ) *COS (GAMA-lH! -SCAT1 ( 1 , - 1 j 

RTLI (3,4) =- (SCATR ( 1, -1 ) *SIN (GAMA-lH) +SCAT1 i 1, -1. ) *COS (GAMA- 1 ti j j 

END ; 

ELSE 

DO : /* RTL(2,4) WHEN GAMA(-l! IS 1MfiGINAF;Y */ 

RTLR (2.4) =-EXP (-GAMA-lH) *SCATR ( 1, - 1 ) ; 

RTLI (2,4)=-EXP(-GAMA_lH)*SCATI (l,-l) ; 

DO ; /* RTL(2,4) WHEN GRMA(-l) IS REAL +/ 

SIN (GAMA-lH) ) : 

END ; 

/* TO COMPUTE RTL(3,3) */ 

DO I = 0 TO -1 BY -1 

/* REAL & IMAGINARY PARTS OF S(O,O)**2 $< S(C),-l!*S(-l,O) */ 

AB (5-1 i =SCATH (0, I) *SCATR ( I, 0) -SCAT1 (0, I) +SCAT1 ( I, 0) 

CD (5-1 ) =SCATR (0, I )+SCAT1 ( I ,(I)) +SCAT1 ((3 , I )*SCATR ( I , (1)) : 

j 

END ; /* END OF COMPUTING REAL 8, IMAGINARY PARTS OF SCAT */ /* RTL(3,I) */ 

IF GAMAD(-l) > 0 THEN 

DO : /* RTL(3,3) WHEN GAMA(-l) IS HEAL */ 

RTLR ( 3,3) =-AB ( 5) *COS ( GAMA02H ) +CD (5) *SIN ( GAMA02H) 

RTL I (3,;5) =- (AB (5) *SIN (GAMA02H) +CD (5) *COS i GAMA02H) 

END ; 

ELSE 

DO : /* RTL(3,S) WHEN GAMA(-l) IS IMAGINARY */ 

KTL~(3,3)=l-AB(S)*COS(GAMAO2H)+CD(S)*SIN(GAMAi~~2H) -EXF'(-GAMA_lH 

-AB(6)*COS(GAMAO-lH)+CD(6!*SIN(GAMA(:)_lHI+l j 

+AB(6)*SIN(GAMAO~lH)+CD(6)+COS(GAMA~S~1H) ) ; 

57 
* 

a 

. 

830. 

821. 

822. 

823. 

824. 

82.5. 

827. 

828. 

829. 

83(:, . 831. 

832. 

833. 

834. 

8.35. 

836. 

837. 

838. 

839. 

84(?. 

841. 

842. 

844. 

845. 

846. 

847. 

848. 

849. 

850. 

851. 

852 * 

853. 

854. 

855. 

856. 

857. 

858. 

859. 

860. 

861. 

862. 

863. 

864. 

865. 

866. 

867. 

868. 

868.1 

869. 

898. 

904. 

905. 

906. 

907. 

908. 

909. 

910. 

911. 

311.1 

512. 

914. 

915. 

916. 

717. 

825. 

843. 

- 

RTLI 14,3) =- (AB (9) *SIN (GAMAi:)sH) +CD (9) *COS (GAMA02H) 

END j 

ELSE 

DO ; /* RTL(4,3) WHEN GAMA(-l) IS IMAGINAEY */ 

RTLR (4,3) =-AB (9 ) *COS ( GAMA02H) +CD (9) *SIN (GAMA02H ) 

RTLI (4,3)=-AP!9)*SIN(GAMA(~2H)-CD(9)*COS(GAMA!:~2H) 

+AB ( 1 (? ) *S I N (GAMA(:)- 1 H ) +CD ( 10 ) *.COS < GAMA!:)- 1 H ) ) : 

-EXP (-GAMA-lH) * (AI3 (10) *COS (GAMAOH) -CD (10) *SIN (GAMAOt-4) ) ; 

-EXP(-GAMA-lH)w!AB(l(:))*SIN(GAMA(:)H)+CD( 10j*COS(GRMA(?Ii) ? F 

END j /* TO COMPUTE HTL(4,4) */ /* REAL D IMAGINARY PARTS OF S(-l,(l)*S((l,-l) D S(-1,-1)**2 */ 

AB ( 11-1 ) =SCATR (-1, I) *SCATR ( I, -1 ) -SCAT1 (-1, I 1 *SCAT1 ( I , -1 ) j 

CD ( 1 1-1 ) =SCATR (-1 , I ) *SCAT1 ( I, -1 ) +SCAT1 (-1 , I ) *SCATR ( I, -1 ) ; 

DO I = 0 TO -1 BY -1 ; 

END j 

IF GAMAD(-l) > (1) THEN 

DO j /* HTL(4,4) WHEN GAMA(-l) IS REAL */ 

RTLH~4,4~~-~AB~l1~*COS~GAMA~~~~1H~-CD~ll~*SIN~G

RTLI ~4,4~~-~AB~1l~*SIN~GAMA~~~~lH~+CD~ll~*COS~G1H) 

/* RTL(4,4) */ 

+AB (12) *COS (GAMA-lZH) -CD (12) *SIN (GAMA-12H) ) +I: 

+AB (12) *SIN (GAMA-~~H)+cD (121 *COS (GAMA-I~H) : 

END ; 

ELSE 

DO ; /* RTL(4,4) WHEN GAMA(-l) IS IMAGINARY */ 

RTLR(4,4)=1-EXP(-GAMA-lH) * (AB (1 1) *COS (GAMAOH) -CD(11) * 

SIN(GAMAOH))-EXP(-GAMA-l2H)*AB(12) : 

RTLI (4,4)=-EXP (-GAMA-lH) it (AB (11) *SIN (GAMAOH) +CD (1 1) * 

COS (GAMAOH) ) -EXP (-GAMA_12H) *CD ( 12) j 

END j /* END OF COMPUTING COEFFICIENT MATRIX RTL FOR SOLVING */ /* R AND RPRIME, T(O) AND T(-l) IN A SYSTEM OF EQUATIONS */ /* TO PRINT ELEMENTS OF RTL MATRIX */ 

PUT SKIP j 

PUT SKIP EDIT(REPEAT( '*' ,551 ) (X (3) ,A) j 

PUT SKIP EDIT('PR1NT COEFFICIENT MATRIX RTL') (X(lO),A); 

PUT SKIP EDIT(REPEAT('*',55)) (X(3),A) ; 

PUT SKIP ; 

DOI*lTO4; 

PUT SKIP: 

PUT SKIP EDIT('I=',I) (X(10) ,A,F(2,0)); 

PUT SKIP EDIT(REPEAT('*',lO)) (X(5) ,A) ; 

PUT SKIP ; 

DO J = 1 TO 4 ; 

PUT SKIP EDIT('RTLR(',I,', ',J,')=',RTLR(I,J), 'RTLI(',I, ', ', 

J,')=',RTLI(I,J)) 

(X (5) ,2 (X (2) ,A,F(2,(3) ,A,F(2,0) ,A,E( 12,s) ) ) j 

END ; /* END OF PRINTING RTL(1,J) FOR I,J = 1 TO 4 */ 

END RTL ; 

MM=8; 

NN= 8; 

PIVOT= 1; 

END ; /* END OF PRINTING RTL(1,J) FOR I,J = 1 TO 4 */ 

/* TO COMPUTE R,RPRIME,T(O) ,AND T(-1) */ 

RS=1 ; 

AAII1: DO 1=1 TO NN/2: 

AAJJl: DO J= 1 TO MM/2; 

AA(1, J) =RTLR (I , J) : 

END AAJJ1 : 

END AAIl'l; 

AAII2: DO I= NN/2+1 TO NN j 

AAJJ2: DO J= 1 TO MM/2 ; 

PtA (I, J+NN/2) ~-RTL-I ( I, J! : 

AA(I,J)PRTLI(I-NN/~,J); 

AA(I.J+NN/2)=RTLR(I-NN/Z.J): 

58 
0 

T 

918. 

919. 

920. 

921. 

922. 

923. 

924. 

925. 

926. 

927. 

929. 

930. 

1042. 

1044. 

1045. 

1046. 

1047. 

1048. 

1049. 

1050. 

1051. 

1052. 

1053. 

1054. 

1055. 

1056. 

1057. 

1058. 

1059. 

1060. 

1061. 

1062. 

1063. 

1064. 

1065. 

1066. 

1067. 

1068. 

1069. 

1070. 

1071. 

1072. 

1073. 

1074. 

1075. 

1076. 

1077. 

1078. 

1079. 

1080. 

11 10. 

1111. 

11 12. 

11 13. 

1 114. 

11 1s. 

11 16. 

928. 

. 

*/ 

I .. END RAJJ2; 

END aclrrz; 

/* RIHGT HAND VECTOR FOR THE SYSTEM OF EQUATIONS +/ 

/* TO SOLVE FOR R, RPRIME, T(O),Pr Ti-1) * / 

AA (1 ,MM+1) =SCATR ( 1,l) ; 

fiA (2,MM+1) =O j 

AA(6,MM+l)=O : 

AA ( I, MM+1) =SCATR (3-1,l) : 

AA ( I+NN/2, MM+ 1 1 =SCAT I (3-1 , 1 ) 

GA(S,MM+l)=SCATI(l,l); 

RHSCAT: DO I = 3 TO NN/2; 

; 

END RHSCGT ; 

CALL GAUSS(AA,MM,NN,PIVOT,RS,X); /* 

RREAL=X(l,l); /* REAL PART OF R */ 

RIMAG=X(l,S) : /* IMAGINARY PART OF R */ 

RPRIMER=X(l,2); /* REAL PART OF R PRIME */ 

RPRIMEI=X(1,6); /* IMAGINARY PART OF R PRIME +/ 

TOR = X(1,3): /* REAL PART OF (1,O) */ 

TOI= X(1,7) /* IMAGINARY PART OF (1,O) */ 

TMINUSlR=X(1,4); /* REAL PART OF (l,-l) +/ 

TMINUSlI=X(1,8); /* IMAGINARY PART OF (l,-l! */ 

RAbS=SQRT(RREf4L**2+HIMAG**2) j 

RPHASE=ATAND(RIMAG/RREAL); /* PHASE ANGLE OF R */ 

RPHAbS=SQRT(RPRIMEH**2~RFRIMEI**2) : 

RPRFHASE=ATAND(RPRIMER/RPRINEI~; /*PHASE OF RPRIME */ 

PUT SKIP; PUT SKIP j 

PUT SKIP EDIT('REAL OF R =',RREAL,'IMAGINARY OF R=',RIMGG) 

/* TO PRINT R */ 

(X (3) ,2 (X (2) ,A,E(12,5) ) ) : 

PUT SKIP : 

PUT SKIP EDIT ( 'ABS (R) = ' , RABS, 'PHASE (R) = ' , RFHASE) 

(X (3) ,2 (X (2) ,A,E(l2,5) ) ); /* TO PRINT R PRIME */ 

PUT SKIP; 

PUT SKIP EDIT("REAL PART OF R PRIME=',RPRIMER, 

'IMAGINARY PART OF R PRIME=',HPHIMEI) 

(X (3) ,2 (X (2) ,A,E(l2,5) )) ; 

PUT SKIP; 

PUT SKIP EDIT('ABS(RPRIME)=',RPRA~S,'PH~SE~~~~~ME)='

PUT SKIP; 

PUT SKIP: 

PUT SKIP EDIT('REAL PART OF T(O)=',TOR, 

RPRPHASE) (X(3) ,2 (X(2) ,A,E(l2,5) ) f: 

/* TO PTINT T(0) */ 

'IMAGINARY PART OF T(O)=',TOI) 

(X (3) ,2 (X(2) ,A,E(12,5) )) ; 

PUT SKIP ; 

PUT SKIP EDIT('REAL PART OF T(-l)=',TMINUSlR, 

'IMAGINARY PCIRT OF T(-i)=',TMINUSlI) 

(X(3) ,2 (X(2) ,A,E(12,5)) 1; 

END TEMODE; /* 

//GO.SYSIN DD * 

.54,. 27,. 27,2.56, 1.44,l. OO(:l, 45,l. 1 

-6, -5, -4, -3, -2, - 1 ,O , 1 ,2, .3,4,5 /* 

//