WAVEGIDE—An interactive waveguide program

WAVEGIDE

An interactive waveguide program

GRANVILLE SEWELL University o f Texas at El Paso, Mathematics Department, El Paso, Texas 79968, USA

SRBA CVETKOVIC University o f Surrey, Department o f Electronic and Electrical Engineering, Guildford, Surrey GU2 5XH, UK

WAVEGIDE is an interactive program which, when used in conjunction with IMSL's P D E / P R O T R A N , can solve waveguide problems on arbitrary cross-sections, with variable permittivity and permeability. The user only has to answer a series of basic questions about his problem, posed in engineering terms, and all propagating modes in a specified interval are calculated, and the corresponding electric and magnetic fields are plotted.

be piecewise constants, as required when the waveguide is a composite of several materials. The partial differential equations satisfied by the electric (Ex, Ey, Ez) and magnetic (Hx, Hy, Hz) components of a propagating wave in such a (non-homogeneous) waveguide structure are (see z Section 8.2i:

OEy OE~ -- + j~Hz ax Oy aHy

aH~

Ox

Oy

= 0

(1)

j~eEz = 0

(2)

where 1. INTRODUCTION

Ex= - j 13 -~x + CO~

P D E / P R O T R A N t, the IMSL partial differential equation solver (developed by the first author), is a finite element program which solves a wide range of PDEs in general two-dimensional regions. Because of its general applicability, it can be used to solve many difficult waveguide problems. Although P D E / P R O T R A N allows the user to input the problem description in a high-level, keyword-oriented, language, it still requires a substantial effort, and considerable experience in the use of P D E / P R O T R A N , to solve many waveguide problems using this package. In order to make the impressive power of P D E / P R O T R A N easily accessible to engineers who want to solve waveguide problems, the authors have developed an interactive driver, called WAVEGIDE, which facilitates the solution of a broad class of waveguide problems. W A V E G I D E asks the user a few simple questions about his/her problem, and generates the appropriate P D E / P R O T R A N input program based on the user's responses.

E~=j(-/3 0E~+'~ ~ff-x~)/(°~c- ~)

2. THE WAVEGUIDE PROBLEM WAVEGIDE solves a cylindrical waveguide problem of arbitrary cross-section, with conducting boundary. The permittivity, e(x, y), and the permeability, #(x, y), are allowed to vary across the cross-section, and may even Accepted April 1989. Discussion closes April 1990.

:~ 1989 Computational Mechanics Publications

-

/(~2p, e - / 3 2 )

ay

H~=j

Hy= - j

-/3~xx +~oc

/(~o2.e-~2)

[~-~y +we

/(w2#8-/32 )

and co = 2~rf ( f = frequency) a n d / 3 j is the propagation constant (which is purely imaginary for propagating waves). When the expressions for Ex, Ey, Hx, Hy are inserted into equations (1) and (2), there result two partial differential equations for the two unknown functions Ez and Hz. The boundary conditions on the conducting boundary are

Ez = 0

and

OHz/On= 0

(3)

For a given ~0, the problem is to determine values of/3 for which nonzero solutions (eigenfunctions) of (1) and (2), with boundary conditions (3), exist. When e and # are constants, equations (1) and (2) reduce to the simple forms: 02Hz + 02Hz OX2 ~ + (co2/zC-/32)Hz = 0

(4)

OZEz + 02E-'---~z+ (w2/~6" - / 3 2 ) E z = 0

(5)

OX 2

Oy 2

Adv. Eng. Software, 1989, Vol. 11, No. 4

169

In this case, the eigenvalues {3 separate out into the eigenvalues of (4, with OHz/On = 0 on the boundary), called transverse electric (TE) modes, and those of (5, with E z = 0 on the boundary), called transverse magnetic (TM) modes. However, in the general case where e ( x , y ) and /x(x,y) are variable (such as is usually the case for optical waveguides), we must solve (1) and (2) simultaneously, and the eigenvalues are called hybrid modes. Since P D E / P R O T R A N uses a Galerkin finite element method, it can handle composite media, where e and/or t* are discontinuous at the interfaces between different materials, without the need for explicit interface conditions.

nonzero right-hand side, returning a function value. g(5), which represents the inverse of a norm of the solution. W A V E G I D E tries to locate all minima of g(/3) in (BMIN, BMAX) using the following procedure: g(~) is evaluated at a user-specified number (NINTER) of uniformly-spaced values t5i between BMIN and BMAX. When an index i is found such that g(/Si) is less than both g([~i-1) and g(15i+1), it searches for a minimum of g in (/3i-1,~i+1), using a public domain subroutine F M I N 3 to locate the minimum accurately. In our experience minima which do not correspond to eigenvalues are very rare, and they betray themselves by their shallowness (g(5) does not go to zero near these minima). Thus W A V E G I D E is usually able to discard these pseudo-eigenvalues. The user may, alternatively, supply a value of/3 and limits on f. In this case W A V E G I D E will find values of the frequency which correspond to propagating modes. If a value of /3 = 0 is supplied, the frequencies thus found will be the cut-off frequencies.

3. SOLVING THE N O N L I N E A R EIGENVALUE PROBLEM P D E / P R O T R A N can solve eigenvalue problems such as (4) or (5) by reducing them to matrix eigenvalue problems of the form A x = XBx, and by using the inverse power method to solve the discretized matrix problem. In the more general case (1-2) the eigenvalue/3 appears nonlinearly, and the problem cannot be reduced to a linear matrix eigenvalue problem. However, eigenvalues can be found using the following idea. If we introduce a " r a n d o m " non-zero righthand side to (1) and (2) and vary t3 over a range of values (BMIN, BMAX), as /3 approaches an eigenvalue the solution becomes infinitely large in magnitude. Using this fact, we can attempt to locate all eigenvalues in a given t3-range. If we are very unlucky in our choice of right hand sides, we could conceivably miss some eigenvalues, but the probability of such an unlucky choice is zero, and this has not proven to be a problem. Thus W A V E G I D E asks the user to supply a frequency, f , (co = 2~rf) and limits (BMIN, BMAX) on/3. For a given value of 5, W A V E G I D E solves (1-2) with a

****

Welcome

to

4. EXAMPLES If the user desires to work through a prepared example before trying a problem of his own, W A V E G I D E will provide answers, corresponding to one of two prepared examples, for each of its questions. The first example is a simple circular tunnel waveguide problem, while the second is a more complex optical waveguide problem, involving two dielectrics. Figure 1 shows a simulated W A V E G I D E interactive session, in which the user asks to work through the first example. Lines beginning with a ? are user responses. Notice that the user is only asked to supply basic physical parameters, and may not even be aware that he is solving a system of partial differential equations. After finishing the interactive session, the user runs the P D E / P R O T R A N program generated by W A V E G I D E ,

WAVEGIDE

****

T h e W A V E G I D E p r o g r a m p e r f o r m s a modal a n a l y s i s of a w a v e g u i d e w i t h arbitrary cross-section R, a n d b o u n d e d w i t h p e r f e c t l y c o n d u c t i n g metallic walls. The permittivity (EPS) and p e r m e a b i l i t y (AMU) of the g u i d i n g m e d i a a r e real s c a l a r q u a n t i t i e s w h i c h c a n v a r y o v e r the cross-section. The user may specify the frequency (F), and al 1 propagating modes with propagation constants (BETA) in a specified BETA-range will (hopefully) be found. Alternatively, the user may specify the propagation constant (BETA), and all propagating modes with frequencies (F) in a specified F-range will be found. In both cases, a full profile of each mode is also displayed graphically. Each question asked by WAVEGIDE has a The long form is recommended for first to see the long form of each question? ***** Enter YES o r NO ? YES

If t h i s through example 0 1 2 ***** Enter 71

170

is an and if if if an

short time

form and a long form. users. Would you like

y o u r f i r s t t i m e u s i n g W A V E G I D E , you m a y w a n t to w o r k e x a m p l e p r o b l e m b e f o r e t r y i n g o n e of y o u r own. A simple a more complex one have been prepared. Enter you do not want to see either example you want to work through the simple example problem you want to work through the more complex example integer value in the range 0 to 2

Adv. Eng. Software, 1989, Vol. 11, No. 4

EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ The example problem is to find the cutoff frequency for the TEll mode in a hollow circular tunnel with perfectly conducting boundary and of radius 4.8 meters. The TEll mode has the lowest cutoff frequency, and it can be shown analytically that its cutoff frequency is 1 8 . 2 9 E6 c y c l e s / s e c (18.29 MHz).

$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

$ $ $ $ $

$

$ After each question, the correct answer for this example will be $ indicated. $ $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ If you specify a value for the frequency (F), values of the propagation constant (BETA) for a] ] propagating modes in a specified r a n g e B M I N < BETA < BMAX w i l l be calcu|ated.

$ $ $ $ $

$

$ $ $

If you specify a value for the propagation constant (BETA), values of the frequency (F) for all modes in a specified range F M I N < F < FMAX w i l l be calculated. Thus if you want the cutoff frequencies, specify BETA=O. Notice that if the waveguide is once the cutoff frequency FC i s propagation constant (BETA) and functions of frequency F: if

F>FC

if

F
BETA ALPHA BETA ALPHA

= = = =

homogeneous (EPS a n d AMU c o n s t a n t ) , known (for a given mode), the attenuation (ALPHA) are known as

2*PI*SQRT(AMU*EPS)*SQRT(F**2-FC**2) 0 0 2*PI*SQRT(AMU*EPS)*SQRT(FC**2-F**2).

In the case of a homogeneous waveguide with losses (EPS = E R - E I * I , where I=SQRT(-1)), the attenuation and propagation constants are still known as functions o f F , o n c e FC i s k n o w n . If FC i s t h e cutoff frequency found w h e n EPS i s s e t t o ER, t h e n ALPHA a n d BETA --for complex EPS--can be calculated from: ALPHA

+ BETA*I

= 2*PI*SQRT(AMU*ER)*SQRT(FC**2-F**2

+ F**2*EI/ER*I)

Do y o u w a n t t o s p e c i f y a frequency (F)? $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ enter: NO $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ ***** Enter YES o r NO ? NO Then, enter the the frequencies

propagation calculated

constant (BETA). wil] be the cutoff

Suggested units f o r BETA a r e $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ $ enter: 0.0 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ ***** Enter constant o r FORTRAN ?

If you set frequencies.

$

BETA=O,

1/METERS. EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ expression *********************

0.0

The interval correspond Enter

to

a value

(FMIN,FMAX) propagating for

FMIN.

will be waves. FMIN

searched

must

be

for

positive

values if

of

BETA

F which is

zero.

Suggested units for F M I N a r e HERTZ ( = C Y C L E S / S E C O N D ) . $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ Since we a l r e a d y know that the true cutoff frequency for the $ TEll mode is 18.29 E6, let us search the interval (15.E6,20.E6). $ enter: 1 5 . E6 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ ***** Enter constant o r FORTRAN e x p r e s s i o n ********************* ? 1 5 . E6 Now e n t e r

a value

for

$ $ $

FMAX.

Suggested units for FMAX a r e HERTZ ( = C Y C L E S / S E C O N D ) . $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ enter: 2 0 . E6 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ ***** Enter constant o r FORTRAN e x p r e s s i o n ********************* ? 2 0 . E6

$

Adv. Eng. Software, 1989, Vol. 11, No. 4

171

NINTER v a l u e s of F will be c h e c k e d in the interval (FMIN,FMAX). The larger NINTER i s c h o s e n , the more ]ikely the program is to find all modes in this interval, but also the more computer time will be r e q u i r e d . Enter a value for NINTER. I f no i d e a , enter 20. $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ enter: 10 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ ***** Enter constant o r FORTRAN e x p r e s s i o n ********************* ? 10 Do y o u w a n t t o g e n e r a t e plots of the Ez,Hz,(Ex,Ey) and (Hx,Hy) fields? It is important to look at the eigenfunction plots before accepting the corresponding eigenvalue as valid. $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ enter: NO $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ ***** E n t e r YES o r NO ? NO

$

$

(To conserve space, we o m i t h e r e t h e p o r t i o n of the interactive session which is used to determine the initial and final triangulations, and the boundary of the cross-section.)

Now y o u m u s t w r i t e a FORTRAN s e g m e n t t o d e f i n e the permittivity (EPS) and the permeability (AMU). T h e s e may b e d e f i n e d in terms of X a n d Y. T h e y may a l s o reference the constants EPSO ( = 8 . 8 5 4 E-12) a n d AMUO ( = 1 2 . 5 6 6 E-7), which hold the permittivity and permeability in vacuum, in units o f FARADS/METER a n d H E N R I E S / M E T E R . Suggested Suggested

units units

for for

EPS a r e AMU a r e

EPS,AMU may a l s o b e d e f i n e d KTRI, if the non-rectangular triangulation.

FARADS/METER. HENRIES/METER. in terms of option was

the initial triangle used to define the

number initial

T h e FORTRAN s t a t e m e n t s must obey the usual rules of FORTRAN--for example, statements must begin i n c o l u m n 7. $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ Since the tunnel contains air, enter, when p r o m p t e d : $ AMU = AMUO $ EPS = EPSO $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ ***** Enter first line o f FORTRAN * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ? AMU = AMUO ***** ?

Enter another EPS = EPSO

***** Enter ? [RETURN] The for and

another

line

(press

'return'

if

no

more)

*******************

line

(press

'return"

if

no

more)

*******************

solution is s a v e d postprocessing by n y in t h e f o r m :

$ $ $

o n a n nx by ny r e c t a n g u l a r g r i d of p o i n t s , the plotting routines. E n t e r v a l u e s f o r nx

nx,ny nx and ny must e a c h be b e t w e e n 3 and 51, i n c l u s i v e . If you have no i d e a , enter 21,21. $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ enter: 21,21 $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ ***** Enter answer as specified ******************************** ?

21,21

$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ Solution: You should find one cutoff frequency in t h e r a n g e $ (FMIN,FMAX), namely, F = 1 8 . 2 9 E6. $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ EXAMPLE $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

***** *****

Input program is completed, ***** and should be e x e c u t e d using***** I M S L ' s PDE/PROTRAN 1 . 0 . *****

***** ***************************************

Figure 1 172

$

Adv. Eng. Software, 1989, Vol. 11, No. 4

$ $

and the result is that one cutoff frequency is found. f = 18.30 E6. The second example involves a 2/~m by 2 #m square optical waveguide (glass), of refractive index n = 1.55,

•

ET.

U3

FOIl F: O . 3 5 2 6 9 E ~ ' I S . B U T A =

embedded in an infinite substrate (silica) o f refractive index n = 1.50. The wavelength in a vacuum is given to be X = 0.85/~m, so the frequency is f = c/X = 3.5269 E l 4 cycles/second. Since the permeability is generally taken

0.11339£,08

c~

KEY

-2.07E-25 -1 .61E-25

_=o

-.1. 1 5 E - 2 5

co

-6.9lE-26

•

'5" I(

-2

z

• 2 30E-26

)"

30E-26

• 6 9lC-26

o t:3

.1

c~

1SE-2S

I

• 1 61E-25

• 2 CITE-2S

o b3

I

-5.00

-3.00

-~.o0

~.oo

5.0(3

3.00

Figure 2A FOR F= 0 . 3 5 2 6 9 E * l ~ . B E T £ =

O.t1339E*OB

~EY

0

-8.28E-28

1

-6,qLIE-2B

2

- q . BoE-28

3

-2.76E-28

q

-9.20E-29

5

.9.20E-29

6

.2.76E-28

7

*LI. 6 0 E - 2 B

8

*6.qqE-28

9

*8.28E-2B

I 0

Co

I

I

-5,00

-3.00 -1.00 X (,,lO~o -6)

I .00

3.00

S.O0

Figure 2B Adv. Eng. Software, 1989, Vol. 11, No. 4

173

qzJ

Ln

(CX,EY3 .,F= 0 . 3 5 2 6 9 E * l S . B E l A =

i

0.I'~339Ce08

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

'

.

.

.

.

.

.

~

.

i

f

:

1

.

.

.

K/// I.Z/z. i ........

I

.

.

.

•

'

,

•

L,

. . . . . . .

~.SSE-Zq SCALE

I" I ' ' ' "

c~

"7

.

I

.

.

.

.

.

.

.

.

.

.

.

,

.

.

.

.

.

.

.

,

.

.

.

.

,

.

.

.

.

.

.

.

.

.

i

.

~7 -S.O0

-3.00

-1.00 X

(-10--

1.00

3.00

S,O0

-6)

Figure 2C

tMX,h~);F=

O,35269E+lS,BETN=

0.11339E,08

! .....................................................................

!

.

.

.

.

.

.

.

.

i

.

.

.

.

.

.

.

.

i

.

.

.

.

.

.

.

.

.

,

.

.

•

,

.

,

P= ,4 .

.

.

.

.

.

.

.

.

.

.

.

.

.

•

.

.

.

.

.

.

.

.

.

.

.

.

.

.

o

. . . . . .

, ~ \ \ \

. . . . . . .

,

. . . . . . . . .

.

.

.

.

.

.

.

. . . . . . . .

,,

.....

,,,,il........ \ ~ \ \

,

"

. . . . .

\\~

I.

OSE-~6

,, \ \ \

SCALE

I c3

i

i

-S.O0

-3.00

-1.00

X

(-lO-*

1.00

-6)

Figure 2D

174

Adv. Eng. Software, 1989, Vol. 11, No. 4

3.00

-~ . 0 0

CHARACTER'I L 1 N E ( 8 0 ) OPEN ( 2 , F I L E = ' M E S S O L D ' ) OPEN ( 3 , F I L E = ' M E S S A G E ' , A C C E S S = ' O ] R E C T ' , R E C L = 8 0 , S T A T L I S = FORH='FORMATTED') DO 10 ] R E C = I , I O 0 0 READ ( 2 , 5 , E N D = 9 9 ) LINE FORMAT ( 8 0 A l ) 5 WRITE ( 3 . 5 , R E C = ] R E C ) LINE ]0 CONTINUE 9 9 STOP END

'NEW',

Figure 3 to be/+o for optical materials, the permittivity is proportional to n 2, e = con 2, and thus has different values in the two dielectrics. This only requires that when W A V E G I D E asks for e and t~, we write a short F O R T R A N 7 7 segment which defines e to be eol.552 when (x, y) is in the glass guide, and col .502 when (x, y) is in the silica substrate. We ask W A V E G I D E to find all propagating modes [3 which correspond to "effective refraction" indices (t3=27rfne~]c) in the range 1.530~
modifications to the O P E N statements and logical unit specifications, which are clearly marked near the beginning of the program. 2. Set the characteristics of the MESSAGE file so that it can be read by W A V E G I D E . F 7 7 in direct access mode, with record length = 80. The easiest way to do this varies from computer to computer, but one technique which will work in any environment is to read the original MESSAGE file in sequential access m o d e and write out a new copy in direct access mode. On an IBM 4341 computer the F O R T R A N 7 7 program shown in Fig. 3 accomplished this task. The original MESSAGE file was renamed to MESSOLD before execution of this program. The last two steps are performed each time a problem is solved using W A V E G I D E : 3. Load W A V E G I D E . O B and execute the resulting load module. After the user has answered all of the interactive questions, a P D E / P R O T R A N input program will have been created, in the file P D E P R O . The user m a y want to edit this file to correct any errors made during the interactive session. 4. Execute P D E / P R O T R A N on the input program PDEPRO.

REFERENCES 1 Sewell, G. Analysis of a Finite Element Method: PDE/PROTRAN, Springer-Verlag, 1985 2 Ramo, S., Whinnery, .1. R. and Van Duzer, T. Fields and Waves in Communication Electronics, John Wiley and Sons (second edition) 1984 3 Forsythe, G., Malcolm, M. and Moler C. Computer Methods for Mathematical Computations, Prentice-Hall, 1977

Adv. Eng. Software, 1989, I.Iol. 11, No. 4

175