X-component of vector of intensity of electric field [V/m]
In this chapter we shall focus on the needs that many researchers, scientists and even students have very often. When using commercial simulation software for numerical simulation of electromagnetic field we frequently encounter many insufficiencies which those software products have. Usually, main aim of computational software developers is to optimize and refine so called core of these programmes – EM field solver. After that CAD (Computer Assisted Design) and post processing parts of EM simulators are dealt with. Mainly this can be an issue with newer, short-time in development products but one’s own post processing using Matlab can be greatly beneficial even when using well established simulators of EM field. This is largely due to its flexibility which cannot be overcome by any EM field simulator.
Throughout this chapter we will show you many ways of post processing which we use in our research of EM field in industrial and medical applications (Vrba et al., 2008; Vydra et al., 2011) and in primary research of EM field around living cells and structures inside cell bodies (Cifra et al., 2011; Havelka et al., 2011).
We hope that this chapter will aid many researchers and students in the vast field of EM research. Here, we present our knowledge and tips which we have gathered through our studies and our research.
1.1. Technical introduction
Generally, rough results we obtain using simulators of electromagnetic field - or from analytical solution of systems described by discrete elements - are in the form of complex vector components of intensity of electric and magnetic field (i.e. time dependent – periodical – components in the directions of coordinate axes). We process these results using Matlab and interpret them to draw conclusions. In this chapter we would like to present basic processing of rough data, calculation of specific absorption rate and other parameters in particular regions of simulation domain, visualization of results in many ways (pcolor, slices, histograms, multiple iso-surface, surf interpretation on various shapes according to specific task etc.). We will provide detailed examples with practical applications and explanation of advantages provided by presented solutions.
2. Rough results from EM field simulator
As mentioned above, in this chapter we suppose that we have obtained rough data from any numerical simulator of EM field and now we want to interpret them. First of all we should look at how the structure of this data looks like. To get the full understanding we shall briefly go through some EM field basics.
2.1. EM field basics
Electromagnetic Field can be described using well known Maxwell’s equations (for more information on Maxwell’s equations please refer to any book dealing with EM field theory).
Simply by solving those equations EM field can be completely described at all points of space and time. This leads us to complete description of EM field using only phasors of intensity of electric and magnetic field E and H (or D and B where D = εE and B = µH). This means that output of conventional commercial simulator is in the form of time dependent vectors that have components in axis x, y and z. These vectors are defined for each part of computational domain (e.g. when using FDTD (Thomas et. al., 1994), vectors are defined for each voxel – block discretizing computational domain).
We can see that this type of data can be extracted in form of matrices (multi-dimensional, e.g. 4D). Now, we shall look closer at those matrices.
2.2. Data structure
As we mentioned in previous chapter, results from simulators of EM field are represented as matrices, which directly predestines them to be processed in Matlab, which is the perfect tool for matrix operations.
There is a sample of data obtained from simulation in the following table. It depicts x-component of vector of intensity of electric field [V/m] in y-axis section in a part of some model.
Following graphical representation can help us shed some more light on the structure of data we obtained. These data are represented as four dimensional matrices (for phasors E and H separately) depicting whole computational domain and they are time dependent.
Generally we can describe phasors as follows.
Note: It may be necessary to convert data to suitable matrix form (e.g. rough data are in the form of a row vector with axial information for each element). We will look into it in the chapter 3.
Now that we know what our data source looks like we can simply process it to view the results and highlight some of their aspects according to our needs (see Table 2. for axial information).
3. Viewing the results
In this section we are going to show some examples of how obtained data can be viewed, how to interpret those results, what type of projection should we use etc. We shall illustrate this on some practical examples of EM field applications.
3.1. Basic transformation of rough data
As mentioned above we might obtain rough data in the form of a row vector. Let us illustrate this in this simple example. Our computational domain is 2 by 2 by 2 thus obtained row vector (x-component, apmlitude) has 8 elements. See Table 2.
|Vector of values||8||6||5||2||5||6||8||9|
From this we can extract axis. As long as we do not know the length of each axis we need to utilize this method to find actual axial data (see Figure 2.).
Note that actual axis arrangement can be different in your case (e.g. x and y axis arrangement may be commuted). Thus it is vital to get familiar with axial arrangement in your exporter from EM field simulator.
Furthermore, actualAxis vectors are underlined because they are growing on every loop iteration. Since axial vectors are not usually very long, this poses only mild concern. They cannot be preallocated because we generally do not know their actual length. If you are expecting very long axial vectors you may consider preallocating them safely longer than your expectations and then just using part of them which is non-zero (you may select this part of a vector using find – please refer to Matlab documentation).
Now that we know the length of each axial vector we may sort the vector of exported values and transform it to a matrix which will better represent three dimensional nature of our computational domain and will allow us to plot data with axial information. In this case we can very well utilize reshape which is built in Matlab.
MATRIX = reshape(vector_of_values,lenght(X),lenght(Y),lenght(Z))
Thanks to this process we now have every component of each vector (i.e. E and H) represented as three dimensional matrix and we can utilize it further.
3.2. Basic plotting of data
First of all, we need to bear in mind that we have time-dependent data. The most basic process is to plot actual situation (distribution of intensity of electric or magnetic) at a given time, or amplitude of vector. (In some applications we may need to plot just one component of this vector. This is even simpler because then we can disregard following method.)
Phasor of intensity of electric or magnetic field can be represented by modulus and phase or real and imaginary part. We need to merge all the components of the vector and obtain real and imaginary part. This can be simply done (i.e. vector adding component matrices together). Then we have one matrix of complex numbers. We can choose specific time in which we need EM field to be plotted simply by adding <0,2π> to the phase of each vector and then we can plot real and imaginary modulus of the vector in specified time. Or we take just the amplitude of vectors and plot them.
It is very usual to plot RMS (i.e. Root Mean Square) value of vectors which is defined as follows.
This can be again obtained very simply from amplitude of intensity of electric field. (Note that this same procedure can be used also in the case of intensity of magnetic field H, usually in applications involving heating and/or drying we deal only with intensity of electric field E, because it is the source of heat generation in exposed samples.)
Now that we have three dimensional matrix of values of RMS|E| we can plot it to see what our results look like. In the following section there is an example we prepared to illustrate how the results can be viewed and interpreted.
3.2.1. Example of basic data plotting
In this section we shall extract data from a simulation which has setup according to the Figure 3. Please note that this is just an example without any practical use, it serves only as an illustration.
We simulated simple section of a waveguide (inner dimensions 100x50x200 mm) with one side shorted (there is an excitation probe in form of a cylinder in the distance of 17 mm from the shorted end) and the other side open (absorbing boundary condition – absorbs 99.9% of incident power). We extracted the data and plotted them using Matlab.
As you can see (Fig 4ab.) we used several colormaps which enable us to highlight different aspects in our results interpretation. Sometimes it is needed to have contrast colormap (jet, lines – for more information see Product Help of Matlab), at other circumstances you may need to use fine and moderate colormaps (hot, gray, bone, pink – for more information see Product Help of Matlab). In the fourth graph in Fig. 4 we used different shading – faceted. This enables us to highlight structure of computational grid. In many commercial simulators of EM field parts of a model need to be meshed finer than others (e.g. in our case the excitation probe needs to be meshed four times more than the rest of the waveguide to be voxeled sufficiently). In other examples we used shading interp to get more clear view.
We can also utilize custom colormaps. This can be exceptionally beneficial in applications where we need to find out where values are at some critical level or higher. We illustrated this feature in the first image in Fig. 4a. In Fig. 5. there is an Colormap Editor which can be accessed through: Figure – Edit – Figure Properties – Colormap pull-down menu – Custom.
In our example we set segment in the middle to black colour and segment next to it to white colour. This resulted in the graph as seen in Fig. 4a. For more information on colormaps please refer to the Product Help of Matlab.
Furthermore, we illustrated how real modulus of vector of intensity of electric field at phase 0° is interpreted using Matlab (second image in Fig. 4a.). This is the most basic interpretation of obtained data we can do.
Note that this kind of results interpretation is much more flexible than the interpretation allowed by post processing tools in commercial EM simulators. In the following example we shall show how to work with time dependency of phasors. Since the results of EM field simulator are extracted when the steady state is reached time dependency is reduced to angle of phasors depicting the field of vectors. Through the following method we can alter phase of those phasors and show real part and imaginary part through one period. The results can be seen in Fig. 6. Figure 7. shows example of data processing to achieve this.
Note: In many EM field simulators you may encounter various errors. Pay special attention to the data structure of your exported data since it may not be useful in the way we have shown here (e.g. real and imaginary parts are exported as absolute values so the vital information about phase is lost).
Note: In the part of the script (Fig. 7.) where lowerThan variable is used we are changing the range of values. Since there are parts of model where values of intensity of electric or magnetic field are very near to zero, minimum of these values in dB would be around -400 dB. This renders produced images useless (value range is huge but most of the relevant values are in the region <-50,0>). Using function find we identify indexes of elements with values lower than -50 dB and we replace those elements with value -50 dB. For more information on function find please refer to the Product Help of Matlab.
Note: In this example we use function subaxis which has similar usage as function subplot but allows users to set the layout of plots in the figure more accurately (options Padding, Spacing, Margin etc.). For more information on subaxis please see internet documentation.
3.2.2. Treatment efficiency analysis
In the following example we will go through one of many useful applications of EM simulators today – evaluation of hyperthermia cancer treatment. Generally, hyperthermia is a method through which tissue is overheated (usually using microwave energy) and cells die (principle of this method is that energy is focused into the cancerous tissue which is less perfused and thus it is more heated – temperature in treated area rises above levels that trigger cell apoptosis). For more information on microwave hyperthermia see for example (Vrba & Oppl., 2008).
Model of simulated experiment can be seen in Figure 8. In this example we use waveguide applicator which is fed by coaxial line ending in protruding inner wire which is located near the shorted end of the waveguide section. There is a horn aperture which helps focus microwave energy to the desired area.
Relatively complex and complicated structure of human body is replaced with so called “phantom” that represents simple muscle tissue. In the model there is a tumor located 1 cm below the surface of phantom. This tumor has the same dielectric parameters as the surrounding muscle tissue (usually the only difference in simulations between muscle tissue and tumorous tissue is in their perfusion, heat transfer rate and heat generation rate – generated heat in tumorous tissue is usually transferred slower than in physiological surrounding tissue).
Additionally there is a water bolus which serves as a coolant body protecting the surface tissue of patient and moving the maximum of temperature to the lower layers. In this example we use source at 434 MHz and the applicator is filled with water (required dimensions of the applicator are effectively lower and impedance matching between waveguide-bolus and body are much better). For more information on waveguide hyperthermia applicators please refer to (Vera et al., 2006).
In this simulation we again extract rough data from EM simulator and we process them further using Matlab. We again have intensity of electric and magnetic field defined in every element of the model. For effective analysis of the treatment we need to evaluate SAR (Specific Absorption Ratio) which describes how much power is absorbed in a weight unit [W/kg]. For more information on SAR see http://www.ets-lindgren.com/pdf/sar_lo.pdf.
To determine how SAR is distributed we need to use this formula.
As we see, in this case we can very well utilize RMS|E|. From formula (3) we can say that SAR is defined by following equation.
SAR is thus depending on RMS value of intensity of electric field, on conductivity of a material and on its density. We need to obtain those values somehow. In our case whole phantom has homogeneous density and electric conductivity thus we can only mask matrix of RMS|E| and multiply each element by coefficient produced by ratio of σ and ρ.
Note: There are other options how to do that. For example some EM simulator allow users to export matrix of dielectric parameters of a model (this will produce similar matrix to our masking matrix).
Now we should look into how to produce masking matrix efficiently. From our model we know that every element of the voxelized model is phantom if its value of y axis is higher than 120 mm (we can find out this in CAD part of EM simulator). We can prepare our three dimensional matrix of RMS|E| according to previous sections. We also have actual axis of this matrix which we obtained from exported file as well. We simply need to identify which element is at the 120 mm value and then we shall mask the elements with lower value of y axis.
An elegant way to find the index at y axis of element which has the value nearest to 120 mm (note that the value is not usually exactly 120 so if we looked for the exact match we would probably fail) is to use function min. Example of this method is as follows.
In this example, we can illustrate one other useful feature of Matlab processing of EM simulator results. We can instantly normalize obtained results. Usually it is possible to extract actual power in the source which was simulated (e.g. 0.002496 W) and we can utilize this to obtain normalization coefficient. Power normalization coefficient is given by the following formula.
Since this is the coefficient of power ratio we need to use coefficient of intensity of electric field ratio which is the square root of coefficient of power (see the equation below).
Simply by multiplying the matrix of data by this coefficient we obtain normalized data to the desired value of delivered power.
Now we can show several interpretations of results we obtained (Z section in the middle of the model). Thanks to Matlab’s flexibility and versatility we can highlight the results in the way that is far more efficient than by using post processing of the EM field simulator. In the following figures we show intensity of electric field and SAR in the treated area.
As you can see in figure 11. each colour scheme highlights different thing (as for the first three). In the fourth graph we used custom colormap to show some kind of a critical zone. In this example we expected that higher values of SAR than 1000 W/kg are considered to be dangerous in some regions (note that these critical limits differ significantly in every application – many factors contribute, e.g. vital or sensitive organs near the treated area etc.) and thus we highlighted the zone with higher values using bright red colour. This may be very important but commercial simulation software generally omits such option.
It may be also very useful to show obtained data only in the region of the tumour or to show values in form of several slices, semi-transparent layers or iso-surface view. Now we shall show how to prepare masking matrix for some simple geometrical objects representing tumours. For more information on the mentioned advanced techniques of results representation please see section 4. Advanced Viewing Techniques.
3.2.3. Preparation of masking matrix
Very simple masking (i.e. when the region that needs to be masked is in the form of a cuboid) was shown in previous section. Now we should look into masking regions round in shape (i.e. spheres, ellipsoids and other common shapes that parts of model can be represented by).
Spheres are case of ellipsoid and we shall treat them as such. So our main aim now is to mask region generally in ellipsoidal form. Equation defining an ellipsoid (with its origin at xc, yc, zc) in the three dimensional coordinate system is shown below.
The graphical representation of a general ellipsoid can be seen in figure 12. There we can see that a, b and c are semi-principal axes of ellipsoid and define its dimensions.
Albeit there is a built-in function designed to generate ellipsoid (see the Product Help of Matlab) we shall show you an easy way which allows you to generate desired masking matrix (e.g. with values 0 or 1).
As you can see in Fig. 13. we simply utilized the formula representing ellipsoid and all elements whose coordinates meet the given restrictions are filled with ones. The other elements remain zero. Constants elementsX, elementsY or elementsZ define the length of axis of computational domain (i.e. actualAxisX and so on). Variables xc, yc and zc define centre of our ellipsoid.
Now we can simply use generated masking matrix and multiply every element of result matrix by according element of masking matrix. Then we can look how our results are interpreted in detailed view in the tumour.
From masked matrix of SAR we can also easily extract total power radiated to the tumorous region. Total power lost in the healthy tissue is than given by the difference between total power (in the case of perfectly matched source) and power lost in the tumour.
Note: In some cases it might be beneficial to set the colorbar range the same in all pictures (as shown in Fig. 14.). The simplest way to do this is to set value of [1,1] element of displayed matrix (i.e. a slice) to the maximum value of the whole three dimensional matrix (if needed, the minimum value should be set to element [max,max] of displayed matrix).
Additional information that we can obtain from masked result matrix is total absorbed power in the tumour. Since we know the dimensions of the tumour we can easily calculate its volume and use following formula.
The volume of standard ellipsoid is given as follows.
For example in our case tumour has volume only 8.37766e-6 m3. Information on total absorbed power may be very efficient way of preliminary evaluation for the treatment planning. Through this method we can determine power lost in any part of a simulated model. There may be some intricate volumes which are not as easily described as ellipsoid. Then we need to export their volume or masking matrix from the EM simulator (if possible).
We can determine the volume of such complicated shapes by summing elementary volumes of voxels representing those shapes. This may be unnecessary in some more advanced EM field simulators since they allow us to export model data in many suitable forms for Matlab processing. But the method presented in the following text is universal and can serve for better understanding of the model, its grid and working principals of EM simulations.
First of all we need to define which voxels are occupied by the model we want to evaluate. For this example we shall use our previously generated ellipsoid. As you can see in Fig. 14. generally grid of a simulation does not have to be symmetrical (i.e. voxels are not cubes but they are in the form of general cuboids). This means that each element may be representing different volume. We have our matrix of zeros and ones which we generated in Fig. 13. Now we need to apply this matrix to the actual coordinate system.
We can use built-in function meshgrid to produce three dimensional matrices which contain coordinates for specified element (for more information on function meshgrid please refer to the Product Help of Matlab). Now we can easily determine exact coordinates of elements occupied by our model (see Fig. 15.).
At this moment we need to get better understanding of how models are defined. Since the values (in this case ones) are located in the nodes of the grid not in the centers we need to determine how the model actually looks like. See the following figure for better understanding (2D-plane is used to illustrate).
As you can see there will be an error caused by this node to cell transformation. When the simulation grid is defined appropriately (i.e. it is fine enough) the error will be only marginal. If such error is unacceptable more advanced techniques of node to cell transformation should be used. But for purposes of this example this method is more than sufficient. For example, volume of previously defined ellipsoid determined voxel by voxel is 8.5302e-6 m3. This means that the relative error is 1.8 % (and this error includes error caused by voxelization itself – this means that node to cell transformation really brings only marginal error in simulations with fine grid).
In the Fig. 16. you can see that node to cell transformation can be done in a few (precisely 8) ways. We can easily determine which way is the most precise (i.e. [(A+1) – A] or [(A-1) – A] and its combinations with B and C). Since our model is voxelized symmetrically there is almost no difference between volumes computed accordingly to various node to cell transformations (ranging between 8.5302e-6 and 8.5304e-6).
4. Advanced viewing techniques
This section deals with some examples of advanced viewing of results that we have found useful during our work. As you have noticed, up to this point we have represented 3D results only in 2D graphs. This is because usually these 2D representations are clearer and easier to interpret. But for overviewing or demonstrational purposes we might need to show whole situation in 3D. For doing this we shall use multiple sections of a computational domain (slice or surf) and iso-surfaces.
4.1. Multiple layer viewing
We can display results in multiple 2D sections with variable transparency. This method of results visualization can give us much better overview of how a situation looks like.
It can be also beneficial to use surface view as semi-transparent layers to view results in semi-masked way. Through this we can highlight some regions of model without completely blocking visual output from other regions. For example, in our case, we can use one surface view of RMS|E| in agar phantom and semi-mask it with alphamap of tumour (values 1) and the rest of tissue (values 0.5).
Figure 19. shows us how this graph is plotted and how transparencyMatrix is prepared.
Note that when using this type of transparency mapping maximum and minimum values must be present in transparencyMatrix (i.e. one edge of transparencyMatrix is set to 0 to denote the minimum value of transparency).
Through this method we can produce much more intricate visualizations of results which are unseen in commercial EM simulators. For example we can mask results in uninteresting regions (immediate vicinity of power source etc.) almost entirely, semi-mask results in exposed tissue and left results in tumour unmasked, vital organs and other key regions of simulation domain. Note that masking results with transparency matrix does not alter presented values.
This may pose a problem when values around a power source are extreme and render the rest of results unclear. For this purpose you may consider utilizing simmilar approach as presented in Figure 7. Instead of finding values lower than some value you may find values higher than some reasonable value and lower all higher values to this level.
Many more possibilities are offered thanks to post processing of EM results in Matlab. As mentioned before, the greatest advantage is that it is extremely flexible and can meet very specific requirements which can arrise with various applications of EM field. In the following section we shall demonstrate some results obtained during our primary research of EM field of microtubules (nanostructures in living cells which serve as a crude frame of a cell and have other important roles in life cycle of a cell).
4.2. Iso-surface viewing
In this section we are going to show several results of EM field around microtubule. This structure is generally consisting of protofilaments which are polymerized tubulin heterodimers. Thirteen protofilaments bound together form a microtubule structure which resembles a long hollow tube (see Fig. 20).
Tubulin heterodimers (i.e. basic building block of microtubules) are highly polar structures and provided some form of external energy (movement of microtubules, mechanical vibrations etc.) they may produce EM field around themselves (Pohl, 1981; Fröhlich, 1978).
In our work EM field simulations of microtubule model were entirely conducted using Matlab. Tubulin heterodimers were represented as vibrating elementary electrical dipoles (EED) and EM field was determined for each of these EEDs. Combining the results led to unravelling the EM field produced by these complicated structures at whole (Havelka, 2009).
For the purpose of this text we shall look at visualization of these results we used to present obtained EM field. Because microtubules are symmetrical structures we have found out that representing its field by iso-surface view is very clear and easy to interpret.
In the figure 23. you can see part of the code used to generate such visualization. More information on isosurface can be found in the Product Help of Matlab. In this particular example we wanted to view results in the form of several iso-surfaces. First of all we determined the range of data obtained through Matlab analysis of EM field around our sample microtubule.
Then we need to choose which values we want to be visualized as iso-surfaces (in this case, values range is extremely wide therefore we choose only exponents – 7, 5, 2). Then we need to find which positions at colour scale are lower than our actually viewed value (we generate the colour scale using function jet(length of colour scale)). Then we simply use function patch (to build 3D wire model of locations with desired values – i.e. 1e2, 1e5, 1e7) and we set its colour which we choose accordingly from our generated colour scheme.
Additionally you can choose lighting (particularly useful in this case is gouraud lighting which does not produce any glances on multiple iso-surfaces). It is also beneficial to use alpha which is lowered with each loop (to retain clarity of the visualization).
Last part is generating tiny triangle which is concealed in the middle of viewed data and contains the minimum and the maximum values which allows us to show colorbar appropriately for all the iso-surfaces.
In this work we have presented basic methods of how to process obtained rough data from commercial EM simulators or even process data from our own simulations done in Matlab. We show how Matlab processing can be greatly beneficial in highlighting results in many ways that are unseen in commercial EM simulators.
We present very simple way to modify data in the form suitable for further processing and then we illustrate how to view these data in ways highlighting specific aspects (i.e. values in specific regions (tumorous tissue), evaluation of treatment efficiency, utilization of Matlab in primary research, 3D viewing etc.). These innovative ways of combination of specialized software with researcher’s versatile tool, such as Matlab is, yield in very productive and efficient way of scientific exploration of the vast field of electromagnetism.
The results that are obtained in our research of EM field around living cell help us understand crucial facts about this part of our lives which has to be truly discovered yet. Matlab in this instance allows us to visualize result so we can support or on the other hand disapprove many theories (e.g. transportation of particles around microtubule via EM field).
This research is supported by the Grant Agency of the Czech Republic by project 102/08/H081: "Non Standard Applications of Physical Fields" and by CTU project SGS12/071/OHK3/1T/13: “Advanced techniques in the processing of industrial materials and biopolymers using electromagnetic field - multi-frequency processing “. Further it is supported also by research program MSM6840770012: “Transdisciplinary Research in the Area of Biomedical Engineering II” of the CTU in Prague, sponsored by the Ministry of Education, Youth and Sports of the Czech Republic.