PECSxxx.PRO: PEC sphere

3D MMP, Electrodynamics, simple axisymmetric model

Ch. Hafner, Computational Optics Group, IFH, ETH, 8092 Zurich, Switzerland

Back to MaX-1, MaX-1 Examples Overview

Download ZIP file containing the MaX-1 projects described here

A perfectly conducting sphere is the simplest 3D object for MMP. Therefore, it is reasonably to study it carefully. In the following, some important hints are given, but it is assumed that the reader is already familiar with the MaX-1 software. Note that Version 2 is required and that an OpenGL Graphics window should be used for easily inspecting the models and results. Make sure that you have downloaded the mos recent version of MaX-1 as well as the GLUT32.DLL. When you start MaX-1, answer the question "Display OpenGL window?" in the affirmative. When the corresponding window pops up, its background color should be white, otherwise, your graphic driver might not correctly support OpenGL.

Defining 3D objects

In general, a 3D object consists of a boundary in 3D space and some appropriate 3D expansions. In most cases, y boundary in 3D space can be generated from a 2D boundary in the xy plane with an appropriate transformation. The most useful transformations are defined in the 3D Objects dialog. For generating a sphere, you may use the transformation that rotates a 2D boundary around the y axis. Since this transformation generates a torus in general, it is called Torus in the 3D Objects dialog.

The simplest way to generate a sphere is to rotate a circle with center at x=0, y=0 around the y axis with an angle of 180 degrees. When you want to take symmetry planes into account, you only discretize a quarter of the sphere. For doing this, rotating an entire circle is not reasonable. Therefore, we first generate a semi-circle in the right half of the xy plane with center at x=0, y=0 and radius 0.5m. Remember that this is a C-polygon with 4 corners. Once you have defined a semi-circular 2D boundary, you also specify the left and right domain numbers (left is the inner material of the sphere, i.e., domain number 0, right is free-space, i.e., domain number 1). You may use default values for all other parameters of this 2D boundary.

Once you have the 2D boundary that shall be used for generating the sphere, you can open the 3D Objects dialog. Select Torus and specify Min. boundary 1 as well as Max. boundary 1 (You only have defined one single boundary. For more complicated objects you might start with a sequence of 2D boundaries that are used all together for defining the object!). Note that you also may specify Min. expansion and Max. expansion. These parameters are used for generating sets of 3D expansions from sets of 2D expansions. We will manually place the expansions in the following. Therefore, you may set the Min. expansion and Max. expansion values 0. Since you have specified a semi circle, you must rotate it around the y axis with an angle of 360 degrees: Set Minimum angle 0 and Sector angle 360.

When MaX-1 generates 3D matching points on the surface of an object, it first generates the 2D matching points along the corresponding 2D boundaries. Afterwards, it applies the same transformation (rotation around the y axis) to each 2D matching point. Finally, it associates rectangles to all 3D matching points. These rectangles are tangential to the surface. The aspect ratio of the sides of the rectangles are specified in the Mat.Pt.aspect ratio box. Select the ratio 1 for obtaining "square" matching points.

When you want to make sure that MaX-1 generates the desired object, you best select the following settings in the 3D Objects dialog: Specify the desired Object # (0 for all objects), check the Matching pts box, don't check the Field on grid box, and select Transparent representation (this ist quickest). Close the dialog and open the OpenGL dialog by right clicking the OpenGL Graphics window. Select draw 3D objects. You now should see your sphere and the 3 axes of the global coordinate system. Take time to play around with the OpenGL visualization of your 3D object. Try different settings in the 3D Objects dialog. Use OpenGL to move, rotate, zoom, resize the window, etc. Select Tools -> OpenGL to open the OpenGL window dialog, where you manually can define different settings of the OpenGL Graphics window. Note that the electromagnetic field can also be visualized in this window (as soon as it has been computed). The outfit of the field is affected by the settings in the Field dialog.

Note that the rotations an translations in the OpenGL Graphics window do not modify the location of the 3D objects in the 3D space. What is modified is to point from where you look at the objects. When you want to move a 3D object to another place, press the Location… button in the 3D Objects dialog and modify the data in the dialog that will pop up.

Plane wave Excitation

Assume that a plane wave is incident on the PEC sphere. You specify the plane wave in the Expansions dialog as you specify a plane wave for 2D computations, but now, you must specify its location and orientation in 3D space. For doing this, press the 3D location… button in the Expansions dialog. Since the Origin of the plane wave has no influence on the results, set it at (0,0,0). The X direction is a vector that specifies the direction of the electric field (when the E wave filter is set in the Expansion dialog), whereas the Y direction is a vector that specifies the direction of the magnetic field. Let (1,0,0) be the X direction and (0,0,-1) the Y direction. The Z direction is then (0,1,0), i.e., the plane wave propagates along the y axis.

Brute-force solution: PECS000.PRO

You already have defined the boundary of the PEC sphere and the excitation. What is missing, is an appropriate expansion for the scattered field in domain 1, i.e., outside the sphere. As in Mie scattering, you best set a single 3D multipole expansion in the center of the sphere. You can easily do this manually in the Expansions dialog. Note that the plane wave excitation should be the last expansion. Therefore, you insert an 3D multipole as first expansion. As for the plane wave, you press the 3D Location… button for defining the location and orientation of the multipole. It is reasonable to use the same X, Y, and Z direction vectors as for the plane waves although an arbitrary orientation might be used for the brute-force computation.

The main problem is now to define the orders and degrees of the multipole. Since we have defined the Z direction of the multipole by (0,1,0) in the global y direction, the degree of the multipole indicates the variation of the field around the y axis (local phi direction of the spherical coordinates). When the frequency is, for example, 600MHz, the wavelength is approximately 0.5m. In this case, the circumference of the sphere is approximately 6.3 wavelengths. Thus, a Maximum order and Maximum degree of 20 should be sufficient for obtaining good results (especially since the surface is very smooth and because the incident wave is simple). Specify these values in the corresponding box. As for 2D multipoles, you can set Minimum order 0. For 3D, you also may specify Minimum degree 0. When you select expansion number 2 (the plane wave) and then expansion number 1 (the multipole), you will see that MaX-1 has replaced Minimum order 0 by Minimum order 1, because order 0 multipoles do not exist in 3D electrodynamics. For given minimum and maximum orders and degree, depending on whether electric (E), magnetic (H), or both (HE) multipoles are used you obtain a certain number of parameters for each multipoles. You should enter this number in the Max box of the Parameters group of the Expansion dialog. Since you probably don't want to evaluate this number, you place a relatively large number, for example, 1000, in the Max box. When you now press the All0 button, MaX-1 will reduce this value to the correct one (880). If you would have entered a smaller value, the value would remain unchanged and MaX-1 would not use all of the desired degrees and orders. Because of memory allocation problems, you should not use extremely high numbers in the Max box. Furthermore, when the number of parameters is high, the computation time will be long.

Once you know the number of parameters to be computed by MMP, you also can estimate the number of matching points that is required. Since the 3D matching points are generated from the 2D matching points, the values in the Matching point definition group of the MMP dialog affect the number of the 3D matching point. Even when you know how many 2D matching points you have along the 2D boundary, it is difficult to estimate the number of 3D matching points. Therefore, you best press the Points!! button in the MMP dialog. MaX-1 will then generate the 2D and 3D matching points and display the corresponding numbers in the MMP dialog. On a PC surface, you have 2 independent boundary conditions for each matching points. Therefore, you should have at least 440 3D matching points. If you have not enough or much too many matching points, modiy the values in the Matching point definition group of the MMP dialog and press Points!! again.

When you have finished the model, proceed as in 2D computations. Press the Sol+err!! button in the MMP dialog. This will take some time when your machine is not very fast! Inspect the output in the Info window. The results should be very accurate. You now can compute the field in the xy plane or in an arbitrary plane as in 2D and visualize it. For inspecting the results, you can also compute the field on the surface of the sphere (select Derived field in the 2D Objects dialog) and use the OpenGL Graphics window for drawing. Figure 1 shows a typical OpenGL representation.

Figure 1: Time average of the Poynting vector field. Brute-force computation.

Symmetry planes: PECS001.PRO

Since you have obtained highly accurate results, you may reduce the computation time by reducing the maximum orders and degrees of the multipole. Because of the fast convergence, this will also strongly reduce the accuracy. An alternative way to reduce the computation time (and memory) without any reduction of the accuracy is taking symmetries into account. Obviously, the incident wave has XY plane symmetry number 1 and YZ plane symmetry number 2. When you specify these numbers in the Project dialog, MaX-1 will automatically reduce the number of parameters of the multipole from 880 to one quarter, i.e., 220. You now can also reduce the discretization of the surface of the sphere to one quarter by reducing the Sector angle in the 3D Objects dialog form 360 to 90 degrees. Thus the MMP matrix size will be reduced by a factor of 16 and the computation time of this matrix should be approximately 64 times shorter! After this, you can proceed as before.

As you can see from Figure 2, MaX-1 now only evaluates a quarter of the sphere, but you have all information of the field because of the symmetry.

Figure 2: Time average of the Poynting vector field. Brute-force computation with 2 symmetry planes, only one quarter is modeled.

Reduce multipole degrees: PECS002.PRO

When you inspect the values of the parameters of the multipole in the previous project, you see that many of them are very small, i.e., numerically zero. Since the incident plane wave travels in y direction, its angular dependence in phi direction (around the y axis) is essentially cos(phi). When you have oriented the Z direction of the 3D multipole also along the y axis, this means that only the multipole degree 1 is excited because only this degree has cos(phi) dependence. Therefore, you can reduce Maximum degree from 20 to 1. As soon as you have done this, you have reduced the number of parameters from 220 to 40. This obviously also drastically reduced both memory and computation time without reducing the accuracy of the result.

Take rotational symmetry into account, 2 slices: PECS003.PRO

When you reduced the maximum degree of the multipole, you implicitly took the axisymmetry of the problem into account. It is now natural to do this also for the matching points. When you know that you body is axisymmetric, it is sufficient to discretize an infinitesimal slice in the xy plane. You therefore can define a 3D object with a very small Sector angle in the 3D Objects dialog. If the sector is small enough, MaX-1 will generate only one 3D matching point for each 2D matching point, i.e., the problem is almost 2D. By properly setting the Minimum angle (for example, –0.5) and the Sector angle (for example, 1.0) you obtain all 3D matching points at those positions, where also the 2D matching points are, i.e., in the xy plane. When you do this, the following problem is encountered: Because of the polarization of the incident plane wave, the Ez components of the field are zero on all 3D matching points. In order to match also these components, you should introduce a second slice of 3D matching points in the yz plane. For doing this, you add a second Torus in the 3D Objects dialog with the same parameters as the first one, but with a different Minimum angle (for example, 89.5).

When you proceed as before, you should obtain the solution shown in Figure 3. Note that the memory and computation time have been reduced once more. The computation time should now be so short that the display of the progress in the Info and movie directives box takes more time than the MMP solution itself. In order to save computation time, you can turn the display off.

Figure 3: Time average of the Poynting vector field. Computation with two thin slices.

2 slices improved, show the field on the entire sphere and on an additional plane: PECS004.PRO

Since the Ez boundary condition on the first slice are automatically fulfilled, you may turn these boundary conditions off by selecting Specific boundary conditions instead of the Usual boundary conditions in the Boundary dialog. Check only the Et box. Since you have used the 2D boundary also for generating the second slice, you would also impose the same boundary condition on the corresponding matching points in the yz plane. Because of the symmetry of the incident wave, Et is zero in the yz plane. You have introduced the second slice for matching the Ez component that is zero on the first slice. Therefore, you now should copy the 2D boundary in the Boundary dialog and replace Et boundary condition by Ez boundary condition. Then, you should set Min. and Max. boundary of the second slice, i.e., of the second object in the 3D Objects dialog equal to 2. By doing this, you reduce the number of equations for each 3D matching point from three to one. This will reduce the memory and computation time once more. Provided that you have enough matching points you should obtain almost the same accuracy as before.

When you want to compare the field you now have on the sphere with the brute-force solution, you want that MaX-1 computes the field on the entire sphere rather than only on the two thin slices. Therefore, you could start with the two slices for computing the parameters of the multipoles. As soon as these parameters are known, you can remove the second slice and increase the Sector angle of the first slice to 360 degrees. After this, you let MaX-1 draw the field on the first object, which is now an entire sphere. Although this provides the desired output, it is inconvenient when you want to modify some part of the model and repeat the entire procedure. Therefore, you best define an entire sphere as a third 3D object in addition to the two slices. Since you don't want that the corresponding matching points are used for the computation of the parameters of the multipole, you should do the following: 1) In the Boundary dialog, you copy boundary number two. 2) Set the Weight of the new, third boundary equal to 0. This makes sure that the corresponding 2D matching points and the resulting 3D matching points will be ignored during the MMP matrix computation. 3) Copy the second object in the 3D Objects dialog and modify the data of the new, third object in such a way that that it is an entire sphere produced with the data of the third 2D boundary.

Maybe, you want to visualize the field not only in the xy, but also in the yz (or in some other) plane. Of course, you can modify the field plane in the Field dialog as usual, but then, you will only see the field in the yz plane. When you want to visualize the field in both the xy and the yz plane at the same time, you can generate a rectangular section of the yz plane as a new, forth object in the 3D Objects dialog. Similarly to the third object, this object shall not be used for the MMP matrix computation, i.e., it shall be based on a 2D boundary with Weight 0. The simples way is to define a Rectangle in the 3D Objects dialog. Note that you must also provide a 2D boundary for generating a rectangle although all geometric data of a rectangle are defined in the 3D Objects dialog. The reason for this is that additional data (type of the boundary conditions, weight, domain numbers on both sides of the rectangle, etc.) are required. Such data are defined for 2D boundaries. Since the third 2D boundary has the weight 0, you might use it for generating the fourth, rectangular object. When you do this, the domain numbers on both sides of the object will be used for computing the domain numbers of the field points that are sufficiently close. As a consequence, you will obtain wrong domain numbers in some field points. To avoid this, you must define a fourth, dummy 2D boundary with Weight 0 and domain numbers –1 on both the left and right hand side. When you define a straight line along the y axis as the fourth 2D boundary, you may also us it for generating a Cylinder as the fourth object in the 3D Object dialog.

If you wish to visualize the field also in the xz plane, you can copy the fourth object and modify the location and orientation of the resulting fifth object by pressing the Location… button in the 3D Objects dialog.

Take rotational symmetry into account, 1 slice: PECS005.PRO

We have modeled the boundary of the sphere with two slices in the two symmetry planes because not some components of the field vanish in these planes. As an alternative, we can also define a single slice between these planes (for example, a Torus with Minimum angle 45 degrees) and impose both Ez and Et boundary conditions along the corresponding 2D boundary. Try and compare the results!

Dipole Excitation

When you replace the plane wave excitation by a dipole in front of the sphere, the line from the dipole to the center of the sphere is a symmetry axis. Let this be the y axis, i.e., the natural axis for generating a Torus with MaX-1. Assume that the dipole is at some position yD on the y axis. It then may have two principal orientations: 1) perpendicular to the y axis and 2) parallel to the y axis.

Orientation perpendicular to axis, sufficient distance: PECS014.PRO

When the dipole is perpendicular to the y axis, you can select its direction in such a way that you obtain the same symmetry of the field as for the plane wave excitation we have considered before. In this case, you can essentially use the same models as before. For example, start with PECS004.PRO. All you have to do is to replace the second expansion, i.e., the plane wave by a dipole, i.e., a multipole with order 1 and degree 1 (or 0 – depending on its orientation in the 3D space). Place the dipole at a sufficient distance from the sphere, for example, at y=1m, when the radius of the sphere is 0.5m. Note that MaX-1 checks the data of the expansion data when you press one of the buttons in the Expansion dialog or when you change the number of the Current expansion. When you have a wrong orientation of the dipole with respect to the symmetry planes defined in the Project dialog, the Expansion check box will pop up and tell you that the expansion has no parameters. If this happens, modify the orientation (or the degree) of the dipole. When the dipole is properly defined, you can proceed exactly as for the plane wave excitation.

You will notice that the accuracy of the results is considerably reduced (compared with the plane wave excitation). This effect gets stronger when you move the dipole closer to the sphere. The reason is simple: the illumination of the sphere is less uniform than in the plane wave case. However, when the distance of the dipole is not short, you obtain accurate results.

Figure 4: Time average of the Poynting vector field. Dipole excitation with dipole at y=-2R , dipole perpendicular to y axis.

Orientation perpendicular to axis, short distance: PEC015.PRO

When the distance of the dipole from the sphere is short enough, you will obtain inaccurate results when you model the scattered field with a multipole in the center of the sphere even when you use quite high orders. In order to improve the model, remember that a PEC sphere acts like a mirror. When you observe a light source in front of a spherical mirror, you see it in the mirror point. The mirror point is at yM=R*R/yS, when R denotes the radius of the sphere and yS the location of the source on the y axis. Therefore, it is reasonable to insert a second multipole in the mirror point (in addition to the first multipole in the center of the sphere). Depending on the maximum orders of the two multipoles you can drastically improve the accuracy of the results without increasing the total number of parameters. After this simple modification, you can proceed as before. Note that you now have two multipoles instead of a single one that is sufficient from the theoretical point of view, i.e., you have a Multiple multipole expansion.

The single multipole expansion with the multipole in the center has an important benefit: The corresponding basis functions are orthogonal on the sphere. Therefore, it would be possible to compute the parameters without solving a matrix equation. This can be done in specialized Mie codes, but not in a general purpose code such as MaX-1. Since the basis of the brute force solution is orthogonal, the condition number of the corresponding MMP matrix is small (As soon as you have computed the rectangular MMP matrix either with the GUR or CG matrix solver, you may let MaX-1 compute the condition number by pressing the Get condition!! Button in the MMP dialog. Note that this computation is both time- and memory consuming for big matrices because a singular value decomposition is done.). As a consequence, the CG matrix solver of MaX-1 will rapidly converge when the square scaling (Matrix scaling 2 in the Matrix solver group of the MMP dialog) is used, i.e., you can save computation time by replacing the standard Givens matrix solver by the CG solver. However, with the appropriate axisymmetric modeling, the MMP matrix computation is fast anyway.

When you add a second multipole, the condition number of the matrix is drastically increased, but with an appropriate matrix solver, you may obtain more accurate results. Try the GUR, GUT, and CG solver! The reason for the high condition number is the fact that the two multipoles are not independent.

When the dipole and the multipole in the mirror point are close to the boundary, it might be necessary to increase the matching point density near the dipole. In 2D MMP computations, this is automatically done when the Min. Overdetermination factor in the MMP dialog is sufficiently big (usually 1.0 or more). Since the generation of the 2D matching points takes the maximum order of the 2D multipoles into account, you may add a dummy 2D multipole at the location of the 3D multipoles. To make sure that these 2D multipoles are dummy, i.e., are not used for the MMP matrix computation, you specify a negative value in the Object # box in the Expansion dialog. Note that the dummy 2D multipoles should have the same maximum orders as the corresponding 3D multipoles. Incidentally, you might also use such dummy 2D multipoles for automatically generating appropriate 3D expansions.

Figure 5: Time average of the Poynting vector field. Dipole excitation with dipole at y=-1.2R , dipole perpendicular to y axis.

Orientation along Axis: PEC016.PRO

When you rotate the orientation of the dipole in such a way that it is parallel to the y axis, the symmetry of the field is different. To take this into account, you must modify the YZ plane symmetry number in the Project dialog. The correct symmetry number now is 1. Furthermore, the angular dependence of the field is now constant, i.e., cos(0*phi) rather than cos (1*phi). Therefore, the minimum and maximum degree of the multipoles should now be 0 rather than 1. Modify this in the Expansion dialog for both multipoles. After these three modifications, you may proceed as before.

Figure 6 shows the resulting field. For getting experience, try different positions. Note that you should also adapt the location of the 3D multipole and the 2D dummy multipole in the mirror point.

Figure 6: Time average of the Poynting vector field. Dipole excitation with dipole at y=-1.2R , dipole along the y axis.